SUBROUTINE DISPL(P,U,I1) C C ROUTINE FOR THE COMPUTATION OF DISPLACEMENT VECTOR OF GENERATED C WAVES AT AN INTERFACE. IT WORKS EVEN FOR A COMPLEX REFRACTION VECTOR. C ROUTINES CHRM AND POLAR WORK ONLY FOR REAL REFRACTION VECTORS C C INPUT PARAMETERS: C P(1-3)... SLOWNESS VECTOR C I1... TYPE OF THE WAVE WITH THE SLOWNESS VECTOR P C I1=1... REFLECTED S WAVE (SV COMPONENT IN ISOTROPIC CASE) C I1=2... REFLECTED S WAVE (SH COMPONENT IN ISOTROPIC CASE) C I1=3... REFLECTED P WAVE C I1=4... TRANSMITTED S WAVE (SV COMPONENT IN ISOTROPIC CASE) C I1=5... TRANSMITTED S WAVE (SH COMPONENT IN ISOTROPIC CASE) C I1=6... TRANSMITTED P WAVE C C OUTPUT PARAMETERS: C U(1-3)... CORRESPONDING DISPLACEMENT VECTOR C complex uxn(3,3) COMPLEX U,U1,U2,U3,P,P1,P2,P3,P2P3,P1P2,P1P3,P1P1,P2P2,P3P3,C11, 1 C12,C13,C22,C23,C33,C11N,C22N,C33N,C23SQ,C13SQ,C12SQ,CD11,CD12, 2 CD13,CD22,CD23,CD33,CZ01,CZ02,CZ03,CDTR,cp1,cp2,cp3,cab DIMENSION U(3),P(3),VSQ(3),Y(6),UN(3),xx(3) COMMON /APROX/ A11,A12,A13,A14,A15,A16,A22,A23,A24,A25,A26,A33, 1 A34,A35,A36,A44,A45,A46,A55,A56,A66, 1 DXA11,DXA12,DXA13,DXA14,DXA15,DXA16,DXA22,DXA23, 1 DXA24,DXA25,DXA26,DXA33,DXA34,DXA35,DXA36,DXA44, 1 DXA45,DXA46,DXA55,DXA56,DXA66, 1 DYA11,DYA12,DYA13,DYA14,DYA15,DYA16,DYA22,DYA23, 1 DYA24,DYA25,DYA26,DYA33,DYA34,DYA35,DYA36,DYA44, 1 DYA45,DYA46,DYA55,DYA56,DYA66, 1 DZA11,DZA12,DZA13,DZA14,DZA15,DZA16,DZA22,DZA23, 1 DZA24,DZA25,DZA26,DZA33,DZA34,DZA35,DZA36,DZA44, 1 DZA45,DZA46,DZA55,DZA56,DZA66, 1 A2546,A1266,A1355,A1456,A3645,A2344 COMMON /AUXI/ IANI(20),INTR,INT1,IOUT,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOUT, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori COMPLEX PS COMMON /RAY/ AY(28,400),DS(20,20),KINT(20),HHH(3,3),tmax, 1 PS(3,7,20),IS(8,20),DINC,DTOLD,N,IREF,IND,IND1 C IRF=IREF IF(I1.LE.3)THEN MC=IS(5,IRF) ELSE MC=IS(6,IRF) END IF ISW=0 C C THIS PART CORRESPONDS TO CHRM ,BUT FOR COMPLEX REFRACTION VECTORS C 900 P1=P(1) P2=P(2) P3=P(3) 910 continue IF(ICOEF.GT.0)WRITE(LOUT,'(A,3(/2F12.6))')'DISPL: START P', 1p1,p2,p3 P2P3=P2*P3 P1P2=P1*P2 P1P3=P1*P3 P1P1=P1*P1 P2P2=P2*P2 P3P3=P3*P3 C11=P1P1*A11+P2P2*A66+P3P3*A55 1+2.*P2P3*A56+2.*P1P3*A15+2.*P1P2*A16 C22=P1P1*A66+P2P2*A22+P3P3*A44 1+2.*P2P3*A24+2.*P1P3*A46+2.*P1P2*A26 C33=P1P1*A55+P2P2*A44+P3P3*A33 1+2.*P2P3*A34+2.*P1P3*A35+2.*P1P2*A45 C23=P1P1*A56+P2P2*A24+P3P3*A34 1 +P2P3*A2344+P1P3*A3645+P1P2*A2546 C13=P1P1*A15+P2P2*A46+P3P3*A35 1 +P2P3*A3645+P1P3*A1355+P1P2*A1456 C12=P1P1*A16+P2P2*A26+P3P3*A45 1 +P2P3*A2546+P1P3*A1456+P1P2*A1266 C11N=C11-1. C22N=C22-1. C33N=C33-1. C23SQ=C23*C23 C13SQ=C13*C13 C12SQ=C12*C12 CD11=C22N*C33N-C23SQ CD22=C11N*C33N-C13SQ CD33=C11N*C22N-C12SQ CD12=C13*C23-C12*C33N CD13=C12*C23-C13*C22N CD23=C12*C13-C23*C11N CDTR=CD11+CD22+CD33 IF(ICOEF.GT.0) 1WRITE(LOUT,'(A,4(/4F12.6))')'DISPL: CDIJ' 1,CD11,CD12,CD13,CD22,CD23,CD33,CDTR if(abs(cdtr).lt..0000001)write(lout,'(A)')'DISPL: SHEAR WAVE 1 SINGULARITY IN CALCULATION OF R/T COEFFICIENT' C C THIS PART CORRESPONDS TO ROUTINE POLAR BUT FOR COMPLEX POLARIZATION C VECTORS C U1=CD11 U2=CD12 U3=CD13 uxn(1,1)=u1 uxn(1,2)=u2 uxn(1,3)=u3 XN=real(U1*CONJG(U1)+U2*CONJG(U2)+U3*CONJG(U3)) xx(1)=xn IF(ICOEF.GT.0)WRITE(LOUT,'(A,F14.7)')'DISPL: XN',XN U1=CD12 U2=CD22 U3=CD23 uxn(2,1)=u1 uxn(2,2)=u2 uxn(2,3)=u3 XN=real(U1*CONJG(U1)+U2*CONJG(U2)+U3*CONJG(U3)) xx(2)=xn IF(ICOEF.GT.0)WRITE(LOUT,'(A,F14.7)')'DISPL: XN',XN U1=CD13 U2=CD23 U3=CD33 uxn(3,1)=u1 uxn(3,2)=u2 uxn(3,3)=u3 XN=real(U1*CONJG(U1)+U2*CONJG(U2)+U3*CONJG(U3)) xx(3)=xn IF(ICOEF.GT.0)WRITE(LOUT,'(A,F14.7)')'DISPL: XN',XN xn=0. do 920 l=1,3 if(xn.ge.xx(l))go to 920 xn=xx(l) nx=l 920 continue xn=1./sqrt(xn) u1=uxn(nx,1)*xn u2=uxn(nx,2)*xn u3=uxn(nx,3)*xn U(1)=U1 U(2)=U2 U(3)=U3 IF(ICOEF.GT.0)WRITE(LOUT,'(A,3(/2F12.6))')'DISPL: U',U C C CHECK OF PRECISION C CZ01=C11N*u1+C12*u2+C13*u3 CZ02=C12*u1+C22N*u2+C23*u3 CZ03=C13*u1+C23*u2+C33N*u3 IF(ICOEF.GT.0) 1WRITE(LOUT,'(A,3(/2F14.7))')'DISPL: PRECISSION OF DISPL' 1,CZ01,CZ02,CZ03 IF(ISW.GT.0)THEN GOTO 1100 END IF RETURN C C COMPUTATION OF DISPLACEMENT VECTORS AT AN INTERFACE IF A SHEAR WAVE C SINGULARITY OCCURS. C FIRST: THE DISPLACEMENT VECTOR OF THE QP-WAVE IS COMPUTED FOR THE C NORMAL DIRECTION OF THE S-WAVE C SECOND: TWO DISPLACEMENT VECTORS ARE CONSTRUCTED SO THAT THEY ARE C PERPENDICULAR TO THE QP-WAVE DISPLACEMENT VECTOR C (THEY CORRESPOND TO SV AND SH COMPONENT) C 1000 ISW=1 do 950 k=1,3 Y(k)=AY(k+1,N) UN(k)=DS(k,IRF) 950 continue Y(4)=REAL(P1) Y(5)=REAL(P2) Y(6)=REAL(P3) IF(ICOEF.GT.0) 1WRITE(LOUT,'(A,3I5,6F12.6)')'DISPL: I1,IRF,MC,Y',I1,IRF,MC,Y UUN=Y(4)*Y(4)+Y(5)*Y(5)+Y(6)*Y(6) UUN=1./SQRT(UUN) Y(4)=Y(4)*UUN Y(5)=Y(5)*UUN Y(6)=Y(6)*UUN CALL INIT(Y,VSQ) V=AMIN1(VSQ(1),VSQ(2),VSQ(3)) V=1./SQRT(V) IF(ICOEF.GT.0)WRITE(LOUT,'(A,3F12.4)')'DISPL: VSQ(1-3)=',VSQ IF(ICOEF.GT.0)WRITE(LOUT,'(A,F12.4)')'DISPL: V=',V p1=y(4)/v p2=y(5)/v p3=y(6)/v go to 910 1100 if(i1.eq.1.or.i1.eq.4)then cp1=U1*U3 cp2=U2*U3 cp3=-(U1*U1+U2*U2) cab=sqrt(real(cp1*conjg(cp1)+cp2*conjg(cp2)+cp3*conjg(cp3))) U(1)=CP1/cab U(2)=CP2/cab U(3)=cp3/cab IF(ICOEF.GT.0)WRITE(LOUT,'(A,7(/2F12.6))')'DISPL: U,P',U,P RETURN end if if(i1.eq.2.or.i1.eq.5)then cp1=-U2 cp2=U1 cp3=cmplx(0.,0.) cab=sqrt(real(CP1*conjg(CP1)+CP2*conjg(CP2)+CP3*conjg(CP3))) U(1)=CP1/cab U(2)=CP2/cab U(3)=cmplx(0.,0.) IF(ICOEF.GT.0)WRITE(LOUT,'(A,7(/2F12.6))')' DISPL U,P',U,P RETURN end if END