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: 7.20 C Date: 2015, April 8 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,LUEIG) C C---------------------------------------------------------------------- INTEGER LU1,LU2,LU3,LUQI,LUEIG 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 LUEIG.. Number of the logical unit connected to the output file C EIGLST with the eigenvectors of Christoffel matrix. 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,WAANGL,WAVECR,WAPROJ, *WAPERT,WASUM,WASUM3,WASUM4,WASUM5, *ERROR,RSEP1,RSEP3T,RSEP3R,RSEP3I,FORM2,LENGTH,EIGEN, *HIVD2,BLOCK,VELOC,PARM3,AP00,AP21,TRANS,HDER INTEGER LENGTH REAL WAANGL C WAPTS,WAREAD,WACHRI,WAMAT,WAANGL,WAVECR,WAPROJ,WAPERT,WASUM, C WASUM3,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,VELOC ... File C model.for. C PARM3 ... File C parm.for. C AP00,AP21 ... File C ap.for. C TRANS ... File C trans.for. C HDER... File hder.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,SINGLF COMMON/WASEP/CRTANI,DEGREE,SINGLF SAVE /WASEP/ C CRTANI ... Switch between isotropic and anisotropic ray tracing. C DEGREE ... Degree of the homogeneous Hamiltonian. C SINGLF ... Frequency of the prevailing-frequency approximation. C....................................................................... C C Auxiliary storage locations: INTEGER MPTS REAL PIGRA(2,2),PIGIA(2,2),PIGR(2,2),PIGI(2,2) REAL PIGRB(2,2),PIGIB(2,2) REAL DPIGRA(2,2),DPIGIA(2,2),DPIGR(2,2),DPIGI(2,2) REAL AR11,AI11,AR21,AI21,AR12,AI12,AR22,AI22 REAL G11,G21,G12,G22,C,S,DPIG REAL PI PARAMETER (PI=3.1415926) REAL GREENA(32) REAL TTCOR,FREQ,GAMA,AUX0,AUX1,AUX2,AUX3,AUX4 INTEGER NPTS INTEGER I,I1,I2,I3,J,IRT,ISINGL,NSINGL INTEGER NFFT,NF,KQIPV,KQITT,KQIRAY,NF1 SAVE NFFT,NF,KQIPV,KQITT,KQIRAY REAL DT,FMIN,FMAX,DF,OF SAVE DT,FMIN,FMAX,DF,OF INTEGER IPTTMP,NYFTMP,ISRFFT REAL YI1TMP,YF1TMP,STEP,TISO,DTLIN,DTAU(7),QICOR(16),EIGOUT(9) LOGICAL LPWAVE,LSWAVE CHARACTER*135 FORMQI CHARACTER*79 FORMEI 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,' FORMEI(1:7)='(1I6,A,' CALL RSEP3I('QIRAY',KQIRAY,0) CALL RSEP3R('CRTANI',CRTANI, 0.) CALL RSEP3R('DEGREE',DEGREE,-1.) CALL RSEP3R('SINGLF',SINGLF, 0.) 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 C Replacing the reference polarization vectors in common block C /POINTC/ by the eigenvectors of the Christoffel matrix: DO 9 I=1,9 PVI(I)=RPTS(24+I,1) PV (I)=RPTS(24+I,NPTS) 9 CONTINUE C C Output to the file EIGLST: IF (LUEIG.NE.0) THEN WRITE(LUEIG,'(A,1I6,A)') 'IREC =',IREC,' /' DO 10, I1=1,NPTS EIGOUT(1)=RPTS( 2,I1) EIGOUT(2)=RPTS(25,I1) EIGOUT(3)=RPTS(26,I1) EIGOUT(4)=RPTS(27,I1) EIGOUT(5)=RPTS(28,I1) EIGOUT(6)=RPTS(29,I1) EIGOUT(7)=RPTS(30,I1) EIGOUT(8)=RPTS(14,I1) EIGOUT(9)=RPTS(15,I1) CALL FORM2(9,EIGOUT,EIGOUT,FORMEI(8:79)) WRITE(LUEIG,FORMEI) IPTS(1,I1),(' ',EIGOUT(I),I=1,9),' /' 10 CONTINUE 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. C Updating common block /POINTC/ for subroutine AP21: NPVI=3 NPV=3 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. C Updating common block /POINTC/ for subroutine AP21: NPVI=6 NPV=6 ENDIF 30 CONTINUE C IF (LPWAVE) THEN C P wave along the whole ray: Y(1)=TTCOR+YI(1) CALL AP21(KQIPV,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: NF1=NF-1 IF (SINGLF.NE.0.) THEN NF1=0 ENDIF DO 90, I2=0,NF1 FREQ=OF+I2*DF IF (SINGLF.NE.0.) THEN FREQ=SINGLF DPIGR(1,1)=0. DPIGR(2,1)=0. DPIGR(1,2)=0. DPIGR(2,2)=0. DPIGI(1,1)=0. DPIGI(2,1)=0. DPIGI(1,2)=0. DPIGI(2,2)=0. ENDIF 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: C WAN-17 CALL ERROR('WAN-17: Ray crossed an interface.') C Not yet coded. * IRT=IRT+1 * IF (IRT.GT.NRT) THEN C WAN-XX * CALL ERROR('WAN-XX: 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. AUX3=-1./3. ELSE AUX2=SIN(AUX0)/AUX0 AUX3=(COS(AUX0)-AUX2)/AUX0**2 ENDIF C Matrix for this step along the ray: C R22, Eq.(49): 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) IF (SINGLF.NE.0.) THEN C R22, Eq.(56): AUX4=GAMA/(2.*PI*FREQ) DPIGRA(1,1)=-AUX2*GAMA*AUX4 DPIGRA(2,1)=-AUX3*RPTS(20,I1)*GAMA*AUX4 DPIGRA(1,2)=-DPIGRA(2,1) DPIGRA(2,2)= DPIGRA(1,1) DPIGIA(1,1)=-(AUX3*GAMA**2+AUX2)*AUX4 DPIGIA(2,1)= 0. DPIGIA(1,2)= 0. DPIGIA(2,2)=-DPIGIA(1,1) ENDIF ELSE C P wave: GOTO 39 ENDIF ENDIF C Matrix for all the steps along the ray to current point: IF (SINGLF.NE.0.) THEN C R22, Eq.(53): CALL WAMAT(PIGRA,PIGIA,DPIGR,DPIGI) PIGRB(1,1)=PIGR(1,1) PIGRB(2,1)=PIGR(2,1) PIGRB(1,2)=PIGR(1,2) PIGRB(2,2)=PIGR(2,2) PIGIB(1,1)=PIGI(1,1) PIGIB(2,1)=PIGI(2,1) PIGIB(1,2)=PIGI(1,2) PIGIB(2,2)=PIGI(2,2) CALL WAMAT(DPIGRA,DPIGIA,PIGRB,PIGIB) DPIGR(1,1)=DPIGR(1,1)+PIGRB(1,1) DPIGR(2,1)=DPIGR(2,1)+PIGRB(2,1) DPIGR(1,2)=DPIGR(1,2)+PIGRB(1,2) DPIGR(2,2)=DPIGR(2,2)+PIGRB(2,2) DPIGI(1,1)=DPIGI(1,1)+PIGIB(1,1) DPIGI(2,1)=DPIGI(2,1)+PIGIB(2,1) DPIGI(1,2)=DPIGI(1,2)+PIGIB(1,2) DPIGI(2,2)=DPIGI(2,2)+PIGIB(2,2) ENDIF C R22, Eq.(47): CALL WAMAT(PIGRA,PIGIA,PIGR,PIGI) 39 CONTINUE 40 CONTINUE C C Prevailing-frequency approximation NSINGL=0 IF (SINGLF.NE.0.) THEN C R22, Eq.(21): DPIG=DPIGR(1,1)**2+DPIGI(1,1)**2+DPIGR(2,1)**2+DPIGI(2,1)**2 DPIG=SQRT(DPIG) PIGIA(1,1)=PIGI(1,1) PIGRA(1,1)=PIGR(1,1) PIGIA(2,1)=PIGI(2,1) PIGRA(2,1)=PIGR(2,1) PIGIA(1,2)=PIGI(1,2) PIGRA(1,2)=PIGR(1,2) PIGIA(2,2)=PIGI(2,2) PIGRA(2,2)=PIGR(2,2) NSINGL=32 ENDIF C C Computing the Green function: DO 80 ISINGL=0,NSINGL,32 C Transforming the propagator matrix from eigenvectors C to polarization vectors: Y(1)=TTCOR+YI(1) IF (SINGLF.NE.0.) THEN IF (ISINGL.EQ.0) THEN C R22, Eq.(41): Y(1)=TTCOR+YI(1)-DPIG C R22, Eq.(22): IF (DPIG.NE.0) THEN PIGIB(1,1)=(PIGIA(1,1)+DPIGR(1,1)/DPIG)/2. PIGRB(1,1)=(PIGRA(1,1)-DPIGI(1,1)/DPIG)/2. PIGIB(2,1)=(PIGIA(2,1)+DPIGR(2,1)/DPIG)/2. PIGRB(2,1)=(PIGRA(2,1)-DPIGI(2,1)/DPIG)/2. PIGIB(1,2)=(PIGIA(1,2)+DPIGR(1,2)/DPIG)/2. PIGRB(1,2)=(PIGRA(1,2)-DPIGI(1,2)/DPIG)/2. PIGIB(2,2)=(PIGIA(2,2)+DPIGR(2,2)/DPIG)/2. PIGRB(2,2)=(PIGRA(2,2)-DPIGI(2,2)/DPIG)/2. ELSE PIGIB(1,1)=PIGIA(1,1) PIGRB(1,1)=PIGRA(1,1) PIGIB(2,1)=PIGIA(2,1) PIGRB(2,1)=PIGRA(2,1) PIGIB(1,2)=0. PIGRB(1,2)=0. PIGIB(2,2)=0. PIGRB(2,2)=0. ENDIF C=COS(2.*PI*SINGLF*DPIG) S=SIN(2.*PI*SINGLF*DPIG) ELSE C R22, Eq.(42): Y(1)=TTCOR+YI(1)+DPIG C R22, Eq.(23): IF (DPIG.NE.0) THEN PIGIB(1,1)=(PIGIA(1,1)-DPIGR(1,1)/DPIG)/2. PIGRB(1,1)=(PIGRA(1,1)+DPIGI(1,1)/DPIG)/2. PIGIB(2,1)=(PIGIA(2,1)-DPIGR(2,1)/DPIG)/2. PIGRB(2,1)=(PIGRA(2,1)+DPIGI(2,1)/DPIG)/2. PIGIB(1,2)=(PIGIA(1,2)-DPIGR(1,2)/DPIG)/2. PIGRB(1,2)=(PIGRA(1,2)+DPIGI(1,2)/DPIG)/2. PIGIB(2,2)=(PIGIA(2,2)-DPIGR(2,2)/DPIG)/2. PIGRB(2,2)=(PIGRA(2,2)+DPIGI(2,2)/DPIG)/2. ELSE PIGIB(1,1)=0. PIGRB(1,1)=0. PIGIB(2,1)=0. PIGRB(2,1)=0. PIGIB(1,2)=PIGIA(1,2) PIGRB(1,2)=PIGRA(1,2) PIGIB(2,2)=PIGIA(2,2) PIGRB(2,2)=PIGRA(2,2) ENDIF S=-S ENDIF C R22, Eq.(43): PIGI(1,1)=C*PIGIB(1,1)+S*PIGRB(1,1) PIGR(1,1)=C*PIGRB(1,1)-S*PIGIB(1,1) PIGI(2,1)=C*PIGIB(2,1)+S*PIGRB(2,1) PIGR(2,1)=C*PIGRB(2,1)-S*PIGIB(2,1) PIGI(1,2)=C*PIGIB(1,2)+S*PIGRB(1,2) PIGR(1,2)=C*PIGRB(1,2)-S*PIGIB(1,2) PIGI(2,2)=C*PIGIB(2,2)+S*PIGRB(2,2) PIGR(2,2)=C*PIGRB(2,2)-S*PIGIB(2,2) ENDIF C 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,GREENA) IF (I2.EQ.0) THEN DO 50 I=1,14 GREEN(ISINGL+I)=GREENA(I) 50 CONTINUE NQ=ISINGL+14 ENDIF J=I2*18 DO 60 I=15,32 GREEN(ISINGL+J+I)=GREENA(I) 60 CONTINUE NQ=NQ+18 80 CONTINUE 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 ........................... C Common block /WASEP/ to store the input SEP parameters used C in subroutine WAN and their subroutines. REAL CRTANI,DEGREE,SINGLF COMMON/WASEP/CRTANI,DEGREE,SINGLF SAVE /WASEP/ C CRTANI ... Switch between isotropic and anisotropic ray tracing. C DEGREE ... Degree of the homogeneous Hamiltonian. C SINGLF ... Frequency of the prevailing-frequency approximation. C....................................................................... C External functions required: EXTERNAL WAVECR,WAANGL REAL WAANGL C....................................................................... C C Auxiliary storage locations: INTEGER NPTSE,LPTISO INTEGER KQIDT SAVE 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 SAVE ERRWAN PARAMETER (PI=3.1415926) INTEGER I1,I2 INTEGER MQANT,MNEWP,NNEWP PARAMETER (MNEWP=50) 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 LPTISO..LPTISO=1 in case that all already checked points were C isotropic. 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. LPTISO=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)) LPTISO=0 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: Point along ray missing') C Here we assume that we are checking the interval of points C NPTSE to NPTSE+1, and thus NPTS should be at least NPTSE+1 C or more. This error should not appear. 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 (LPTISO.EQ.1) THEN C New point is the first anisotropic point on the ray, C angular change is not to be computed: DELTFI=0. LPTISO=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 CALL WAVECR(ROLDP(18),ROLDP(21),RNEWP(18,I1),RNEWP(21,I1), * RNEWP(12,I1),RNEWP(13,I1)) DELTFI=WAANGL(ROLDP(18),ROLDP(21),RNEWP(18,I1),RNEWP(21,I1)) DEFI=1./SQRT(RNEWP(12,I1))-1./SQRT(RNEWP(13,I1)) DEFI=DEFI-1./SQRT(ROLDP(12))+1./SQRT(ROLDP(13)) IF (SINGLF.NE.0.) THEN DEFI=DEFI*PI*SINGLF*RPTS(2,NPTS)/6. ELSE DEFI=DEFI*PI*(OF+FLOAT(NF-1)*DF)*RPTS(2,NPTS)/6. ENDIF 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 (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 CALL WAVECR(ROLDP(18),ROLDP(21), * RNEWP(18,NNEWP),RNEWP(21,NNEWP), * RNEWP(12,NNEWP),RNEWP(13,NNEWP)) DELTFI=WAANGL(ROLDP(18),ROLDP(21), * RNEWP(18,NNEWP),RNEWP(21,NNEWP)) DEFI=1./SQRT(RNEWP(12,NNEWP))-1./SQRT(RNEWP(13,NNEWP)) DEFI=DEFI-1./SQRT(ROLDP(12))+1./SQRT(ROLDP(13)) IF (SINGLF.NE.0.) THEN DEFI=DEFI*PI*SINGLF*RPTS(2,NPTS)/6. ELSE DEFI=DEFI*PI*(OF+FLOAT(NF-1)*DF)*RPTS(2,NPTS)/6. ENDIF 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,SINGLF COMMON/WASEP/CRTANI,DEGREE,SINGLF SAVE /WASEP/ C CRTANI ... Switch between isotropic and anisotropic ray tracing. C DEGREE ... Degree of the homogeneous Hamiltonian. C SINGLF ... Frequency of the prevailing-frequency approximation. C....................................................................... C C Auxiliary storage locations: INTEGER ISRFA,ISRF2,KODNEW,ICBNEW,IEN,ISB2A,ISB2B,ICB2A,ICB2B,I1 INTEGER IYFA(13) 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 IYFA(13)=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 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 C Subroutines and external functions required: EXTERNAL HDER,WASUM,WASUM3,WASUM4,WASUM5 REAL WASUM,WASUM3,WASUM4,WASUM5 C HDER... File hder.for. 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,SINGLF COMMON/WASEP/CRTANI,DEGREE,SINGLF SAVE /WASEP/ C CRTANI ... Switch between isotropic and anisotropic ray tracing. C DEGREE ... Degree of the homogeneous Hamiltonian. C SINGLF ... Frequency of the prevailing-frequency approximation. C C....................................................................... C INTEGER IERR 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 REAL EE(9) 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 EE... Eigenvectors of the Christoffel matrix. Not used here. C 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,EE,IERR) 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