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 subroutine aniso6cg(a,rho,px,py,eval,evec) implicit real*8 (a-h,o-z) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c c Finds e-vals and e-vecs for given layer. c Based on EISPACK routines for finding vertical slowness given c horizontal slownesses by solving the 6X6 matrix problem. c c COMPLEX SYSTEM MATRIX ALLOWED -- this version is for complex c horizontal slowness. c dimension a(3,3,3,3) complex*16 px,py,p(3),eval(6),evec(6,6),xnorm logical iso common/info/lay,id,iso iso=.false. p(1)=px p(2)=py call eigc (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 eigc (den,ae,p,wc,zc) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c c Based on subroutine eig by: c W.S. Guest, October 1988. c c Modified by C J Thomson, February 1995. c This version allows for complex horizontal slowness. c c Sean's comments from the original follow (but do not relate exactly c to the new version) ... 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 c This new version calls isoeigc or Eispack routine cg and the sorting c is done by routine sortslwc. c c C. Thomson February 1995. 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) dimension ct2(3,3,3),ct(3,3),qc(3,3,3,3) dimension amr(6,6),ami(6,6) dimension fv1(6),fv2(6),fv3(6) dimension zi(6,6),wr(6),wi(6),zr(6,6) complex*16 p(3),wc(6),zc(6,6) complex*16 t(3,3),tt(3,3),s(3,3) c Common Block common/info/lay,id,iso c c c write (*,*) 'in eigc',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 (*,*)'isoeigc used to solve "A"' c33=c(3,3,3,3) c44=c(2,3,2,3) c write (*,*) c33,c44 call isoeigc(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,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 amr(i,k) = dreal(tt(i,k)) ami(i,k) = dimag(tt(i,k)) amr(i,l) = ci(i,k) ami(i,l) = 0.d0 amr(j,k) = dreal(s(i,k)) ami(j,k) = dimag(s(i,k)) amr(j,l) = dreal(t(i,k)) ami(j,l) = dimag(t(i,k)) 140 continue 135 continue c c Find the eigenvalues and eigenvectors c n6=6 call cg(n6,n6,amr,ami,wr,wi,1,zr,zi,fv1,fv2,fv3,ierr) if(ierr.ne.0)print *,'**IERR.ne.0 afer CG -- in eigc' c call sortslwc(wr,wi,zr,zi,wc,zc) 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 sortslwc(wr,wi,zr,zi,wc,zc) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. implicit real*8 (a-h,o-z) c c Sorts eigenvalues/vectors into separate up/down sets. c c This version sorts real and imaginary parts into positive and c negative groups. Does not sort within a group (unecessary?). c c C. Thomson, February 1995. c dimension wr(6),wi(6),zr(6,6),zi(6,6) dimension wwrkr(6),wwrki(6),zwrkr(6,6),zwrki(6,6) complex*16 wc(6),zc(6,6) c c Sort according to size of p3 (vertical slowness) c ipos=0 ineg=0 acc=1.d-12 do 100 l=1,6 c test=0.d0 if(wr(l).eq.0.d0)then if(wi(l).eq.0.d0)then print *,'*** error 1 in sortslwc' stop endif test=2.d0*acc else if(wi(l).ne.0.d0)then test=dabs(wi(l)/wr(l)) endif endif if(test.gt.acc)then if(wi(l).gt.0.d0)then ipos=ipos+1 wwrkr(ipos)=wr(l) wwrki(ipos)=wi(l) call trnsfv6(zr(1,l),zi(1,l), . zwrkr(1,ipos),zwrki(1,ipos)) else ineg=ineg+1 wwrkr(ineg+3)=wr(l) wwrki(ineg+3)=wi(l) call trnsfv6(zr(1,l),zi(1,l), . zwrkr(1,ineg+3),zwrki(1,ineg+3)) endif else if(wr(l).gt.0.d0)then ipos=ipos+1 wwrkr(ipos)=wr(l) wwrki(ipos)=wi(l) call trnsfv6(zr(1,l),zi(1,l), . zwrkr(1,ipos),zwrki(1,ipos)) else ineg=ineg+1 wwrkr(ineg+3)=wr(l) wwrki(ineg+3)=wi(l) call trnsfv6(zr(1,l),zi(1,l), . zwrkr(1,ineg+3),zwrki(1,ineg+3)) endif endif 100 continue if(ipos.ne.3.or.ineg.ne.3)print *,'**error 2 in sortslwc' do 101 i=1,6 wc(i)=dcmplx(wwrkr(i),wwrki(i)) do 102 l=1,6 zc(l,i)=dcmplx(zwrkr(l,i),zwrki(l,i)) 102 continue 101 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine trnsfv6(vin1,vin2,vout1,vout2) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. implicit real*8 (a-h,o-z) dimension vin1(6),vin2(6),vout1(6),vout2(6) do 100 l=1,6 vout1(l)=vin1(l) vout2(l)=vin2(l) 100 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c