;$Id: pvar2.pro,v 1.2 2014/11/27 11:23:58 brandenb Exp $
;
;  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)))
END