common cdat,x,y,z,nx,ny,nz,nw,ntmax,date0,time0
;
restore,'data/VAR3_16' & restore,'data/eVAR3_16' & t=150.
;restore,'data/eVAR1_16' & t=50.
;restore,'data/eVAR2_16' & t=100.
s=size(bb) & nx=s[1] & ny=s[2] & nz=s[3]
x=2.*!pi*findgen(nx)/nx
y=2.*!pi*findgen(ny)/ny
z=2.*!pi*findgen(nz)/nz
;
;  Compute <A^2> in the intermediate eigenvector gauge
;
bb=reform(bb,1L*nx*ny*nz,3)
uij=reform(uij,1L*nx*ny*nz,3,3)
sij=uij
;
;  compute sij
;
for i=0,2 do begin
for j=0,2 do begin
  sij(*,i,j)=.5*(uij(*,i,j)+uij(*,j,i))
end
end
;
;  the last row, vecs(*,*,3), contains the eigenvalues
;
b=1./sqrt(dot2(bb))
vecs=eigvec3_arr(sij)
cosb1=b*dot(bb,vecs(*,*,0)) & pdf,fmin=0,fmax=1,n=500,abs(cosb1),xb1,yb1
cosb2=b*dot(bb,vecs(*,*,1)) & pdf,fmin=0,fmax=1,n=500,abs(cosb2),xb2,yb2
cosb3=b*dot(bb,vecs(*,*,2)) & pdf,fmin=0,fmax=1,n=500,abs(cosb3),xb3,yb3
save,file='align_tmp.sav',xb1,yb1,xb2,yb2,xb3,yb3
;
plot,xb3,yb3,xr=[0,1],yr=[0,1.6]
oplot,xb2,yb2,col=122
oplot,xb1,yb1,col=55
;
;  pdfs of eigenvalues
;
pdf,n=500,reform(vecs(*,0,3)),xe1,ye1
pdf,n=500,reform(vecs(*,1,3)),xe2,ye2
pdf,n=500,reform(vecs(*,2,3)),xe3,ye3
save,file='eigenvals.sav',xe1,ye1,xe2,ye2,xe3,ye3
;
bb=reform(bb,nx,ny,nz,3)
sij=reform(sij,nx,ny,nz,3,3)
vecs=reform(vecs,nx,ny,nz,3,4)
nn=float(reform(vecs(*,*,*,*,1)))
;
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)
;
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)
;
k2=dot2(kk)
Ak=complexarr(nx,ny,nz,3)
for j=0,2 do Ak(*,*,*,j)=fft(aa(*,*,*,j),-1)
nk=dot(nn,kk)
term1=dot(kk,Ak)-nk*dot(nn,Ak)
term2=k2-nk^2
bad=where(term2 eq 0.)
term2(bad)=1.
ratio=term1/term2
ratio(bad)=0.
ii=complex(0.,1.)
LLam=fltarr(nx,ny,nz,3)
for j=0,2 do LLam(*,*,*,j)=fft(-kk(*,*,*,j)*ratio,1)
bbz=fft(ii*kk(*,*,*,0)*Ak(*,*,*,1)-ii*kk(*,*,*,1)*Ak(*,*,*,0),1)
;
b2m=mean(dot2(bb))
print,mean(dot2(aa))
print,mean(dot2(aa+LLam))
;
fo='(f7.1,2e10.2,2f8.4)'
openw,1,'a2.tmp'
printf,1,t,mean(dot2(aa)),mean(dot2(aa+LLam)),mean(dot2(aa+LLam))/b2m,sqrt(b2m),fo=fo
close,1
;
spawn,'cat a2.tmp'
spawn,'cat a2.tmp >> a2.dat'
END