c Zachary H. Levine 22 March 2001 - 4 June 2001 c zlevine@nist.gov c This code is in the public domain. Do not remove comments. c This code is a product of the National Institute of Standards c and Technology. As a courtesy, please acknowledge if used in c a publication or report and inform the author. Thank you. c This routine calculates the eigenvalues and eigenvectors c of the correction to the dielectric matrix the anisotropic c term in a cubic crystal. The interaction is c eps_ij = T_ijkl q_k q_l interaction c where eps is the dielectric matrix, T is a fourth rank tensor c See JH Burnett, ZH Levine, and EL Shirley, unpublished c for more details; T here is alpha in the reference and tt in c the code. (Expect publication second half 2001.) c This subroutine is a translation of part of anEig.m c a Mathematica routine. NIST does not recommend or endorse c the use of particular commercial products; commercial c products identified here are not necessarily the best available c for the purpose. c In eigSys, the input variables are t for "theta" and p for c "phi" which have the usual definition from spherical c co-ordinates. The direction of propagation q-hat is in this c direction. It is also called "r" because that is the name c in spherical co-ordinates. "th" and "ph" are the other unit c vectors theta-hat and phi-hat of spherical co-ordinates. c Trr is the contraction of T with q-hat (twice). c proj is a 2D subspace orthogonal to q. The polarization c direction is assumed to lie in this plane. The matrix Trr, c is projected into this subspace. There, it is diagonalized c to give the eigenvalues and eigenvectors. The eigenvectors c are projected back into the 3D space. They are normalized c to unit vectors and sorted so that the directions are continuous c functions of the input (except at directions q which have of c eigenvalue degeneracy); also, there may be discontinuities c of an overall sign. c The user (optionally) may use ttSet to obtain the tensor T. c A companion test program, aneigt.f, is available. subroutine aneig (theta, phi, eigVal, eigVec) implicit none double precision theta, phi, eigVal(2), eigVec(3,2) c Input c theta - spherical co-ordinate theta value for direction of c propagation q (in radians) c phi - spherical co-ordinate phi value for direction of c propagation q (in radians) c Output c eigVal - 2 eigenvalues of polarization (sorted by family) c eigVec - eigenvector of polarization c index1: xyz component index2: eigenvector index c Internal c proj - index1: xyz component, index2: thetaHat,phiHat c eig2D - index1: thetaHat,phiHat component, index2: eigenvector index c tt - fourth rank anisotropic tensor from ttSet double precision tiny parameter (tiny=1.d-10) integer dir, i, j, k, l double precision tt(3,3,3,3), ct, st, cp, sp, $ eig2D(3,2), Trr(3,3), proj(3,2), $ r(3), th(3), ph(3), TProj(2,2), $ temp(3,2), temp2, aa, bb, dd, disc, norm logical first data first/.true./ save first, tt if (first) then call ttSet(tt) first = .false. endif ct = Cos(theta) st = Sin(theta) cp = Cos(phi) sp = Sin(phi) th(1) = ct * cp th(2) = ct * sp th(3) = -st ph(1) = -sp ph(2) = cp ph(3) = 0 r(1) = st * cp r(2) = st * sp r(3) = ct do 1 i = 1,3 do 1 j = 1,3 Trr(i,j) = 0.d0 do 1 k = 1,3 do 1 l = 1,3 1 Trr(i,j) = Trr(i,j) + tt(i,j,k,l) * r(k) * r(l) proj(1,1) = th(1) proj(2,1) = th(2) proj(3,1) = th(3) proj(1,2) = ph(1) proj(2,2) = ph(2) proj(3,2) = ph(3) do 2 i = 1,3 do 2 j = 1,2 temp(i,j) = 0.d0 do 2 k = 1,3 2 temp(i,j) = temp(i,j) + Trr(i,k) * proj(k,j) do 3 i = 1,2 do 3 j = 1,2 TProj(i,j) = 0.d0 do 3 k = 1,3 3 TProj(i,j) = TProj(i,j) + proj(k,i) * temp(k,j) aa = TProj(1,1) bb = TProj(1,2) + TProj(2,1) dd = TProj(2,2) if (abs(bb).gt.tiny) then disc = sqrt( (aa-dd)**2 + bb**2 ) eigVal(1) = 0.5d0*(aa + dd - disc) eigVal(2) = 0.5d0*(aa + dd + disc) eig2D(1,1) = aa-dd-disc eig2D(2,1) = bb eig2D(1,2) = aa-dd+disc eig2D(2,2) = bb else eigVal(1) = aa eigVal(2) = dd eig2D(1,1) = 1.d0 eig2D(2,1) = 0.d0 eig2D(1,2) = 0.d0 eig2D(2,2) = 1.d0 endif do 4 i = 1,3 do 4 j = 1,2 eigVec(i,j) = 0 do 4 k = 1,2 4 eigVec(i,j) = eigVec(i,j) + proj(i,k)*eig2D(k,j) do 5 i = 1,2 norm = sqrt(eigVec(1,i)**2 + eigVec(2,i)**2 + eigVec(3,i)**2) eigVec(1,i) = eigVec(1,i)/norm eigVec(2,i) = eigVec(2,i)/norm 5 eigVec(3,i) = eigVec(3,i)/norm if (abs(r(1)) .ge. abs(r(2)) .and. abs(r(1)) .ge. abs(r(3))) then dir = 1 else if (abs(r(2)) .ge. abs(r(3))) then dir = 2 else dir = 3 endif if (abs(eigVec(dir,1)) .lt. abs(eigVec(dir,2))) then temp2 = eigVal(1) eigVal(1) = eigVal(2) eigVal(2) = temp2 do 6 i = 1, 3 temp2 = eigVec(i,1) eigVec(i,1) = eigVec(i,2) 6 eigVec(i,2) = temp2 endif return end subroutine ttSet(tt) implicit none c Set the independent tensor component which exists in a cubic c system but not the isotropic system is made below with the c eigenvalue splitting along 110 scaled to eps_(1-10)-eps(001)=1 c Output double precision tt(3,3,3,3) integer i, j, k, l do 1 i = 1,3 do 1 j = 1,3 do 1 k = 1,3 do 1 l = 1,3 1 tt(i,j,k,l) = 0.d0 tt(1,1,1,1) = 0.8d0 tt(2,2,2,2) = 0.8d0 tt(3,3,3,3) = 0.8d0 tt(1,1,2,2) = -0.4d0 tt(1,1,3,3) = -0.4d0 tt(2,2,3,3) = -0.4d0 tt(2,2,1,1) = -0.4d0 tt(3,3,1,1) = -0.4d0 tt(3,3,2,2) = -0.4d0 tt(2,3,2,3) = -0.4d0 tt(2,3,3,2) = -0.4d0 tt(3,2,2,3) = -0.4d0 tt(3,2,3,2) = -0.4d0 tt(1,3,1,3) = -0.4d0 tt(1,3,3,1) = -0.4d0 tt(3,1,1,3) = -0.4d0 tt(3,1,3,1) = -0.4d0 tt(1,2,1,2) = -0.4d0 tt(1,2,2,1) = -0.4d0 tt(2,1,1,2) = -0.4d0 tt(2,1,2,1) = -0.4d0 return end