C
C Subroutine file 'wan.for' to compute quantities along a ray necessary C for computation of the Green function by means of coupling ray-theory C in weakly anisotropic models with interfaces, and to calculate second- C order travel time perturbations in models without interfaces. C C Version: 6.00 C Date: 2006, June 5 C C Coded by: Petr Bulant C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C http://sw3d.cz/staff/bulant.htm C C======================================================================= C SUBROUTINE WAN(LU1,LU2,LU3,GREEN,MGREEN,RPTS,IPTS,MAXPTS,NQ, * ERR,LUQI) C C---------------------------------------------------------------------- INTEGER LU1,LU2,LU3,LUQI INTEGER MGREEN,MAXPTS,NQ INTEGER NQPTS PARAMETER (NQPTS=61) REAL GREEN(MGREEN),RPTS(NQPTS,MAXPTS/NQPTS),ERR(2) INTEGER IPTS(NQPTS,MAXPTS/NQPTS) C ------------------------------------------------------------------ C Input: C LU1 ... Number of the logical unit connected to the CRT output C file with quantities along rays. C LU2 ... Number of the logical unit connected to the CRT output C file with quantities at the storing surface. C LU3 ... Number of the logical unit connected to the CRT output C file with quantities at the initial points of rays. C MGREEN..Dimension of an output array GREEN. C MAXPTS..Maximum number of records in arrays RPTS and IPTS. C ERR... Output of previous invocations of WAN, zero for the first C invocation. C LUQI... Number of the logical unit connected to the output file C QILST with the second-order perturbations of anisotropic C travel times of S waves. C Zero if the file is not to be generated. C C Output: C GREEN...Array containing the Green function for the given ray C and for all the frequencies: C GREEN(1)... Travel time between receiver and source. C GREEN(2)... Imaginary part of the complex-valued travel time C between receiver and source due to attenuation. C GREEN(3:8)... Coordinates of the receiver and coordinates C of the source. C GREEN(9:14)... Derivatives of the travel time with respect C to the coordinates of the receiver and coordinates of the C source. C GREEN(15:(14+NF*18)) ... C for P wave once, for S wave NF times following 18 numbers, C specifying 1 000 000 times enlarged amplitude of the C Green function: contravariant components of the complex- C valued 3*3 matrix Gij in model coordinates, where the C first subscript corresponds to the receiver and the second C subscript corresponds to the source. The components are C ordered as C Re(G11),Im(G11),Re(G21),Im(G21),Re(G31),Im(G31), C Re(G12),Im(G12),Re(G22),Im(G22),Re(G32),Im(G32), C Re(G13),Im(G13),Re(G23),Im(G23),Re(G33),Im(G33). C NQ ... Number of records stored in the array GREEN. C IPTS,RPTS... Auxiliary arrays, redimensioned in each invocation. C ERR... Maximum of calculated errors of the quasi-isotropic C projection of the polarization vectors. C ERR(1): at the source, C ERR(2): at the receiver. C C....................................................................... C Subroutines and external functions required: EXTERNAL WAPTS,WAREAD,WACHRI,WAMAT,WACHAN,WAPROJ, *WAPERT,WASUM,WASUM3,WASUM4,WASUM5, *ERROR,RSEP1,RSEP3T,RSEP3R,RSEP3I,FORM2,LENGTH,EIGEN, *HIVD2,BLOCK,PARM3,AP00,AP21,TRANS INTEGER LENGTH REAL WACHAN C WAPTS,WAREAD,WACHRI,WAMAT,WACHAN,WAPROJ,WAPERT,WASUM,WASUM3, C WASUM4,WASUM5 ... This file. C ERROR ... File C error.for. C RSEP1,RSEP3I,RSEP3T,RSEP3R ... File C sep.for. C FORM2 ... File C forms.for. C LENGTH... File C length.for. C EIGEN ... File C eigen.for. C HIVD2 ... File C means.for. C BLOCK ... File C model.for. C PARM3 ... File C parm.for. C AP00,AP21 ... File C ap.for. C TRANS ... File C trans.for. C C Common block /POINTC/ to store the results of complete ray tracing: INCLUDE 'pointc.inc' C ........................... C Common block /RTMAT/ to store the matrices of reflection-transmittion C coefficients. INTEGER MRT,NRT PARAMETER (MRT=100) REAL PIRTR(2,2,MRT),PIRTI(2,2,MRT) COMMON/RTMAT/NRT,PIRTR,PIRTI SAVE /RTMAT/ C MRT ... Dimension of arrays. C NRT ... Number of stored R-T matrices. C PIRTR .. Real part of R-T matrices. C PIRTI .. Imaginary part of R-T matrices. C The matrices are expressed in terms of polarization vectors. C ........................... C Common block /WASEP/ to store the input SEP parameters used C in subroutine WAN and their subroutines. REAL CRTANI,DEGREE COMMON/WASEP/CRTANI,DEGREE SAVE /WASEP/ C CRTANI ... Switch between isotropic and anisotropic ray tracing. C DEGREE ... Degree of the homogeneous Hamiltonian. C....................................................................... C C Auxiliary storage locations: INTEGER MPTS REAL PIGRA(2,2),PIGIA(2,2),PIGR(2,2),PIGI(2,2) REAL AR11,AI11,AR21,AI21,AR12,AI12,AR22,AI22 REAL G11,G21,G12,G22 REAL PI PARAMETER (PI=3.1415926) REAL GREENA(32) REAL TTCOR,FREQ,GAMA,AUX0,AUX1,AUX2 INTEGER NPTS INTEGER I,I1,I2,I3,J,IRT INTEGER NFFT,NF,KQIPV,KQITT,KQIRAY REAL DT,FMIN,FMAX,DF,OF INTEGER IPTTMP,NYFTMP,ISRFFT REAL YI1TMP,YF1TMP,STEP,TISO,DTLIN,DTAU(7),QICOR(16) LOGICAL LPWAVE,LSWAVE CHARACTER*135 FORMQI CHARACTER*256 TXTERR DATA NFFT/0/ C C IPTS,RPTS... Quantities in the points on the ray: C IPTS(1,I)... Index of the I-th point, zero for points C added to the ray by interpolation. C RPTS(2,I)... Travel time in I-th point. C RPTS(3-5,I).. Coordinates of the point. C RPTS(6-8,I).. Slowness vector in the point. C RPTS(9-11,I) Polarization vector in the point. C IPTS(12,I)... Index of complex block. C RPTS(13,I)... Eigenvalue of Christoffel matrix in I-th point C corresponding to qP wave. C RPTS(14,I)... Eigenvalue of Christoffel matrix in I-th point C corresponding to qS1 (faster) wave. C RPTS(15,I)... Eigenvalue of Christoffel matrix in I-th point C corresponding to qS2 (slower) wave. C RPTS(16-19,I) Eigenvectors projected to the plane defined by C polarization vectors. C RPTS(20,I)... Angular difference of eigenvectors of Christoffel C matrix corresponding to qS1 wave in I-th and C in (I-1)-st point. C RPTS(21,I)... Time integral of the difference between qS1 and C qS2 eigenvalues along the ray from the (I-1)-st C to the I-th point. C RPTS(22-24,I) Derivatives of the velocity. C RPTS(25-33,I) Eigenvectors of qS1, qS2 and P wave. C RPTS(34,I)... S-wave velocity. C RPTS(35-38,I) Matrix of geometrical spreading. C RPTS(39-42,I) Transformation matrix P. C IPTS(43,I)... Index of surface, zero inside a complex block. C RPTS(44-52,I) Covariant components of the basis vectors C of the ray-centred coordinate system. C RPTS(53-61,I) Contravariant components of the basis vectors C of the ray-centred coordinate system. C C----------------------------------------------------------------------- C IF (NFFT.EQ.0) THEN CALL RSEP3I('NFFT',NFFT,1) CALL RSEP3R('DT ',DT ,1.) CALL RSEP3R('FMIN',FMIN,0.) CALL RSEP3R('FMAX',FMAX,0.5/DT) CALL RSEP3R('DF ',DF ,1./(FLOAT(NFFT)*DT)) CALL RSEP3R('OF ',OF ,DF*NINT(FMIN/DF)) CALL RSEP3I('NF ',NF ,NINT((FMAX-OF)/DF)+1) CALL RSEP3I('QIPV',KQIPV,0) CALL RSEP3I('QITT',KQITT,0) FORMQI(1:7)='(1I6,A,' CALL RSEP3I('QIRAY',KQIRAY,0) CALL RSEP3R('CRTANI',CRTANI, 0.) CALL RSEP3R('DEGREE',DEGREE,-1.) ENDIF C C Reading in the quantities stored in individual points on C a ray, computing all the auxiliary quantities in the points, C and, if necessary, adding new points on the ray by interpolation: MPTS=MAXPTS/NQPTS NPTS=0 CALL WAPTS(LU1,LU2,LU3,NF,OF,DF,RPTS,IPTS,NQPTS,MPTS,NPTS) IF (NPTS.EQ.0) THEN C End of rays: NQ=0 RETURN ENDIF C TISO=Y(1) C C Computing the values of travel time corrections along the ray: LPWAVE=.TRUE. LSWAVE=.TRUE. TTCOR=0. DO 30, I1=1,NPTS-1 I2=I1+1 IF (IPTS(12,I2).GT.0) THEN C P wave: TTCOR=TTCOR+(1./SQRT(RPTS(13,I2))+1./SQRT(RPTS(13,I1)))*0.5* * (RPTS(2,I2)-RPTS(2,I1)) LSWAVE=.FALSE. ELSE C S wave: IF(KQITT.LE.0) THEN TTCOR=TTCOR+(0.25/SQRT(RPTS(14,I2))+0.25/SQRT(RPTS(15,I2)) * +0.25/SQRT(RPTS(14,I1))+0.25/SQRT(RPTS(15,I1))) * *(RPTS(2,I2)-RPTS(2,I1)) ELSE TTCOR=TTCOR+(1.5-0.125*(RPTS(14,I2)+RPTS(15,I2) * +RPTS(14,I1)+RPTS(15,I1))) * *(RPTS(2,I2)-RPTS(2,I1)) END IF LPWAVE=.FALSE. ENDIF 30 CONTINUE C IF (LPWAVE) THEN C P wave along the whole ray: Y(1)=TTCOR+YI(1) CALL AP21(KQIPV,RPTS(25,1),RPTS(25,NPTS),GREEN) NQ=32 GOTO 91 ENDIF C IF (LSWAVE.AND.((LUQI.NE.0).OR.(KQIRAY.NE.0))) THEN C S wave along the whole ray, C calculating second-order travel-time perturbations: DTAU(1)=0. DTAU(3)=0. IPTTMP=IPT YI1TMP=YI(1) NYFTMP=NYF ISRFFT=ISRFF YF1TMP=YF(1) YI(1) =RPTS(2,1) NYF =1 DTLIN=0. DO 95, I1=1,NPTS IF (IPTS(1,I1).EQ.2) THEN STEP=RPTS(2,I1)-RPTS(2,1) ENDIF IF (IPTS(1,I1).EQ.0) THEN STEP=0. GOTO 96 ENDIF 95 CONTINUE 96 CONTINUE DO 100, I1=1,NPTS DTLIN=DTLIN+RPTS(21,I1) IPT =IPTS(1,I1)-1 ISRFF =IPTS(43,I1) YF(1) =RPTS(2,I1) IF (KQIRAY.NE.0) THEN RPTS(21,I1)=RPTS(21,I1)+DTAU(1)-DTAU(3) ENDIF CALL WAPERT(RPTS(1,I1),IPTS(1,I1),NQPTS,STEP,DTAU) IF (KQIRAY.NE.0) THEN RPTS(21,I1)=RPTS(21,I1)+DTAU(3)-DTAU(1) ENDIF 100 CONTINUE C IF (LUQI.NE.0) THEN DTLIN=DTLIN/2. QICOR(1)=TISO QICOR(2)=TTCOR-DTLIN-TISO QICOR(3)=TTCOR+DTLIN-TISO QICOR(4)=DTAU(1) QICOR(5)=DTAU(3) QICOR(6)=TTCOR-DTLIN+DTAU(1) QICOR(7)=QICOR(6) QICOR(8)=TTCOR+DTLIN+DTAU(3) QICOR(9)=QICOR(8) IF(CRTANI.EQ.0.) THEN QICOR(10)=TTCOR-TISO QICOR(11)=0.25*(DTAU(1)+DTAU(2)+DTAU(3)) QICOR(12)=TISO+QICOR(10)+QICOR(11) QICOR(13)=QICOR(12) QICOR(14)=0.25*(DTAU(1)-DTAU(2)+DTAU(3)) QICOR(15)=QICOR(14) QICOR(16)=DTAU(1)-DTAU(2)+DTAU(3) ELSE QICOR(10)=DTAU(7)-TISO QICOR(11)=DTAU(6) QICOR(12)=DTAU(7)+DTAU(6) QICOR(13)=QICOR(12) QICOR(14)=DTAU(1)-DTAU(4)+DTAU(6) QICOR(15)=DTAU(3)-DTAU(5)+DTAU(6) QICOR(16)=DTAU(1)-DTAU(2)+DTAU(3) END IF CALL FORM2(16,QICOR,QICOR,FORMQI(8:135)) WRITE(LUQI,FORMQI) IREC,(' ',QICOR(I),I=1,16),' /' ENDIF C IPT =IPTTMP YI(1) =YI1TMP NYF =NYFTMP ISRFF =ISRFFT YF(1) =YF1TMP C IF (KQIRAY.NE.0) THEN TTCOR=TTCOR+0.5*(DTAU(3)+DTAU(1)) ENDIF ELSEIF ((.NOT.LSWAVE).AND.((LUQI.NE.0).OR.(KQIRAY.NE.0))) THEN C WAN-20 WRITE(TXTERR,'(2A,1I6,A)') * 'WAN-20: Second-order perturbations of anisotropic travel ', * 'times may not be calculated for the ray ',IRAY, * ', the ray crossed an interface.' CALL WARN(TXTERR(1:LENGTH(TXTERR))) C The perturbations may be calculated only along the rays C which do not cross any interface. ENDIF C C S wave in some part of the ray: C Loop over the frequencies: DO 90, I2=0,NF-1 FREQ=OF+I2*DF PIGR(1,1)=1. PIGR(2,1)=0. PIGR(1,2)=0. PIGR(2,2)=1. PIGI(1,1)=0. PIGI(2,1)=0. PIGI(1,2)=0. PIGI(2,2)=0. IRT=0 C Computing the propagator matrix PiGe along the ray: C Loop along points on the ray: DO 40, I1=2,NPTS I3=I1-1 IF (RPTS(2,I1).EQ.RPTS(2,I3)) THEN C The ray crossed an interface: IRT=IRT+1 IF (IRT.GT.NRT) THEN C WAN-17 CALL ERROR('WAN-17: Disorder in number of R-T.') C This error should not appear. Number of reflections- C -transmittions should be the same in WAREAD and here. ENDIF C Transforming R-T matrix from polarization vectors to C eigenvectors: G11=RPTS(16,I3) G21=RPTS(17,I3) G12=RPTS(18,I3) G22=RPTS(19,I3) AR11=PIRTR(1,1,IRT)*G11+PIRTR(1,2,IRT)*G21 AI11=PIRTI(1,1,IRT)*G11+PIRTI(1,2,IRT)*G21 AR21=PIRTR(2,1,IRT)*G11+PIRTR(2,2,IRT)*G21 AI21=PIRTI(2,1,IRT)*G11+PIRTI(2,2,IRT)*G21 AR12=PIRTR(1,1,IRT)*G12+PIRTR(1,2,IRT)*G22 AI12=PIRTI(1,1,IRT)*G12+PIRTI(1,2,IRT)*G22 AR22=PIRTR(2,1,IRT)*G12+PIRTR(2,2,IRT)*G22 AI22=PIRTI(2,1,IRT)*G12+PIRTI(2,2,IRT)*G22 G11=RPTS(16,I1) G21=RPTS(18,I1) G12=RPTS(17,I1) G22=RPTS(19,I1) PIGRA(1,1)=G11*AR11+G12*AR21 PIGIA(1,1)=G11*AI11+G12*AI21 PIGRA(2,1)=G21*AR11+G22*AR21 PIGIA(2,1)=G21*AI11+G22*AI21 PIGRA(1,2)=G11*AR12+G12*AR22 PIGIA(1,2)=G11*AI12+G12*AI22 PIGRA(2,2)=G21*AR12+G22*AR22 PIGIA(2,2)=G21*AI12+G22*AI22 ELSE C Element of the ray inside a complex block: IF (IPTS(12,I1).LT.0) THEN C S wave: GAMA=PI*FREQ*RPTS(21,I1) AUX0=SQRT(RPTS(20,I1)**2 + GAMA**2) AUX1=COS(AUX0) IF (AUX0.EQ.0.) THEN AUX2=1. ELSE AUX2=SIN(AUX0)/AUX0 ENDIF C Matrix for this step along the ray: PIGRA(1,1)= AUX1 PIGRA(2,1)=-RPTS(20,I1)*AUX2 PIGRA(1,2)=-PIGRA(2,1) PIGRA(2,2)= AUX1 PIGIA(1,1)=-GAMA*AUX2 PIGIA(2,1)= 0. PIGIA(1,2)= 0. PIGIA(2,2)=-PIGIA(1,1) ELSE C P wave: GOTO 39 ENDIF ENDIF C Matrix for all the steps along the ray to current point: CALL WAMAT(PIGRA,PIGIA,PIGR,PIGI) 39 CONTINUE 40 CONTINUE C C Computing the Green function: C Transforming the propagator matrix from eigenvectors C to polarization vectors: Y(1)=TTCOR+YI(1) IF(KQIPV.EQ.0) THEN Y(28)=PIGR(1,1) Y(29)=PIGI(1,1) Y(30)=PIGR(2,1) Y(31)=PIGI(2,1) Y(32)=PIGR(1,2) Y(33)=PIGI(1,2) Y(34)=PIGR(2,2) Y(35)=PIGI(2,2) ELSE G11=RPTS(16,1) G12=RPTS(17,1) G21=RPTS(18,1) G22=RPTS(19,1) AR11=PIGR(1,1)*G11+PIGR(1,2)*G21 AI11=PIGI(1,1)*G11+PIGI(1,2)*G21 AR21=PIGR(2,1)*G11+PIGR(2,2)*G21 AI21=PIGI(2,1)*G11+PIGI(2,2)*G21 AR12=PIGR(1,1)*G12+PIGR(1,2)*G22 AI12=PIGI(1,1)*G12+PIGI(1,2)*G22 AR22=PIGR(2,1)*G12+PIGR(2,2)*G22 AI22=PIGI(2,1)*G12+PIGI(2,2)*G22 G11=RPTS(16,NPTS) G21=RPTS(17,NPTS) G12=RPTS(18,NPTS) G22=RPTS(19,NPTS) Y(28)=G11*AR11+G12*AR21 Y(29)=G11*AI11+G12*AI21 Y(30)=G21*AR11+G22*AR21 Y(31)=G21*AI11+G22*AI21 Y(32)=G11*AR12+G12*AR22 Y(33)=G11*AI12+G12*AI22 Y(34)=G21*AR12+G22*AR22 Y(35)=G21*AI12+G22*AI22 ENDIF C CALL AP21(KQIPV,RPTS(25,1),RPTS(25,NPTS),GREENA) IF (I2.EQ.0) THEN DO 50 I=1,14 GREEN(I)=GREENA(I) 50 CONTINUE NQ=14 ENDIF J=I2*18 DO 60 I=15,32 GREEN(J+I)=GREENA(I) 60 CONTINUE NQ=NQ+18 C 90 CONTINUE 91 CONTINUE C IF (KQIPV.EQ.1) THEN C Calculating the error of the QI projection of polariz. vectors G11=RPTS(16,1) G21=RPTS(17,1) G12=RPTS(18,1) G22=RPTS(19,1) ERR(1)=AMAX1(SQRT(2.*ABS(1.-ABS(G11*G22-G12*G21))),ERR(1)) G11=RPTS(16,NPTS) G21=RPTS(17,NPTS) G12=RPTS(18,NPTS) G22=RPTS(19,NPTS) ERR(2)=AMAX1(SQRT(2.*ABS(1.-ABS(G11*G22-G12*G21))),ERR(2)) ENDIF C C RETURN END C C======================================================================= C SUBROUTINE WAPTS(LU1,LU2,LU3,NF,OF,DF,RPTS,IPTS,NQPTS,MPTS,NPTS) C C----------------------------------------------------------------------- C Subroutine to read in the quantities stored in individual points C on the ray, to compute all the auxiliary quantities in the points, C and, if necessary, to add new points on the rays by interpolation. C Reading the files with points on the rays is done by a simple C invocation of subroutine WAREAD. C Computation of elastic parameters is completed by invocation C of subroutine PARM3 of file 'parm.for'. C Then Christoffel matrix is evaluated and its eigenvalues and C eigenvectors are computed by subroutine EIGEN of file 'eigen.for'. C In the next step the angular difference DELTFI is computed for C each subinterval along the ray. C If the value of DELTFI is greater than prescribed limit, new points C are added on the ray. Coordinates, slowness vector and polarization C vector are interpolated by subroutine HIVD2 of the file 'means.for'. C Material parameters, Christoffel matrix and their eigenvalues and C eigenvectors are calculated in the new point. Matrix of geometrical C spreading, transformation matrix P and covariant and contravariant C components of the basis vectors of the RCCS are interpolated linearly. C INTEGER LU1,LU2,LU3,NF REAL OF,DF INTEGER NPTS,NQPTS,MPTS REAL RPTS(NQPTS,MPTS) INTEGER IPTS(NQPTS,MPTS) C Input: C LU1 ... Number of the logical unit connected to the CRT output C file with quantities along rays. C LU2 ... Number of the logical unit connected to the CRT output C file with quantities at the initial points of rays. C LU3 ... Number of the logical unit connected to the CRT output C file with quantities at the storing surface. C NF ... Number of frequencies. C OF ... First frequency. C DF ... Frequency increment. C NQPTS.. Dimension of arrays RPTS and IPTS. C MPTS... Dimension of arrays RPTS and IPTS. C Output: C RPTS,IPTS ... The arrays are filled with all the quantities for C single two-point ray during one invocation of this C subroutine: C IPTS(1,I) ... Index of the I-th point, zero for points C added to the ray by interpolation. C RPTS(2,I) ... Travel time in I-th point. C RPTS(3-5,I) . Coordinates of the point. C RPTS(6-8,I) . Slowness vector in the point. C RPTS(9-11,I) Polarization vector in the point. C IPTS(12,I) .. Index of complex block. C RPTS(13,I)... Eigenvalue of Christoffel matrix in I-th point C corresponding to qP wave. C RPTS(14,I)... Eigenvalue of Christoffel matrix in I-th point C corresponding to qS1 (faster) wave. C RPTS(15,I)... Eigenvalue of Christoffel matrix in I-th point C corresponding to qS2 (slower) wave. C RPTS(16-19,I) Eigenvectors projected to the plane defined by C polarization vectors. C RPTS(20,I) .. Angular difference of eigenvectors of Christoffel C matrix corresponding to qS1 wave in I-th and C in (I-1)-st point. C RPTS(21,I) .. Time integral of the difference between qS1 and C qS2 eigenvalues along the ray from the (I-1)-st C to the I-th point. C RPTS(22-24,I) Derivatives of the velocity. C RPTS(25-33,I) Eigenvectors of qS1, qS2 and P wave. C RPTS(34,I)... S-wave velocity. C RPTS(35-38,I) Matrix of geometrical spreading. C RPTS(39-42,I) Transformation matrix P. C IPTS(43,I)... Index of surface, zero inside a complex block. C RPTS(44-52,I) Covariant components of the basis vectors C of the ray-centred coordinate system. C RPTS(53-61,I) Contravariant components of the basis vectors C of the ray-centred coordinate system. C NPTS ... Number of points stored in RPTS (IPTS). C C External functions required: EXTERNAL WACHAN REAL WACHAN C....................................................................... C C Auxiliary storage locations: INTEGER NPTSE INTEGER KQIDT REAL UP(10),US(10),RHO,QP,QS,VP,VS,VD(10),QL,AA(10,21),QQ(21) REAL EE(9),DER(9) REAL ERRWAN,DEFI,DELTFI,PI PARAMETER (PI=3.1415926) INTEGER I1,I2 INTEGER MQANT,MNEWP,NNEWP PARAMETER (MNEWP=5) PARAMETER (MQANT=56) REAL ROLDP(MQANT),RNEWP(MQANT,MNEWP) INTEGER KNEWP(MNEWP) REAL TTOLD,XOLD(3),DXOLD(3),POLD(3),DPOLD(3),EOLD(3),DEOLD(3) REAL TTNEW,XNEW(3),DXNEW(3),PNEW(3),DPNEW(3),ENEW(3),DENEW(3),AUX C C ROLDP(I),RNEWP(I,J) C I=1 ... Travel time. C 2-4 ... Coordinates. C 5-7 ... Slowness vector. C 8-10 ... Polarization vector E1. C 11 ... qP eigenvalue of Christoffel matrix. C 12 ... qS1 eigenvalue of Christoffel matrix. C 13 ... qS2 eigenvalue of Christoffel matrix. C 14 ... First component of the qS1 eigenvector C projected to the plane perpendicular to the ray, C defined by two basis vectors of ray-centered C coordinate system. C 15 ... Second component of the qS1 eigenvector. C 16 ... First component of the qS2 eigenvector. C 17 ... Second component of the qS2 eigenvector. C 18-26... Eigenvectors of qS1, qS2 and P wave. C 27 ... S-wave velocity. C 28-30... Derivatives of S-wave velocity. C 31-34... Matrix of geometrical spreading. C 34-38... Transformation matrix P. C 39-56... Covariant and contravariant components of C the basis vectors of the RCCS. C KNEWP...Division index of points in RNEWP. C AA.. Values, first and second partial derivatives of real C parts of 21 reduced (divided by the density) elastic C parameters. The order of the value, first and second C partial derivatives of each parameter Aij is: C Aij,Aij1,Aij2,Aij3,Aij11,Aij12,Aij22,Aij13,Aij23,Aij33. C The order of parameters (second array index) is: C A11,A12,A22,A13,A23,A33,A14,A24,A34,A44,A15,A25,A35,A45, C A55,A16,A26,A36,A46,A56,A66. C RHO... Density at the given point. C QQ ... Imaginary parts of 21 reduced elastic parameters at the C given point, ordered as C Q11,Q12,Q22,Q13,Q23,Q33,Q14,Q24,Q34,Q44,Q15,Q25,Q35,Q45, C Q55,Q16,Q26,Q36,Q46,Q56,Q66. C EE ... Eigenvectors of the Christoffel matrix. C DER ... Derivatives dx/dt (anisotropic). C ERRWAN..Maximum error in propagator matrix along whole ray. C DEFI .. Maximum allowed angular change for eigenvectors of C Christoffel matrix between neighboring points on the ray. C TTOLD...Travel time in the old point stored for interpolation. C XOLD ...Coordinates of the old point stored for interpolation. C DXOLD ..Derivatives in the old point stored for interpolation. C POLD ...Slowness vector in the old point stored for interpolation. C DPOLD ..Derivatives of the slowness vector in the old point. C EOLD ...Polarization vector in the old point. C DEOLD ..Derivatives of the polarization vector in the old point. C TTNEW,XNEW,DXNEW,PNEW,DPNEW,ENEW,DENEW ... The same quantities C in the new point. Here old point and new point mean two C consecutive points in which the output of CRT was C recorded. C NPTSE...Number of points along the ray, where all the quantities C have been checked. C C----------------------------------------------------------------------- C CALL RSEP3R('ERRWAN',ERRWAN,0.001) CALL RSEP3I('QIDT',KQIDT,0) C C Reading the quantities computed by CRT for one ray: CALL WAREAD(LU1,LU2,LU3,RPTS,IPTS,NQPTS,MPTS,NPTS) IF (NPTS.EQ.0) C End of rays: * RETURN IF (NPTS.LT.2) THEN C WAN-01 CALL ERROR('WAN-01: A ray formed by single point') C This error should not appear. C Each ray should be represented by at least two points. ENDIF C C C Reading the material parameters in the initial point: CALL PARM3(IPTS(12,1),RPTS(3,1),AA,RHO,QQ) C Computing eigenvectors and eigenvalues of C the Christoffel matrix in the initial point: CALL WACHRI(RPTS(6,1),RPTS(7,1),RPTS(8,1), * AA(1,1),AA(1,2),AA(1,3),AA(1,4),AA(1,5),AA(1,6),AA(1,7), * AA(1,8),AA(1,9),AA(1,10),AA(1,11),AA(1,12),AA(1,13),AA(1,14), * AA(1,15),AA(1,16),AA(1,17),AA(1,18),AA(1,19),AA(1,20),AA(1,21), * RPTS(13,1),RPTS(14,1),RPTS(15,1),EE,DER) C Storing the eigenvectors of the Christoffel matrix RPTS(25,1)=EE(4) RPTS(26,1)=EE(5) RPTS(27,1)=EE(6) RPTS(28,1)=EE(7) RPTS(29,1)=EE(8) RPTS(30,1)=EE(9) RPTS(31,1)=EE(1) RPTS(32,1)=EE(2) RPTS(33,1)=EE(3) IF (ABS(RPTS(15,1)-RPTS(14,1)).LT.0.00001) THEN C Isotropic case, no projection of eigenvectors: RPTS(16,1)=1. RPTS(17,1)=0. RPTS(18,1)=0. RPTS(19,1)=1. ELSE C Projecting the qS eigenvectors to the plane C perpendicular to the ray: CALL WAPROJ(RPTS(6,1),RPTS(7,1),RPTS(8,1), * RPTS(9,1),RPTS(10,1),RPTS(11,1), * EE(4),EE(5),EE(6),EE(7),EE(8),EE(9), * RPTS(16,1),RPTS(17,1),RPTS(18,1),RPTS(19,1)) ENDIF RPTS(20,1)=0. RPTS(21,1)=0. C C Quantities for future possible interpolation: TTOLD=RPTS(2,1) C Coordinates: XOLD(1)=RPTS(3,1) XOLD(2)=RPTS(4,1) XOLD(3)=RPTS(5,1) C dx/dt=p*v**2=p/sum(pi pi): AUX=(RPTS(6,1)**2+RPTS(7,1)**2+RPTS(8,1)**2) DXOLD(1)=RPTS(6,1)/AUX DXOLD(2)=RPTS(7,1)/AUX DXOLD(3)=RPTS(8,1)/AUX C Slowness: POLD(1)=RPTS(6,1) POLD(2)=RPTS(7,1) POLD(3)=RPTS(8,1) C dpi/dt=-1/v*dv/dxi: AUX=SQRT(AUX) DPOLD(1)=-RPTS(22,1)/AUX DPOLD(2)=-RPTS(23,1)/AUX DPOLD(3)=-RPTS(24,1)/AUX C Polarization: EOLD(1)=RPTS( 9,1) EOLD(2)=RPTS(10,1) EOLD(3)=RPTS(11,1) C dei/dt=v*sum(dv/dxk ek)*pi AUX=AUX* *(RPTS(22,1)*RPTS(9,1)+RPTS(23,1)*RPTS(10,1)+RPTS(24,1)*RPTS(11,1)) DEOLD(1)=AUX*RPTS(6,1) DEOLD(2)=AUX*RPTS(7,1) DEOLD(3)=AUX*RPTS(8,1) C ROLDP( 1)=RPTS( 2,1) ROLDP( 2)=RPTS( 3,1) ROLDP( 3)=RPTS( 4,1) ROLDP( 4)=RPTS( 5,1) ROLDP( 5)=RPTS( 6,1) ROLDP( 6)=RPTS( 7,1) ROLDP( 7)=RPTS( 8,1) ROLDP( 8)=RPTS( 9,1) ROLDP( 9)=RPTS(10,1) ROLDP(10)=RPTS(11,1) ROLDP(11)=RPTS(13,1) ROLDP(12)=RPTS(14,1) ROLDP(13)=RPTS(15,1) ROLDP(14)=RPTS(16,1) ROLDP(15)=RPTS(17,1) ROLDP(16)=RPTS(18,1) ROLDP(17)=RPTS(19,1) ROLDP(18)=RPTS(25,1) ROLDP(19)=RPTS(26,1) ROLDP(20)=RPTS(27,1) ROLDP(21)=RPTS(28,1) ROLDP(22)=RPTS(29,1) ROLDP(23)=RPTS(30,1) ROLDP(24)=RPTS(31,1) ROLDP(25)=RPTS(32,1) ROLDP(26)=RPTS(33,1) ROLDP(27)=RPTS(34,1) ROLDP(28)=RPTS(22,1) ROLDP(29)=RPTS(23,1) ROLDP(30)=RPTS(24,1) ROLDP(31)=RPTS(35,1) ROLDP(32)=RPTS(36,1) ROLDP(33)=RPTS(37,1) ROLDP(34)=RPTS(38,1) ROLDP(35)=RPTS(39,1) ROLDP(36)=RPTS(40,1) ROLDP(37)=RPTS(41,1) ROLDP(38)=RPTS(42,1) ROLDP(39)=RPTS(44,1) ROLDP(40)=RPTS(45,1) ROLDP(41)=RPTS(46,1) ROLDP(42)=RPTS(47,1) ROLDP(43)=RPTS(48,1) ROLDP(44)=RPTS(49,1) ROLDP(45)=RPTS(50,1) ROLDP(46)=RPTS(51,1) ROLDP(47)=RPTS(52,1) ROLDP(48)=RPTS(53,1) ROLDP(49)=RPTS(54,1) ROLDP(50)=RPTS(55,1) ROLDP(51)=RPTS(56,1) ROLDP(52)=RPTS(57,1) ROLDP(53)=RPTS(58,1) ROLDP(54)=RPTS(59,1) ROLDP(55)=RPTS(60,1) ROLDP(56)=RPTS(61,1) NNEWP=0 KNEWP(1)=0 NPTSE=1 C C C Loop along the ray: 5 CONTINUE IF (NNEWP.LE.0) THEN C Reading the material parameters in a new point on the ray: IF (NNEWP.LT.0) THEN C WAN-02 CALL ERROR('WAN-02: Negative number of points') C This error should not appear. C The number of new points should be zero or positive integer. ENDIF IF (NPTSE+1.GT.NPTS) THEN C WAN-03 CALL ERROR('WAN-03: Array RPTS small') C The dimension of the array RPTS is given by the dimension MRAM C in the include file C ram.inc. ENDIF I1=NPTSE+1 RNEWP( 1,1)=RPTS( 2,I1) RNEWP( 2,1)=RPTS( 3,I1) RNEWP( 3,1)=RPTS( 4,I1) RNEWP( 4,1)=RPTS( 5,I1) RNEWP( 5,1)=RPTS( 6,I1) RNEWP( 6,1)=RPTS( 7,I1) RNEWP( 7,1)=RPTS( 8,I1) RNEWP( 8,1)=RPTS( 9,I1) RNEWP( 9,1)=RPTS(10,I1) RNEWP(10,1)=RPTS(11,I1) RNEWP(27,1)=RPTS(34,I1) RNEWP(28,1)=RPTS(22,I1) RNEWP(29,1)=RPTS(23,I1) RNEWP(30,1)=RPTS(24,I1) RNEWP(31,1)=RPTS(35,I1) RNEWP(32,1)=RPTS(36,I1) RNEWP(33,1)=RPTS(37,I1) RNEWP(34,1)=RPTS(38,I1) RNEWP(35,1)=RPTS(39,I1) RNEWP(36,1)=RPTS(40,I1) RNEWP(37,1)=RPTS(41,I1) RNEWP(38,1)=RPTS(42,I1) RNEWP(39,1)=RPTS(44,I1) RNEWP(40,1)=RPTS(45,I1) RNEWP(41,1)=RPTS(46,I1) RNEWP(42,1)=RPTS(47,I1) RNEWP(43,1)=RPTS(48,I1) RNEWP(44,1)=RPTS(49,I1) RNEWP(45,1)=RPTS(50,I1) RNEWP(46,1)=RPTS(51,I1) RNEWP(47,1)=RPTS(52,I1) RNEWP(48,1)=RPTS(53,I1) RNEWP(49,1)=RPTS(54,I1) RNEWP(50,1)=RPTS(55,I1) RNEWP(51,1)=RPTS(56,I1) RNEWP(52,1)=RPTS(57,I1) RNEWP(53,1)=RPTS(58,I1) RNEWP(54,1)=RPTS(59,I1) RNEWP(55,1)=RPTS(60,I1) RNEWP(56,1)=RPTS(61,I1) NNEWP=1 CALL PARM3(IPTS(12,I1),RNEWP(2,1),AA,RHO,QQ) C Computing the Christoffel matrix in the new point: CALL WACHRI(RNEWP(5,1),RNEWP(6,1),RNEWP(7,1), * AA(1,1),AA(1,2),AA(1,3),AA(1,4),AA(1,5),AA(1,6),AA(1,7), * AA(1,8),AA(1,9),AA(1,10),AA(1,11),AA(1,12),AA(1,13),AA(1,14), * AA(1,15),AA(1,16),AA(1,17),AA(1,18),AA(1,19),AA(1,20), * AA(1,21),RNEWP(11,1),RNEWP(12,1),RNEWP(13,1),EE,DER) C Storing the eigenvectors of the Christoffel matrix RNEWP(18,1)=EE(4) RNEWP(19,1)=EE(5) RNEWP(20,1)=EE(6) RNEWP(21,1)=EE(7) RNEWP(22,1)=EE(8) RNEWP(23,1)=EE(9) RNEWP(24,1)=EE(1) RNEWP(25,1)=EE(2) RNEWP(26,1)=EE(3) C Projecting the qS eigenvectors to the plane C perpendicular to the ray: CALL WAPROJ(RNEWP(5,1),RNEWP(6,1),RNEWP(7,1),RNEWP(8,1), * RNEWP(9,1),RNEWP(10,1),EE(4),EE(5),EE(6),EE(7),EE(8),EE(9), * RNEWP(14,1),RNEWP(15,1),RNEWP(16,1),RNEWP(17,1)) C C Quantities for future possible interpolation: TTNEW=RPTS(2,I1) C Coordinates: XNEW(1)=RPTS(3,I1) XNEW(2)=RPTS(4,I1) XNEW(3)=RPTS(5,I1) C dx/dt=p*v**2=p/sum(pi pi): AUX=(RPTS(6,I1)**2+RPTS(7,I1)**2+RPTS(8,I1)**2) DXNEW(1)=RPTS(6,I1)/AUX DXNEW(2)=RPTS(7,I1)/AUX DXNEW(3)=RPTS(8,I1)/AUX C Slowness: PNEW(1)=RPTS(6,I1) PNEW(2)=RPTS(7,I1) PNEW(3)=RPTS(8,I1) C dpi/dt=-1/v*dv/dxi: AUX=SQRT(AUX) DPNEW(1)=-RPTS(22,I1)/AUX DPNEW(2)=-RPTS(23,I1)/AUX DPNEW(3)=-RPTS(24,I1)/AUX C Polarization: ENEW(1)=RPTS( 9,I1) ENEW(2)=RPTS(10,I1) ENEW(3)=RPTS(11,I1) C dei/dt=v*sum(dv/dxk ek)*pi AUX=AUX*(RPTS(22,I1)*RPTS(9,I1)+RPTS(23,I1)*RPTS(10,I1)+ * RPTS(24,I1)*RPTS(11,I1)) DENEW(1)=AUX*RPTS(6,I1) DENEW(2)=AUX*RPTS(7,I1) DENEW(3)=AUX*RPTS(8,I1) C IF (ABS(RNEWP(12,1)-RNEWP(13,1)).LT.0.00001) THEN C New point is isotropic, qS eigenvalues are the same. C No change in eigenvalues and eigenvectors: RNEWP(14,1)=ROLDP(14) RNEWP(15,1)=ROLDP(15) RNEWP(16,1)=ROLDP(16) RNEWP(17,1)=ROLDP(17) RNEWP(18,1)=ROLDP(18) RNEWP(19,1)=ROLDP(19) RNEWP(20,1)=ROLDP(20) RNEWP(21,1)=ROLDP(21) RNEWP(22,1)=ROLDP(22) RNEWP(23,1)=ROLDP(23) RNEWP(24,1)=ROLDP(24) RNEWP(25,1)=ROLDP(25) RNEWP(26,1)=ROLDP(26) DELTFI=0. GOTO 20 ENDIF IF (ABS(ROLDP(12)-ROLDP(13)).LT.0.00001) THEN C Old point is isotropic, new point is anisotropic: IF ((ROLDP(14).EQ.1.).AND.(ROLDP(15).EQ.0.).AND. * (ROLDP(16).EQ.0.).AND.(ROLDP(17).EQ.1.)) THEN C New point is the first anisotropic point on the ray, C angular change is not to be computed: DELTFI=0. GOTO 20 ENDIF ENDIF IF (RPTS(2,NPTSE).EQ.RPTS(2,NPTSE+1)) THEN C The ray is crossing an interface. DELTFI=0. GOTO 20 ENDIF ENDIF C C Computing the angular change in eigenvectors, adding new points C on the ray if necessary: DO 10, I1=1,NNEWP DELTFI=WACHAN(ROLDP(14),RNEWP(11,I1)) DEFI=1./SQRT(RNEWP(12,I1))-1./SQRT(RNEWP(13,I1)) DEFI=DEFI-1./SQRT(ROLDP(12))+1./SQRT(ROLDP(13)) DEFI=DEFI*PI*(OF+FLOAT(NF-1)*DF)*RPTS(2,NPTS)/6. DEFI=ABS(DEFI) IF (DEFI*PI/8..GT.ERRWAN) THEN DEFI=ERRWAN/DEFI ELSE DEFI=PI/8. END IF IF (ABS(DELTFI).LE.DEFI) THEN C Angular change is less than prescribed limit, the I1-th point C of array RNEWP will be used as the next point on the ray: NNEWP=I1 GOTO 20 ENDIF 10 CONTINUE C Angular change is greater than prescribed limit for all the points C of array RNEWP. 15 CONTINUE C Loop for adding new points on the ray until the angular change is C small enough: IF (KNEWP(NNEWP).GE.MNEWP-2) THEN C WAN-05 CALL ERROR('WAN-05: Maximum number of divisions exceeded') C The angular change of eigenvectors in two consecutive points C is too big. More divisions than fits into the memory is C needed to keep the change under the prescribed limit. Try to C decrease the parameter C STORE, C or to increase the dimension MNEWP or the angular change DEFI. ENDIF C Adding new point to the ray: KNEWP(NNEWP)=KNEWP(NNEWP)+1 NNEWP=NNEWP+1 IF (NNEWP.GT.MNEWP) THEN C WAN-06 CALL ERROR('WAN-06: Array KNEWP too small') C This error should not appear, error05 should appear instead. ENDIF C Travel time: RNEWP(1,NNEWP)=(ROLDP(1)+RNEWP(1,NNEWP-1))*0.5 C Coordinates: CALL HIVD2(3,TTOLD,XOLD,DXOLD,TTNEW,XNEW,DXNEW, * RNEWP(1,NNEWP),RNEWP(2,NNEWP),DER) C Slowness vector: CALL HIVD2(3,TTOLD,POLD,DPOLD,TTNEW,PNEW,DPNEW, * RNEWP(1,NNEWP),RNEWP(5,NNEWP),DER) C Polarization vector: CALL HIVD2(3,TTOLD,EOLD,DEOLD,TTNEW,ENEW,DENEW, * RNEWP(1,NNEWP),RNEWP(8,NNEWP),DER) C Matrices Q and P, basis vectors of RCCS: DO 17, I1=31,56 RNEWP(I1,NNEWP)=(ROLDP(I1)+RNEWP(I1,NNEWP-1))*0.5 17 CONTINUE C Material parameters: CALL PARM2(IPTS(12,NPTSE+1),RNEWP(2,NNEWP),UP,US,RHO,QP,QS) CALL VELOC(IPTS(12,NPTSE+1),UP,US,QP,QS,VP,VS,VD,QL) RNEWP(27,NNEWP)=VD(1) RNEWP(28,NNEWP)=VD(2) RNEWP(29,NNEWP)=VD(3) RNEWP(30,NNEWP)=VD(4) CALL PARM3(IPTS(12,NPTSE+1),RNEWP(2,NNEWP),AA,RHO,QQ) C Christoffel matrix and eigenvalues: CALL WACHRI(RNEWP(5,NNEWP),RNEWP(6,NNEWP),RNEWP(7,NNEWP), * AA(1,1),AA(1,2),AA(1,3),AA(1,4),AA(1,5),AA(1,6),AA(1,7), * AA(1,8),AA(1,9),AA(1,10),AA(1,11),AA(1,12),AA(1,13),AA(1,14), * AA(1,15),AA(1,16),AA(1,17),AA(1,18),AA(1,19),AA(1,20),AA(1,21), * RNEWP(11,NNEWP),RNEWP(12,NNEWP),RNEWP(13,NNEWP),EE,DER) C Storing the eigenvectors of the Christoffel matrix RNEWP(18,NNEWP)=EE(4) RNEWP(19,NNEWP)=EE(5) RNEWP(20,NNEWP)=EE(6) RNEWP(21,NNEWP)=EE(7) RNEWP(22,NNEWP)=EE(8) RNEWP(23,NNEWP)=EE(9) RNEWP(24,NNEWP)=EE(1) RNEWP(25,NNEWP)=EE(2) RNEWP(26,NNEWP)=EE(3) C Projection of the qS eigenvectors to the plane C perpendicular to the ray: CALL WAPROJ(RNEWP(5,NNEWP),RNEWP(6,NNEWP),RNEWP(7,NNEWP), * RNEWP(8,NNEWP),RNEWP(9,NNEWP),RNEWP(10,NNEWP), * EE(4),EE(5),EE(6),EE(7),EE(8),EE(9), * RNEWP(14,NNEWP),RNEWP(15,NNEWP),RNEWP(16,NNEWP),RNEWP(17,NNEWP)) C Index of the division: KNEWP(NNEWP)=KNEWP(NNEWP-1) C IF (ABS(RNEWP(12,NNEWP)-RNEWP(13,NNEWP)).LT.0.00001) THEN C Isotropic case, qS eigenvalues are the same, C no change in eigenvectors: RNEWP(14,NNEWP)=ROLDP(14) RNEWP(15,NNEWP)=ROLDP(15) RNEWP(16,NNEWP)=ROLDP(16) RNEWP(17,NNEWP)=ROLDP(17) RNEWP(18,NNEWP)=ROLDP(18) RNEWP(19,NNEWP)=ROLDP(19) RNEWP(20,NNEWP)=ROLDP(20) RNEWP(21,NNEWP)=ROLDP(21) RNEWP(22,NNEWP)=ROLDP(22) RNEWP(23,NNEWP)=ROLDP(23) RNEWP(24,NNEWP)=ROLDP(24) RNEWP(25,NNEWP)=ROLDP(25) RNEWP(26,NNEWP)=ROLDP(26) DELTFI=0. ELSE DELTFI=WACHAN(ROLDP(14),RNEWP(11,NNEWP)) DEFI=1./SQRT(RNEWP(12,NNEWP))-1./SQRT(RNEWP(13,NNEWP)) DEFI=DEFI-1./SQRT(ROLDP(12))+1./SQRT(ROLDP(13)) DEFI=DEFI*PI*(OF+FLOAT(NF-1)*DF)*RPTS(2,NPTS)/6. DEFI=ABS(DEFI) IF (DEFI*PI/8..GT.ERRWAN) THEN DEFI=ERRWAN/DEFI ELSE DEFI=PI/8. END IF ENDIF IF (ABS(DELTFI).LE.DEFI) THEN C Angular change is less than prescribed limit, this point C of array RNEWP will be used as the next point on the ray: GOTO 20 ELSE C Angular change is greater than prescribed limit, adding C a new point to the ray: GOTO 15 ENDIF C End of the loop. C 20 CONTINUE C Angular change DELTFI for points ROLDP, RNEWP(i,NNEWP) is less C than prescribed limit. Recording the computed quantities. NPTSE=NPTSE+1 IF (NNEWP.NE.1) THEN C The new point was computed by interpolation. C Shifting the array RPTS: NPTS=NPTS+1 IF (NPTS.GT.MPTS) THEN C WAN-07 CALL ERROR('WAN-07: Array RPTS small') C The dimension of the array RPTS is given by the dimension MRAM C in the include file C ram.inc. ENDIF DO 31, I1=NPTS,NPTSE+1,-1 DO 30, I2=1,NQPTS IF ((I2.EQ.1).OR.(I2.EQ.12).OR.(I2.EQ.43)) THEN IPTS(I2,I1)=IPTS(I2,I1-1) ELSE RPTS(I2,I1)=RPTS(I2,I1-1) ENDIF 30 CONTINUE 31 CONTINUE C Recording interpolated quantities: IPTS( 1,NPTSE)=0 RPTS( 2,NPTSE)=RNEWP( 1,NNEWP) RPTS( 3,NPTSE)=RNEWP( 2,NNEWP) RPTS( 4,NPTSE)=RNEWP( 3,NNEWP) RPTS( 5,NPTSE)=RNEWP( 4,NNEWP) RPTS( 6,NPTSE)=RNEWP( 5,NNEWP) RPTS( 7,NPTSE)=RNEWP( 6,NNEWP) RPTS( 8,NPTSE)=RNEWP( 7,NNEWP) RPTS( 9,NPTSE)=RNEWP( 8,NNEWP) RPTS(10,NPTSE)=RNEWP( 9,NNEWP) RPTS(11,NPTSE)=RNEWP(10,NNEWP) IPTS(12,NPTSE)=IPTS(12,NPTSE-1) RPTS(34,NPTSE)=RNEWP(27,NNEWP) RPTS(22,NPTSE)=RNEWP(28,NNEWP) RPTS(23,NPTSE)=RNEWP(29,NNEWP) RPTS(24,NPTSE)=RNEWP(30,NNEWP) RPTS(35,NPTSE)=RNEWP(31,NNEWP) RPTS(36,NPTSE)=RNEWP(32,NNEWP) RPTS(37,NPTSE)=RNEWP(33,NNEWP) RPTS(38,NPTSE)=RNEWP(34,NNEWP) RPTS(39,NPTSE)=RNEWP(35,NNEWP) RPTS(40,NPTSE)=RNEWP(36,NNEWP) RPTS(41,NPTSE)=RNEWP(37,NNEWP) RPTS(42,NPTSE)=RNEWP(38,NNEWP) IPTS(43,NPTSE)=0 RPTS(44,NPTSE)=RNEWP(39,NNEWP) RPTS(45,NPTSE)=RNEWP(40,NNEWP) RPTS(46,NPTSE)=RNEWP(41,NNEWP) RPTS(47,NPTSE)=RNEWP(42,NNEWP) RPTS(48,NPTSE)=RNEWP(43,NNEWP) RPTS(49,NPTSE)=RNEWP(44,NNEWP) RPTS(50,NPTSE)=RNEWP(45,NNEWP) RPTS(51,NPTSE)=RNEWP(46,NNEWP) RPTS(52,NPTSE)=RNEWP(47,NNEWP) RPTS(53,NPTSE)=RNEWP(48,NNEWP) RPTS(54,NPTSE)=RNEWP(49,NNEWP) RPTS(55,NPTSE)=RNEWP(50,NNEWP) RPTS(56,NPTSE)=RNEWP(51,NNEWP) RPTS(57,NPTSE)=RNEWP(52,NNEWP) RPTS(58,NPTSE)=RNEWP(53,NNEWP) RPTS(59,NPTSE)=RNEWP(54,NNEWP) RPTS(60,NPTSE)=RNEWP(55,NNEWP) RPTS(61,NPTSE)=RNEWP(56,NNEWP) ENDIF C Recording quantities for computation of anisotropic corrections: RPTS(13,NPTSE)=RNEWP(11,NNEWP) RPTS(14,NPTSE)=RNEWP(12,NNEWP) RPTS(15,NPTSE)=RNEWP(13,NNEWP) RPTS(16,NPTSE)=RNEWP(14,NNEWP) RPTS(17,NPTSE)=RNEWP(15,NNEWP) RPTS(18,NPTSE)=RNEWP(16,NNEWP) RPTS(19,NPTSE)=RNEWP(17,NNEWP) RPTS(20,NPTSE)=DELTFI IF (KQIDT.LE.0) THEN RPTS(21,NPTSE)= * (0.5/SQRT(RNEWP(13,NNEWP))-0.5/SQRT(RNEWP(12,NNEWP)) * +0.5/SQRT(ROLDP(13)) -0.5/SQRT(ROLDP(12))) * *(RNEWP(1,NNEWP)-ROLDP(1)) ELSE RPTS(21,NPTSE)=0.25*(RNEWP(12,NNEWP)-RNEWP(13,NNEWP) * +ROLDP(12) -ROLDP(13)) * *(RNEWP(1,NNEWP)-ROLDP(1)) END IF RPTS(25,NPTSE)=RNEWP(18,NNEWP) RPTS(26,NPTSE)=RNEWP(19,NNEWP) RPTS(27,NPTSE)=RNEWP(20,NNEWP) RPTS(28,NPTSE)=RNEWP(21,NNEWP) RPTS(29,NPTSE)=RNEWP(22,NNEWP) RPTS(30,NPTSE)=RNEWP(23,NNEWP) RPTS(31,NPTSE)=RNEWP(24,NNEWP) RPTS(32,NPTSE)=RNEWP(25,NNEWP) RPTS(33,NPTSE)=RNEWP(26,NNEWP) C IF (NPTSE.LT.NPTS) THEN C Continuing with the next point on the ray: DO 100, I1=1,MQANT ROLDP(I1)=RNEWP(I1,NNEWP) 100 CONTINUE KNEWP(NNEWP)=0 NNEWP=NNEWP-1 IF (NNEWP.LE.0) THEN TTOLD=TTNEW DO 101, I1=1,3 XOLD(I1)=XNEW(I1) DXOLD(I1)=DXNEW(I1) POLD(I1)=PNEW(I1) DPOLD(I1)=DPNEW(I1) EOLD(I1)=ENEW(I1) DEOLD(I1)=DENEW(I1) 101 CONTINUE ENDIF GOTO 5 ENDIF RETURN END C C======================================================================= C SUBROUTINE WAREAD(LU1,LU2,LU3,RPTS,IPTS,NQPTS,MPTS,NPTS) C C---------------------------------------------------------------------- C Subroutine reads the unformatted output of program CRT. C Reading the output files is completed by a simple invocation of C subroutine AP00 of file 'ap.for'. Subroutine reads the quantities C stored in individual points on the ray into array RPTS. C INTEGER LU1,LU2,LU3 INTEGER NPTS,NQPTS,MPTS REAL RPTS(NQPTS,MPTS) INTEGER IPTS(NQPTS,MPTS) C Input: C LU1 ... Number of the logical unit connected to the CRT output C file with quantities along rays. C LU2 ... Number of the logical unit connected to the CRT output C file with quantities at the initial points of rays. C LU3 ... Number of the logical unit connected to the CRT output C file with quantities at the storing surface. C NQPTS.. Dimension of arrays RPTS and IPTS. C MPTS... Dimension of arrays RPTS and IPTS. C Output: C RPTS,IPTS ... The arrays are filled with the quantities stored in C CRT output files for single two-point ray during one C invocation of this subroutine: C IPTS(1,I) ... Index of the I-th point, zero for points C added to the ray by interpolation. C RPTS(2,I) ... Travel time in I-th point. C RPTS(3-5,I) . Coordinates of the point. C RPTS(6-8,I) . Slowness vector in the point. C RPTS(9-11,I) Polarization vector in the point. C IPTS(12,I) .. Index of complex block. C RPTS(22-24,I) Derivatives of the velocity. C RPTS(34,I)... S-wave velocity. C RPTS(35-38,I) Matrix of geometrical spreading. C RPTS(39-42,I) Transformation matrix P. C IPTS(43,I)... Index of surface, zero inside a complex block. C RPTS(44-52,I) Covariant components of the basis vectors C of the ray-centred coordinate system. C RPTS(53-61,I) Contravariant components of the basis vectors C of the ray-centred coordinate system. C NPTS ... Number of points stored in RPTS (IPTS). C C ........................... C Common block /POINTC/ to store the results of complete ray tracing: INCLUDE 'pointc.inc' C ........................... C Common block /RTMAT/ to store the matrices of reflection-transmittion C coefficients. INTEGER MRT,NRT PARAMETER (MRT=100) REAL PIRTR(2,2,MRT),PIRTI(2,2,MRT) COMMON/RTMAT/NRT,PIRTR,PIRTI SAVE /RTMAT/ C MRT ... Dimension of arrays. C NRT ... Number of stored R-T matrices. C PIRTR .. Real part of R-T matrices. C PIRTI .. Imaginary part of R-T matrices. C The matrices are defined here in terms of polarization vectors. C ........................... C Common block /WASEP/ to store the input SEP parameters used C in subroutine WAN and their subroutines. REAL CRTANI,DEGREE COMMON/WASEP/CRTANI,DEGREE SAVE /WASEP/ C CRTANI ... Switch between isotropic and anisotropic ray tracing. C DEGREE ... Degree of the homogeneous Hamiltonian. C....................................................................... C C Auxiliary storage locations: INTEGER ISRFA,ISRF2,KODNEW,ICBNEW,IEN,ISB2A,ISB2B,ICB2A,ICB2B,I1 INTEGER IYFA(12) REAL UP(10),US(10),RHO,QP,QS,VP,VS,VD(10),QL REAL YLFA(6),YFA(35),YYFA(5) REAL HI(18),H(18),HUI(9) C C----------------------------------------------------------------------- C ISRFA=0 NRT=0 5 CONTINUE CALL AP00(LU1,LU2,LU3) IF (IWAVE.LT.1) THEN C End of rays: NPTS=0 RETURN ENDIF IF (IREC.LE.0) C Only two-point rays are to be written into RPTS (IPTS). * GOTO 5 IF ((IPT.NE.0).AND.(IPT.NE.1)) C End of previous ray behind the reference surface: * GOTO 5 C C Initial point of a new two-point ray: NPTS=1 C Index of a point: IPTS(1,NPTS)=NPTS C Travel time: RPTS(2,NPTS)=YI(1) C Coordinates: RPTS(3,NPTS)=YI(3) RPTS(4,NPTS)=YI(4) RPTS(5,NPTS)=YI(5) C Slowness vector: RPTS(6,NPTS)=YI(6) RPTS(7,NPTS)=YI(7) RPTS(8,NPTS)=YI(8) C Polarization vector: RPTS(9,NPTS)= YI(9) RPTS(10,NPTS)=YI(10) RPTS(11,NPTS)=YI(11) C Index of complex block: IPTS(12,NPTS)=ICB1I C S-wave velocity and their derivatives: IF (CRTANI.EQ.0) THEN RPTS(34,NPTS)=YLI(2) RPTS(22,NPTS)=YLI(4) RPTS(23,NPTS)=YLI(5) RPTS(24,NPTS)=YLI(6) ELSE CALL PARM2(IPTS(12,NPTS),RPTS(3,NPTS),UP,US,RHO,QP,QS) CALL VELOC(IPTS(12,NPTS),UP,US,QP,QS,VP,VS,VD,QL) RPTS(34,NPTS)=VD(1) RPTS(22,NPTS)=VD(2) RPTS(23,NPTS)=VD(3) RPTS(24,NPTS)=VD(4) ENDIF C Matrix of geometrical spreading: RPTS(35,NPTS)=YI(12) RPTS(36,NPTS)=YI(13) RPTS(37,NPTS)=YI(16) RPTS(38,NPTS)=YI(17) C Transformation matrix P: RPTS(39,NPTS)=YI(14) RPTS(40,NPTS)=YI(15) RPTS(41,NPTS)=YI(18) RPTS(42,NPTS)=YI(19) C Index of the surface: IPTS(43,NPTS)=0 C Covariant and contravariant components of the basis vectors C of the RCCS: CALL AP03(0,RPTS(44,NPTS),H,HUI) IF (CRTANI.NE.0) THEN CALL AN03(ICB1I,YI,RPTS(44,NPTS)) ENDIF C C 20 CONTINUE C New point: IF (NYF.EQ.0) THEN C WAN-04 CALL ERROR('WAN-04: Undefined quantities at point O/F.') C Quantities at the point O/F are not defined. This may happen C if the CRT output quantities are stored at the specified surface C only, and not along rays. ENDIF IF (YF(1).LT.Y(1)) THEN C The point along the ray is before the reference surface, C recording the point along the ray: NPTS=NPTS+1 IF (NPTS.GT.MPTS) THEN C WAN-08 CALL ERROR('WAN-08: Array RPTS small') C The dimension of the array RPTS is given by the dimension MRAM C in the include file C ram.inc. ENDIF C Index of a point: IPTS(1,NPTS)=NPTS C Travel time: RPTS(2,NPTS)=YF(1) C Coordinates: RPTS(3,NPTS)=YF(3) RPTS(4,NPTS)=YF(4) RPTS(5,NPTS)=YF(5) C Slowness vector: RPTS(6,NPTS)=YF(6) RPTS(7,NPTS)=YF(7) RPTS(8,NPTS)=YF(8) C Polarization vector: RPTS(9,NPTS)= YF(9) RPTS(10,NPTS)=YF(10) RPTS(11,NPTS)=YF(11) C Index of complex block: IPTS(12,NPTS)=ICB1F C S-wave velocity and their derivatives: IF (CRTANI.EQ.0) THEN RPTS(34,NPTS)=YLF(2) RPTS(22,NPTS)=YLF(4) RPTS(23,NPTS)=YLF(5) RPTS(24,NPTS)=YLF(6) ELSE CALL PARM2(IPTS(12,NPTS),RPTS(3,NPTS),UP,US,RHO,QP,QS) CALL VELOC(IPTS(12,NPTS),UP,US,QP,QS,VP,VS,VD,QL) RPTS(34,NPTS)=VD(1) RPTS(22,NPTS)=VD(2) RPTS(23,NPTS)=VD(3) RPTS(24,NPTS)=VD(4) ENDIF C Matrix of geometrical spreading: RPTS(35,NPTS)=YF(12)*YI(12)+YF(16)*YI(13)+ * YF(20)*YI(14)+YF(24)*YI(15) RPTS(36,NPTS)=YF(13)*YI(12)+YF(17)*YI(13)+ * YF(21)*YI(14)+YF(25)*YI(15) RPTS(37,NPTS)=YF(12)*YI(16)+YF(16)*YI(17)+ * YF(20)*YI(18)+YF(24)*YI(19) RPTS(38,NPTS)=YF(13)*YI(16)+YF(17)*YI(17)+ * YF(21)*YI(18)+YF(25)*YI(19) C Transformation matrix P: RPTS(39,NPTS)=YF(14)*YI(12)+YF(18)*YI(13)+ * YF(22)*YI(14)+YF(26)*YI(15) RPTS(40,NPTS)=YF(15)*YI(12)+YF(19)*YI(13)+ * YF(23)*YI(14)+YF(27)*YI(15) RPTS(41,NPTS)=YF(14)*YI(16)+YF(18)*YI(17)+ * YF(22)*YI(18)+YF(26)*YI(19) RPTS(42,NPTS)=YF(15)*YI(16)+YF(19)*YI(17)+ * YF(23)*YI(18)+YF(27)*YI(19) C Index of the surface: IPTS(43,NPTS)=ISRFF C Covariant and contravariant components of the basis vectors C of the RCCS: CALL APYYF CALL AP03(0,HI,RPTS(44,NPTS),HUI) CALL APYYF IF (CRTANI.NE.0) THEN CALL AN03(ICB1F,YF,RPTS(44,NPTS)) ENDIF C IF ((ISRFF.NE.0).AND.(ISRFA.EQ.0)) THEN C Crossing an interface - first point at an interface: CALL BLOCK(YF(3),0,0,ISRF2,ISB2A,ICB2A) CALL BLOCK(YF(3),ISRFF,ISB2A,ISRF2,ISB2B,ICB2B) IF (ICB2A.NE.IABS(ICB1F)) THEN I1=ICB2B ICB2B=ICB2A ICB2A=I1 I1=ISB2B ISB2B=ISB2A ISB2A=I1 ENDIF IF (ICB2A.NE.IABS(ICB1F)) THEN C WAN-15 CALL ERROR('WAN-15: Wrong indices of complex blocks.') C This error should not appear, at least one of the indices C ICB2A, ICB2B should equal IABS(ICB1F). ENDIF DO 28, I1=1,5 YYFA(I1)=0. YLFA(I1)=YLF(I1) 28 CONTINUE YLFA(6)=YLF(6) DO 30, I1=1,27 YFA(I1)=YF(I1) 30 CONTINUE YFA(28)=1. YFA(29)=0. YFA(30)=0. YFA(31)=0. YFA(32)=0. YFA(33)=0. YFA(34)=1. YFA(35)=0. IF (ICB1F.GE.0) THEN IYFA(1)=31 ELSE IYFA(1)=35 ENDIF IYFA(2)=0 IYFA(3)=0 IYFA(4)=ISB2A IYFA(5)=ICB1F IYFA(6)=ISRFF IYFA(7)=ISB2B IYFA(8)=ICB2B IYFA(9)=0 IYFA(10)=0 IYFA(11)=0 IYFA(12)=0 KODNEW=0 ISRFA=ISRFF ELSEIF ((ISRFF.NE.0).AND.(ISRFA.NE.0)) THEN C Crossing an interface - second point at an interface: ICBNEW=ICB1F CALL TRANS(YLFA,YFA,YYFA,IYFA,KODNEW,ICBNEW,IEN) NRT=NRT+1 IF (NRT.GT.MRT) THEN C WAN-16 CALL ERROR('WAN-16: Small dimension of PIRTR and PIRTI.') C The arrays store R-T coefficients along the ray, the C dimension MRT should exceed the number of interfaces C passed by the ray. ENDIF PIRTR(1,1,NRT)=YFA(28) PIRTI(1,1,NRT)=YFA(29) IF (IPTS(12,NPTS-1).GE.0) THEN C Incident P wave: PIRTR(2,1,NRT)=YFA(30) PIRTI(2,1,NRT)=YFA(31) PIRTR(1,2,NRT)=0. PIRTI(1,2,NRT)=0. PIRTR(2,2,NRT)=0. PIRTI(2,2,NRT)=0. ELSE C Incident S wave: PIRTR(1,2,NRT)=YFA(30) PIRTI(1,2,NRT)=YFA(31) PIRTR(2,1,NRT)=YFA(32) PIRTI(2,1,NRT)=YFA(33) PIRTR(2,2,NRT)=YFA(34) PIRTI(2,2,NRT)=YFA(35) ENDIF RPTS(2,NPTS)=RPTS(2,NPTS-1) ISRFA=0 ENDIF C C Reading the results of the complete ray tracing: CALL AP00(LU1,LU2,LU3) IF ((IWAVE.LT.1).OR.(IPT.LE.1)) THEN C This should not happen, the ray must reach the C reference surface. C WAN-09 CALL ERROR('WAN-09: The ray missed the reference surface') C As only the two-point rays are considered by the subroutine C "WAN", each of the rays should pass the reference surface. C Check, whether you have specified right names C of the input files with points along rays, C points at their initial points and points at the reference C surface in file C CRTOUT, C or whether you have correctly specified its name C CRTOUT. ENDIF C GOTO 20 ENDIF C C The point along the ray is at or above the reference surface, C recording the point at the reference surface: NPTS=NPTS+1 IF (NPTS.GT.MPTS) THEN C WAN-10 CALL ERROR('WAN-10: Array RPTS small') C The dimension of the array RPTS is given by the dimension MRAM C of in the include file C ram.inc. ENDIF C Index of a point: IPTS(1,NPTS)=NPTS C Travel time: RPTS(2,NPTS)=Y(1) C Coordinates: RPTS(3,NPTS)=Y(3) RPTS(4,NPTS)=Y(4) RPTS(5,NPTS)=Y(5) C Slowness vector: RPTS(6,NPTS)=Y(6) RPTS(7,NPTS)=Y(7) RPTS(8,NPTS)=Y(8) C Polarization vector: RPTS(9,NPTS)= Y(9) RPTS(10,NPTS)=Y(10) RPTS(11,NPTS)=Y(11) C Index of complex block: IPTS(12,NPTS)=ICB1 C S-wave velocity and their derivatives: IF (CRTANI.EQ.0) THEN RPTS(34,NPTS)=YL(2) RPTS(22,NPTS)=YL(4) RPTS(23,NPTS)=YL(5) RPTS(24,NPTS)=YL(6) ELSE CALL PARM2(IPTS(12,NPTS),RPTS(3,NPTS),UP,US,RHO,QP,QS) CALL VELOC(IPTS(12,NPTS),UP,US,QP,QS,VP,VS,VD,QL) RPTS(34,NPTS)=VD(1) RPTS(22,NPTS)=VD(2) RPTS(23,NPTS)=VD(3) RPTS(24,NPTS)=VD(4) ENDIF C Matrix of geometrical spreading: CALL AP05(0,HUI, * RPTS(35,NPTS),RPTS(36,NPTS),RPTS(37,NPTS),RPTS(38,NPTS)) C Transformation matrix P: CALL AP06(0,HUI, * RPTS(39,NPTS),RPTS(40,NPTS),RPTS(41,NPTS),RPTS(42,NPTS)) C Index of the surface: IPTS(43,NPTS)=ISRF C Covariant and contravariant components of the basis vectors C of the RCCS: CALL AP03(0,HI,RPTS(44,NPTS),HUI) IF (CRTANI.NE.0) THEN CALL AN03(ICB1,Y,RPTS(44,NPTS)) ENDIF C RETURN END C C======================================================================= C SUBROUTINE WACHRI(P1,P2,P3,B11,B12,B22,B13,B23,B33, * B14,B24,B34,B44,B15,B25,B35,B45,B55, * B16,B26,B36,B46,B56,B66, * G1,G2,G3,EE,DER) C C---------------------------------------------------------------------- C Subroutine to compute the Christoffel matrix, its eigenvalues C and eigenvectors. REAL P1,P2,P3 REAL B11,B12,B22,B13,B23,B33,B14,B24,B34,B44 REAL B15,B25,B35,B45,B55,B16,B26,B36,B46,B56,B66 REAL G1,G2,G3,EE(9),DER(9) C C Input: C P1,P2,P3... Slowness vector. C Bii ... Values of real parts of 21 reduced C (divided by the density) elastic parameters. C Output: C G1,G2,G3 ... Eigenvalues of the Christoffel matrix. C EE ... Eigenvectors of the Christoffel matrix. C DER ... Derivatives dx/dt=dH/dp=Aijkl Ej Ek pl stored columnwise C for qP, qS1 and qS2 waves. C C----------------------------------------------------------------------- C REAL A11,A12,A13,A14,A21,A22,A23,A24,A31,A32,A33,A34,A44 REAL A15,A25,A35,A45,A55,A16,A26,A36,A46,A56,A66 REAL A1111,A2111,A3111,A1211,A2211,A3211,A1311,A2311,A3311 REAL A1121,A2121,A3121,A1221,A2221,A3221,A1321,A2321,A3321 REAL A1131,A2131,A3131,A1231,A2231,A3231,A1331,A2331,A3331 REAL A1112,A2112,A3112,A1212,A2212,A3212,A1312,A2312,A3312 REAL A1122,A2122,A3122,A1222,A2222,A3222,A1322,A2322,A3322 REAL A1132,A2132,A3132,A1232,A2232,A3232,A1332,A2332,A3332 REAL A1113,A2113,A3113,A1213,A2213,A3213,A1313,A2313,A3313 REAL A1123,A2123,A3123,A1223,A2223,A3223,A1323,A2323,A3323 REAL A1133,A2133,A3133,A1233,A2233,A3233,A1333,A2333,A3333 EQUIVALENCE (A11,A1111) EQUIVALENCE (A22,A2222) EQUIVALENCE (A33,A3333) EQUIVALENCE (A16,A1112,A1121,A1211,A2111) EQUIVALENCE (A26,A2221,A2212,A2122,A1222) EQUIVALENCE (A15,A1113,A1131,A1311,A3111) EQUIVALENCE (A35,A3331,A3313,A3133,A1333) EQUIVALENCE (A24,A2223,A2232,A2322,A3222) EQUIVALENCE (A34,A3332,A3323,A3233,A2333) EQUIVALENCE (A23,A2233,A3322) EQUIVALENCE (A13,A1133,A3311) EQUIVALENCE (A12,A1122,A2211) EQUIVALENCE (A44, A2323,A3232,A2332,A3223) EQUIVALENCE (A55, A1313,A3131,A1331,A3113) EQUIVALENCE (A66, A1212,A2121,A1221,A2112) EQUIVALENCE (A14,A1123,A1132,A2311,A3211) EQUIVALENCE (A25,A2213,A2231,A1322,A3122) EQUIVALENCE (A36,A3312,A3321,A1233,A2133) EQUIVALENCE (A56,A1321,A3112,A2113,A1231,A1213,A2131,A1312,A3121) EQUIVALENCE (A46,A2312,A3221,A1223,A2132,A2123,A1232,A2321,A3212) EQUIVALENCE (A45,A3213,A2331,A1332,A3123,A3132,A1323,A3231,A2313) REAL GAMMA(6),E11,E21,E31,E12,E22,E32,E13,E23,E33 C EQUIVALENCE (GAMMA(1),G1),(GAMMA(3),G2),(GAMMA(6),G3) C EQUIVALENCE (EE(1),E11),(EE(4),E12),(EE(7),E13) C EQUIVALENCE (EE(2),E21),(EE(5),E22),(EE(8),E23) C EQUIVALENCE (EE(3),E31),(EE(6),E32),(EE(9),E33) REAL A111,A112,A121,A122,A113,A123,A131,A132,A133 REAL A211,A212,A221,A222,A213,A223,A231,A232,A233 REAL A311,A312,A322,A313,A321,A323,A331,A332,A333 REAL AUX INTEGER KQICHM DATA KQICHM/-1/ C C GAMMA,G1,G2,G3...Christoffel matrix, later its eigenvalues. C (E11,E12,E13) C EE=(E21,E22,E23)... Eigenvectors of the christoffel matrix. C (E31,E32,E33) C A111,A211,A311,A112,A212,A312,A122,A222,A322,A113,A213,A313,A123, C A223,A323,A133,A233,A333... A(I,J,K,L)*P(L) summed over L. C A11,A21,A31,A12,A22,A32,A13,A23,A33 ... Aijk*Ek C C....................................................................... C A11=B11 A22=B22 A33=B33 A16=B16 A26=B26 A15=B15 A35=B35 A24=B24 A34=B34 A23=B23 A13=B13 A12=B12 A44=B44 A55=B55 A66=B66 A14=B14 A25=B25 A36=B36 A56=B56 A46=B46 A45=B45 C Christoffel matrix: A111=A1111*P1+A1112*P2+A1113*P3 A112=A1121*P1+A1122*P2+A1123*P3 A121=A1211*P1+A1212*P2+A1213*P3 A122=A1221*P1+A1222*P2+A1223*P3 A113=A1131*P1+A1132*P2+A1133*P3 A123=A1231*P1+A1232*P2+A1233*P3 A131=A1311*P1+A1312*P2+A1313*P3 A132=A1321*P1+A1322*P2+A1323*P3 A133=A1331*P1+A1332*P2+A1333*P3 A211=A2111*P1+A2112*P2+A2113*P3 A212=A2121*P1+A2122*P2+A2123*P3 A221=A2211*P1+A2212*P2+A2213*P3 A222=A2221*P1+A2222*P2+A2223*P3 A213=A2131*P1+A2132*P2+A2133*P3 A223=A2231*P1+A2232*P2+A2233*P3 A231=A2311*P1+A2312*P2+A2313*P3 A232=A2321*P1+A2322*P2+A2323*P3 A233=A2331*P1+A2332*P2+A2333*P3 A311=A3111*P1+A3112*P2+A3113*P3 A312=A3121*P1+A3122*P2+A3123*P3 A322=A3221*P1+A3222*P2+A3223*P3 A313=A3131*P1+A3132*P2+A3133*P3 A321=A3211*P1+A3212*P2+A3213*P3 A323=A3231*P1+A3232*P2+A3233*P3 A331=A3311*P1+A3312*P2+A3313*P3 A332=A3321*P1+A3322*P2+A3323*P3 A333=A3331*P1+A3332*P2+A3333*P3 GAMMA(1)=P1*A111+P2*A211+P3*A311 GAMMA(2)=P1*A112+P2*A212+P3*A312 GAMMA(3)=P1*A122+P2*A222+P3*A322 GAMMA(4)=P1*A113+P2*A213+P3*A313 GAMMA(5)=P1*A123+P2*A223+P3*A323 GAMMA(6)=P1*A133+P2*A233+P3*A333 C C Quasi-isotropic modification of the Christoffel matrix: IF(KQICHM.LT.0) THEN CALL RSEP3I('QICHM',KQICHM,0) END IF IF(KQICHM.GT.0) THEN AUX=SQRT(P1*P1+P2*P2+P3*P3) E11=P1/AUX E21=P2/AUX E31=P3/AUX E12=GAMMA(1)*E11+GAMMA(2)*E21+GAMMA(4)*E31 E22=GAMMA(2)*E11+GAMMA(3)*E21+GAMMA(5)*E31 E32=GAMMA(4)*E11+GAMMA(5)*E21+GAMMA(6)*E31 AUX=2.*(E11*E12+E21*E22+E31*E32) GAMMA(1)=GAMMA(1)-E11*E12-E12*E11+AUX*E11*E11 GAMMA(2)=GAMMA(2)-E11*E22-E12*E21+AUX*E11*E21 GAMMA(3)=GAMMA(3)-E21*E22-E22*E21+AUX*E21*E21 GAMMA(4)=GAMMA(4)-E11*E32-E12*E31+AUX*E11*E31 GAMMA(5)=GAMMA(5)-E21*E32-E22*E31+AUX*E21*E31 GAMMA(6)=GAMMA(6)-E31*E32-E32*E31+AUX*E31*E31 END IF C C Selecting eigenvalue and eigenvector of the Christoffel matrix: CALL EIGEN(GAMMA,EE,3,0) G1=GAMMA(1) G2=GAMMA(3) G3=GAMMA(6) E11=EE(1) E21=EE(2) E31=EE(3) E12=EE(4) E22=EE(5) E32=EE(6) E13=EE(7) E23=EE(8) E33=EE(9) IF (G3.LT.0.) THEN C WAN-11 CALL ERROR('WAN-11: Negative eigenvalue of Christoffel matrix') C This error should not appear. END IF AUX=E11*E11+E21*E21+E31*E31 IF (ABS(AUX-1.).GT.0.00001) THEN C WAN-12 CALL ERROR('WAN-12: Eigenvector is not normalized') C This error should not appear. ENDIF AUX=E12*E12+E22*E22+E32*E32 IF (ABS(AUX-1.).GT.0.00001) THEN C WAN-13 CALL ERROR('WAN-13: Eigenvector is not normalized') C This error should not appear. ENDIF AUX=E13*E13+E23*E23+E33*E33 IF (ABS(AUX-1.).GT.0.00001) THEN C WAN-14 CALL ERROR('WAN-14: Eigenvector is not normalized') C This error should not appear. ENDIF C C Computation of derivatives dx/dt: A11= A111*E11+A112*E21+A113*E31 A21= A211*E11+A212*E21+A213*E31 A31= A311*E11+A312*E21+A313*E31 A12= A121*E11+A122*E21+A123*E31 A22= A221*E11+A222*E21+A223*E31 A32= A321*E11+A322*E21+A323*E31 A13= A131*E11+A132*E21+A133*E31 A23= A231*E11+A232*E21+A233*E31 A33= A331*E11+A332*E21+A333*E31 DER(1)=A11*E11+ A12*E21+ A13*E31 DER(2)=A21*E11+ A22*E21+ A23*E31 DER(3)=A31*E11+ A32*E21+ A33*E31 A11= A111*E12+A112*E22+A113*E32 A21= A211*E12+A212*E22+A213*E32 A31= A311*E12+A312*E22+A313*E32 A12= A121*E12+A122*E22+A123*E32 A22= A221*E12+A222*E22+A223*E32 A32= A321*E12+A322*E22+A323*E32 A13= A131*E12+A132*E22+A133*E32 A23= A231*E12+A232*E22+A233*E32 A33= A331*E12+A332*E22+A333*E32 DER(4)=A11*E12+ A12*E22+ A13*E32 DER(5)=A21*E12+ A22*E22+ A23*E32 DER(6)=A31*E12+ A32*E22+ A33*E32 A11= A111*E13+A112*E23+A113*E33 A21= A211*E13+A212*E23+A213*E33 A31= A311*E13+A312*E23+A313*E33 A12= A121*E13+A122*E23+A123*E33 A22= A221*E13+A222*E23+A223*E33 A32= A321*E13+A322*E23+A323*E33 A13= A131*E13+A132*E23+A133*E33 A23= A231*E13+A232*E23+A233*E33 A33= A331*E13+A332*E23+A333*E33 DER(7)=A11*E13+ A12*E23+ A13*E33 DER(8)=A21*E13+ A22*E23+ A23*E33 DER(9)=A31*E13+ A32*E23+ A33*E33 C RETURN END C C======================================================================= C SUBROUTINE WAMAT(A,B,C,D) C C---------------------------------------------------------------------- C Subroutine to compute the product of two 2x2 complex matrices. C The second matrix (C+iD) is destroyed in the computation. REAL A(2,2),B(2,2),C(2,2),D(2,2) C Input: C A,B,C,D ... Real and imaginary parts of the two matrices. C Output: C C,D ... Real and imaginary parts of resulting matrix. C C....................................................................... C Auxiliary storage locations: REAL CR11,CR21,CR12,CR22,CI11,CI21,CI12,CI22 C....................................................................... CR11=A(1,1)*C(1,1)-B(1,1)*D(1,1)+A(1,2)*C(2,1)-B(1,2)*D(2,1) CR21=A(2,1)*C(1,1)-B(2,1)*D(1,1)+A(2,2)*C(2,1)-B(2,2)*D(2,1) CR12=A(1,1)*C(1,2)-B(1,1)*D(1,2)+A(1,2)*C(2,2)-B(1,2)*D(2,2) CR22=A(2,1)*C(1,2)-B(2,1)*D(1,2)+A(2,2)*C(2,2)-B(2,2)*D(2,2) C CI11=A(1,1)*D(1,1)+B(1,1)*C(1,1)+A(1,2)*D(2,1)+B(1,2)*C(2,1) CI21=A(2,1)*D(1,1)+B(2,1)*C(1,1)+A(2,2)*D(2,1)+B(2,2)*C(2,1) CI12=A(1,1)*D(1,2)+B(1,1)*C(1,2)+A(1,2)*D(2,2)+B(1,2)*C(2,2) CI22=A(2,1)*D(1,2)+B(2,1)*C(1,2)+A(2,2)*D(2,2)+B(2,2)*C(2,2) C C(1,1)=CR11 C(2,1)=CR21 C(1,2)=CR12 C(2,2)=CR22 D(1,1)=CI11 D(2,1)=CI21 D(1,2)=CI12 D(2,2)=CI22 C RETURN END C C======================================================================= C REAL FUNCTION WACHAN(G,H) C C---------------------------------------------------------------------- C Subroutine to compute the smallest angle between the two-dimensional C vector G and vectors U,V,-U,-V. C The subroutine reorganizes the vectors U and V in such way, C that the pair U,V is equal to pair G,H rotated with angle WACHAN. C The real numbers R1 (associated with U) and R2 (associated with V), C and the vectors PU (associated with U) and PV (associated with V), C are reorganized in the same way. C REAL G(4),H(13) C G(1:4)=G1,G2,H1,H2 C H(1:7)=R1,R2,R3,U1,U2,V1,V2,PU1,PU2,PU3,PV1,PV2,PV3 C Input: C G1,G2,H1,H2 ... A pair of two-dimensional orthonormal vectors. C U1,U2,V1,V2 ... A pair of two-dimensional orthonormal vectors. C PU1,PU2,PU3,PV1,PV2,PV3 ... A pair of three-dimensional vectors C associated to U and V. C R1,R2 ... Real numbers associated to U and V. C Output: C WACHAN ... The smallest one from the angles between vector G C and vectors U,V,-U,-V. C U1,U2,V1,V2 ... Selection from U,V,-U,-V in such way, that C the pair U,V is equal to pair G,H rotated C with angle WACHAN. C PU1,PU2,PU3,PV1,PV2,PV3 ... A pair of three-dimensional vectors C associated to U and V. C R1,R2 ... Real numbers associeted to U and V. C C....................................................................... C C Auxiliary storage locations: REAL SP1,SP2,SP3,SP4,A1,A2,A3,AUX REAL G1,G2,H1,H2,U1,U2,V1,V2,PU1,PU2,PU3,PV1,PV2,PV3,R1,R2 C C----------------------------------------------------------------------- G1=G(1) G2=G(2) H1=G(3) H2=G(4) R1=H(2) R2=H(3) U1=H(4) U2=H(5) V1=H(6) V2=H(7) PU1=H(8) PU2=H(9) PU3=H(10) PV1=H(11) PV2=H(12) PV3=H(13) C SP1=ABS(1 -( G1*U1+G2*U2)) SP2=ABS(1 -(-G1*U1-G2*U2)) SP3=ABS(1 -( G1*V1+G2*V2)) SP4=ABS(1 -(-G1*V1-G2*V2)) AUX=AMIN1(SP1,SP2,SP3,SP4) IF (AUX.EQ.SP1) THEN C No action. ELSEIF (AUX.EQ.SP2) THEN U1=-U1 U2=-U2 PU1=-PU1 PU2=-PU2 PU3=-PU3 ELSEIF (AUX.EQ.SP3) THEN A1=U1 A2=U2 U1=V1 U2=V2 V1=A1 V2=A2 A1=PU1 A2=PU2 A3=PU3 PU1=PV1 PU2=PV2 PU3=PV3 PV1=A1 PV2=A2 PV3=A3 AUX=R1 R1=R2 R2=AUX ELSEIF (AUX.EQ.SP4) THEN A1=U1 A2=U2 U1=-V1 U2=-V2 V1=A1 V2=A2 A1=PU1 A2=PU2 A3=PU3 PU1=-PV1 PU2=-PV2 PU3=-PV3 PV1=A1 PV2=A2 PV3=A3 AUX=R1 R1=R2 R2=AUX ENDIF SP1=ABS(1 - ( H1*V1+H2*V2)) SP2=ABS(1 - (-H1*V1-H2*V2)) AUX=AMIN1(SP1,SP2) IF (AUX.EQ.SP1) THEN C No action. ELSEIF (AUX.EQ.SP2) THEN V1=-V1 V2=-V2 PV1=-PV1 PV2=-PV2 PV3=-PV3 ENDIF WACHAN=ASIN(0.5*((G1+U1)*(H1-V1)+(G2+U2)*(H2-V2))) H(2)=R1 H(3)=R2 H(4)=U1 H(5)=U2 H(6)=V1 H(7)=V2 H(8)=PU1 H(9)=PU2 H(10)=PU3 H(11)=PV1 H(12)=PV2 H(13)=PV3 RETURN END C C C======================================================================= C SUBROUTINE WAPROJ(P1,P2,P3,E1,E2,E3,G1,G2,G3,H1,H2,H3, * GE1,GE2,HE1,HE2) C C---------------------------------------------------------------------- C Subroutine to project vectors G and H to the plane defined by vector E C and vector PxE. C REAL P1,P2,P3,E1,E2,E3,G1,G2,G3,H1,H2,H3,GE1,GE2,HE1,HE2 C Input: C P1,P2,P3,E1,E2,E3 ... Vectors defining the plane. C G1,G2,G3,H1,H2,H3 ... Vectors to be projected. C Output: C GE1,GE2,HE1,HE2 ... Projected vectors. C C....................................................................... C C Auxiliary storage locations: REAL F1,F2,F3,AUX C C----------------------------------------------------------------------- C Second vector defining the plane: F1=P2*E3-E2*P3 F2=P3*E1-E3*P1 F3=P1*E2-E1*P2 AUX=SQRT(F1*F1+F2*F2+F3*F3) F1=F1/AUX F2=F2/AUX F3=F3/AUX C Projecting vectors G and H to the plane defined by E and F: GE1=E1*G1+E2*G2+E3*G3 GE2=F1*G1+F2*G2+F3*G3 HE1=E1*H1+E2*H2+E3*H3 HE2=F1*H1+F2*H2+F3*H3 C RETURN END C C======================================================================= C SUBROUTINE WAPERT(RPTS,IPTS,NQPTS,STEP,DTAU) C C---------------------------------------------------------------------- C Subroutine to compute second-order travel-time corrections. INTEGER NQPTS REAL RPTS(NQPTS),STEP,DTAU(7) INTEGER IPTS(NQPTS) C C Input: C IPTS,RPTS,NQPTS... Quantities in the point on the ray. C See the description above. C STEP... Step of the independent variable along the ray, see C the CRT input parameter C STORE. C Output: C DTAU... Second-order travel-time corrections. C C....................................................................... C Subroutines and external functions required: EXTERNAL WASUM,WASUM3,WASUM4,WASUM5 REAL WASUM,WASUM3,WASUM4,WASUM5 C WASUM,WASUM3,WASUM4,WASUM5 ... This file. C ........................... C Common block /WASEP/ to store the input SEP parameters used C in subroutine WAN and their subroutines. REAL CRTANI,DEGREE COMMON/WASEP/CRTANI,DEGREE SAVE /WASEP/ C CRTANI ... Switch between isotropic and anisotropic ray tracing. C DEGREE ... Degree of the homogeneous Hamiltonian. C....................................................................... C REAL AA(10,21),RHO,QQ(21), * AAA(3,3,3,3),AAA1(3,3,3,3),AAA2(3,3,3,3),AAA3(3,3,3,3) REAL V0,VV0(3),P(3),E1(3),E2(3),H1(3),H2(3),H3(3), * HC1(3),HC2(3),HC3(3), * HHC0(3),HHC0SP(3,3),HHC0S(3,3),HHTC0S(2,2), * QQT(4),PPT(4),QIT(4),GIJT(2,2), * HH0(3),HHT0(3),HH1,HH2,HH3, * HH1L(3),HH2L(3),HH3L(3),HH1U(3),HH2U(3),HH3U(3), * HHT1L(2),HHT2L(2),HHT3L(2),HHT1U(2),HHT2U(2),HHT3U(2), * QTKT(7),TT(7),TAUT(4),TOUT(2,3),VK(6),T(6) REAL HDEP,HDES,HDE,HDE1,HDE2,HDE3,HDE4,HDE5,HDE6, * HDE11,HDE12,HDE22,HDE13,HDE23,HDE33,HDE14,HDE24,HDE34, * HDE44,HDE15,HDE25,HDE35,HDE45,HDE55,HDE16,HDE26,HDE36, * HDE46,HDE56,HDE66 INTEGER ISRF1T,IFUN1T(7),IFUN2T(7),ISRF1K,IFUN1K(6),IFUN2K(6), * NSUM1,NSUM2,IX,NDER,NFUN11,NFUN12,NFUN21,NFUN22 REAL AUX,PIPI,X1T,FUN1T(7),X1K,FUN1K(6) SAVE NSUM1,IX,NDER,X1T,ISRF1T,NFUN11,NFUN21,IFUN1T,FUN1T,IFUN2T, * NSUM2,X1K,ISRF1K,NFUN12,NFUN22,IFUN1K,FUN1K,IFUN2K C AA ... Values, first and second partial derivatives of real C parts of 21 reduced (divided by the density) elastic C parameters. The order of the value, first and second C partial derivatives of each parameter Aij is: C Aij,Aij1,Aij2,Aij3,Aij11,Aij12,Aij22,Aij13,Aij23,Aij33. C The order of parameters (second array index) is: C A11,A12,A22,A13,A23,A33,A14,A24,A34,A44,A15,A25,A35,A45, C A55,A16,A26,A36,A46,A56,A66. C AAA... 3*3*3*3 matrix of reduced elastic parameters. C AAA1... 3*3*3*3 matrix of first partial derivatives of reduced C elastic parameters according to X1. C AAA2... 3*3*3*3 matrix of first partial derivatives of reduced C elastic parameters according to X2. C AAA3... 3*3*3*3 matrix of first partial derivatives of reduced C elastic parameters according to X3. C----------------------------------------------------------------------- C IF (CRTANI.EQ.0.) THEN V0=RPTS(34) VV0(1)=RPTS(22) VV0(2)=RPTS(23) VV0(3)=RPTS(24) P(1)=RPTS(6) P(2)=RPTS(7) P(3)=RPTS(8) E1(1)=RPTS(25) E1(2)=RPTS(26) E1(3)=RPTS(27) E2(1)=RPTS(28) E2(2)=RPTS(29) E2(3)=RPTS(30) QQT(1)=RPTS(35) QQT(2)=RPTS(36) QQT(3)=RPTS(37) QQT(4)=RPTS(38) PPT(1)=RPTS(39) PPT(2)=RPTS(40) PPT(3)=RPTS(41) PPT(4)=RPTS(42) CALL PARM3(IPTS(12),RPTS(3),AA,RHO,QQ) CALL WAIJKL(AA,AAA,1) CALL WAIJKL(AA,AAA1,2) CALL WAIJKL(AA,AAA2,3) CALL WAIJKL(AA,AAA3,4) H1(1)=RPTS(44) H1(2)=RPTS(45) H1(3)=RPTS(46) H2(1)=RPTS(47) H2(2)=RPTS(48) H2(3)=RPTS(49) H3(1)=RPTS(50) H3(2)=RPTS(51) H3(3)=RPTS(52) HH0(1)=VV0(1)/V0 HH0(2)=VV0(2)/V0 HH0(3)=VV0(3)/V0 HHT0(1)=WASUM(HH0,H1) HHT0(2)=WASUM(HH0,H2) HHT0(3)=WASUM(HH0,H3) HH1=-1./SQRT(WASUM5(AAA,E1,P,E1,P)) HH2=-1./SQRT(WASUM5(AAA,E2,P,E2,P)) HH1L(1)=-0.5*(HH1**3)*WASUM5(AAA1,E1,P,E1,P) HH1L(2)=-0.5*(HH1**3)*WASUM5(AAA2,E1,P,E1,P) HH1L(3)=-0.5*(HH1**3)*WASUM5(AAA3,E1,P,E1,P) HH2L(1)=-0.5*(HH2**3)*WASUM5(AAA1,E2,P,E2,P) HH2L(2)=-0.5*(HH2**3)*WASUM5(AAA2,E2,P,E2,P) HH2L(3)=-0.5*(HH2**3)*WASUM5(AAA3,E2,P,E2,P) HH1U(1)=-(HH1**3)*WASUM4(AAA,E1,E1,P,1) HH1U(2)=-(HH1**3)*WASUM4(AAA,E1,E1,P,2) HH1U(3)=-(HH1**3)*WASUM4(AAA,E1,E1,P,3) HH2U(1)=-(HH2**3)*WASUM4(AAA,E2,E2,P,1) HH2U(2)=-(HH2**3)*WASUM4(AAA,E2,E2,P,2) HH2U(3)=-(HH2**3)*WASUM4(AAA,E2,E2,P,3) HHT1L(1)=WASUM(HH1L,H1) HHT1L(2)=WASUM(HH1L,H2) HHT2L(1)=WASUM(HH2L,H1) HHT2L(2)=WASUM(HH2L,H2) HHT1U(1)=WASUM(HH1U,H1) HHT1U(2)=WASUM(HH1U,H2) HHT2U(1)=WASUM(HH2U,H1) HHT2U(2)=WASUM(HH2U,H2) QTKT(1)=-((HHT1L(1)+HH1*HHT0(1))*QQT(1)+ * (HHT1L(2)+HH1*HHT0(2))*QQT(2)+ * HHT1U(1)*PPT(1)+HHT1U(2)*PPT(2)) QTKT(2)=-((HHT2L(1)+HH2*HHT0(1))*QQT(1)+ * (HHT2L(2)+HH2*HHT0(2))*QQT(2)+ * HHT2U(1)*PPT(1)+HHT2U(2)*PPT(2)) QTKT(3)=-((HHT1L(1)+HH1*HHT0(1))*QQT(3)+ * (HHT1L(2)+HH1*HHT0(2))*QQT(4)+ * HHT1U(1)*PPT(3)+HHT1U(2)*PPT(4)) QTKT(4)=-((HHT2L(1)+HH2*HHT0(1))*QQT(3)+ * (HHT2L(2)+HH2*HHT0(2))*QQT(4)+ * HHT2U(1)*PPT(3)+HHT2U(2)*PPT(4)) IF (IPTS(1).EQ.1) THEN C Initial point of the ray: C Initiating first integration: NSUM1=4 IX=1 NDER=1 NFUN11=4 NFUN21=4 IFUN1T(1)=1 IFUN1T(2)=3 IFUN1T(3)=2 IFUN1T(4)=4 FUN1T(1)=QTKT(1) FUN1T(2)=QTKT(2) FUN1T(3)=QTKT(3) FUN1T(4)=QTKT(4) IFUN2T(1)=1 IFUN2T(2)=3 IFUN2T(3)=2 IFUN2T(4)=4 C Initiating second integration: NSUM2=3 NFUN12=3 NFUN22=3 IFUN1K(1)=1 IFUN1K(2)=2 IFUN1K(3)=3 AUX=V0**2 TAUT(1)=-HHT1U(1)/AUX TAUT(2)=-HHT1U(2)/AUX TAUT(3)=-HHT2U(1)/AUX TAUT(4)=-HHT2U(2)/AUX FUN1K(1)=-(2.*(HHT1U(1)*TAUT(1)+HHT1U(2)*TAUT(2)) + * (TAUT(1)*TAUT(1)+TAUT(2)*TAUT(2))*AUX) FUN1K(2)=-(2.*(HHT2U(1)*TAUT(3)+HHT2U(2)*TAUT(4)) + * (TAUT(3)*TAUT(3)+TAUT(4)*TAUT(4))*AUX) FUN1K(3)=-(HHT1U(1)*TAUT(3)+HHT1U(2)*TAUT(4) + * HHT2U(1)*TAUT(1)+HHT2U(2)*TAUT(2) + * (TAUT(1)*TAUT(3)+TAUT(2)*TAUT(4))*AUX) IFUN2K(1)=1 IFUN2K(2)=2 IFUN2K(3)=3 C Output values: DTAU(1)=0. DTAU(2)=0. DTAU(3)=0. DTAU(4)=0. DTAU(5)=0. DTAU(6)=0. DTAU(7)=0. RETURN ENDIF CALL AP28(NSUM1,TT,IX,NDER,STEP,X1T,ISRF1T,NFUN11,IFUN1T,FUN1T, * NFUN21,IFUN2T,QTKT) AUX=QQT(1)*QQT(4)-QQT(2)*QQT(3) IF (AUX.EQ.0) THEN C WAN-18 CALL ERROR('WAN-18: Zero determinant of geom. spread.') C This error should not appear. ENDIF QIT(1)= QQT(4)/AUX QIT(4)= QQT(1)/AUX QIT(2)=-QQT(2)/AUX QIT(3)=-QQT(3)/AUX TAUT(1)=QIT(1)*TT(1)+QIT(2)*TT(2) TAUT(2)=QIT(3)*TT(1)+QIT(4)*TT(2) TAUT(3)=QIT(1)*TT(3)+QIT(2)*TT(4) TAUT(4)=QIT(3)*TT(3)+QIT(4)*TT(4) VK(1)=-(2.*(HHT1U(1)*TAUT(1)+HHT1U(2)*TAUT(2)) + * (TAUT(1)*TAUT(1)+TAUT(2)*TAUT(2))*V0**2) VK(2)=-(2.*(HHT2U(1)*TAUT(3)+HHT2U(2)*TAUT(4)) + * (TAUT(3)*TAUT(3)+TAUT(4)*TAUT(4))*V0**2) VK(3)=-(HHT1U(1)*TAUT(3)+HHT1U(2)*TAUT(4) + * HHT2U(1)*TAUT(1)+HHT2U(2)*TAUT(2) + * (TAUT(1)*TAUT(3)+TAUT(2)*TAUT(4))*V0**2) CALL AP28(NSUM2,T,IX,NDER,STEP,X1K,ISRF1K,NFUN12,IFUN1K,FUN1K, * NFUN22,IFUN2K,VK) C Second order perturbations towards anisotropic rays: DTAU(1)=T(1)/2. DTAU(3)=T(2)/2. C Estimation of the two equal second order perturbations from C anisotropic common ray towards anisotropic rays: DTAU(2)=T(3) DTAU(4)=0. DTAU(5)=0. DTAU(6)=0. DTAU(7)=0. RETURN ENDIF C IF (CRTANI.NE.0.) THEN V0=RPTS(34) VV0(1)=RPTS(22) VV0(2)=RPTS(23) VV0(3)=RPTS(24) P(1)=RPTS(6) P(2)=RPTS(7) P(3)=RPTS(8) E1(1)=RPTS(25) E1(2)=RPTS(26) E1(3)=RPTS(27) E2(1)=RPTS(28) E2(2)=RPTS(29) E2(3)=RPTS(30) QQT(1)=RPTS(35) QQT(2)=RPTS(36) QQT(3)=RPTS(37) QQT(4)=RPTS(38) PPT(1)=RPTS(39) PPT(2)=RPTS(40) PPT(3)=RPTS(41) PPT(4)=RPTS(42) CALL PARM3(IPTS(12),RPTS(3),AA,RHO,QQ) CALL WAIJKL(AA,AAA,1) CALL WAIJKL(AA,AAA1,2) CALL WAIJKL(AA,AAA2,3) CALL WAIJKL(AA,AAA3,4) H1(1)=RPTS(44) H1(2)=RPTS(45) H1(3)=RPTS(46) H2(1)=RPTS(47) H2(2)=RPTS(48) H2(3)=RPTS(49) H3(1)=RPTS(50) H3(2)=RPTS(51) H3(3)=RPTS(52) HC1(1)=RPTS(53) HC1(2)=RPTS(54) HC1(3)=RPTS(55) HC2(1)=RPTS(56) HC2(2)=RPTS(57) HC2(3)=RPTS(58) HC3(1)=RPTS(59) HC3(2)=RPTS(60) HC3(3)=RPTS(61) CALL HDER(-1,AA,P(1),P(2),P(3),HDEP,HDES,HDE, * HDE1,HDE2,HDE3,HDE4,HDE5,HDE6, * HDE11,HDE12,HDE22,HDE13,HDE23,HDE33,HDE14,HDE24,HDE34, * HDE44,HDE15,HDE25,HDE35,HDE45,HDE55,HDE16,HDE26,HDE36, * HDE46,HDE56,HDE66) HH0(1)=HDE1 HH0(2)=HDE2 HH0(3)=HDE3 HHC0(1)=HDE4 HHC0(2)=HDE5 HHC0(3)=HDE6 HHC0SP(1,1)=HDE44 HHC0SP(2,1)=HDE45 HHC0SP(3,1)=HDE46 HHC0SP(1,2)=HDE45 HHC0SP(2,2)=HDE55 HHC0SP(3,2)=HDE56 HHC0SP(1,3)=HDE46 HHC0SP(2,3)=HDE56 HHC0SP(3,3)=HDE66 C Eq. 10: HHC0S(1,1)=HHC0SP(1,1)-3.*HHC0(1)*HHC0(1) HHC0S(2,1)=HHC0SP(2,1)-3.*HHC0(2)*HHC0(1) HHC0S(3,1)=HHC0SP(3,1)-3.*HHC0(3)*HHC0(1) HHC0S(1,2)=HHC0SP(1,2)-3.*HHC0(1)*HHC0(2) HHC0S(2,2)=HHC0SP(2,2)-3.*HHC0(2)*HHC0(2) HHC0S(3,2)=HHC0SP(3,2)-3.*HHC0(3)*HHC0(2) HHC0S(1,3)=HHC0SP(1,3)-3.*HHC0(1)*HHC0(3) HHC0S(2,3)=HHC0SP(2,3)-3.*HHC0(2)*HHC0(3) HHC0S(3,3)=HHC0SP(3,3)-3.*HHC0(3)*HHC0(3) C Eq. 39: HHT0(1)=WASUM(HH0,HC1) HHT0(2)=WASUM(HH0,HC2) HHT0(3)=WASUM(HH0,HC3) C Eq. 41: HHTC0S(1,1)=WASUM3(HHC0S,H1,H1) HHTC0S(1,2)=WASUM3(HHC0S,H1,H2) HHTC0S(2,1)=WASUM3(HHC0S,H2,H1) HHTC0S(2,2)=WASUM3(HHC0S,H2,H2) C Eq. 11: HH1=-1./SQRT(WASUM5(AAA,E1,P,E1,P)) HH2=-1./SQRT(WASUM5(AAA,E2,P,E2,P)) C Eq. 15: PIPI=(P(1)**2+P(2)**2+P(3)**2) HH3=-1./V0/SQRT(PIPI) C Eq. 12: HH1L(1)=-0.5*(HH1**3)*WASUM5(AAA1,E1,P,E1,P) HH1L(2)=-0.5*(HH1**3)*WASUM5(AAA2,E1,P,E1,P) HH1L(3)=-0.5*(HH1**3)*WASUM5(AAA3,E1,P,E1,P) HH2L(1)=-0.5*(HH2**3)*WASUM5(AAA1,E2,P,E2,P) HH2L(2)=-0.5*(HH2**3)*WASUM5(AAA2,E2,P,E2,P) HH2L(3)=-0.5*(HH2**3)*WASUM5(AAA3,E2,P,E2,P) C Eq. 16: HH3L(1)=VV0(1)/V0**2/SQRT(PIPI) HH3L(2)=VV0(2)/V0**2/SQRT(PIPI) HH3L(3)=VV0(3)/V0**2/SQRT(PIPI) C Eq. 13: HH1U(1)=-(HH1**3)*WASUM4(AAA,E1,E1,P,1) HH1U(2)=-(HH1**3)*WASUM4(AAA,E1,E1,P,2) HH1U(3)=-(HH1**3)*WASUM4(AAA,E1,E1,P,3) HH2U(1)=-(HH2**3)*WASUM4(AAA,E2,E2,P,1) HH2U(2)=-(HH2**3)*WASUM4(AAA,E2,E2,P,2) HH2U(3)=-(HH2**3)*WASUM4(AAA,E2,E2,P,3) C Eq. 17: HH3U(1)=1./V0/SQRT(PIPI**3)*P(1) HH3U(2)=1./V0/SQRT(PIPI**3)*P(2) HH3U(3)=1./V0/SQRT(PIPI**3)*P(3) C Eq. 42: HHT1L(1)=WASUM(HH1L,HC1) HHT1L(2)=WASUM(HH1L,HC2) HHT2L(1)=WASUM(HH2L,HC1) HHT2L(2)=WASUM(HH2L,HC2) HHT3L(1)=WASUM(HH3L,HC1) HHT3L(2)=WASUM(HH3L,HC2) C Eq. 43: HHT1U(1)=WASUM(HH1U,H1) HHT1U(2)=WASUM(HH1U,H2) HHT2U(1)=WASUM(HH2U,H1) HHT2U(2)=WASUM(HH2U,H2) HHT3U(1)=WASUM(HH3U,H1) HHT3U(2)=WASUM(HH3U,H2) C Eq. 55: QTKT(1)=-((HHT1L(1)+HH1*HHT0(1))*QQT(1)+ * (HHT1L(2)+HH1*HHT0(2))*QQT(2)+ * HHT1U(1)*PPT(1)+HHT1U(2)*PPT(2)) QTKT(2)=-((HHT2L(1)+HH2*HHT0(1))*QQT(1)+ * (HHT2L(2)+HH2*HHT0(2))*QQT(2)+ * HHT2U(1)*PPT(1)+HHT2U(2)*PPT(2)) QTKT(3)=-((HHT3L(1)+HH3*HHT0(1))*QQT(1)+ * (HHT3L(2)+HH3*HHT0(2))*QQT(2)+ * HHT3U(1)*PPT(1)+HHT3U(2)*PPT(2)) QTKT(4)=-((HHT1L(1)+HH1*HHT0(1))*QQT(3)+ * (HHT1L(2)+HH1*HHT0(2))*QQT(4)+ * HHT1U(1)*PPT(3)+HHT1U(2)*PPT(4)) QTKT(5)=-((HHT2L(1)+HH2*HHT0(1))*QQT(3)+ * (HHT2L(2)+HH2*HHT0(2))*QQT(4)+ * HHT2U(1)*PPT(3)+HHT2U(2)*PPT(4)) QTKT(6)=-((HHT3L(1)+HH3*HHT0(1))*QQT(3)+ * (HHT3L(2)+HH3*HHT0(2))*QQT(4)+ * HHT3U(1)*PPT(3)+HHT3U(2)*PPT(4)) QTKT(7)=-HH3 IF (IPTS(1).EQ.1) THEN C Initial point of the ray: C Initiating first integration: NSUM1=7 IX=1 NDER=1 NFUN11=7 NFUN21=7 IFUN1T(1)=1 IFUN1T(2)=3 IFUN1T(3)=5 IFUN1T(4)=2 IFUN1T(5)=4 IFUN1T(6)=6 IFUN1T(7)=7 FUN1T(1)=QTKT(1) FUN1T(2)=QTKT(2) FUN1T(3)=QTKT(3) FUN1T(4)=QTKT(4) FUN1T(5)=QTKT(5) FUN1T(6)=QTKT(6) FUN1T(7)=QTKT(7) IFUN2T(1)=1 IFUN2T(2)=3 IFUN2T(3)=5 IFUN2T(4)=2 IFUN2T(5)=4 IFUN2T(6)=6 IFUN2T(7)=7 C Initiating second integration: NSUM2=6 NFUN12=6 NFUN22=6 IFUN1K(1)=1 IFUN1K(2)=2 IFUN1K(3)=3 IFUN1K(4)=4 IFUN1K(5)=5 IFUN1K(6)=6 AUX=HHTC0S(1,1)*HHTC0S(2,2)-HHTC0S(2,1)*HHTC0S(1,2) IF (AUX.EQ.0) THEN C WAN-19 CALL ERROR * ('WAN-19: Zero determinant of second der. of Ham.') C This error should not appear. ENDIF GIJT(1,1)= HHTC0S(2,2)/AUX GIJT(2,2)= HHTC0S(1,1)/AUX GIJT(2,1)=-HHTC0S(2,1)/AUX GIJT(1,2)=-HHTC0S(1,2)/AUX TOUT(1,1)=-GIJT(1,1)*HHT1U(1)-GIJT(1,2)*HHT1U(2) TOUT(2,1)=-GIJT(2,1)*HHT1U(1)-GIJT(2,2)*HHT1U(2) TOUT(1,2)=-GIJT(1,1)*HHT2U(1)-GIJT(1,2)*HHT2U(2) TOUT(2,2)=-GIJT(2,1)*HHT2U(1)-GIJT(2,2)*HHT2U(2) TOUT(1,3)=-GIJT(1,1)*HHT3U(1)-GIJT(1,2)*HHT3U(2) TOUT(2,3)=-GIJT(2,1)*HHT3U(1)-GIJT(2,2)*HHT3U(2) FUN1K(1)=-(HHT1U(1)*TOUT(1,1)+HHT1U(2)*TOUT(2,1) + * HHT1U(1)*TOUT(1,1)+HHT1U(2)*TOUT(2,1) + * HHTC0S(1,1)*TOUT(1,1)*TOUT(1,1)+ * HHTC0S(2,1)*TOUT(2,1)*TOUT(1,1)+ * HHTC0S(1,2)*TOUT(1,1)*TOUT(2,1)+ * HHTC0S(2,2)*TOUT(2,1)*TOUT(2,1)) FUN1K(2)=-(HHT1U(1)*TOUT(1,2)+HHT1U(2)*TOUT(2,2) + * HHT2U(1)*TOUT(1,1)+HHT2U(2)*TOUT(2,1) + * HHTC0S(1,1)*TOUT(1,1)*TOUT(1,2)+ * HHTC0S(2,1)*TOUT(2,1)*TOUT(1,2)+ * HHTC0S(1,2)*TOUT(1,1)*TOUT(2,2)+ * HHTC0S(2,2)*TOUT(2,1)*TOUT(2,2)) FUN1K(3)=-(HHT2U(1)*TOUT(1,2)+HHT2U(2)*TOUT(2,2) + * HHT2U(1)*TOUT(1,2)+HHT2U(2)*TOUT(2,2) + * HHTC0S(1,1)*TOUT(1,2)*TOUT(1,2)+ * HHTC0S(2,1)*TOUT(2,2)*TOUT(1,2)+ * HHTC0S(1,2)*TOUT(1,2)*TOUT(2,2)+ * HHTC0S(2,2)*TOUT(2,2)*TOUT(2,2)) FUN1K(4)=-(HHT1U(1)*TOUT(1,3)+HHT1U(2)*TOUT(2,3) + * HHT3U(1)*TOUT(1,1)+HHT3U(2)*TOUT(2,1) + * HHTC0S(1,1)*TOUT(1,1)*TOUT(1,3)+ * HHTC0S(2,1)*TOUT(2,1)*TOUT(1,3)+ * HHTC0S(1,2)*TOUT(1,1)*TOUT(2,3)+ * HHTC0S(2,2)*TOUT(2,1)*TOUT(2,3)) FUN1K(5)=-(HHT2U(1)*TOUT(1,3)+HHT2U(2)*TOUT(2,3) + * HHT3U(1)*TOUT(1,2)+HHT3U(2)*TOUT(2,2) + * HHTC0S(1,1)*TOUT(1,2)*TOUT(1,3)+ * HHTC0S(2,1)*TOUT(2,2)*TOUT(1,3)+ * HHTC0S(1,2)*TOUT(1,2)*TOUT(2,3)+ * HHTC0S(2,2)*TOUT(2,2)*TOUT(2,3)) FUN1K(6)=-(HHT3U(1)*TOUT(1,3)+HHT3U(2)*TOUT(2,3) + * HHT3U(1)*TOUT(1,3)+HHT3U(2)*TOUT(2,3) + * HHTC0S(1,1)*TOUT(1,3)*TOUT(1,3)+ * HHTC0S(2,1)*TOUT(2,3)*TOUT(1,3)+ * HHTC0S(1,2)*TOUT(1,3)*TOUT(2,3)+ * HHTC0S(2,2)*TOUT(2,3)*TOUT(2,3)) IFUN2K(1)=1 IFUN2K(2)=2 IFUN2K(3)=3 IFUN2K(4)=4 IFUN2K(5)=5 IFUN2K(6)=6 C Output values: DTAU(1)=0. DTAU(2)=0. DTAU(3)=0. DTAU(4)=0. DTAU(5)=0. DTAU(6)=0. DTAU(7)=0. RETURN ENDIF C Eq. 51: CALL AP28(NSUM1,TT,IX,NDER,STEP,X1T,ISRF1T,NFUN11,IFUN1T,FUN1T, * NFUN21,IFUN2T,QTKT) C Eq. 57: AUX=QQT(1)*QQT(4)-QQT(2)*QQT(3) IF (AUX.EQ.0) THEN C WAN-21 CALL ERROR('WAN-21: Zero determinant of geom. spread.') C This error should not appear. ENDIF QIT(1)= QQT(4)/AUX QIT(4)= QQT(1)/AUX QIT(2)=-QQT(2)/AUX QIT(3)=-QQT(3)/AUX TOUT(1,1)=QIT(1)*TT(1)+QIT(2)*TT(2) TOUT(2,1)=QIT(3)*TT(1)+QIT(4)*TT(2) TOUT(1,2)=QIT(1)*TT(3)+QIT(2)*TT(4) TOUT(2,2)=QIT(3)*TT(3)+QIT(4)*TT(4) TOUT(1,3)=QIT(1)*TT(5)+QIT(2)*TT(6) TOUT(2,3)=QIT(3)*TT(5)+QIT(4)*TT(6) C Eq. 62: VK(1)=-(HHT1U(1)*TOUT(1,1)+HHT1U(2)*TOUT(2,1) + * HHT1U(1)*TOUT(1,1)+HHT1U(2)*TOUT(2,1) + * HHTC0S(1,1)*TOUT(1,1)*TOUT(1,1)+HHTC0S(2,1)*TOUT(2,1)*TOUT(1,1)+ * HHTC0S(1,2)*TOUT(1,1)*TOUT(2,1)+HHTC0S(2,2)*TOUT(2,1)*TOUT(2,1)) VK(2)=-(HHT1U(1)*TOUT(1,2)+HHT1U(2)*TOUT(2,2) + * HHT2U(1)*TOUT(1,1)+HHT2U(2)*TOUT(2,1) + * HHTC0S(1,1)*TOUT(1,1)*TOUT(1,2)+HHTC0S(2,1)*TOUT(2,1)*TOUT(1,2)+ * HHTC0S(1,2)*TOUT(1,1)*TOUT(2,2)+HHTC0S(2,2)*TOUT(2,1)*TOUT(2,2)) VK(3)=-(HHT2U(1)*TOUT(1,2)+HHT2U(2)*TOUT(2,2) + * HHT2U(1)*TOUT(1,2)+HHT2U(2)*TOUT(2,2) + * HHTC0S(1,1)*TOUT(1,2)*TOUT(1,2)+HHTC0S(2,1)*TOUT(2,2)*TOUT(1,2)+ * HHTC0S(1,2)*TOUT(1,2)*TOUT(2,2)+HHTC0S(2,2)*TOUT(2,2)*TOUT(2,2)) VK(4)=-(HHT1U(1)*TOUT(1,3)+HHT1U(2)*TOUT(2,3) + * HHT3U(1)*TOUT(1,1)+HHT3U(2)*TOUT(2,1) + * HHTC0S(1,1)*TOUT(1,1)*TOUT(1,3)+HHTC0S(2,1)*TOUT(2,1)*TOUT(1,3)+ * HHTC0S(1,2)*TOUT(1,1)*TOUT(2,3)+HHTC0S(2,2)*TOUT(2,1)*TOUT(2,3)) VK(5)=-(HHT2U(1)*TOUT(1,3)+HHT2U(2)*TOUT(2,3) + * HHT3U(1)*TOUT(1,2)+HHT3U(2)*TOUT(2,2) + * HHTC0S(1,1)*TOUT(1,2)*TOUT(1,3)+HHTC0S(2,1)*TOUT(2,2)*TOUT(1,3)+ * HHTC0S(1,2)*TOUT(1,2)*TOUT(2,3)+HHTC0S(2,2)*TOUT(2,2)*TOUT(2,3)) VK(6)=-(HHT3U(1)*TOUT(1,3)+HHT3U(2)*TOUT(2,3) + * HHT3U(1)*TOUT(1,3)+HHT3U(2)*TOUT(2,3) + * HHTC0S(1,1)*TOUT(1,3)*TOUT(1,3)+HHTC0S(2,1)*TOUT(2,3)*TOUT(1,3)+ * HHTC0S(1,2)*TOUT(1,3)*TOUT(2,3)+HHTC0S(2,2)*TOUT(2,3)*TOUT(2,3)) C Eq. 59: CALL AP28(NSUM2,T,IX,NDER,STEP,X1K,ISRF1K,NFUN12,IFUN1K,FUN1K, * NFUN22,IFUN2K,VK) C Second order perturbations towards anisotropic rays: DTAU(1)=T(1)/2. DTAU(2)=T(2) DTAU(3)=T(3)/2. C Second order perturbation towards isotropic ray: DTAU(4)=T(4) DTAU(5)=T(5) DTAU(6)=T(6)/2. C Estimation of the isotropic ray travel time: DTAU(7)=TT(7) RETURN ENDIF END C C======================================================================= C SUBROUTINE WAIJKL(A,B,I) C C---------------------------------------------------------------------- REAL A(10,21),B(3,3,3,3) INTEGER I B(1,1,1,1)=A(I,1) B(1,1,2,2)=A(I,2) B(2,2,1,1)=A(I,2) B(2,2,2,2)=A(I,3) B(1,1,3,3)=A(I,4) B(3,3,1,1)=A(I,4) B(2,2,3,3)=A(I,5) B(3,3,2,2)=A(I,5) B(3,3,3,3)=A(I,6) B(1,1,2,3)=A(I,7) B(1,1,3,2)=A(I,7) B(2,3,1,1)=A(I,7) B(3,2,1,1)=A(I,7) B(2,2,2,3)=A(I,8) B(2,2,3,2)=A(I,8) B(2,3,2,2)=A(I,8) B(3,2,2,2)=A(I,8) B(3,3,3,2)=A(I,9) B(3,3,2,3)=A(I,9) B(3,2,3,3)=A(I,9) B(2,3,3,3)=A(I,9) B(2,3,2,3)=A(I,10) B(3,2,3,2)=A(I,10) B(2,3,3,2)=A(I,10) B(3,2,2,3)=A(I,10) B(1,1,1,3)=A(I,11) B(1,1,3,1)=A(I,11) B(1,3,1,1)=A(I,11) B(3,1,1,1)=A(I,11) B(2,2,1,3)=A(I,12) B(2,2,3,1)=A(I,12) B(1,3,2,2)=A(I,12) B(3,1,2,2)=A(I,12) B(3,3,3,1)=A(I,13) B(3,3,1,3)=A(I,13) B(3,1,3,3)=A(I,13) B(1,3,3,3)=A(I,13) B(3,2,1,3)=A(I,14) B(2,3,3,1)=A(I,14) B(1,3,3,2)=A(I,14) B(3,1,2,3)=A(I,14) B(3,1,3,2)=A(I,14) B(1,3,2,3)=A(I,14) B(3,2,3,1)=A(I,14) B(2,3,1,3)=A(I,14) B(1,3,1,3)=A(I,15) B(3,1,3,1)=A(I,15) B(1,3,3,1)=A(I,15) B(3,1,1,3)=A(I,15) B(1,1,1,2)=A(I,16) B(1,1,2,1)=A(I,16) B(1,2,1,1)=A(I,16) B(2,1,1,1)=A(I,16) B(2,2,2,1)=A(I,17) B(2,2,1,2)=A(I,17) B(2,1,2,2)=A(I,17) B(1,2,2,2)=A(I,17) B(3,3,1,2)=A(I,18) B(3,3,2,1)=A(I,18) B(1,2,3,3)=A(I,18) B(2,1,3,3)=A(I,18) B(2,3,1,2)=A(I,19) B(3,2,2,1)=A(I,19) B(1,2,2,3)=A(I,19) B(2,1,3,2)=A(I,19) B(2,1,2,3)=A(I,19) B(1,2,3,2)=A(I,19) B(2,3,2,1)=A(I,19) B(3,2,1,2)=A(I,19) B(1,3,2,1)=A(I,20) B(3,1,1,2)=A(I,20) B(2,1,1,3)=A(I,20) B(1,2,3,1)=A(I,20) B(1,2,1,3)=A(I,20) B(2,1,3,1)=A(I,20) B(1,3,1,2)=A(I,20) B(3,1,2,1)=A(I,20) B(1,2,1,2)=A(I,21) B(2,1,2,1)=A(I,21) B(1,2,2,1)=A(I,21) B(2,1,1,2)=A(I,21) RETURN END C C======================================================================= C REAL FUNCTION WASUM(A,B) C C---------------------------------------------------------------------- REAL A(3),B(3) INTEGER I WASUM=0. DO 10, I=1,3 WASUM=WASUM+A(I)*B(I) 10 CONTINUE RETURN END C C======================================================================= C REAL FUNCTION WASUM3(AA,B,C) C C---------------------------------------------------------------------- REAL AA(3,3),B(3),C(3) INTEGER I,J WASUM3=0. DO 20, I=1,3 DO 10, J=1,3 WASUM3=WASUM3+AA(I,J)*B(I)*C(J) 10 CONTINUE 20 CONTINUE RETURN END C C======================================================================= C REAL FUNCTION WASUM5(A,B,C,D,E) C C---------------------------------------------------------------------- REAL A(3,3,3,3),B(3),C(3),D(3),E(3) INTEGER I,J,K,L WASUM5=0. DO 13, I=1,3 DO 12, J=1,3 DO 11, K=1,3 DO 10, L=1,3 WASUM5=WASUM5+A(I,J,K,L)*B(I)*C(J)*D(K)*E(L) 10 CONTINUE 11 CONTINUE 12 CONTINUE 13 CONTINUE RETURN END C C======================================================================= C REAL FUNCTION WASUM4(A,B,C,D,J) C C---------------------------------------------------------------------- REAL A(3,3,3,3),B(3),C(3),D(3) INTEGER I,J,K,L WASUM4=0. DO 13, I=1,3 DO 11, K=1,3 DO 10, L=1,3 WASUM4=WASUM4+A(I,J,K,L)*B(I)*C(K)*D(L) 10 CONTINUE 11 CONTINUE 13 CONTINUE RETURN END C C======================================================================= C