C
c c........|.........|.........|.........|.........|.........|.........|.. c c Subroutines used by program rmatrix1.f, which c computes r/t matrices for a stack of homogeneous c layers. Uses matrix algorithm of B. L. N. Kennett. c General anisotropy permitted. c c C. J. Thomson. c c This version allows complex px,py and omega. February 1995. c c See subroutine evset (next) for important comments regarding c eigenvalue/eigenvector calculations. 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 subroutine evset(px,py,omega,ilay,eval,evec) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c Sets up the e-values and e-vectors for a given layer. c c Calls routine isotroc f the layer is isotropic. c c For anisotropic layers variable itype=1,2,3 controls the method c used to find the e-vals/e-vecs. c c itype=1 : calls Eispack routine cg for a complex general system c matrix A -- i.e. for the case of complex slowness or c elastic parameters (note complex a(i,j,k,l) not allowed c at present but code easily extended). c c itype=2 : for the case of a real system matrix A; calls aniso6c1 c which in turn calls Sean Guest's routine eig which uses c public domain e-value/e-vector routines that look c similar to Eispack routines. c c itype=3 : again for real A matrices; this time the vertical c slownesses are found by a 3x3 matrix method rather c like ray theory. Calls aniso3c which uses Georg c R"umpker's routine eigp3 which uses a public-domain c routine for the roots of a sextic. c c See routines aniso6cg, aniso6c1 and polarc to check normalisation c of the eigenvectors (amplitude of displacement part set to c unity). c c IMPORTANT NOTE CONCERNING EIGENVECTOR ORDERING: c The first three e-vals/e-vecs correspond to downgoing or c downwardly decaying waves and the last three to upgoing or c upwardly decaying waves. This is all that the R/T matrix c calculation really requires. In the isotropic case there is c easy sorting within these two groups into P, SV and SH. c In the anisotropic case the sorting is generally from largest c to smallest vertical slowness. Polarisations in anisotropic c media can behave in interesting ways and so it is assumed the c user will take care to understand which wave function should c be called, say, qSH and which column it occupies (simply c by looking at the eigenvectors). c implicit real*8(a-h,o-z) complex*16 px,py,omega complex*16 eval(6), evec(6,6) c c COMMON BLOCKS to be INCLUDEd c c max number of layers --- 1 = half-space above or top layer c in the stack just below free surface c nlay = half-space below c parameter (nlaymx=30) c c Earth model c common /elast/isoflg(nlaymx),thickn(nlaymx),a(3,3,3,3,nlaymx), . rho(nlaymx) c complex*16 evals,evecs common /evprob/evals(6,nlaymx),evecs(6,6,nlaymx) complex*16 q11i,q12i,q21i,q22i complex*16 tui,rui,tdi,rdi common /intfce/q11i(3,3,nlaymx),q12i(3,3,nlaymx), . q21i(3,3,nlaymx),q22i(3,3,nlaymx), . tui(3,3,nlaymx),rui(3,3,nlaymx), . tdi(3,3,nlaymx),rdi(3,3,nlaymx) complex*16 tusi,rusi,tdsi,rdsi common /stackrt/tusi(3,3,nlaymx),rusi(3,3,nlaymx), . tdsi(3,3,nlaymx),rdsi(3,3,nlaymx) c c symmetry testing data c complex*16 diag1,j2cst1,j2cstn common /symts/diag1(6,nlaymx) common /symts2/j2cst1(3,3),j2cstn(3,3) c c END OF COMMON BLOCKS c if(isoflg(ilay).eq.1)then call isotroc(a(1,1,1,1,ilay),rho(ilay), . px,py,eval,evec) c-LK call write1(ilay,isoflg(ilay),px,py,eval,evec) else itype=1 if(itype.eq.1)then c find e-vals + e-vecs of complex general matrix ... call aniso6cg(a(1,1,1,1,ilay),rho(ilay), . px,py,eval,evec) c-LK call write1(ilay,isoflg(ilay),px,py,eval,evec) elseif(itype.eq.2)then c find e-vals + e-vecs for case of real horizontal slownesses ... call aniso6c1(a(1,1,1,1,ilay),rho(ilay), . px,py,eval,evec) c-LK call write1(ilay,isoflg(ilay),px,py,eval,evec) elseif(itype.eq.3)then c find e-vals + e-vecs for case of real horizontal slownesses c using 3x3 matrix method (like ray theory) ... call aniso3c(a(1,1,1,1,ilay),rho(ilay), . px,py,eval,evec) c-LK call write1(ilay,isoflg(ilay),px,py,eval,evec) endif endif return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine estore(eval,evec,evals,evecs) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. complex*16 eval(6),evec(6,6),evals(6),evecs(6,6) do 100 i=1,6 evals(i)=eval(i) do 200 j=1,6 evecs(i,j)=evec(i,j) 200 continue 100 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine qstore(q11,q12,q21,q22,q11i,q12i,q21i,q22i) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. complex*16 q11(3,3),q12(3,3),q21(3,3),q22(3,3) complex*16 q11i(3,3),q12i(3,3),q21i(3,3),q22i(3,3) do 100 i=1,3 do 200 j=1,3 q11i(i,j)=q11(i,j) q12i(i,j)=q12(i,j) q21i(i,j)=q21(i,j) q22i(i,j)=q22(i,j) 200 continue 100 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine rtstore(tu,ru,td,rd,tui,rui,tdi,rdi) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. complex*16 tu(3,3),ru(3,3),td(3,3),rd(3,3) complex*16 tui(3,3),rui(3,3),tdi(3,3),rdi(3,3) do 100 i=1,3 do 200 j=1,3 tui(i,j)=tu(i,j) rui(i,j)=ru(i,j) rdi(i,j)=rd(i,j) tdi(i,j)=td(i,j) 200 continue 100 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine rtget(tu,ru,td,rd,tui,rui,tdi,rdi) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. complex*16 tu(3,3),ru(3,3),td(3,3),rd(3,3) complex*16 tui(3,3),rui(3,3),tdi(3,3),rdi(3,3) do 100 i=1,3 do 200 j=1,3 tu(i,j)=tui(i,j) ru(i,j)=rui(i,j) rd(i,j)=rdi(i,j) td(i,j)=tdi(i,j) 200 continue 100 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine xmult6(a,b,c) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c c complex matrix multiplication a*b=c c dimension 6 ... allows overwriting c complex*16 a(6,6), b(6,6), c(6,6), wrk(6,6) do 11 i=1,6 do 13 j=1,6 wrk(i,j)=(0.d0,0.d0) do 15 k=1,6 wrk(i,j)=wrk(i,j) + a(i,k)*b(k,j) 15 continue 13 continue 11 continue do 17 i=1,6 do 19 j=1,6 c(i,j)=wrk(i,j) 19 continue 17 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine xeveci(nlay,ilay,evec,eveci) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. implicit real*8(a-h,o-z) c c Inverts the eigenvector matrix, and sets up the diagonal 'constant' c of J-primed symmetry in the process (see theory notes, section 6, c equations (6.5,6.6)). c complex*16 evec(6,6), eveci(6,6), wrk(6,6) c c COMMON BLOCKS to be INCLUDEd c c max number of layers --- 1 = half-space above or top layer c in the stack just below free surface c nlay = half-space below c parameter (nlaymx=30) c c Earth model c common /elast/isoflg(nlaymx),thickn(nlaymx),a(3,3,3,3,nlaymx), . rho(nlaymx) complex*16 evals,evecs common /evprob/evals(6,nlaymx),evecs(6,6,nlaymx) complex*16 q11i,q12i,q21i,q22i complex*16 tui,rui,tdi,rdi common /intfce/q11i(3,3,nlaymx),q12i(3,3,nlaymx), . q21i(3,3,nlaymx),q22i(3,3,nlaymx), . tui(3,3,nlaymx),rui(3,3,nlaymx), . tdi(3,3,nlaymx),rdi(3,3,nlaymx) complex*16 tusi,rusi,tdsi,rdsi common /stackrt/tusi(3,3,nlaymx),rusi(3,3,nlaymx), . tdsi(3,3,nlaymx),rdsi(3,3,nlaymx) c c symmetry testing data c complex*16 diag1,j2cst1,j2cstn common /symts/diag1(6,nlaymx) common /symts2/j2cst1(3,3),j2cstn(3,3) c c END OF COMMON BLOCKS c do 11 i=1,3 do 15 j=1,3 eveci(i,j)=evec(j+3,i) eveci(i,j+3)=evec(j,i) eveci(i+3,j)=evec(j+3,i+3) eveci(i+3,j+3)=evec(j,i+3) 15 continue 11 continue call xmult6(eveci,evec,wrk) iwrite=0 if(iwrite.eq.1)then write(11,*) write(11,*)'------------------------------------------' write(11,*) write(11,*)'in xeveci...diag constant, layer ',ilay,'...' write(11,*)(wrk(i,i),i=1,6) endif do 17 i=1,6 diag1(i,nlay+1-ilay)=wrk(i,i) do 19 j=1,6 eveci(i,j)=eveci(i,j)/wrk(i,i) 19 continue 17 continue if(iwrite.eq.1)then write(11,*) write(11,*)'------------------------------------------' write(11,*) endif return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine trnsf6(e1,e2,x1,x2) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. complex*16 e1(6),e2(6),x1(6,6),x2(6,6) do 100 i=1,6 e1(i)=e2(i) do 200 j=1,6 x1(i,j)=x2(i,j) 200 continue 100 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine retra(q11,q12,q21,q22,tu,ru,td,rd) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c Extracts the R/T coefficients from the Q-matrix for a given single c interface (see theory notes, equation (3.4)). implicit real*8(a-h,o-z) complex*16 q11(3,3),q12(3,3),q21(3,3),q22(3,3) complex*16 tu(3,3),ru(3,3),td(3,3),rd(3,3) complex*16 tuneg(3,3),xhelp(3,3) call xinv3(q22,tu) call xmult3(q12,tu,ru) call xneg3(tu,tuneg) call xmult3(tuneg,q21,rd) call xmult3(q12,rd,xhelp) call xadd3(q11,xhelp,td) return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine qmatrix(eveci1,evec2,q11,q12,q21,q22) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c Sets up the Q matrix for a single interface from the e-vectors c defined on either side (see theory notes, equation (3.3)). implicit real*8(a-h,o-z) complex*16 eveci1(6,6),evec2(6,6),qq(6,6) complex*16 q11(3,3),q12(3,3),q21(3,3),q22(3,3) call xmult6(eveci1,evec2,qq) do 17 i=1,3 do 19 j=1,3 q11(i,j)=qq(i,j) q12(i,j)=qq(i,j+3) q21(i,j)=qq(i+3,j) q22(i,j)=qq(i+3,j+3) 19 continue 17 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine xmult3(a,b,c) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c c Complex matrix multiplication a*b=c c dimension 3 ... allows overwriting c complex*16 a(3,3), b(3,3), c(3,3), wrk(3,3) do 11 i=1,3 do 13 j=1,3 wrk(i,j)=(0.d0,0.d0) do 15 k=1,3 wrk(i,j)=wrk(i,j) + a(i,k)*b(k,j) 15 continue 13 continue 11 continue do 17 i=1,3 do 19 j=1,3 c(i,j)=wrk(i,j) 19 continue 17 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine xneg3(a,aneg) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c c Complex matrix multiplication with (-1) c complex*16 a(3,3), aneg(3,3) do 11 i=1,3 do 13 j=1,3 aneg(i,j)=-a(i,j) 13 continue 11 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine xadd3(a,b,c) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c c Complex matrix addition a+b=c c complex*16 a(3,3), b(3,3), c(3,3) do 11 i=1,3 do 13 j=1,3 c(i,j)=a(i,j)+b(i,j) 13 continue 11 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine xinv3(a,ainv) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. complex*16 a(3,3),ainv(3,3) complex*16 co(3,3), deta, detb, detc c c Inversion of 3x3 matrix c co(1,1)=(a(2,2)*a(3,3)-a(2,3)*a(3,2)) co(1,2)=-(a(2,1)*a(3,3)-a(2,3)*a(3,1)) co(1,3)=(a(2,1)*a(3,2)-a(2,2)*a(3,1)) co(2,1)=-(a(1,2)*a(3,3)-a(1,3)*a(3,2)) co(2,2)=(a(1,1)*a(3,3)-a(1,3)*a(3,1)) co(2,3)=-(a(1,1)*a(3,2)-a(1,2)*a(3,1)) co(3,1)=(a(1,2)*a(2,3)-a(1,3)*a(2,2)) co(3,2)=-(a(1,1)*a(2,3)-a(1,3)*a(2,1)) co(3,3)=(a(1,1)*a(2,2)-a(1,2)*a(2,1)) deta=a(1,1)*co(1,1)+a(1,2)*co(1,2)+a(1,3)*co(1,3) detb=a(2,1)*co(2,1)+a(2,2)*co(2,2)+a(2,3)*co(2,3) detc=a(3,1)*co(3,1)+a(3,2)*co(3,2)+a(3,3)*co(3,3) c write(11,*) c write(11,*)'det check in xinv3 ... deta, detb, detc =' c write(11,*)deta,detb,detc do 10 i=1,3 do 11 j=1,3 ainv(i,j)=co(j,i)/deta 11 continue 10 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine addit(nlay,ilay,tus,rus,tds,rds, . tu,ru,td,rd) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c c Addition/recursion formulas (see theory notes, section 4, c equations (4.8)). c complex*16 tus(3,3),rus(3,3),tds(3,3),rds(3,3) complex*16 tu(3,3),ru(3,3),td(3,3),rd(3,3) c workspaces complex*16 tusw(3,3),rusw(3,3),tdsw(3,3),rdsw(3,3) complex*16 ximat(3,3),reverb(3,3),reverbi(3,3) save ximat if(ilay.eq.nlay-1)then c first interface -- trivial call trnsf3(tus,tu) call trnsf3(rus,ru) call trnsf3(tds,td) call trnsf3(rds,rd) do 100 i=1,3 do 101 j=1,3 ximat(i,j)=(0.d0,0.d0) 101 continue ximat(i,i)=(1.d0,0.d0) 100 continue return endif iwrite=0 if(iwrite.eq.1)then write(11,*) write(11,*)'--------------------------------------------------' write(11,*) write(11,*)'inside addit on recursion' do 200 i=1,3 write(11,1000)(ximat(i,j),j=1,3) 200 continue endif 1000 format(3(' (',d10.4,',',d10.4,')')) c symbolic notation... reverb = ximat - rds * ru call xmult3(rds,ru,reverb) call xneg3(reverb,reverb) call xadd3(ximat,reverb,reverb) call xinv3(reverb,reverbi) if(iwrite.eq.1)then write(11,*) write(11,*)'reverb and reverbi...' do 201 i=1,3 write(11,1000)(reverb(i,j),j=1,3) 201 continue do 202 i=1,3 write(11,1000)(reverbi(i,j),j=1,3) 202 continue endif c addition rules... with symbolic notation comments c tusw = tu * reverbi * tus call xmult3(reverbi,tus,tusw) call xmult3(tu,tusw,tusw) c rdsw = rd + tu * reverbi * rds * td call xmult3(rds,td,rdsw) call xmult3(reverbi,rdsw,rdsw) call xmult3(tu,rdsw,rdsw) call xadd3(rd,rdsw,rdsw) c rusw = rus + tds * ru * reverbi * tus call xmult3(reverbi,tus,rusw) call xmult3(ru,rusw,rusw) call xmult3(tds,rusw,rusw) call xadd3(rus,rusw,rusw) c tdsw = tds * [ ximat + ru*reverbi*rds ] * td call xmult3(reverbi,rds,tdsw) call xmult3(ru,tdsw,tdsw) call xadd3(ximat,tdsw,tdsw) call xmult3(tdsw,td,tdsw) call xmult3(tds,tdsw,tdsw) c copy over call trnsf3(tus,tusw) call trnsf3(rus,rusw) call trnsf3(tds,tdsw) call trnsf3(rds,rdsw) return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine layer(ilay,omega,xh,eval,tu,ru,td,rd) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c c Includes the phase shift across a layer (theory notes, section 4, c equations (4.6)). c implicit real*8(a-h,o-z) complex*16 omega complex*16 xi complex*16 eval(6) complex*16 tu(3,3),ru(3,3),td(3,3),rd(3,3) complex*16 ed(3,3),eui(3,3) iwrite=0 if(iwrite.eq.1)then write(11,*) write(11,*)'.................................................' write(11,*) write(11,*)'crossing layer ',ilay,' omega, xh = ',omega,xh write(11,*) endif xi=(0.d0,1.d0) do 11 i=1,3 do 13 j=1,3 ed(i,j)=(0.d0,0.d0) eui(i,j)=(0.d0,0.d0) 13 continue 11 continue c c elements of layer matrix for downgoing waves c ed(1,1)=cdexp(xi*omega*eval(1)*xh) if(cdabs(ed(1,1)).gt.1.d0+1.d-12) .print *,'**error ed11 in layer |ed11|=',cdabs(ed(1,1)) ed(2,2)=cdexp(xi*omega*eval(2)*xh) if(cdabs(ed(2,2)).gt.1.d0+1.d-12) .print *,'**error ed22 in layer |ed22|=',cdabs(ed(2,2)) ed(3,3)=cdexp(xi*omega*eval(3)*xh) if(cdabs(ed(3,3)).gt.1.d0+1.d-12) .print *,'**error ed33 in layer |ed33|=',cdabs(ed(3,3)) c c elements of inverse layer matrix for upgoing waves c eui(1,1)=cdexp(-xi*omega*eval(4)*xh) if(cdabs(eui(1,1)).gt.1.d0+1.d-12) .print *,'**error eui11 in layer |eui11|=',cdabs(eui(1,1)) eui(2,2)=cdexp(-xi*omega*eval(5)*xh) if(cdabs(eui(2,2)).gt.1.d0+1.d-12) .print *,'**error eui22 in layer |eui22|=',cdabs(eui(2,2)) eui(3,3)=cdexp(-xi*omega*eval(6)*xh) if(cdabs(eui(3,3)).gt.1.d0+1.d-12) .print *,'**error eui33 in layer |eui33|=',cdabs(eui(3,3)) call xmult3(td,ed,td) call xmult3(rd,ed,rd) call xmult3(eui,rd,rd) call xmult3(eui,tu,tu) if(iwrite.eq.1)then write(11,*) write(11,*)'.................................................' write(11,*) endif return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine trnsf3(x1,x2) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. complex*16 x1(3,3),x2(3,3) do 100 i=1,3 do 200 j=1,3 x1(i,j)=x2(i,j) 200 continue 100 continue return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine j1symtst(nlay,ilay,tus,rus,tds,rds) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c Performs the J_1 symmetry test. See section 6 of the c theory notes (equation (6.21)). Symmetries here satisfied c by perfectly elastic medium, real slownesses/frequency c (propagating waves) and general anisotropy. implicit real*8(a-h,o-z) complex*16 tus(3,3),rus(3,3),tds(3,3),rds(3,3) complex*16 wrk1(3,3),wrk2(3,3),wrk3(3,3),wrk4(3,3) complex*16 wrk5(3,3),wrk6(3,3),wrk7(3,3),wrk8(3,3) c c COMMON BLOCKS to be INCLUDEd c c max number of layers --- 1 = half-space above or top layer c in the stack just below free surface c nlay = half-space below c parameter (nlaymx=30) c c Earth model c common /elast/isoflg(nlaymx),thickn(nlaymx),a(3,3,3,3,nlaymx), . rho(nlaymx) complex*16 evals,evecs common /evprob/evals(6,nlaymx),evecs(6,6,nlaymx) complex*16 q11i,q12i,q21i,q22i complex*16 tui,rui,tdi,rdi common /intfce/q11i(3,3,nlaymx),q12i(3,3,nlaymx), . q21i(3,3,nlaymx),q22i(3,3,nlaymx), . tui(3,3,nlaymx),rui(3,3,nlaymx), . tdi(3,3,nlaymx),rdi(3,3,nlaymx) complex*16 tusi,rusi,tdsi,rdsi common /stackrt/tusi(3,3,nlaymx),rusi(3,3,nlaymx), . tdsi(3,3,nlaymx),rdsi(3,3,nlaymx) c c symmetry testing data c complex*16 diag1,j2cst1,j2cstn common /symts/diag1(6,nlaymx) common /symts2/j2cst1(3,3),j2cstn(3,3) c c END OF COMMON BLOCKS c do 100 i=1,3 do 101 j=1,3 wrk1(i,j)= diag1(i,1)*rus(i,j) wrk2(i,j)= diag1(i+3,nlay+1-ilay)*tus(i,j) wrk3(i,j)= diag1(i,1)*tds(i,j) wrk4(i,j)= diag1(i+3,nlay+1-ilay)*rds(i,j) 101 continue 100 continue tst1=(0.d0,0.d0) tst2=(0.d0,0.d0) do 200 i=1,3 do 201 j=1,3 wrk5(i,j)=(0.d0,0.d0) wrk6(i,j)=(0.d0,0.d0) wrk7(i,j)=(0.d0,0.d0) wrk8(i,j)=(0.d0,0.d0) do 202 k=1,3 wrk5(i,j)= wrk5(i,j)+dconjg(rus(k,i))*wrk1(k,j) . - dconjg(tus(k,i))*wrk2(k,j) wrk6(i,j)= wrk6(i,j)+dconjg(rus(k,i))*wrk3(k,j) . - dconjg(tus(k,i))*wrk4(k,j) wrk7(i,j)= wrk7(i,j)+dconjg(tds(k,i))*wrk1(k,j) . - dconjg(rds(k,i))*wrk2(k,j) wrk8(i,j)= wrk8(i,j)+dconjg(tds(k,i))*wrk3(k,j) . - dconjg(rds(k,i))*wrk4(k,j) 202 continue c form norms for accuracy/confidence measure tst1=tst1+cdabs(wrk6(i,j))+cdabs(wrk7(i,j)) if(i.ne.j)then tst1=tst1+cdabs(wrk5(i,j))+cdabs(wrk8(i,j)) endif if(i.eq.j)then tst2=tst2+cdabs(wrk5(i,j))+cdabs(wrk8(i,j)) endif 201 continue 200 continue write(11,*) write(11,*)'------------------------------------------' write(11,*) write(11,*)'test of stack r/t j1 symmetries...' write(11,*)'(propagating waves only)...' write(11,*)'diagonal constants check...' write(11,*)(diag1(i+3,1),i=1,3) write(11,*)(diag1(i,nlay+1-ilay),i=1,3) c1000 format(3(' (',d9.3,',',d9.3,')')) write(11,*) write(11,*)'11 block...' write(11,*)wrk5(1,1) do 300 i=1,3 write(11,1001)(wrk5(i,j),j=1,3) 300 continue 1001 format(3(' (',d9.3,',',d9.3,')')) write(11,*) write(11,*)'12 block...' write(11,*)wrk6(1,1) do 301 i=1,3 write(11,1001)(wrk6(i,j),j=1,3) 301 continue write(11,*) write(11,*)'21 block...' write(11,*)wrk7(1,1) do 302 i=1,3 write(11,1001)(wrk7(i,j),j=1,3) 302 continue write(11,*) write(11,*)'22 block...' write(11,*)wrk8(1,1) do 303 i=1,3 write(11,1001)(wrk8(i,j),j=1,3) 303 continue write(11,*) write(11,*)'in F format...' write(11,*)'11 block...' write(11,*)wrk5(1,1) do 400 i=1,3 write(11,1002)(wrk5(i,j),j=1,3) 400 continue 1002 format(3(' (',f8.1,',',f8.1,')')) write(11,*) write(11,*)'12 block...' write(11,*)wrk6(1,1) do 401 i=1,3 write(11,1002)(wrk6(i,j),j=1,3) 401 continue write(11,*) write(11,*)'21 block...' write(11,*)wrk7(1,1) do 402 i=1,3 write(11,1002)(wrk7(i,j),j=1,3) 402 continue write(11,*) write(11,*)'22 block...' write(11,*)wrk8(1,1) do 403 i=1,3 write(11,1002)(wrk8(i,j),j=1,3) 403 continue symacc=tst1/tst2 write(11,*) write(11,*)'measure of j1 accuracy... symacc = ',symacc write(11,*) write(11,*)'------------------------------------------' write(11,*) return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine j2sym(nlay,ilay,evec,j2cnst) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c Sets up 'constants' for use in J_2 symmetry test c (horizontal mirror plane media, see theory section (6)). implicit real*8(a-h,o-z) complex*16 evec(6,6),j2cnst(3,3) complex*16 wrk1(6,6),wrk2(6,6) iwrite=0 if(iwrite.eq.1)then write(11,*) write(11,*)'------------------------------------------' write(11,*) write(11,*)'in j2sym, nlay=',nlay,' ilay=',ilay endif do 100 i=1,3 do 101 j=1,3 if(i.lt.3)then wrk1(i,j)= evec(i+3,j) wrk1(i+3,j)= -evec(i,j) wrk1(i,j+3)= evec(i+3,j+3) wrk1(i+3,j+3)= -evec(i,j+3) else wrk1(i,j)= -evec(i+3,j) wrk1(i+3,j)= evec(i,j) wrk1(i,j+3)= -evec(i+3,j+3) wrk1(i+3,j+3)= evec(i,j+3) endif 101 continue 100 continue do 200 i=1,6 do 201 j=1,6 wrk2(i,j)=(0.d0,0.d0) do 202 k=1,6 wrk2(i,j)=wrk2(i,j)+evec(k,i)*wrk1(k,j) 202 continue 201 continue 200 continue if(iwrite.eq.1)then write(11,*) write(11,*)'value of j2 constant...' do 300 i=1,6 write(11,3000)(wrk2(i,j),j=1,6) 300 continue write(11,*)'and in F format...' do 301 i=1,6 write(11,3001)(wrk2(i,j),j=1,6) 301 continue endif 3000 format(6(' (',d10.4,',',d10.4,')')) 3001 format(6(' (',f8.1,',',f8.1,')')) do 400 i=1,3 do 401 j=1,3 j2cnst(i,j)= wrk2(i,j+3) 401 continue 400 continue if(iwrite.eq.1)then write(11,*) write(11,*)'------------------------------------------' write(11,*) endif return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine j2symtst(nlay,ilay,tus,rus,tds,rds) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c Performs J_2 symmtery test using previously set up c 'constants' (theory section (6)). Media with horizontal c mirror plane, but general anisotropy, anelasticity, allowed. implicit real*8(a-h,o-z) complex*16 tus(3,3),rus(3,3),tds(3,3),rds(3,3) complex*16 wrk1(3,3),wrk2(3,3),wrk3(3,3),wrk4(3,3) c c COMMON BLOCKS to be INCLUDEd c c max number of layers --- 1 = half-space above or top layer c in the stack just below free surface c nlay = half-space below c parameter (nlaymx=30) c c Earth model c common /elast/isoflg(nlaymx),thickn(nlaymx),a(3,3,3,3,nlaymx), . rho(nlaymx) complex*16 evals,evecs common /evprob/evals(6,nlaymx),evecs(6,6,nlaymx) complex*16 q11i,q12i,q21i,q22i complex*16 tui,rui,tdi,rdi common /intfce/q11i(3,3,nlaymx),q12i(3,3,nlaymx), . q21i(3,3,nlaymx),q22i(3,3,nlaymx), . tui(3,3,nlaymx),rui(3,3,nlaymx), . tdi(3,3,nlaymx),rdi(3,3,nlaymx) complex*16 tusi,rusi,tdsi,rdsi common /stackrt/tusi(3,3,nlaymx),rusi(3,3,nlaymx), . tdsi(3,3,nlaymx),rdsi(3,3,nlaymx) c c symmetry testing data c complex*16 diag1,j2cst1,j2cstn common /symts/diag1(6,nlaymx) common /symts2/j2cst1(3,3),j2cstn(3,3) c c END OF COMMON BLOCKS c do 100 i=1,3 do 101 j=1,3 wrk1(i,j)=(0.d0,0.d0) wrk2(i,j)=(0.d0,0.d0) wrk3(i,j)=(0.d0,0.d0) wrk4(i,j)=(0.d0,0.d0) do 102 k=1,3 wrk1(i,j)= wrk1(i,j)+rus(k,i)*j2cstn(k,j) wrk2(i,j)= wrk2(i,j)+rds(k,i)*j2cst1(j,k) wrk3(i,j)= wrk3(i,j)+j2cstn(k,i)*tds(k,j) wrk4(i,j)= wrk4(i,j)+tus(k,i)*j2cst1(j,k) 102 continue 101 continue 100 continue tst1a=0.d0 tst1b=0.d0 tst2a=0.d0 tst2b=0.d0 tst3a=0.d0 tst3b=0.d0 do 103 i=1,3 do 104 j=1,3 tst1a=tst1a+cdabs(wrk1(i,j)) tst1b=tst1b+cdabs(wrk1(i,j)-wrk1(j,i)) tst2a=tst2a+cdabs(wrk2(i,j)) tst2b=tst2b+cdabs(wrk2(i,j)-wrk2(j,i)) tst3a=tst3a+cdabs(wrk3(i,j))+cdabs(wrk4(i,j)) tst3b=tst3b+cdabs(wrk3(i,j)-wrk4(i,j)) 104 continue 103 continue write(11,*) write(11,*)'------------------------------------------' write(11,*) write(11,*)'test of stack r/t j2 symmetries...' write(11,*)'(TI media only)...' c......................................................................| write(11,*)'wrk1...' do 201 i=1,3 write(11,3000)(wrk1(i,j),j=1,3) 201 continue write(11,*)'wrk2...' do 202 i=1,3 write(11,3000)(wrk2(i,j),j=1,3) 202 continue write(11,*)'wrk3...' do 203 i=1,3 write(11,3000)(wrk3(i,j),j=1,3) 203 continue write(11,*)'wrk4...' do 204 i=1,3 write(11,3000)(wrk4(i,j),j=1,3) 204 continue write(11,*) write(11,*)'and in F format...' write(11,*)'wrk1...' do 301 i=1,3 write(11,3001)(wrk1(i,j),j=1,3) 301 continue write(11,*)'wrk2...' do 302 i=1,3 write(11,3001)(wrk2(i,j),j=1,3) 302 continue write(11,*)'wrk3...' do 303 i=1,3 write(11,3001)(wrk3(i,j),j=1,3) 303 continue write(11,*)'wrk4...' do 304 i=1,3 write(11,3001)(wrk4(i,j),j=1,3) 304 continue 3000 format(3(' (',d10.4,',',d10.4,')')) 3001 format(3(' (',f8.1,',',f8.1,')')) symac1=tst1b/tst1a symac2=tst2b/tst2a symac3=2.d0*tst3b/tst3a write(11,*) write(11,*)'measure of j2 accuracy... symac1 = ',symac1 write(11,*)'measure of j2 accuracy... symac2 = ',symac2 write(11,*)'measure of j2 accuracy... symac3 = ',symac3 write(11,*) write(11,*)'------------------------------------------' write(11,*) return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine jprtst(evalsf,evecsf,nlay, . tus,rus,tds,rds,tuspr,ruspr,tdspr,rdspr) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c Performs the J-primed symmetry test (theory section (6), eqn. (6.16)). c This is the most general symmetry. See comments below regarding the c effect of column ordering in the e-vector matrices. implicit real*8(a-h,o-z) complex*16 tus(3,3),rus(3,3),tds(3,3),rds(3,3) complex*16 ttus(3,3),trus(3,3),ttds(3,3),trds(3,3) complex*16 tuspr(3,3),ruspr(3,3),tdspr(3,3),rdspr(3,3) complex*16 ttuspr(3,3),truspr(3,3),ttdspr(3,3),trdspr(3,3) complex*16 evalsf(6,2),evecsf(6,6,2) complex*16 diagt(6,6),diagb(6,6) complex*16 tstru1,tsttd1,tsttu1,tstrd1 c c COMMON BLOCKS to be INCLUDEd c c max number of layers --- 1 = half-space above or top layer c in the stack just below free surface c nlay = half-space below c parameter (nlaymx=30) c c Earth model c common /elast/isoflg(nlaymx),thickn(nlaymx),a(3,3,3,3,nlaymx), . rho(nlaymx) c complex*16 evals,evecs common /evprob/evals(6,nlaymx),evecs(6,6,nlaymx) complex*16 q11i,q12i,q21i,q22i complex*16 tui,rui,tdi,rdi common /intfce/q11i(3,3,nlaymx),q12i(3,3,nlaymx), . q21i(3,3,nlaymx),q22i(3,3,nlaymx), . tui(3,3,nlaymx),rui(3,3,nlaymx), . tdi(3,3,nlaymx),rdi(3,3,nlaymx) complex*16 tusi,rusi,tdsi,rdsi common /stackrt/tusi(3,3,nlaymx),rusi(3,3,nlaymx), . tdsi(3,3,nlaymx),rdsi(3,3,nlaymx) c c symmetry testing data c complex*16 diag1,j2cst1,j2cstn common /symts/diag1(6,nlaymx) common /symts2/j2cst1(3,3),j2cstn(3,3) c c END OF COMMON BLOCKS c iwrite=1 if(iwrite.eq.1)then write(11,*) write(11,*)'-----------------------------------------' write(11,*) write(11,*)'in jprtst -- testing reversal symmetries' endif c set up matrices which should come out to be diagonal -- or .... c Important Note. The up/down columns of the e-vector matrix c are interchanged for the reversed system, at least in the c theoretical description (papers, notes). However, in the code c the columns of the e-vector matrix are always down first, then c up (set eg in subroutine vertsl2) irrespective of the sign of the c slowness. Therefore, in this code the so-called diagonal c matrices actually come out to be trailing diagonal. This causes c no particular problem so long as one is aware of it. c c Furthermore, the sorting in, eg, vertsl2 does not ensure that c the columns within the up/down blocks are always in the order c P, S1, S2. They may be interchanged and this means that the c matrix d-primed in the notes may be block diagonal (with 1's and c zeros) not truly diagonal. do 100 i=1,3 do 200 j=1,3 diagt(i,j)=dcmplx(0.d0,0.d0) diagt(i,j+3)=dcmplx(0.d0,0.d0) diagt(i+3,j)=dcmplx(0.d0,0.d0) diagt(i+3,j+3)=dcmplx(0.d0,0.d0) diagb(i,j)=dcmplx(0.d0,0.d0) diagb(i,j+3)=dcmplx(0.d0,0.d0) diagb(i+3,j)=dcmplx(0.d0,0.d0) diagb(i+3,j+3)=dcmplx(0.d0,0.d0) do 300 l=1,3 diagt(i,j) = diagt(i,j) . +evecs(l,i,1)*evecsf(l+3,j,1) . -evecs(l+3,i,1)*evecsf(l,j,1) diagt(i,j+3) = diagt(i,j+3) . +evecs(l,i,1)*evecsf(l+3,j+3,1) . -evecs(l+3,i,1)*evecsf(l,j+3,1) diagt(i+3,j) = diagt(i+3,j) . +evecs(l,i+3,1)*evecsf(l+3,j,1) . -evecs(l+3,i+3,1)*evecsf(l,j,1) diagt(i+3,j+3) = diagt(i+3,j+3) . +evecs(l,i+3,1)*evecsf(l+3,j+3,1) . -evecs(l+3,i+3,1)*evecsf(l,j+3,1) diagb(i,j) = diagb(i,j) . +evecs(l,i,nlay)*evecsf(l+3,j,2) . -evecs(l+3,i,nlay)*evecsf(l,j,2) diagb(i,j+3) = diagb(i,j+3) . +evecs(l,i,nlay)*evecsf(l+3,j+3,2) . -evecs(l+3,i,nlay)*evecsf(l,j+3,2) diagb(i+3,j) = diagb(i+3,j) . +evecs(l,i+3,nlay)*evecsf(l+3,j,2) . -evecs(l+3,i+3,nlay)*evecsf(l,j,2) diagb(i+3,j+3) = diagb(i+3,j+3) . +evecs(l,i+3,nlay)*evecsf(l+3,j+3,2) . -evecs(l+3,i+3,nlay)*evecsf(l,j+3,2) 300 continue 200 continue 100 continue if(iwrite.eq.1)then write(11,*)'top diagonal(?) constant' do 400 i=1,6 write(11,3000)(diagt(i,j),j=1,6) 3000 format(3('(',d10.4,',',d10.4,')'),/, . 5X,3('(',d10.4,',',d10.4,')')) 400 continue write(11,*) write(11,*)'bottom diagonal(?) constant' do 500 i=1,6 write(11,3000)(diagb(i,j),j=1,6) 500 continue endif c r/t coefficient symmetries tstru1=dcmplx(0.d0,0.d0) tsttd1=dcmplx(0.d0,0.d0) tsttu1=dcmplx(0.d0,0.d0) tstrd1=dcmplx(0.d0,0.d0) tstru2=0.d0 tsttd2=0.d0 tsttu2=0.d0 tstrd2=0.d0 do 600 i=1,3 do 700 j=1,3 trus(i,j) = dcmplx(0.d0,0.d0) ttds(i,j) = dcmplx(0.d0,0.d0) ttus(i,j) = dcmplx(0.d0,0.d0) trds(i,j) = dcmplx(0.d0,0.d0) truspr(i,j) = dcmplx(0.d0,0.d0) ttdspr(i,j) = dcmplx(0.d0,0.d0) ttuspr(i,j) = dcmplx(0.d0,0.d0) trdspr(i,j) = dcmplx(0.d0,0.d0) do 701 k=1,3 trus(i,j) = trus(i,j) + diagb(i+3,k)*rus(k,j) ttds(i,j) = ttds(i,j) + diagb(i+3,k)*tds(k,j) ttus(i,j) = ttus(i,j) + diagt(i,k+3)*tus(k,j) trds(i,j) = trds(i,j) + diagt(i,k+3)*rds(k,j) truspr(i,j) = truspr(i,j) + ruspr(k,i)*diagb(k,j+3) ttdspr(i,j) = ttdspr(i,j) + tuspr(k,i)*diagt(k+3,j) ttuspr(i,j) = ttuspr(i,j) + tdspr(k,i)*diagb(k,j+3) trdspr(i,j) = trdspr(i,j) + rdspr(k,i)*diagt(k+3,j) 701 continue tstru1= tstru1+truspr(i,j)+trus(i,j) tstru2= tstru2+cdabs(truspr(i,j))+cdabs(trus(i,j)) tsttd1= tsttd1+ttdspr(i,j)-ttds(i,j) tsttd2= tsttd2+cdabs(ttdspr(i,j))+cdabs(ttds(i,j)) tsttu1= tsttu1+ttuspr(i,j)-ttus(i,j) tsttu2= tsttu2+cdabs(ttuspr(i,j))+cdabs(ttus(i,j)) tstrd1= tstrd1+trdspr(i,j)+trds(i,j) tstrd2= tstrd2+cdabs(trdspr(i,j))+cdabs(trds(i,j)) 700 continue 600 continue if(iwrite.eq.1)then write(11,*)'reversed and scaled r/t matrices' write(11,*)'ru-primed' write(11,3001)((truspr(i,j),j=1,3),i=1,3) 3001 format(3('(',d10.4,',',d10.4,')')) write(11,*)'ru' write(11,3001)((trus(i,j),j=1,3),i=1,3) write(11,*)'td-primed' write(11,3001)((ttdspr(i,j),j=1,3),i=1,3) write(11,*)'td' write(11,3001)((ttds(i,j),j=1,3),i=1,3) write(11,*)'tu-primed' write(11,3001)((ttuspr(i,j),j=1,3),i=1,3) write(11,*)'tu' write(11,3001)((ttus(i,j),j=1,3),i=1,3) write(11,*)'rd-primed' write(11,3001)((trdspr(i,j),j=1,3),i=1,3) write(11,*)'rd' write(11,3001)((trds(i,j),j=1,3),i=1,3) accru=cdabs(tstru1)/tstru2 acctd=cdabs(tsttd1)/tsttd2 acctu=cdabs(tsttu1)/tsttu2 accrd=cdabs(tstrd1)/tstrd2 write(11,*) write(11,*)'accuracy measures:' write(11,*)'acc-ru =',accru write(11,*)'acc-td =',acctd write(11,*)'acc-tu =',acctu write(11,*)'acc-rd =',accrd endif return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine frees(evec,tds,rds,det1,det2) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. c Computes the determinants associated with the free--surface c calculation (see theory notes section (5)). Not a required c subroutine. implicit real*8(a-h,o-z) complex*16 evec(6,6) complex*16 tds(3,3),rds(3,3) complex*16 wrk1(3,3),det1,det2 do 100 i=1,3 do 200 j=1,3 wrk1(i,j)=(0.d0,0.d0) do 300 k=1,3 wrk1(i,j)=wrk1(i,j)+evec(i+3,k+3)*rds(k,j) 300 continue wrk1(i,j)=evec(i+3,j)+wrk1(i,j) 200 continue 100 continue det1=wrk1(1,1)*(wrk1(2,2)*wrk1(3,3)-wrk1(2,3)*wrk1(3,2)) . -wrk1(1,2)*(wrk1(2,1)*wrk1(3,3)-wrk1(2,3)*wrk1(3,1)) . +wrk1(1,3)*(wrk1(2,1)*wrk1(3,2)-wrk1(2,2)*wrk1(3,1)) det2=tds(1,1)*(tds(2,2)*tds(3,3)-tds(2,3)*tds(3,2)) . -tds(1,2)*(tds(2,1)*tds(3,3)-tds(2,3)*tds(3,1)) . +tds(1,3)*(tds(2,1)*tds(3,2)-tds(2,2)*tds(3,1)) return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine write1(ilay,isoflg,px,py,eval,evec) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. implicit real*8(a-h,o-z) complex*16 px,py complex*16 eval(6),evec(6,6) write(11,*) write(11,*)'---------------------------------------------' write(11,*) write(11,*) write(11,*)'layer ',ilay,' px=',px,' py=',py if(isoflg.eq.0)write(11,*)'anisotropic layer' if(isoflg.eq.1)write(11,*)'isotropic layer' write(11,*) write(11,*)'eigenvalues or vertical slownesses...' write(11,*)eval write(11,*)'and eigenvectors...' do 1000 j=1,6 write(11,*)(evec(i,j),i=1,6) write(11,*) 1000 continue write(11,*) write(11,*)'---------------------------------------------' write(11,*) return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine write2(ilay,evec,eveci) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. implicit real*8 (a-h,o-z) complex*16 evec(6,6),eveci(6,6),testm(6,6) call xmult6(evec,eveci,testm) write(11,*) write(11,*)'---------------------------------------------' write(11,*) write(11,*) write(11,*)'test of inverse e-vector matrix, layer ',ilay do 55 i=1,6 write(11,3000)(testm(i,j),j=1,6) 55 continue write(11,*)'and in F format...' do 56 i=1,6 write(11,3001)(testm(i,j),j=1,6) 56 continue call xmult6(eveci,evec,testm) write(11,*) write(11,*)'second test of inverse...' do 57 i=1,6 write(11,3000)(testm(i,j),j=1,6) 57 continue write(11,*)'and in F format...' do 58 i=1,6 write(11,3001)(testm(i,j),j=1,6) 58 continue 3000 format(6(' (',d10.4,',',d10.4,')')) 3001 format(6(' (',f5.2,',',f5.2,')')) write(11,*) write(11,*)'---------------------------------------------' write(11,*) return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine write3(ilay,tu,ru,td,rd,evec1,evec2) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. implicit real*8(a-h,o-z) complex*16 tu(3,3),ru(3,3),td(3,3),rd(3,3) complex*16 testm1(6,6),testm2(6,6),evec1(6,6),evec2(6,6) write(11,*) write(11,*)'---------------------------------------------' write(11,*) write(11,*) write(11,*)'r/t matrices ... interface ',ilay write(11,*)'tu' do 1021 i=1,3 write(11,3100)(tu(i,j),j=1,3) 1021 continue write(11,*)'ru' do 1022 i=1,3 write(11,3100)(ru(i,j),j=1,3) 1022 continue write(11,*)'td' do 1023 i=1,3 write(11,3100)(td(i,j),j=1,3) 1023 continue write(11,*)'rd' do 1024 i=1,3 write(11,3100)(rd(i,j),j=1,3) 1024 continue 3100 format(3(' (',d9.3,',',d9.3,')')) write(11,*) write(11,*)'testing b.c.s ... interface',ilay do 3201 i=1,6 do 3202 j=1,6 testm1(i,j)=(0.d0,0.d0) 3202 continue 3201 continue do 3200 i=1,3 do 3205 j=1,3 do 3210 k=1,3 testm1(i,j)=testm1(i,j)+evec1(i,k)*ru(k,j) testm1(i,j+3)=testm1(i,j+3)+evec1(i,k)*td(k,j) testm1(i+3,j)=testm1(i+3,j)+evec1(i+3,k)*ru(k,j) testm1(i+3,j+3)=testm1(i+3,j+3)+evec1(i+3,k)*td(k,j) 3210 continue testm1(i,j)=testm1(i,j) + evec1(i,j+3) testm1(i+3,j)=testm1(i+3,j) + evec1(i+3,j+3) 3205 continue 3200 continue write(11,*)testm1(1,1),testm1(2,2),testm1(6,6) write(11,*)'lower stress-displacement matrix...' do 3215 i=1,6 write(11,3000)(testm1(i,j),j=1,6) 3215 continue 3000 format(6(' (',d10.4,',',d10.4,')')) do 3203 i=1,6 do 3204 j=1,6 testm2(i,j)=(0.d0,0.d0) 3204 continue 3203 continue do 3300 i=1,3 do 3305 j=1,3 do 3310 k=1,3 testm2(i,j)=testm2(i,j)+evec2(i,k+3)*tu(k,j) testm2(i,j+3)=testm2(i,j+3)+evec2(i,k+3)*rd(k,j) testm2(i+3,j)=testm2(i+3,j)+evec2(i+3,k+3)*tu(k,j) testm2(i+3,j+3)=testm2(i+3,j+3)+evec2(i+3,k+3)*rd(k,j) 3310 continue testm2(i,j+3)=testm2(i,j+3) + evec2(i,j) testm2(i+3,j+3)=testm2(i+3,j+3) + evec2(i+3,j) 3305 continue 3300 continue write(11,*)'upper stress-displacement matrix...' write(11,*)testm2(1,1),testm2(2,2),testm2(6,6) do 3315 i=1,6 write(11,3000)(testm2(i,j),j=1,6) 3315 continue bcacc1=0.d0 bcacc2=0.d0 do 4000 i=1,6 do 4001 j=1,6 bcacc1=bcacc1+cdabs(testm1(i,j)-testm2(i,j)) bcacc2=bcacc2+cdabs(testm1(i,j))+cdabs(testm2(i,j)) 4001 continue 4000 continue bcacc=bcacc1/bcacc2 write(11,*) write(11,*)'overall bc accuracy measure...diff, sum, bcacc =' write(11,*)bcacc1,bcacc2,bcacc write(11,*) write(11,*)'---------------------------------------------' write(11,*) return end c c........|.........|.........|.........|.........|.........|.........|.. c subroutine write4(tus,rus,tds,rds) c Copyright (c) 1996 C. J. Thomson. c All rights reserved by the author. complex*16 tus(3,3),rus(3,3),tds(3,3),rds(3,3) write(11,*) write(11,*)'---------------------------------------------' write(11,*) write(11,*)'final r/t matrices ...' write(11,*)'tus' do 9021 i=1,3 write(11,3100)(tus(i,j),j=1,3) 9021 continue write(11,*)'rus' do 9022 i=1,3 write(11,3100)(rus(i,j),j=1,3) 9022 continue write(11,*)'tds' do 9023 i=1,3 write(11,3100)(tds(i,j),j=1,3) 9023 continue write(11,*)'rds' do 9024 i=1,3 write(11,3100)(rds(i,j),j=1,3) 9024 continue 3100 format(1(' (',d15.9,',',d15.9,')')) return end c c........|.........|.........|.........|.........|.........|.........|.. c C