;$Id: heldens_normtest.pro,v 1.3 2006/03/25 10:56:59 brandenb Exp $
;
;  calculate B-field from vector potential,
;  exclude ghost zones, and remove horizontal average (hfluctv)
;
bbb_total=pc_noghost(curl(aa))
bbb=hfluctv(bbb_total)
;
;  determine mesh size of the current array
;
s=size(bbb) & nx=s[1] & ny=s[2] & nz=s[3]
if s[0] eq 3 then nvar=1 else nvar=s[4]
;
;  calculate mesh spacing in x and k space
;
dx=x(1)-x(0) & dkx=2.*!pi/(nx*dx)
dy=y(1)-y(0) & dky=2.*!pi/(ny*dy)
dz=z(1)-z(0) & dkz=2.*!pi/(nz*dz)
;
;  calculate components of wave vector
;
kk=fltarr(nx,ny,nz,3)
kk(*,*,*,0)=shift(rebin(reform(dkx*(findgen(nx)-(nx-1)/2),nx,1,1),nx,ny,nz),-(nx-1)/2,-(ny-1)/2,-(nz-1)/2)
kk(*,*,*,1)=shift(rebin(reform(dky*(findgen(ny)-(ny-1)/2),1,ny,1),nx,ny,nz),-(nx-1)/2,-(ny-1)/2,-(nz-1)/2)
kk(*,*,*,2)=shift(rebin(reform(dkz*(findgen(nz)-(nz-1)/2),1,1,nz),nx,ny,nz),-(nx-1)/2,-(ny-1)/2,-(nz-1)/2)
;
;  calculate 1/k^2
;
;k2=dot2(kk)
;k2(0,0,0)=1.
;k21=1./k2 & k21(0,0,0)=0.
;
;  calculate b(k) and b*(k), denoted here by hk.
;
bk=complexarr(nx,ny,nz,3)
for j=0,2 do bk(*,*,*,j)=fft(bbb(*,*,*,j),-1)
hk=conj(bk)
;
;  calculate b(k+1/2), etc
;  p/m stands for +/- and h stands for 1/2
;
kkp1=shift(kk,0,0,-1,0)
bkp1=shift(bk,0,0,-1,0)
hkp1=shift(hk,0,0,-1,0)
kkph=.5*(kkp1+kk)
bkph=.5*(bkp1+bk)
hkph=.5*(hkp1+hk)
;
;  calculate b(k-1/2), etc
;
kkm1=shift(kk,0,0,+1,0)
bkm1=shift(bk,0,0,+1,0)
hkm1=shift(hk,0,0,+1,0)
kkmh=.5*(kkm1+kk)
bkmh=.5*(bkm1+bk)
hkmh=.5*(hkm1+hk)
;
;  calculate b^2
;
ii=complex(0.,1.)
mag0=dot(bk,hk)
mag1=dot(bkph,hkmh)
mag2=dot(bkp1,hkm1)
mag0m=total(mag0)
mag1m=total(mag1)
mag2m=total(mag2)
print,mag0m,mag1m,mag2m
;
;  do the same for Coulomb gauged expression
;
;helC=dot(bkph,cross(hkmh,ii*kkph))/dot2(kkph)
;helCm=total(helC)
;print,helCm
;
;  calculate exp(iKR) factor, for K=(0,0,1)
;
zzz=z(n1:n2)
;expiKR=exp(ii*zzz)
exp0iKR=exp(0.*ii*zzz)
exp1iKR=exp(1.*ii*zzz)
exp2iKR=exp(2.*ii*zzz)
plot,zzz,float(mag0m*exp0iKR+2.*mag1m*exp1iKR+2.*mag2m*exp2iKR)
oplot,zzz,float(mag0m*exp0iKR+2.*mag2m*exp2iKR),li=2
;
;  save coefficients for post-processing
;
save,file=varfile+'.magm',zzz,t,exp0iKR,exp1iKR,exp2iKR,mag0m,mag1m,mag2m
;
end