FUNCTION eigvec3_arr,m
;
;  calculate eigenvector of 3x3 matrix, given in the form m(*,3,3)
;  The last index says to which eigenvalue the vector belongs, i.e.
;  ret[ipoint,*,0] is the eigenvector belonging to ret[ipoint,0,3],
;  ret[ipoint,*,1] is the eigenvector belonging to ret[ipoint,1,3], and
;  ret[ipoint,*,2] is the eigenvector belonging to ret[ipoint,2,3].
;  $Id: eigvec3_arr.pro,v 1.2 2016/12/02 00:26:48 brandenb Exp $
;
res=eigval3_arr(m)
s=size(m) & n=s(1)
j=indgen(n)
;
;  calculate eigenvectors for each eigenvalue
;
vv=complexarr(n,3,3)
for i=0,2 do begin
print,'i=',i
lambda=res(*,i)
mm=complex(m)
for k=0,2 do mm(*,k,k)=mm(*,k,k)-lambda
;
;  m00*x0 + m01*x1 + m02*x2 = 0
;  m10*x0 + m11*x1 + m12*x2 = 0
;
;  eliminate 1st row
;
; (m01*m10 - m11*m00) * x1 + (m02*m10 - m12*m00) * x2 = 0
;  choose x1=1.
;------------------------------------------------------------------------
;
;  eliminate 3rd row
;
;  (m00*m12 - m10*m02) * x0 + (m01*m12 - m11*m02) * x1 = 0
;  choose x0=1.
;
x0=replicate(1.,n)
a0=(mm(*,0,0)*mm(*,1,2)-mm(*,1,0)*mm(*,0,2))
a1=(mm(*,0,1)*mm(*,1,2)-mm(*,1,1)*mm(*,0,2))
x1=-a0*x0/a1
x2=-(mm(*,0,0)*x0+mm(*,0,1)*x1)/mm(*,0,2)
xx=[x0,x1,x2]
xx=reform(xx,n,3)
norm=sqrt(abs(xx(*,0))^2+abs(xx(*,1))^2+abs(xx(*,2))^2)
norm=rebin(reform(1./norm,n,1),n,3)
xx=xx*norm
vv(*,*,i)=xx
endfor
;
;  return 
;
ret=[reform(vv,n*9L),reform(res,n*3L)]
ret=reform(ret,n,3,4)
return,ret
END