Home Manual Packages Global Index Keywords Quick Reference ``` #include "legndr.i" #include "Chris/utls.i" #include "Eric/fft_utils.i" #include "Chris/randfield.i" #include "Eric/random.i" #include "Eric/histo.i" #include "Chris/linterp.i" require, "gamma.i" /*------------------------------------------------ // test alm transform tt0=1+Ylm(3,2,span(0,pi,50)(,-:1:50),span(0,2*pi,50)(-:1:50,)); alms= Falm(tt0,15) tt1= almF(alms,25); INFO,abs(tt1-tt0); // test cl generation cl2= array(0.,50); for(i=1;i<=50;i++){alms=alm_rand(array(1.,50));cls= almCl(alms); cl2+=cls;} plg,cl2/50; // test equal energy per scale PS func f1 (x) { return 1./x/(x+1.);} alms= alm_rand(f1(indgen(25))); plaitof,float(almF(alms,50)); -------------------------------------------------*/ //_YLMMAX =50; // max order in presaved alms func Factorial (n) { require, "gamma.i"; return exp(ln_gamma(n+1)); } //func Factorial(l,m) //{ // require, "gamma.i"; // return exp(ln_gamma(l+1)-ln_gamma(m+1)-ln_gamma(l+1)+ln_gamma(m+1)); //} func Ylm lm(l,m,theta,phi) /* DOCUMENT Ylm spherical harmonic: Ylm(l,m,theta,phi) */ { return ylm_coef(l,m)* legndr(l,m,cos(theta)) * exp(1i*m*phi); } func genYlm (nn) { /* DOCUMENT generates Ylm tables for latter use (dimensions nn*nn). */ lmax=_YLMMAX; theta=span(0,pi,nn)(,-:1:nn); phi=span(0,2*pi,nn)(-:1:nn,); ylms=array(complex,lmax,2*lmax-1,nn,nn); // cant be bigger because of n! for(l=0;l<=lmax-1;l++) { for(m=-l;m<=l;m++) { u=Ylm(l,abs(m),theta,phi); ylms(l+1,m+l+1,,)=((m>0)?u:((-1)^m*conj(u))); } } save,createb("/home1/pichon/Project/flows/data/ylm"+pr1(nn)+".yor"),ylms; return ylms; } func alm_testF (alm,nn) { /* DOCUMENT field=almF(alm,nn) returns the recomposed field (dimensions nn*nn) using alm. field is of the form theta along absiss phi along ordinate */ lmax=dimsof(alm)(2); if (nn==50) upload,"/home1/pichon/Project/flows/data/ylm50.yor",s=1; else if (nn==100) upload,"/home1/pichon/Project/flows/data/ylm100.yor",s=1; else if (nn==150) upload,"/home1/pichon/Project/flows/data/ylm150.yor",s=1; else if (nn==200) upload,"/home1/pichon/Project/flows/data/ylm200.yor",s=1; else if (nn< 300) { write,"you might consider 50,100,150,or 200 for speed"; ylms =genYlm(nn); } res=array(complex,nn,nn); pp =dimsof(alm)(2); res = (alm*ylms(1:pp,1:2*pp+1,,))(sum,sum,,); return res; } func almF (alm,nn) { /* DOCUMENT field=almF(alm,nn) returns the recomposed field (dimensions nn*nn) using alm. field is of the form theta along absiss phi along ordinate */ theta=span(0,pi,nn)(,-:1:nn); phi=span(0,2*pi,nn)(-:1:nn,); res=array(complex,nn,nn); lmax=dimsof(alm)(2); for(l=0;l<=lmax-1;l++) { for(m=-l;m<=l;m++) { u=Ylm(l,abs(m),theta,phi); res+=((m>0)?u:((-1)^m*conj(u)))*alm(l+1,l+1+m); } } return res; } func almF3D (alm,nn, xyz=,rmax=) { /* DOCUMENT field3D=almF3D(alm,nn) returns the 3D recomposed field (dimensions nn*nn*nn) using alm. field is of the form theta along absiss phi along ordinate r along alt ordinate EXAMPLE nn=10; #include "Dom/Density_utils.i" x=span(-1,1,nn)(,-:1:nn)(,,-:1:nn); y=span(-1,1,nn)(-:1:nn,)(,,-:1:nn); z=span(-1,1,nn)(-:1:nn,)(-:1:nn,,); u=1/sqrt((x-1.1)^2+(y-1.1)^2+0.1^2); v=Projec_Halo(u,[0.,0,0],[x,y,z],1.,sz=nn*2); alm=Falm(v,nn); u1=almF3D(alm,nn,xyz=1,rmax=1).re; w=where(x^2+y^2+z^2<1); pl,(u1 -u)(w)(*); pl,(u1)(w)(*),color=-6; */ if (is_void(rmax)) rmax=0.5; if (is_void(xyz)) { theta=span(0,pi,nn)(,-:1:nn); phi=span(0,2*pi,nn)(-:1:nn,); r=span(0.001,1,nn)(-:1:nn,)(-:1:nn,,); } else { x=rmax*span(-1,1,nn)(,-:1:nn)(,,-:1:nn); y=rmax*span(-1,1,nn)(,-:1:nn)(-:1:nn,); z=rmax*span(-1,1,nn)(-:1:nn,)(-:1:nn,,); theta= atan(sqrt(x^2+y^2),z); phi =atan(x,y); r =sqrt(x^2+y^2+z^2); } res=array(complex,nn,nn); lmax=dimsof(alm)(2); for(l=0;l<=lmax-1;l++) { for(m=-l;m<=l;m++) { u=Ylm(l,abs(m),theta,phi); u= u(,,-:1:nn); u *= r^l; res+=((m>0)?u:((-1)^m*conj(u)))*alm(l+1,l+1+m); } } return res; } func _PLM (l,m,x) /* DOCUMENT legendre function which also takes m<0 */ { res=legndr(l,abs(m),x); if (m<0) { res*=double((-1)^m)*Factorial(l+m)/Factorial(l-m); } return res; } func Falm (field,lmax,nn_out=,logg=) { /* DOCUMENT Falm(field,lmax,nn_out=,logg=) returns the field decomposition in Ylm nn_out == >dimension of the interpolated field logg=1 => log(field) instead of field */ dimfield=dimsof(field)(2); alm=array(0.+0.i,[2,lmax,2*lmax+1]); if (logg==1) field=log(10+field); if (is_void(nn_out)) nn_out=dimfield*4; if (nn_out!=dimfield) { tot=nn_out/dimfield; fieldre=cinterp2(double(field),n=tot); fieldim=cinterp2(im(field),n=tot); field=fieldre+1i*fieldim; } TFfield=fft(field(,1:(nn_out-1)),[],[1]); //TFfield=roll(TFfield,[0,dimsof(TFfield)(3)/2]) Vm=fft_indgen(nn_out-1); //Vm=roll(Vm); zero=where(Vm==0); theta=span(0,pi,nn_out); for(l=0;l<=lmax-1;l++) { for(m=-l;m<=l;m++) { PL=_PLM(l,m,cos(theta)); wm=where(m==Vm); temp=TFfield(,wm)*2*pi/double(nn_out-1); temp2=sin(theta)*temp*PL; Pre=integ(temp2(,1).re,theta,theta(nn_out)); Pim=integ(temp2(,1).im,theta,theta(nn_out)); //pf= ylm_coef(l,m); pf=sqrt((2*l+1)/(4*pi)*Factorial(l-m)/Factorial(l+m)); alm(l+1,l+1+m)=(Pre+1i*Pim)*pf; } } return alm; } func almCl (alm) { /* DOCUMENT almCl(alm) returns the power spectrum deduced from alms. */ //nn=100; if(dimsof(alm)(1)==2) { lmax=dimsof(alm)(2); cl=array(double,lmax); cl =(abs(alm)^2)(,sum); cl/=4*pi*((indgen(lmax)-1)*2+1); return cl; } else if(dimsof(alm)(1)==3) { cl=array(double,dimsof(alm)(2),dimsof(alm)(4)); cl=(abs(alm)^2)(,sum,); cl/=4*pi*((indgen(dimsof(alm)(2))-1)*2+1)(,-:1:dimsof(alm)(4)); return cl; } } func alm_cor (alm,nn=) { /*DOCUMENT alm_cor(field,nn=) returns the auto-correlation function of a field using its decomposition in spherical harmonics. nn is the number of points for theta. */ // if (is_void(nn)) nn=100; lmax=dimsof(alm)(2); w=array(double,nn); theta=span(-pi,pi,nn); cls=almCl(alm); for(l=0;l<=lmax-1;l++) w+=double(cls(l+1)*(2*l+1))*legndr(l,0,cos(theta)); return w; } func alm_rand (s) { dims=dimsof(s)(2); l1=span(0,dims-1,dims); almr=array(complex,dims,2*dims+1); s=sqrt(4.*pi*s/(2.-1./double(2*l1+1.))); almr(1,1)=s(1)*random_n(1); for(l=1;l<=dims-1;l++) { temp=s(l+1)*(random_n(l+1)+1i*random_n(l+1)); temp(1)=double(temp(1)); almr(l+1,l+1:2*l+1)=temp; almr(l+1,indgen(l))=(-1)^(indgen(l)(::-1))*conj(almr(l+1,l+1+indgen(l)(::-1))); } return almr; } func PS(x,l) { Ohm=0.001; // return (1-exp(-T*x^2));}; coeff=4*pi/float(2.-1./float(2*l+1)); T=10000 // return 1/float(l*(l+1.))*(0.01+(1/sqrt(l+1)*T*x)^2)^-1.5; return 1/float(l^2)*(0.01+(1/sqrt(l+1)*T*x)^2)^-1.5; //return 1/float(l*(l+1.))*coeff*(1-exp(-x/Ohm)); //return 1/float((l+0.01)^6)*coeff*(1-exp(-x/Ohm)); }; func almT (steps,lmax) { alms=array(complex,lmax+1,2*lmax+1,steps); alms(1,,)=0; for(l=1;l<=lmax;l++) { for(m=0;m<=l;m++) { shared_indx=[l]; arg_func=PS; z=genrandfield1D(steps,fun=foo); zr=fft(z,-1); z=genrandfield1D(steps,fun=foo); zi=fft(z,-1); alms(l+1,m+l+1,)=zr+1i*zi; } alms(l+1,indgen(l),)=(-1)^(indgen(l)(::-1))*conj(alms(l+1,l+1+indgen(l)(::-1),)); alms(l+1,l+1,)= double(alms(l+1,l+1,)); } //vv=array(float,lmax,lmax,nn);for(i=1;i<=nn;i++) vv(,,i)=float(almF(alms(,,i),lmax)); //uu=array(double,lmax+1,steps); //for(i=1;i<=steps;i++) uu(,i)=almCl(alms(,,i))*(indgen(lmax+1)-1)*indgen(lmax+1); //return uu; return alms; } func alm_tst (s,n,logl=,disp=,time=) { /*DOCUMENT alm_tst(s,n) compares Cls with s. n is the number of random realisations for alms. */ ws; clmr=array(double,dimsof(s)(2),n); l=span(0,dimsof(s)(2)-1,dimsof(s)(2)); for (i=1;i<=n-1;i++) { if(time) { almr=almT(1,numberof(s)-1); clmr(,i)=almCl(almr(,,1)); } else { almr=alm_rand(s); clmr(,i)=almCl(almr); } } print,"ok"; //ws,1; //plg,log(clmr(,avg)+0.001),log(l+0.001),color="red"; //plg,log(2*s*s/(4*pi)+0.001),log(l+0.001),color="blue"; cll=clmr(,avg); if(disp==1) { ws,2; if(logl==1) logxy,1,1; plg,l*(l+1)*clmr(,avg),l,color="red",width=5; plg,l*(l+1)*s,l; } return cll(rms); } func Projec (fonc,maxx,maxy,maxz,R200,x=,y=,z=,sz=,rescale=) { /* DOCUMENT Projec(fonc,maxx,maxy,maxz,R200,x=,y=,z=,sz=,rescale=) return "fonc" interpolated on a sphere centered in "maxx" "maxy" "maxz" with a radius of "R200". fonc is defined on x,y,z. sz=size of the projected map. rescale=1 if R200 is given using number of cases (and not in distance unit). */ //print,"Projec"; nn=dimsof(fonc)(2); if (is_void(sz)) sz=50; if(is_void(rescale)) rescale=1; th=span(0,pi,sz)(,-:1:sz); ph= span(0,2*pi,sz)(-:1:sz,); if (is_void(x)) { x=span(-1.,1.,nn)(,-:1:nn)(,,-:1:nn); y=span(-1.,1.,nn)(-:1:nn,)(,,-:1:nn); z=span(-1.,1.,nn)(-:1:nn,)(-:1:nn,,); } r=R200; if (rescale==1) r*=(max(x)-min(x))/nn; x1=r*cos(ph)*sin(th); y1=r*sin(ph)*sin(th); z1=r*cos(th); x2=x1+x(maxx,maxy,maxz); y2=y1+y(maxx,maxy,maxz); z2=z1+z(maxz,maxy,maxz); v=LInterp(fonc,x2,y2,z2,x=x(,1,1),y=y(1,,1),z=z(1,1,)); return v; } ```