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<u(long(candx)+1));
if(is_array(w))
{
grow,xx,candx(w);
n1+=numberof(w);
}
}
while(n1<nn);
return xx(1:nn)/double(pp);
}
func Dgenrandfield2D (nn,&u,fun=,pp=)
/* DOCUMENT
pp= is the discretization of the underlying gaussian random field.
u is the corresponding GRF
EXAMPLE
pp=32; xy=Dgenrandfield2D(pp);
tt=histo3d(xy,indgen(pp),indgen(pp));
tt /= double(sum(tt));
pli,tt;
EXAMPLE
SEE ALSO:
*/
{
require, "Chris/linterp.i";
local xx,yy,zz,pp,uz,u,candx,candy,candz;
if(is_void(pp)) pp=64;
uz=genrandfield2D(pp,fun=fun);
u=double(fft(uz,[-1,-1]));
u-=min(u);um=max(u);
n1=0;
xx=yy=[];
n=nn/10;
do
{
candx=random(n)*pp;
candy=random(n)*pp;
candz=random(n)*um;
u1=LInterp(u,candx,candy);
w=where(candz<u1);
if(is_array(w))
{
grow,xx,candx(w);
grow,yy,candy(w);
n1+=numberof(w);
}
}
while(n1<nn);
return [xx(1:nn),yy(1:nn)]/double(pp);
}
func Dgenrandfield3D (nn,&u,fun=,pp=)
/* DOCUMENT
pp= is the discretization of the underlying gaussian random field.
u is the corresponding GRF
EXAMPLE
pp=32; xyz=Dgenrandfield3D(pp);
tt=histo3d([xx,yy,zz],indgen(pp),indgen(pp),indgen(pp));
tt /= double(sum(tt));
sec3,tt;
SEE ALSO:
*/
{
require, "Chris/linterp.i";
local xx,yy,zz,pp,uz,u;
if(is_void(pp)) pp=32;
uz=genrandfield3D(pp,fun=fun);
u=double(fft(uz,[-1,-1,-1]));
u-=min(u);um=max(u);
n1=0;
xx=yy=zz=[];
n=nn/10;
do
{
candx=random(n)*pp;
candy=random(n)*pp;
candz=random(n)*pp;
candu=random(n)*um;
u1=LInterp(u,candx,candy,candz);
w=where(candu<u1);
if(is_array(w))
{
grow,xx,candx(w);
grow,yy,candy(w);
grow,zz,candz(w);
n1+=numberof(w);
}
}
while(n1<nn);
return [xx(1:nn),yy(1:nn),zz(1:nn)]/double(pp);
}
|