My first job was at Lehman Brothers Fixed Income Research in Andy Morton’s group. I was hard-working and naive about what was required then for vertical ascent. I came with mathematical background but untrained in science and I spent my time coding. My perspective has matured regarding science since, so now I realize the hard work of getting nature to follow the quantitative models now attempting to get a static Einstein universe to produce the CMB anisotropy. What is clear is that the hard part of numerical work remains as much drudgery as I had in the trading floor of Lehman, 13 hours in the neon light punctuated by dark skies. What has not changed when making Octave models for CMB anisotropy (as opposed to the derivative pricing models) is that the annoyance of the coding problems is dwarfed by the difficulty of producing a model that can fit the observed anisotropy. Tonight, I can’t solve the problem: find a distribution of point masses on a three-sphere such that the gravitational redshifts alone may reproduce the anisotropy. So much for the coding problem.

The sample-frequency-concentrated-distribution looks quite good but not in a realistic FREQUENCY SCALE.

pkg load statistics pkg load gsl pkg load image global n=10000 global q0=100 global points=zeros(n,2) global q0 global n global X = mvnrnd( [0,0,0,0],diag([1,1,1,1]),n); for i=1:n, X(i,:)=X(i,:)/norm(X(i,:),2); endfor # stereographic projection + normalize function z=pwrrnd(k,n) u=unifrnd(0,1,n) z=(1-u)^(1/(1-k))*0.001 endfunction global masses = lognrnd(0,4,n,1)*50; function d=direction(x1,x2) d=x2-x1; d=d/norm(d,2); endfunction function d=distance(x1,x2) d=acos(dot(x1,x2)); endfunction function f=gforce( x1,x2, m1,m2) f=m1*m2*direction(x1,x2)/distance(x1,x2)^2; #f=0 endfunction timestep=10 if 0, for a=1:n, force_on_x = 0; x1 = X(a,:); m1 = masses(a); for b=1:n, x2 = X(b,:); d=distance(x1,x2); m2=masses(b); if d>0 && d<pi, force_on_x = force_on_x + gforce(x1,x2,m1,m2); endif endfor x = x+0.5*force_on_x*timestep^2/m1; x = x/norm(x,2); X(a,:)=x; endfor endif global z = zeros(n,3); for i=1:n, y=X(i,:); y1=y(1); y2=y(2); y3=y(3); y4=y(4); zp=[ y1/(1-y4), y2/(1-y4),y3/(1-y4) ]; z(i,:) = zp/norm(zp,2); # select from log-normal #u = lognrnd(0,1); #if u>5, #z(i,:) = [0,0,0]; #endif endfor function t=fitangle(x) t=x; while t>2*pi, t=t-2*pi; endwhile while t<0, t=t+2*pi; endwhile endfunction function b=bin_s2point( z ) if abs(z)<0.001, b=[0,0]; return endif global q0; theta = atan( z(2)/z(1)); theta=fitangle(theta); phi = atan( sqrt(z(1)^2+z(2)^2)/z(3)); phi = fitangle(phi); thetan = floor(theta*q0/(2*pi)); phin = floor(phi*q0/(2*pi)); b=[thetan,phin]; endfunction # create a spherical grid in longitudes/latitudes global g = zeros(q0,q0); global refpoint = [1,0,0,0] function dox() global n; global g; global X; global masses; for i = 1:n, global refpoint; global z; global q0; global points; dg = zeros(q0,q0); cds = bin_s2point(z(i,:)); cds(1)=mod(cds(1),q0); cds(2)=mod(cds(2),q0); disp([cds(1),cds(2)]) points(i,:)=cds; fflush(stdout) #sdir = sign(direction(refpoint,X(i,:)); dist = distance(refpoint,X(i,:)); #g(cds(1)+1,cds(2)+1) = g(cds(1)+1,cds(2)+1)+1; m1=masses(i); dg(cds(1)+1,cds(2)+1) = m1/dist; dg = imsmooth( dg, "Gaussian", ceil(2),ceil(2)); #disp(sum(sum(dg))) #fflush(stdout) dg(1,1)=0; g = g+dg; endfor endfunction dox() sum(sum(g)) #g = (g-min(min(g))); g=g/sum(sum(g)); # value at a spherical harmonic function v=sphCoeff( l,m, g) v = 0; q0=20; for j=1:q0, for k=1:q0, w1 = cos(2*pi*j/q0); w2 = exp(1i*2*pi*k/q0); fnval = gsl_sf_legendre_sphPlm( l, m, w1 )*w2; if ~isnan(fnval*g(j,k)), v = v + fnval*g(j,k); endif endfor endfor disp([l,m,v]) fflush(stdout) endfunction function v=sumSphCoeff( l, g) v=0; for m=0:1 v = v+sphCoeff(l,m,g); endfor v=abs(v)^2 endfunction valsByL = zeros(1,1000); for k=1:1000, valsByL(1,k)=sumSphCoeff(200*k,g); endfor plot( 1:1000, valsByL);

## Leave a Reply