@parameters
;
default,ivar,-1
default,iread,0
if iread eq 0 then begin
  pc_read_param,/param2,obj=param
  if iread eq 0 then begin
    if ivar ge 0 then begin
      pc_read_var,/magic,/trimall,obj=var,variables=['uu','bb','uij','aa'],ivar=ivar
    endif else begin
      pc_read_var,/magic,/trimall,obj=var,variables=['uu','bb','uij','aa']
    endelse
    iread=1
  endif
  xxx=var.x
  yyy=var.y
  zzz=var.z
  iread=1
endif
nz=n_elements(zzz)
;
tvscl,var.bb(*,*,nz-1,2)
;
;  Procedure to compute rate of strain tensor,
;  as well as its eigenvalues and eigenvectors.
;
;  dimensions
;
s=size(var.uij)
nx=s[1]
ny=s[2]
nz=s[3]
;
;aa=reform(var.aa,1L*nx*ny*nz,3)
bb=reform(var.bb,1L*nx*ny*nz,3)
uij=reform(var.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
;
;stop
b=1./sqrt(dot2(var.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.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)))
;
x=var.x
y=var.y
z=var.z
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(var.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(var.bb))
print,mean(dot2(var.aa))
print,mean(dot2(var.aa+LLam))
;
fo='(f7.1,2e10.2,2f8.4)'
openw,1,'a2.tmp'
printf,1,var.t,mean(dot2(var.aa)),mean(dot2(var.aa+LLam)),mean(dot2(var.aa+LLam))/b2m,sqrt(b2m),fo=fo
close,1
;
spawn,'cat a2.tmp'
spawn,'cat a2.tmp >> a2.dat'
END