C
```c
c........|.........|.........|.........|.........|.........|.........|..
c........|.........|.........|.........|.........|.........|.........|..
c
c Rmatrix and associated programs may be used freely by University
c researchers and teachers. Commercial parties interested in using
c the program should contact C. J. Thomson.
c
c........|.........|.........|.........|.........|.........|.........|..
c........|.........|.........|.........|.........|.........|.........|..
c
c
subroutine aniso6c1(a,rho,px,py,eval,evec)
implicit real*8 (a-h,o-z)

c Copyright (c) 1996 C. J. Thomson.

c Finds e-vals and e-vecs for given layer.
c Based on Sean Guest's routines for finding vertical slowness given
c horizontal slownesses by solving the 6X6 matrix problem (using
c public domain library routines - Eispack? They look like it.).
c
dimension a(3,3,3,3),p(3)
complex*16 px,py
complex*16 eval(6),evec(6,6),xnorm
logical iso
common/info/lay,id,iso

iso=.false.

p(1)=dble(px)
p(2)=dble(py)
call eig (rho,a,p,eval,evec)

c Normalise these evecs to unit displacement

do 100 j=1,6
xnorm=cdsqrt(evec(1,j)**2+evec(2,j)**2+evec(3,j)**2)
do 101 i=1,6
evec(i,j)=evec(i,j)/xnorm
101   continue
100 continue

return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine eig (den,ae,p,wc,zc)

c Copyright (c) 1996 W. S. Guest and C. J. Thomson.

c
c          W.S. Guest  October 1988
c
c This subroutine finds the eigenvalues (vertical slowness) and
c eigenvectors (displacement-stress) when given the 81 component
c elastic constants AE, the density DEN, and the horizontal slownesses
c p(1) and p(2). If the medium is isotropic, then a call is made
c to ISOEIG which provide a rapid analytic solution (Fryer & Frazer).
c THe anisotropic formulation is after Woodhouse (1974) and then the
c system matrix is reduced using the `numerical recipes' routine DEIGV.
c The vectors must then be sorted into
c                 Pu,SVu,SHu,Pd,SVd,SHd
c This is done in SORT which bases the S orientation on polarization
c or in SORT_AMP which sorts the eigenvalues from smallest to largest.
c
implicit real*8(a-h,o-z)
c
integer i,j,k,l,m,n,lay,id
logical iso
dimension c(3,3,3,3),ae(3,3,3,3),cm(3,3),ci(3,3),t(3,3),tt(3,3)
dimension ct2(3,3,3),ct(3,3),qc(3,3,3,3),s(3,3),a(6,6),p(3)
dimension zi(6,6),wr(6),wi(6),zr(6,6)
complex*16 wc(6),zc(6,6)
c Common Block
common/info/lay,id,iso
c
c
c     write (*,*) 'ineig',iso
do 10 i=1,3
do 20 j=1,3
do 30 k=1,3
do 40 l=1,3
c(i,j,k,l) = ae(i,j,k,l) * den
40      continue
30     continue
20    continue
10   continue
c
if (iso) then
c        write (*,*)'IsoEig used to solve "A"'
c33=c(3,3,3,3)
c44=c(2,3,2,3)
c        write (*,*) c33,c44
call isoeig(den,p,c33,c44,wc,zc)
goto 999
endif
c
c Calculate the c33 inverse matrix
c
cm(1,1)=c(2,3,2,3)*c(3,3,3,3) - c(2,3,3,3)**2.d0
cm(2,2)=c(1,3,1,3)*c(3,3,3,3) - c(1,3,3,3)**2.d0
cm(3,3)=c(1,3,1,3)*c(2,3,2,3) - c(1,3,2,3)**2.d0
cm(1,2)=-1.d0*(c(2,3,1,3)*c(3,3,3,3) - c(2,3,3,3)*c(3,3,1,3))
cm(1,3)=c(2,3,1,3)*c(3,3,2,3) - c(2,3,2,3)*c(3,3,1,3)
cm(2,3)=-1.d0*(c(1,3,1,3)*c(3,3,2,3) - c(1,3,2,3)*c(3,3,1,3))
cm(2,1)=cm(1,2)
cm(3,1)=cm(1,3)
cm(3,2)=cm(2,3)
det=c(1,3,1,3)*cm(1,1) + c(2,3,1,3)*cm(2,1) + c(3,3,1,3)*cm(3,1)
c
do 50 i=1,3
do 55 j=1,3
ci(i,j)=cm(i,j)/det
55     continue
50   continue
c
c Calculate the 'T' and 'T' transpose matrices
c
do 60 m=1,3
do 65 n=1,3
c1=0.d0
c2=0.d0
do 70 k=1,3
c1 = c1 + c(m,1,k,3)*ci(k,n)
c2 = c2 + c(m,2,k,3)*ci(k,n)
70       continue
tt(n,m) = -1.d0*(p(1)*c1 + p(2)*c2)
t(m,n) = tt(n,m)
65     continue
60   continue
c
c Calculate the 'Q' matrix
c
do 75 i=1,2
do 80 k=1,3
do 85 l=1,3
ct(k,l)=0.d0
do 90 m=1,3
ct(k,l) = ct(k,l) + c(k,i,m,3)*ci(m,l)
90         continue
85       continue
80     continue
c
do 95 j=1,2
do 100 k=1,3
do 105 l=1,3
ct2(j,k,l)=0.d0
do 110 m=1,3
ct2(j,k,l)= ct2(j,k,l) + ct(k,m)*c(m,3,l,j)
110          continue
qc(i,j,k,l)=c(k,i,l,j) - ct2(j,k,l)
105        continue
100      continue
95     continue
75   continue
c
c Calculate the 'S' matrix
c
do 115 k=1,3
do 120 l=1,3
s(k,l)=0.d0
do 125 i=1,2
do 130 j=1,2
s(k,l) = s(k,l) + p(i)*p(j)*qc(i,j,k,l)
130        continue
125      continue
s(k,l) = -1.d0 *s(k,l)
if (k.eq.l) s(k,l)=s(k,l)+den
120    continue
115  continue
c
c Form the complete system matrix 'A'
c
do 135 i=1,3
j=i+3
do 140 k=1,3
l=k+3
a(i,k) =  tt(i,k)
a(i,l) =  ci(i,k)
a(j,k) =  s(i,k)
a(j,l) =  t(i,k)
140    continue
135  continue
c
c Find the eigenvalues and eigenvectors
c
call deigv(6,a,wr,wi,zr,zi)
c
call sort_pol(wr,wi,zr,zi,wc,zc,ae,p)
c     call sort_amp(wr,wi,zr,zi,wc,zc,ae,p)
c
999  continue
c
c      write (19,*)
c      write (19,'(6e12.3)') wc(1),wc(2),wc(3),wc(4),wc(5),wc(6)
c      write (19,*)
c      write (19,'(12d10.2)')  zc(1,1),zc(1,2),zc(1,3),
c    *                       zc(1,4),zc(1,5),zc(1,6)
c      write (19,'(12d10.2)')  zc(2,1),zc(2,2),zc(2,3),
c    *                       zc(2,4),zc(2,5),zc(2,6)
c      write (19,'(12d10.2)')  zc(3,1),zc(3,2),zc(3,3),
c    *                       zc(3,4),zc(3,5),zc(3,6)
c
return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine sort_amp (wr,wi,zr,zi,wc,zc,ae,p)

c Copyright (c) 1996 W. S. Guest and C. J. Thomson.

c
c same as SORT except that the shear eigenvalue separation is based
c soley on magnitudes. If all are real then the eigenvalues are arranged
c as 1, the smallest magnitude, to 3, the largest. This means that 3 will
c be the last to go evanescant. When all three waves are imaginary then the
c magnitude of the imaginary components are monitored and they are arranged
c from 1, the largest, to 3, the smallest. If there are no intersections then
c this will guarentee smooth flowing curves.
c
implicit real*8 (a-h,o-z)
c
integer np(6),neg(3),pos(3),neg_ev(3),pos_ev(3),in_ev,ip_ev,
+        nvmax,pvmax,nvmin,pvmin,in,ip,i,j,k
dimension wr(6),wi(6),zr(6,6),zi(6,6),p(3),v(3),ae(3,3,3,3)
complex*16 wc(6),zc(6,6)
c Common Blocks
common/info/lay,id,iso
c
c
numreal=0
in=0
ip=0
in_ev=0
ip_ev=0
velmaxn=0.d0
velmaxp=0.d0
velminn=10000.d0
velminp=10000.d0
holdn=-10.d0
holdp=10.d0
c
do i=1,6
np(i)=0
enddo
c
c      write (19,'(6e10.3)') wr
c      write (19,'(6e10.3)') wi
c
c
c
do i=1,6
if (dabs(wi(i)) .lt. 1.d-9)then
numreal=numreal+1
p(3)=wr(i)
call vtot(ae,p,v)
velocity=dsqrt(v(1)**2. + v(2)**2. + v(3)**2.)
c          write (19,'(i4,4f10.1)')i,v,velocity
c
if (v(3) .lt. 0.d0)then
in=in+1
neg(in)=i
if (velocity .lt. velminn)then
velminn=velocity
nvmin=i
endif
if (velocity .gt. velmaxn)then
velmaxn=velocity
nvmax=i
endif
else
ip=ip+1
pos(ip)=i
if (velocity .lt. velminp)then
velminp=velocity
pvmin=i
endif
if (velocity .gt. velmaxp)then
velmaxp=velocity
pvmax=i
endif
endif
elseif (wi(i) .lt. 0.d0)then
in=in+1
in_ev=in_ev+1
neg_ev(in_ev)=i
neg(in)=i
else
ip=ip+1
ip_ev=ip_ev+1
pos_ev(ip_ev)=i
pos(ip)=i
endif
enddo
c
c      write (19,*) ip,in,pos,neg
c      write (19,*) 'number of real eigenvalues',numreal
c      write (19,*) pos_ev,neg_ev
c      write (19,*) 'pvmax,pvmin',pvmax,pvmin
c      write (19,*) 'nvmax,nvmin',nvmax,nvmin
c
j=0
k=0
c
if (numreal .ge. 4)then
c
np(3)=pvmin
np(6)=nvmin
c
if (numreal .eq. 6)then
np(1)=pvmax
np(4)=nvmax
else
np(1)=pos_ev(1)
np(4)=neg_ev(1)
endif
c
do i=1,3
if (pos(i).ne.np(1) .and. pos(i).ne.np(3)) np(2)=pos(i)
if (neg(i).ne.np(4) .and. neg(i).ne.np(6)) np(5)=neg(i)
enddo
c
elseif (numreal .ge. 2)then
c
np(3)=pvmin
np(6)=nvmin
c
if (wi(pos_ev(1)) .gt. wi(pos_ev(2)))then
np(1)=pos_ev(1)
np(2)=pos_ev(2)
else
np(2)=pos_ev(1)
np(1)=pos_ev(2)
endif
if (wi(neg_ev(1)) .lt. wi(neg_ev(2)))then
np(4)=neg_ev(1)
np(5)=neg_ev(2)
else
np(5)=neg_ev(1)
np(4)=neg_ev(2)
endif
c
else
c
do 200 i=1,3
if (wi(pos_ev(1)) .gt. wi(pos_ev(2))) then
j=j+1
if (j .eq. 1)then
holdp=wi(pos(i))
np(2)=pos(i)
elseif (wi(pos(i)) .lt. holdp) then
np(3)=np(2)
np(2)=pos(i)
else
np(3)= pos(i)
endif
endif
c
if (neg(i) .ne. np(4)) then
k=k+1
if (j .eq. 1)then
holdn=wi(neg(i))
np(5)=neg(i)
elseif (wi(neg(i)) .gt. holdn) then
np(6)=np(5)
np(5)=neg(i)
else
np(6)= neg(i)
endif
endif
200      continue
c
endif
c
c     write (19,*) 'np',np
c
c Switch the eigenvalues and eigenvectors to the order P,S1,S2
c
do i=1,6
wc(i)=dcmplx(wr(np(i)),wi(np(i)))
do j=1,3
zc(j,i)=dcmplx(zr(j,np(i)),zi(j,np(i)))
c         zc(j+3,i)=dcmplx(zr(j+3,np(i)),zi(j+3,np(i)))*(0.d0,1.d0)
zc(j+3,i)=dcmplx(zr(j+3,np(i)),zi(j+3,np(i)))
enddo
enddo
c
return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine sort_pol (wr,wi,zr,zi,wc,zc,ae,p)

c Copyright (c) 1996 W. S. Guest and C. J. Thomson.

c
c              W.S. Guest   July 1989
c
c This subroutine takes the input of the eigenvalues and eigenvectors
c calculated in eig and orders them into pu,svu,shu,pd,svd,shd. a
c check is first made to see if any of the eigenvalues are complex
c representing evanescant waves. separation of the up and downgoing
c waves is based on th sign of the eigenvalue but a check is made to
c make sure there are three of each sign, if not then the total vel.
c are calculated to get energy flux direction.
c if all of the eigenvalues are real,
c then the p-wave is taken as the smallest eigenvalues, both in the
c up and downgoing direction. the s-waves are not sorted, but
c always come out as sv and then sh. s-waves can be separated based
c on displacement, but problem occur for vertically travelling waves
c
implicit real*8 (a-h,o-z)
c
integer i,j,k,np(6),neg(3),npos(3)
dimension wr(6),wi(6),zr(6,6),zi(6,6),
+          p(3),v(3),ae(3,3,3,3)
complex*16 wc(6),zc(6,6)
c Common Blocks
common/info/lay,id,iso
c Parameter defintions
parameter (pi = 3.1415926535897932d0)
c
j=0
k=0
holdn=-10.d0
holdp=10.d0
numev=0
hldevn=0.d0
hldevp=0.d0
do 5 i=1,6
np(i)=0
5    continue
c
c Check to see is any of the waves are evanesant
c
do 10 i=1,6
if (wi(i) .ne. 0.d0)then
if (wi(i) .lt. 0.d0)then
j=j+1
neg(j)=i
numev=numev+1
if (wi(i) .lt. hldevn) then
np(4)=i
hldevn=wi(i)
endif
else
numev=numev+1
k=k+1
npos(k)=i
if (wi(i) .gt. hldevp) then
np(1)=i
hldevp=wi(i)
endif
endif
else
if (wr(i) .lt. 0.d0)then
j=j+1
neg(j)=i
if (numev .eq. 0) then
if (wr(i) .gt. holdn)then
holdn=wr(i)
np(4)=i
endif
endif
else
k=k+1
npos(k)=i
if (numev .eq. 0) then
if (wr(i) .lt. holdp)then
holdp=wr(i)
np(1)=i
endif
endif
endif
endif
10   continue
c
c     write (19,*) 'numev',numev,'   j',j,'   k',k
c
if(j .gt. 3 .or. k .gt. 3)then
j=0
k=0
do 20 i=1,6
if (dabs(wi(i)) .lt. 1.d-9)then
p(3)=-wr(i)
call vtot(ae,p,v)
if (v(3) .gt. 0.d0)then
j=j+1
neg(j)=i
else
k=k+1
npos(k)=i
endif
elseif (wi(i) .lt. 0.d0)then
j=j+1
neg(j)=i
else
k=k+1
npos(k)=i
endif
20     continue
endif
c
j=0
k=0
if (numev .le. 2)then
do 755 i=1,3
if (npos(i) .ne. np(1)) then
j=j+1
hfrac=dsqrt(dble(zr(1,npos(i))**2.d0)+dble
*    (zr(2,npos(i)))**2.d0)/dsqrt(dble(zr(1,npos(i)))**2.d0
*    +dble(zr(2,npos(i)))**2.d0+dble(zr(3,npos(i)))**2.d0)
if (j .eq. 1) then
hfrac1=hfrac
ihld=i
elseif (dabs(hfrac1) .gt. dabs(hfrac))then
np(3)= npos(ihld)
np(2)= npos(i)
else
np(3)= npos(i)
np(2)= npos(ihld)
endif
endif
if (neg(i) .ne. np(4)) then
k=k+1
hfracn=dsqrt(dble(zr(1,neg(i)))**2.d0+
*    dble(zr(2,neg(i)))**2.d0)/dsqrt(dble(zr(1,neg(i)))**2.d0
*    +dble(zr(2,neg(i)))**2.d0+dble(zr(3,neg(i)))**2.d0)
if (k .eq. 1) then
hfracn1=hfracn
ihld=i
elseif (dabs(hfracn1) .gt. dabs(hfracn))then
np(6)= neg(ihld)
np(5)= neg(i)
else
np(6)= neg(i)
np(5)= neg(ihld)
endif
endif
755    continue
c
elseif (numev .le. 4)then
c
do 756 i=1,3
c       if (npos(i) .ne. np(1)) then
c          j=j+1
c          if (j .eq. 1) then
c            if (wi(npos(i)) .gt. 0.d0)then
c              ihld=i
c            else
c              if (dabs(dble(zr(2,npos(i)))) .gt.
c     *            dabs(dble(zr(3,npos(i)))))then
c                 np(3)=npos(i)
c              else
c                 np(2)=npos(i)
c              endif
c            endif
c          else
c            if (wi(npos(i)) .gt. 0.d0)then
c              if (np(3) .eq. 0)then
c                np(3)=npos(i)
c              else
c                np(2)=npos(i)
c              endif
c            else
c              if (dabs(dble(zr(2,npos(i)))) .gt.
c     *            dabs(dble(zr(3,npos(i)))))then
c                 np(2)=npos(ihld)
c                 np(3)=npos(i)
c              else
c c                np(2)=npos(i)
c                 np(3)=npos(ihld)
c              endif
c            endif
c          endif
c         endif
if (npos(i) .ne. np(1)) then
if (wi(npos(i)) .gt. 0.d0)then
np(5-id)=npos(i)
c              np(2)=neg(i)
else
np(id)=npos(i)
c              np(3)=neg(i)
endif
endif
c         if (neg(i) .ne. np(4)) then
c          k=k+1
c          if (k .eq. 1) then
c            if (wi(neg(i)) .gt. 0.d0)then
c              jhld=i
c            else
c              if (dabs(dble(zr(2,neg(i)))) .gt.
c     *            dabs(dble(zr(3,neg(i)))))then
c                 np(6)=neg(i)
c              else
c c                np(5)=neg(i)
c              endif
c            endif
c          else
c           if (wi(neg(i)) .gt. 0.d0)then
c              if (np(6) .eq. 0)then
c                np(6)=neg(i)
c              else
c                np(5)=neg(i)
c              endif
c           else
c              if (dabs(dble(zr(2,neg(i)))) .gt.
c     *            dabs(dble(zr(3,neg(i)))))then
c                 np(5)=neg(jhld)
c                 np(6)=neg(i)
c              else
c                 np(5)=neg(i)
c                 np(6)=neg(jhld)
c              endif
c            endif
c           endif
c          endif
c
if (neg(i) .ne. np(4)) then
if (wi(neg(i)) .lt. 0.d0)then
np(8-id)=neg(i)
c              np(5)=neg(i)
else
c              np(6)=neg(i)
np(3+id)=neg(i)
endif
endif
756   continue
c
elseif (numev .le. 6)then
c
do 757 i=1,3
if (npos(i) .ne. np(1)) then
if (np(2) .eq. 0.d0)then
np(2)=npos(i)
else
np(3)=npos(i)
endif
endif
if (neg(i) .ne. np(4)) then
if (np(5) .eq. 0.d0)then
np(5)=neg(i)
else
np(6)=neg(i)
endif
endif
757   continue
endif
c
c     write (19,*) 'np',np
c
c Switch the eigenvalues and eigenvectors to the order P,S1,S2
c
do i=1,6
wc(i)=dcmplx(wr(np(i)),wi(np(i)))
do j=1,3
zc(j,i)=dcmplx(zr(j,np(i)),zi(j,np(i)))
c         zc(j+3,i)=dcmplx(zr(j+3,np(i)),zi(j+3,np(i)))*(0.d0,1.d0)
zc(j+3,i)=dcmplx(zr(j+3,np(i)),zi(j+3,np(i)))
enddo
enddo
c
return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
SUBROUTINE VTOT(CC,P,V)
C
IMPLICIT REAL*8 (A-H,O-Z)
integer lay,id
logical iso
DIMENSION CC(3,3,3,3),A(3,3),CP(3,3),DCP(3,3,3),V(3),P(3),DET(3)
COMMON/info/lay,id,iso
C
if (iso)then
c
pmag=p(1)**2.d0+p(2)**2.d0+p(3)**2.d0
do 10 i=1,3
v(i)=p(i)/pmag
10     continue
c
else
c
DO 28 I=1,3
DO 27 L=1,3
PSUM = 0.0D0
DO 26 J=1,3
DO 25 K=1,3
PSUM = PSUM+CC(I,J,K,L)*P(J)*P(K)
25            CONTINUE
26          CONTINUE
CP(I,L) = PSUM
CP(L,I) = PSUM
27       CONTINUE
28     CONTINUE
C
C
DO 34 M=1,3
do 33 i=1,3
DO 32 L=1,3
X = 0.0D0
Y = 0.0D0
DO 30 J=1,3
X = X+CC(I,J,M,L)*P(J)
30           CONTINUE
DO 31 K=1,3
Y = Y+CC(I,M,K,L)*P(K)
31           CONTINUE
DCP(M,I,L) = X+Y
32        CONTINUE
33      CONTINUE
34    CONTINUE
C
C
DO 45 I=1,3
DET(I)=0.0D0
45    CONTINUE
C
C
DO 55 N=1,3
DO 54 M=1,3
DO 51 I=1,3
DO 50 J=1,3
A(I,J)=CP(I,J)
50           CONTINUE
51         CONTINUE
DO 52 I=1,3
A(I,I) = CP(I,I)-1.D0
52         CONTINUE
DO 53 I=1,3
A(I,M) = DCP(N,I,M)
53         CONTINUE
Z1 = A(1,1)*A(2,2)*A(3,3)-A(1,1)*A(2,3)*A(3,2)
Z2 = -A(1,2)*A(2,1)*A(3,3)+A(1,2)*A(2,3)*A(3,1)
Z3 = A(1,3)*A(2,1)*A(3,2)-A(1,3)*A(2,2)*A(3,1)
DET(N) = DET(N)+Z1+Z2+Z3
54      CONTINUE
55    CONTINUE
C
DEN = 0.0D0
DO 65 I=1,3
DEN = DEN + P(I)*DET(I)
65    CONTINUE
C
DO 71 I=1,3
V(I) = DET(I)/DEN
71    CONTINUE
c
endif
C
RETURN
END
c
c........|.........|.........|.........|.........|.........|.........|..
c
subroutine isoeig(den,p,c33,c44,wc,zc)

c Copyright (c) 1996 W. S. Guest and C. J. Thomson.

c
c Original version: W.S. Guest   November 1988
c
c Modified/Corrected: C. Thomson   October 1993
c                     Sean's original lines commented out for
c                     comparison.
c  NOTE: c33 = lambda + 2mu  and  c44 = mu   in Lame notation
c
c
c
c ``This subroutine calculates the displacement-stress eigenvectors
c and the eigenvalues in the case of an isotropic medium. The coding
c follows Fryer and Frazer (1987) with the obvious errors corrected.
c The vectors are choosen to be in a right handed system. P is along
c the slowness vector, SV is pointed upwards, and SH is horizontal
c but its direction depends on the value of p(x) (It goes to the
c left of the p(x) direction). All of the vectors are normalized so
c that the displacements are of unit magnitude.
c
c     Orientations rechecked   07/92     Sean Guest''
c
implicit real*8(a-h,o-z)
c
integer i
dimension p(3)
complex*16 wc(6),zc(6,6),norm,acomp,ccomp,bcomp
c
p2=p(1)**2.d0 +p(2)**2.d0
fac1=den/c33-p2
fac2=den/c44-p2
c     write (19,*) p2,fac1,fac2
c
c Define the eigenvalues
c
wc(1)=cdsqrt(dcmplx(fac1,0.d0))
wc(2)=cdsqrt(dcmplx(fac2,0.d0))
wc(3)=wc(2)
wc(4)=-wc(1)
wc(5)=-wc(2)
wc(6)=-wc(3)
c
c     write (19,'(12e10.3)') wc
c
c Evaluation of the Eigenvectors
c Solution of the Eigenvectors
c
do 10 i=0,3,3
c
zc(1,1+i)=p(1)
zc(2,1+i)=p(2)
zc(3,1+i)=wc(1+i)
zc(4,1+i)=2.d0*p(1)*c44*wc(1+i)
zc(5,1+i)=2.d0*p(2)*c44*wc(1+i)
zc(6,1+i)=den-2.d0*p2*c44
c
c      zc(4,1+i)=-2.d0*p(1)*c44*wc(1+i)
c      zc(5,1+i)=-2.d0*p(2)*c44*wc(1+i)
c      zc(6,1+i)=2.d0*p2*c44-den
c
zc(1,2+i)=p(1)
zc(2,2+i)=p(2)
zc(3,2+i)=-p2/wc(2+i)
zc(4,2+i)=p(1)*zc(6,1+i)/wc(2+i)
zc(5,2+i)=p(2)*zc(6,1+i)/wc(2+i)
zc(6,2+i)=-2.d0*p2*c44
c
c      zc(4,2+i)=p(1)*(2.d0*p2*c44-den)/wc(2+i)
c      zc(5,2+i)=p(2)*(2.d0*p2*c44-den)/wc(2+i)
c      zc(6,2+i)=2.d0*p2*c44
c      if (dble(wc(2+i)) .gt. 0.d0)then
c        do j=1,6
c          zc(j,2+i) = -zc(j,2+i)
c        enddo
c      endif
c
zc(1,3+i)=-p(2)
zc(2,3+i)=p(1)
zc(3,3+i)=0.d0
zc(4,3+i)=-p(2)*c44*wc(3+i)
zc(5,3+i)=p(1)*c44*wc(3+i)
zc(6,3+i)=0.d0
c
10   continue
c
cs      do 40 j=1,6
cs        do 30 i=4,6
cs          zc(i,j)=zc(i,j)*(0.d0,-1.d0)
cs 30    continue
cs 40   continue
c
c Normalize the vectors based on displacement
c
do 50 j=1,6
acomp=zc(1,j)*zc(1,j)
bcomp=zc(2,j)*zc(2,j)
ccomp=zc(3,j)*zc(3,j)
norm=cdsqrt(acomp+bcomp+ccomp)
do 60 i=1,6
zc(i,j)=zc(i,j)/norm
60     continue
50   continue
c
return
end
c
c........|.........|.........|.........|.........|.........|.........|..
c
C                                                                 ```