Home Manual Packages Global Index Keywords Quick Reference ``` /* * randfield.i -- * * Useful routines for Gaussian random field operations in Yorick. * * Copyright (c) 1995-2000 Eric THIEBAUT. * * History: * \$Id\$ * \$Log\$ * *----------------------------------------------------------------------------- */ #include "Eric/random.i" func PS(s){return s=1/(0.001^2+s^(7/3.));} /* // compute the 1D powerspectrum from 3D realisations nr=100; nm=5; uu=array(0.,nr,nm); for(i=1;i<=nm;i++){ z=genrandfield3D(nr,fun=PS); zr=float(fft(z,[-1,-1,-1])); zi=fft(zr,[1,0,0])/nr^2; zi =abs(zi)^2; kk=(indgen(nr)-1)*2*pi/float(nr); zi= reform(zi,[2,nr,nr*nr])(,avg); plg,zi(2:),kk(2:); uu(,i)=zi; logxy,1,1; } plg,uu(2:,avg),kk(2:),width=4,color="red"; plg,PS(kk(2:))*kk(2:)^2/2./pi/(15/3.-1),kk(2:),width=4,color="blue"; print,((log(uu)(nr/6:nr/2,avg))(dif)/(log(kk)(nr/6:nr/2))(dif))(avg) */ /* func PS(s){return s=1/(0.001^2+s^(7/3.));} nr=200; nm=15000; uu=array(0.,nr,nm); zi=array(0.,nr); kk=(indgen(nr)-1)*2*pi/float(nr); for(i=1;i<=nm;i++){ z=genrandfield1D(nr,fun=PS); zi +=abs(z)^2/float(nm); } //plg,zi,kk; logxy,1,1; plg,zi(2:)(dif)/kk(2:)(dif)/kk(2:)(zcen),kk(2:)(zcen); plg,PS(kk(2:)(zcen))/kk(2:)(zcen)^2*(7/3.),kk(2:)(zcen),width=4,color="blue"; */ require, "Eric/fft_utils.i" func GRF_random m(u1, u2, u3, must_be_real=,fun=,param=) { if (is_void(u2 )) s= sqrt(u1*u1); else if(is_void(u3)) s= sqrt(u1*u1+u2*u2); else s= sqrt(u1*u1+u2*u2+u3*u3); if (is_void(fun)) { s=1/(0.001+s^(7/3.)); } else if (is_void(param)) { s= fun(s);} else {s=fun(s,param);} s= sqrt(s); if (must_be_real==1) return (s)*random_normal(numberof(u1),seed=randomize()); else return (s)/sqrt(2.)*(random_normal(numberof(u1),seed=randomize())+ 1i*random_normal(numberof(u1),seed=randomize())); } func genrandfield3D(n,fun=,param=) { /* DOCUMENT genrandfield3D(n); generates Gaussian Random Field of dims [3,n,n,n] which obeys the power spectrum given by fun (default k^-5/3) include, "Eric/histo.i"; ns=64; z=genrandfield3D(ns); zr=float(fft(z,[-1,-1,-1])); //pli,zr(,,5); z2=float(fft(z*conj(z),[-1,-1,-1])); //fft_pli,z2(,,5); x1=(indgen(ns)-1)*2*pi/float(ns); y1=(indgen(ns)-1)*2*pi/float(ns); z1=(indgen(ns)-1)*2*pi/float(ns); r1= abs(x1, y1(-,),z1(-,-,)); zc=array(0.,ns,ns,ns); for(i=1;i<=50;i++){z=genrandfield3D(ns); z2=float(z*conj(z));zc+=z2/50.;}; py= histo2(r1, px, weight=zc, average=1, interp=1,binsize=r1(dif)(min)); plg, py(2:), px(2:), color="red"; //azimuthal average plg, PS(px(2:)),px(2:),color="blue",width=4; logxy,1,1; write," power spectrum slope",(log(py(5:25))(dif)/log(px(5:25))(dif))(2:12:10)(avg) SEE ALSO fft_indgen, fft. */ require,"random.i"; if (dimsof(n)(1)==0) dimlist= [3,n,n,n]; else dimlist=n; z= array(complex, dimlist); // ******* genere un tableau hermitique ****** // 2-tableaux d'indices de frequence : up = indice de la // frequence correspondant a -u au sens de la FFT : u= array(0.,dimlist); u(*)= indgen(numberof(z)); up= fft_symmetric_index(dimlist); n1= dimlist(2); n2= dimlist(3); n3= dimlist(4); u1= 2*pi/n1*fft_indgen(n1)( ,-:1:n2,-:1:n3); u2= 2*pi/n2*fft_indgen(n2)(-:1:n1, ,-:1:n3); u3= 2*pi/n3*fft_indgen(n3)(-:1:n1,-:1:n2, ); if (is_array((i= where(u == up)))) { // cas des elements reels: z(i)=GRF_random(u1(i), u2(i), u3(i), must_be_real=1,fun=fun,param=param); } if (is_array((i= where(u < up)))) { // cas des elements complexes: tmp= GRF_random(u1(i), u2(i), u3(i), must_be_real=0,fun=fun,param=param); z(i)= tmp; z(up(i))= conj(tmp); } z(1,1,1)=0; return z; // sec3,float(fft(z,[-1,-1,-1])) } func genrandfield2D (n,fun=,param=) { /* DOCUMENT genrandfield(n); generates Gaussian Random Field of dims [2,n,n] which obeys the power spectrum given by fun (default k^-5/3) EXAMPLE include, "Eric/histo.i"; ns=128; z=genrandfield2D(ns); zr=float(fft(z,[-1,-1])); //pli,zr; z2=float(fft(z*conj(z),[-1,-1])); // pli,z2; zc=array(0.,ns,ns); for(i=1;i<=50;i++){z=genrandfield2D(ns); z2=z*conj(z);zc+=float(z2)/50.;}; //fft_pli,zc; x1=(indgen(ns)-1)*2*pi/float(ns); y1=(indgen(ns)-1)*2*pi/float(ns); r1= abs(x1, y1(-,)); py= histo2(r1, px, weight=zc, average=1, interp=1,binsize=r1(dif)(min)); plg, py(2:), px(2:), color="red"; //azimuthal average plg, PS(px(2:)),px(2:),color="blue",width=4; logxy,1,1; write," power spectrum slope",(log(py(5:25))(dif)/log(px(5:25))(dif))(2:12:10)(avg) SEE ALSO fft_indgen, fft. */ require,"random.i" if (dimsof(n)(1)==0) dimlist= [2,n,n]; else dimlist=n; dimlist= [2,n,n]; // ca marche en N-dims z= array(complex, dimlist); // ******* genere un tableau hermitique ****** // 2-tableaux d'indices de frequence : up = indice de la // frequence correspondant a -u au sens de la FFT : u= array(0.,dimlist); u(*)= indgen(numberof(z)); up= fft_symmetric_index(dimlist); n1= dimlist(2); n2= dimlist(3); u1= 2*pi/n1*fft_indgen(n1)( ,-:1:n2); u2= 2*pi/n2*fft_indgen(n2)(-:1:n1, ); if (is_array((i= where(u == up)))) { // cas des elements reels: z(i)=GRF_random(u1(i), u2(i), must_be_real=1,fun=fun,param=param); } if (is_array((i= where(u < up)))) { // cas des elements complexes: tmp= GRF_random(u1(i), u2(i), must_be_real=0,fun=fun,param=param); z(i)= tmp; z(up(i))= conj(tmp); } z(1,1)=0; return z; } func genrandfield1D (n,fun=,param=) { /* DOCUMENT genrandfield(n); generates Gaussian Random Field of dims [1,n] which obeys the power spectrum given by fun (default k^-5/3) EXAMPLE compute correlation function in 1D include, "Eric/histo.i"; nm=512; ns=512; z=genrandfield1D(ns,fun=PS); zr=float(fft(z,-1)); z2=float(fft(z*conj(z),-1)); zc=array(0.,ns); for(i=1;i<=nm;i++){z=genrandfield1D(ns,fun=PS); z2=z*conj(z);zc+=z2/float(nm);}; logxy,1,1; kk=2*pi*(indgen(ns)-1)/float(ns); plg,float(zc),kk; plg, PS(kk),kk,color="blue",width=4; write, " power spectrum slope",(log(float(zc)(5:125))(dif)/log(indgen(ns)(5:125))(dif))(avg) SEE ALSO fft_indgen, fft. */ require,"random.i" dimlist= [1,n]; // ca marche en N-dims z= array(complex, dimlist); // ******* genere un tableau hermitique ****** // 2-tableaux d'indices de frequence : up = indice de la // frequence correspondant a -u au sens de la FFT : u= array(0.,dimlist); u(*)= indgen(numberof(z)); up= fft_symmetric_index(dimlist); n1= dimlist(2); u1= 2*pi/n*fft_indgen(n1); if (is_array((i= where(u == up)))) { // cas des elements reels: z(i)=GRF_random(u1(i), must_be_real=1,fun=fun,param=param); } if (is_array((i= where(u < up)))) { // cas des elements complexes: tmp= GRF_random(u1(i), must_be_real=0,fun=fun,param=param); z(i)= tmp; z(up(i))= conj(tmp); } z(1)=0; return z; // z2=float(fft(z*conj(z),[-1,-1])) // pli,float(fft(z,[-1,-1)) } func fft_index (length) { x= indgen(0:(w=length/2)); if ((w+= 1-length)<0) grow, x, indgen(w:-1); return x; } func fft_roll (a) /* DOCUMENT fft_roll, a rolls values of a SEE ALSO fft_pli, roll. */ { if ((dims= dimsof(a))(1)==2) { dim1= dims(2); min1= (max1= dim1/2) - dim1 + 1; dim2= dims(3); min2= (max2= dim2/2) - dim2 + 1; return roll(a, [-min1, -min2]); } else { if ((dims= dimsof(a))(1)==3) dim1= dims(2); min1= (max1= dim1/2) - dim1 + 1; dim2= dims(3); min2= (max2= dim2/2) - dim2 + 1; dim3= dims(3); min3= (max3= dim3/2) - dim3 + 1; return roll(a, [-min1, -min2,-min3]); } error,"requires 2 or 3 dims"; } func fft_sec3 (a, scale=, legend=, hide=, top=, vmax=) /* DOCUMENT fft_sec3, a; Plot 3-D FFT array A as an set of images, taking care of "rolling" A and setting correct world boundaries. Keyword SCALE can be used to indicate the "frequel" scale along both axis (SCALE is a scalar) or along each axis (SCALE is a 2-element vector: SCALE=[XSCALE,YSCALE]); by default, SCALE=[1.0, 1.0]. KEYWORDS legend, hide, top, cmin, cmax. SEE ALSO fft_pli, roll. */ { animate,1; if (is_void(vmax)) vmax=max(a); if(is_void(sec)) sec=3; for(i=1;i<=dimsof(a)(sec+1);i++){ fma; if(sec==3) fft_pli,a(,,i),cmax=vmax, scale=scale, legend=legend, hide=hide, top=top; if(sec==2) fft_pli,a(,i,),cmax=vmax, scale=scale, legend=legend, hide=hide, top=top; if(sec==1) fft_pli,a(i,,),cmax=vmax, scale=scale, legend=legend, hide=hide, top=top; if(!is_void(inter)) { print,i; rdline,prompt="hit RET or Enter to continue";} pause,50;} animate,0; } func exemple_PS (dim,size,N=,plog=,fun=) { require, "Eric/histo.i"; /* DOCUMENT example compute power spectrum in dim DIMENSION*/ if (is_void(plog)) plog=1; if (is_void(size)) size=35; if (is_void(N)) N=1; slope=[]; for (m=1;m<=N;m++) { if (dim==1) { // champ de densite z=genrandfield1D(size,fun=fun); zr=float(fft(z,[-1])); if (N==1) {window,1;ws;plg,zr;pltitle,"champ de densite";} // funtion de correlation par corr_func zc=array(0.,size); for(i=1;i<=50;i++){z=genrandfield1D(size,fun=fun);zr=float(fft(z,[-1])); z2=corr_func(zr,px);zc+=z2(:-1);}; if (N==1) {window,2;ws;plg,zc/50.,px(:-1);pltitle,"corr function";} zc=float(fft(zc/50.,[-1]))/size; u1=fft_indgen(size)*2*pi/size;r1= abs(u1); print,"r1";INFO,r1; pk=px*2*pi/size; print,"pk:";INFO,pk; pycorr= histo2(r1, pk, weight=zc, average=1, interp=1, binsize=r1(dif)(min)); print,"pk:";INFO,pk; if (N==1) {window,3,style="win12.gs";ws;plsys,1;info,pycorr;info,pk; plg, pycorr, pk, color="red",type=2;} zc=array(0.,size); for(i=1;i<=50;i++){z=genrandfield1D(size,fun=fun); z2=float((z*conj(z)));zc+=z2;}; py= histo2(r1, pk, weight=zc/50., average=1, interp=1, binsize=r1(dif)(min)); if (N==1) { if (plog) {window,3,style="boxed.gs";ws; plg, log10(py(2:-2)), log10(pk(2:-2)), color="red"; if (is_void(fun)) plg, log10(PS(pk(2:-2))),log10(pk(2:-2)); else plg, log10(fun(pk(2:-2))),log10(pk(2:-2)); } else {plsys,2; plg, py/50., pk;} //azimuthal average plsys,1;if (plog) pltitle,"spectre de puissance"; else pltitle,"k^3^ P(k)";} } // example compute power spectrum in 2D if (dim==2) {print,"corr_func bizarre pour dim==2"; z=genrandfield2D(size,fun=fun); zr=float(fft(z,[-1,-1])); if (N==1) {window,1;ws;pli,zr;pltitle,"champ de densite";} z2=float(fft(z*conj(z),[-1,-1])); if (N==0) {window,2;ws;pli,z2;} z=genrandfield2D(size,fun=fun);zr=float(fft(z,[-1,-1])); z2=corr_func(zr,px);zc=z2; for(i=1;i<=1;i++){z=genrandfield2D(size,fun=fun);zr=float(fft(z,[-1,-1])); z2=corr_func(zr,px);zc+=z2;}; if (N==1) {window,2;ws;plg,zc/2.,px;pltitle,"corr function";} zc=float(fft(zc/2.,[-1]))/(size^2); u1=(fft_indgen(long(px(0)+1))-1)*2*pi/size;r1= abs(u1); pk=px*2*pi/size; pycorr= histo2(r1, pk, weight=zc, average=1, interp=1,binsize=r1(dif)(min)); if (N==1) {window,3,style="win12.gs";ws;plsys,1; plg, (pycorr(2:-2)), (pk(2:-2)), color="red",type=2;} u1=fft_indgen(size)*2*pi/size;r1= abs(u1, u1(-,)); zc=array(0.,size,size); for(i=1;i<=50;i++){z=genrandfield2D(size,fun=fun); z2=float((z*conj(z)));zc+=z2;}; py= histo2(r1, pk, weight=zc, average=1, interp=1,binsize=r1(dif)(min)); if (N==1) { if (plog) {window,3,style="boxed.gs";ws; plg, log10(py(2:-2)/50.), log10(pk(2:-2)), color="red"; if (is_void(fun)) plg, log10(PS(pk(2:-2))),log10(pk(2:-2)); else plg, log10(fun(pk(2:-2))),log10(pk(2:-2));} else {plsys,2; //plg, (px(2:-2))^3*(py(2:-2)/50.), (px(2:-2)), color="red"; //plg, (px(2:-2))^3*(pycorr(2:-2)), (px(2:-2)), color="red"; // plg, (py(2:-2)/50.), (pk(2:-2)); plsys,2;//plg, (py(2:-2)/50.), (px(2:-2)); } //azimuthal average plsys,1;if (plog) pltitle,"spectre de puissance"; else pltitle,"k^3^ P(k)";} S=array(0.,size); for (l=1;l<=size;l++) S(l)=((zr(l:,)-zr(:1-l,))^2)(*)(avg); if (N==1) {window,4;ws;plg, S, indgen(size); pltitle,"structure function";} } if (dim==3) { // example compute power spectrum in 3D z= genrandfield3D(size,fun=fun); //in 3D if (N==1) {window,1;ws;sec3,float(fft(z,[-1,-1,-1]))(,,::5); pltitle,"champ de densite";} if (N==1) {window,1;ws;zr=float(fft(z,[-1,-1,-1])); pli,zr(,,avg); pltitle,"champ de densite";} z2=float(fft(z*conj(z),[-1,-1])); if (N==0) {window,2;ws;pli,z2(,,avg);pltitle,"corr function";} zc=array(0.,size,size,size); for(i=1;i<=50;i++){z=genrandfield3D(size,fun=fun); z2=float(fft(z*conj(z),[-1,-1,-1]));zc+=z2;}; if (N==1) {window,2;ws;fft_pli,zc(,,avg);} zc=array(0.,size,size,size); for(i=1;i<=50;i++){z=genrandfield3D(size,fun=fun); z2=float((z*conj(z)));zc+=z2;}; u1=fft_indgen(size)*2*pi/size; r1= abs(u1, u1(-,),u1(-,-,)); pk=px*2*pi/size; py= histo2(r1, pk, weight=zc, average=1, interp=1,binsize=r1(dif)(min)); if (N==1) {window,3;ws;if (plog) { plg, log10(py(2:-2)/50), log10(pk(2:-2)), color="red"; if (is_void(fun)) plg, log10(PS(pk(2:-2))),log10(pk(2:-2)); else plg, log10(fun(pk(2:-2))),log10(pk(2:-2)); } else plg, (py(2:-2)/50), (pk(2:-2)), color="red"; pltitle,"spectre de puissance";} } if (plog) { if (is_void(S)) grow,slope,(log(py(5:numberof(pk)-5))(dif)/log(pk(5:numberof(pk)-5))(dif))(avg); else grow,slope,[[(log(py(5:numberof(pk)-5))(dif)/log(pk(5:numberof(pk)-5))(dif))(avg),(log(S(2:size/2-2))(dif)/log(indgen(size)(2:size/2-2))(dif))(avg)]]; } } // write," power spectrum slope",(log(py(5:numberof(px)-5))(dif)/log(2*pi/size.*px(5:numberof(px)-5))(dif))(avg)+2*5./3.; return slope; end; } //============================== func corr_func (delta,&px) { local c,t,r1; dims=dimsof(delta); if (dims(1)==1) { u=indgen(dims(2));r1=abs(u-u(-,)); c=delta(*)*delta(-,*); return histo2(r1,px,weight=c,average=1,interp=1); } if (dims(1)==2) { n=long(sqrt((dims(2)-1)^2+(dims(3)-1)^2))+1; t=array(0.,n+1); z=indgen(dims(2));u=indgen(dims(3)); for (i=1;i<=1;i++) { r1=(sqrt((z-z(-,))^2+(u(-,-,)-i)^2))(*); c=(delta(i,)*delta(-,,))(*); tt=histo2(r1,px,weight=c,average=1,interp=1,binmin=0,binmax=n+1,binsize=1); ttt=0*t;ttt(:numberof(px))=tt; t+=ttt; } return t; } if (dims(1)==3) { print,"this is computed via the powerspectrum"; ns=dimsof(delta)(0); x1=(indgen(ns)-1)*2*pi/float(ns); y1=(indgen(ns)-1)*2*pi/float(ns); z1=(indgen(ns)-1)*2*pi/float(ns); r1= abs(x1, y1(-,),z1(-,-,)); z=fft(delta,[1,1,1]); zc=float(z*conj(z)); py= histo2(r1, px, weight=zc, average=1, interp=1,binsize=r1(dif)(min)); return py; } } func Dgenrandfield1D (nn,&u,fun=,pp=) /* DOCUMENT pp= is the discretization of the underlying gaussian random field. u is the corresponding GRF EXAMPLE SEE ALSO: */ { local xx,yy,zz,pp,uz,u,candx,candy,candz; if(is_void(pp)) pp=64; uz=genrandfield1D(pp,fun=fun); u=double(fft(uz,-1)); u-=min(u);um=max(u); n1=0; xx=[]; n=nn/10; if(n==0) n=nn; do { candx=random(n)*pp; candy=random(n)*um; w=where(candy