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;
}
|