;$Id: phel2.pro,v 1.16 2016/07/26 15:04:13 brandenb Exp $
;
;  This routine reads and works with the FFT result obtained with phel12.
;  This is the version *with* domain doubling in the z direction.
;
restore,'data/bbk.sav'
print,'all read'
;
;  toggle between KK=0 and 1
;
for KK=0,1 do begin
;
b1kp=shift(b1k,0,0,-KK)
b2kp=shift(b2k,0,0,-KK)
b3kp=shift(b3k,0,0,-KK)
print,'all shifted'
;
imax=min([nx,ny,nz]/2)
x1=indgen(nx)-nx/2
y1=indgen(ny)-ny/2
z1=indgen(nz)-nz/2
x1=reform(x1,nx,1,1)
y1=reform(y1,1,ny,1)
z1=reform(z1,1,1,nz)
x1=rebin(x1,nx,ny,nz)
y1=rebin(y1,nx,ny,nz)
z1=rebin(z1,nx,ny,nz)
rad=sqrt(x1^2+y1^2+z1^2)
;
;  unit vector
;
rad1=rad
orig=where(rad1 eq 0.)
rad1(orig)=1.
rad1=1./rad1
rad1(orig)=0.
;
x1=x1*rad1
y1=y1*rad1
z1=z1*rad1
print,'rad & rad1 computed'
;
tEk=fltarr(imax)
tEk1=fltarr(imax)
tEk2=fltarr(imax)
tEk3=fltarr(imax)
kHk=fltarr(imax)
kHk1=fltarr(imax)
kHk2=fltarr(imax)
kHk3=fltarr(imax)
iii=complex(0.,1.)
for i=0,imax-1 do begin
  ii=where((rad ge i-.5) and (rad lt i+.5))
  ;
  tmp1=float(b1kp(ii)*conj(b1k(ii)))
  tmp2=float(b2kp(ii)*conj(b2k(ii)))
  tmp3=float(b3kp(ii)*conj(b3k(ii)))
  tEk1(i)=total(tmp1)
  tEk2(i)=total(tmp2)
  tEk3(i)=total(tmp3)
  tEk(i)=total(tmp1+tmp2+tmp3)
  ;
  tmp1=2*float(-iii*b1kp(ii)*conj(b2k(ii))*z1(ii))
  tmp2=2*float(-iii*b2kp(ii)*conj(b3k(ii))*x1(ii))
  tmp3=2*float(-iii*b3kp(ii)*conj(b1k(ii))*y1(ii))
  kHk1(i)=total(tmp1)
  kHk2(i)=total(tmp2)
  kHk3(i)=total(tmp3)
  kHk(i)=total(tmp1+tmp2+tmp3)
  ;
  print,i
endfor
print,'kHk computed'
;
k=indgen(imax)
plot_oi,k(1:*),kHk(1:*)
plot,k,kHk,xr=[0,20]
;
k1=k
k1[0]=1.
k1=1./k1
k1[0]=0.
;
wfile='single_2Ek_and_kHk'+str(KK)+'.sav'
save,file=wfile,t,k,k1,tEk,tEk1,tEk2,tEk3,kHk,kHk1,kHk2,kHk3
print,'saved to file '+wfile
;
endfor
;
end