;
; strain_eigvec.pro
;
; calculates the eigenvectors of the array st (i.e. strain tensor)
; three eigenvectors at each grid point
;
;default,lguij,1
;if lguij then begin
;  uij=fltarr(nx,ny,3,3)
;  uij[*,*,0,0]=var.guij1
;  uij[*,*,0,1]=var.guij2
;  uij[*,*,0,2]=var.guij3
;  uij[*,*,1,0]=var.guij4
;  uij[*,*,1,1]=var.guij5
;  uij[*,*,1,2]=var.guij6
;  uij[*,*,2,0]=var.guij7
;  uij[*,*,2,1]=var.guij8
;  uij[*,*,2,2]=var.guij9
;endif
;
antisym=fltarr(nx,ny,3,3)
strain=fltarr(nx,ny,3,3)
for j=0,2 do begin
for i=0,2 do begin
  strain[*,*,i,j]=.5*(uij[*,*,i,j]+uij[*,*,j,i])
  antisym[*,*,i,j]=.5*(uij[*,*,i,j]-uij[*,*,j,i])
endfor
endfor
;
s=size(strain) & nx=s(1) & ny=s(2) & nz=s(3)
nz=1
;
help,nx,ny,nz
tmp=reform(strain,nx*ny*nz,3,3)
;help,tmp
;
tmpp=eigvec3_arr(tmp)  
;help,tmpp
;
eigvec_array=reform(tmpp,nx,ny,nz,3,4)
eigvec_array=reform(eigvec_array)
print,'eigvec_array'
;
oee1=dot(o1,eigvec_array[*,*,*,0])
oee2=dot(o1,eigvec_array[*,*,*,1])
oee3=dot(o1,eigvec_array[*,*,*,2])
;
;  norms
;
s2=total(total(strain ^2,4),3)
a2=total(total(antisym^2,4),3)
;
;  flux
;
oxu=cross(var.oo,var.uu)
FS=fltarr(nx,ny,3)
for j=0,2 do begin
  FS[*,*,j]=FS[*,*,j]-dot(var.uu,strain[*,*,*,j])
endfor
;
END