(* Zachary H. Levine    22 March 2001 - 4 June 2001
   zlevine@nist.gov


   This code is a product of the National Institute of Standards
   and Technology.  As a courtesy, please acknowledge if used in
   a publication or report and inform the author.  Thank you.


   This routine calculates the eigenvalues and eigenvectors
   of the correction to the dielectric matrix the anisotropic
   term in a cubic crystal.   The interaction is
                 eps_ij = T_ijkl q_k q_l interaction
   where eps is the dielectric matrix, T is a fourth rank tensor
   The eigenvalues and eigenvectors are calculated in the subspace
   orthogonal to the direction of propagation (as permitted by
   the Maxwell equations).


   See JH Burnett, ZH Levine, and EL Shirley, unpublished
   for more details; T here is alpha in the reference and tt in
   the code. (Expect publication second half 2001.)


   This routine is written for Mathematica 4.0.  NIST does
   not recommend or endorse the use of particular commercial
   products; commercial products identified here are not necessarily
   the best available for the purpose.  The first part is translated
   into Fortran also, called aneig.f.


   To adapt to other uses, make sure T is defined by setAniso.
   (This will occur if the program is loaded with <<anEig.m.)


   The function eigSys takes as input two polar co-ordinates (theta
   and phi) which describe the direction of q.  The output is a
   Mathematica expression of the type returned by Mathematica's
   Eigensystem, namely the first component holds the two eigenvalues
   and the second the two eigenvectors.  Each eigenvector is given
   in Cartesian co-ordinates.  They are sorted into families.
*)


(* The independent tensor component which exists in a cubic system
   but not the isotropic system is made below.  It is scaled so that
   the eigenvalue splitting for q=(1,1,0)/Sqrt(2) is
               eps_ij E_i E_j - eps_ij E'_i E'_j = 1
   where E =(1,-1,0)/Sqrt(2) and E'=(0,0,1).
*)


setAniso[] = 
   Block[{d,T},
d = 3;
T = Table [0, {i,d}, {j,d}, {k,d}, {l,d}];
T[[1,1,1,1]] =  2 ; T[[2,2,2,2]] =  2 ; T[[3,3,3,3]] =  2;
T[[1,1,2,2]] = -1 ; T[[1,1,3,3]] = -1 ; T[[2,2,3,3]] = -1;
T[[2,2,1,1]] = -1 ; T[[3,3,1,1]] = -1 ; T[[3,3,2,2]] = -1;
T[[2,3,2,3]] = -1 ; T[[2,3,3,2]] = -1 ; T[[3,2,2,3]] = -1 ; T[[3,2,3,2]] = -1;
T[[1,3,1,3]] = -1 ; T[[1,3,3,1]] = -1 ; T[[3,1,1,3]] = -1 ; T[[3,1,3,1]] = -1;
T[[1,2,1,2]] = -1 ; T[[1,2,2,1]] = -1 ; T[[2,1,1,2]] = -1 ; T[[2,1,2,1]] = -1;
Return[T*0.4]
]
T = setAniso[]



(* ----------------------------------------- end T set up *)


(* In eigSys, the input variables are t for "theta" and p for
   "phi" which have the usual definition from spherical 
   co-ordinates.  The direction of propagation q-hat is in this
   direction.  It is also called "r" because that is the name
   in spherical co-ordinates.  "th" and "ph" are the other unit
   vectors theta-hat and phi-hat of spherical co-ordinates.
   Trr is the contraction of T with q-hat (twice).


   proj is a 2D subspace orthogonal to q.  The polarization
   direction is assumed to lie in this plane.  The matrix Trr,
   is projected into this subspace.  There, it is diagonalized
   to give the eigenvalues and eigenvectors.  The eigenvectors
   are projected back into the 3D space.  They are normalized
   to unit vectors and sorted so that the directions are continuous
   functions of the input (except at directions q which have of
   eigenvalue degeneracy).
*)
eigSys[theta_,phi_]:=
    Block[{ct,st,cp,sp,th,ph,r,Trr,proj,projTr},
        ct = Cos[theta];
        st = Sin[theta];
        cp = Cos[phi];
        sp = Sin[phi];
        th = {  ct cp,  ct sp, -st };
        ph = { -   sp,     cp,   0 };
        r  = {  st cp,  st sp,  ct };
        Trr = (T.r).r;
        proj = {th,ph};
        projTr = Transpose[proj];
        TProj = proj . Trr . projTr;
        {eigVal, eigVec} = Eigensystem[TProj];
        eigVec3D = eigVec . proj;                (* unproject eigenvector *)
        eigVec3D = {eigVec3D[[1]]/Sqrt[eigVec3D[[1]].eigVec3D[[1]]],
                    eigVec3D[[2]]/Sqrt[eigVec3D[[2]].eigVec3D[[2]]]};
        dir = maxIndex[Abs[r]];                  (* sort eigensystem by dir *)
        If [ Abs[eigVec3D[[1,dir]]] < Abs[eigVec3D[[2,dir]]], 
                eigVec3D = Reverse[eigVec3D];
                eigVal   = Reverse[eigVal] ];
        Return[ {eigVal, eigVec3D} ]
]


(*  Returns the index associated with the maximum vector.
*)


maxIndex[vec_] := Block [{x,y,z,m},
                          x = vec[[1]];
                          y = vec[[2]];
                          z = vec[[3]];
                          m = 1;
                          If [y>x, m=2];
                          If [z>y && z>x, m=3];
                          Return[m] ]


(* The part below is not translated into Fortran.
   It makes a nice plot of the directions of the eigenvectors.
*)


ev[theta_,phi_,i_] := Block[ {evalue,evectr},
                             {evalue,evectr}= eigSys[theta,phi];
                              Return[ evectr[[i]] ] ]


<<Graphics`PlotField3D`  (* From Mathematica 4.0 standard add-on package *)


scale = 0.05


phi0 = -N[3Pi/4]
phi1 =  N[ Pi/4]
dphi =  N[Pi/40]
dth  =  N[Pi/40]
th0  =  dth
th1  =  N[Pi/2]
evTable1 = Flatten [
        Table [ {{Sin[th] Cos[ph], Sin[th] Sin[ph], Cos[th]},
                 ev[th,ph,1]*scale},
                    {ph,phi0,phi1,dphi}, {th,th0,th1,dth} ], 
            1] // Chop
ListPlotVectorField3D[evTable1]
evTable2 = Flatten [
        Table [ {{Sin[th] Cos[ph], Sin[th] Sin[ph], Cos[th]},
                 ev[th,ph,2]*scale},
                    {ph,phi0,phi1,dphi}, {th,th0,th1,dth} ], 
            1] // Chop
ListPlotVectorField3D[evTable2]