C
C Program CRTPFA calculating two prevailing-frequency C coupling-ray-theory S-wave Green tensors along each input reference C ray. The Green tensors are continuous along the reference ray. C The input reference ray is then doubled. The first of the doubled C rays is written with the first Green tensor, the second of the doubled C rays with the second Green tensor. C C Version: 7.0 C Date: 2013, June 13 C C Coded by: Petr Bulant C Department of Geophysics, Charles University in Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C http://sw3d.cz/staff/bulant.htm C C....................................................................... C C C Description of data files: C C Main input data read from external interactive device (*): C The data consist of character strings read by list directed (free C format) input. The strings have thus to be enclosed in C apostrophes. The interactive * external unit may be redirected to C the file containing the data. C (1) 'SEP'/ C 'SEP'...Name of the file with input parameters. C Description of file SEP. C If blank, default values of the corresponding data are C considered. C No default, 'SEP' must be specified and cannot be blank. C C C Input file 'SEP' in the SEP format: C The file has the form of a SEP parameter file. C For the description of the SEP format refer to file 'sep.for'. C Names of input and output files related to ray tracing: C CRTOUT='string'...File with the names of the output files of C program CRT. Only names of the files 'CRT-R', 'CRT-S' C and 'CRT-I' are read. Description of file CRTOUT. C If CRTOUT is blank, default names are considered. C Default: CRTOUT=' ' which means 'CRT-R'='r01.out', C 'CRT-S'=' ' and 'CRT-I'='r01i.out' C (usually sufficient for the first elementary wave) C CRTNEW='string'... File with the names of the output files C to contain the new system of doubled rays. C The files 'CRT-R' and 'CRT-I' and optionally 'CRT-S' are C generated. C The form of file CRTNEW is the same as of file CRTOUT. C Default: CRTNEW=' ' which means 'CRT-R'='r01-pf.out', C 'CRT-S'=' ' and 'CRT-I'='r01i-pf.out' C MODEL='string'... Name of the input formatted file with the input C data for the anisotropic model to be used for the coupling C ray theory. C Description of file MODEL. C No default, MODEL must be specified. C C Prevailing-frequency approximation of the coupling ray theory: C SINGLF=positive real... Prevailing frequency for the C prevailing-frequency approximation of the coupling ray C theory. C No default, SINGLF must be specified. C Auxiliary parameters for anisotropic models, need not be specified: C ERRWAN=real... Maximum error in propagator matrix along whole ray. C Default: ERRWAN=0.001 C Parameters for optional second-order perturbations of anisotropic C travel times of S waves. The perturbations may be calculated only C along the isotropic or anisotropic rays which do not cross any C interface. C The perturbations can not be calculated for P waves. C Furthermore, in more complex models in the vicinity of caustics, C the second-order perturbations may be infinitely large and should C be avoided. C QIRAY=integer ... Optional second-order perturbations C of anisotropic travel times of S waves. C QIRAY=0: No second-order perturbations (default). C QIRAY=1: Second-order perturbations of anisotropic C travel times of S waves are calculated and are taken C into account when calculating the travel-time C corrections. C Default: QIRAY=0 C Parameters used for ray tracing. The values of the C parameters must be the same as used for ray tracing: C CRTANI=real... Switch to anisotropic ray tracing. C Default: CRTANI=0. C Optional parameters specifying the form of the real quantities C written in the output formatted files: C MINDIG,MAXDIG=positive integers ... See the description in file C forms.for. C C C Input formatted file CRTOUT: C For general description of file CRTOUT refer to the file writ.for. C Description specific to this program: C 'CRT-R'.. Must be specified. C 'CRT-S'.. May be specified. C 'CRT-I'.. Must be specified. C 'CRT-T'.. Not used. C Only the first data line is read in the current version. C If 'CRT-R' and 'CRT-I' only are specified, all the rays in file C 'CRT-R' are doubled. This option is suitable for example for C interpolation within ray tubes formed by the coupling rays. C If all the three 'CRT-R', 'CRT-S' and 'CRT-I' are specified, C only the rays common to files 'CRT-R' and 'CRT-S' are doubled. C This option is suitable for example for calculation of Green C function along the coupling rays. C Defaults for the first data line: C 'CRT-R'='r01.out', 'CRT-S'=' ', 'CRT-I'='r01i.out' C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL RSEP1,RSEP3I,RSEP3R,RSEP3T,ERROR,WARN,AP00,AP03,AP05, * AP06,AP21,AP28,APYYF,AN03,APW1,APW2,APW3,APW4, * PARM2,PARM3,MODEL1,VELOC,NSRFC,HIVD2,HDER, * KOOR,WAMAT,WAIJKL,WASUM,WASUM3,WASUM4,WASUM5,WACHRI, * WAPROJ,WAVECR,WAANGL,PFWAN,PFREA1,PFREA2,PFINTP,PFPIGE, * PFSHIF,PFPERT,LENGTH,SRFC2,EIGEN,VAL1,VAL2,CURVB1, * SURFB1,VAL3B1,CURVBD,SURFBD,VAL3BD INTEGER KOOR,NSRFC,LENGTH REAL WASUM,WASUM3,WASUM4,WASUM5,WAANGL C RSEP1,RSEP3I,RSEP3R,RSEP3T ... sep.for. C ERROR,WARN ... File error.for. C AP00,AP03,AP05,AP06,AP21,AP28,APYYF ... File ap.for. C AN03... File an.for. C APW1,APW2,APW3,APW4 ... File apw.for. C PARM2,PARM3 ... File parm.for. C MODEL1,VELOC,NSRFC ... File model.for. C HDER... File hder.for. C HIVD2.. File means.for. C KOOR... File metric.for. C WAMAT,WAIJKL,WASUM,WASUM3,WASUM4,WASUM5,WACHRI,WAPROJ,WAVECR, C WAANGL ... File wanpfa.for. C PFWAN,PFREA1,PFREA2,PFINTP,PFPIGE,PFSHIF,PFPERT ... This file. C LENGTH..File length.for. C EIGEN...File eigen.for. C SRFC2...File srfc.for. C VAL1,VAL2 ... File val.for. C CURVB1,SURFB1,VAL3B1,CURVBD,SURFBD,VAL3BD ... File fit.for. C C----------------------------------------------------------------------- C C Common blocks /PFPTS/ and /PFSEP/ to store quantities at points C of rays and to store the input SEP parameters: INCLUDE 'crtpfa.inc' C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER LU0,LUN INTEGER LU1,LU2,LU3,LU1N,LU2N,LU3N PARAMETER (LU0=1,LUN=2,LU1=3,LU3=4,LU1N=5,LU3N=6) C LU2 can be 0 or 7; LU2N can be 0 or 8 CHARACTER*80 FILSEP,FILMOD,FILCRT,FILNEW CHARACTER*80 FILER,FILES,FILEI,FILOR,FILOS,FILOI CHARACTER*80 FILERN,FILESN,FILEIN,FILORN,FILOSN,FILOIN INTEGER I,IL C C----------------------------------------------------------------------- C C Reading name of SEP file with input data: WRITE(*,'(A)') '+CRTPFA: Enter input filename: ' FILSEP=' ' READ(*,*) FILSEP C C Reading all data from the SEP file into the memory: IF (FILSEP.NE.' ') THEN CALL RSEP1(LU0,FILSEP) ELSE C CRTPFA-01 CALL ERROR('CRTPFA-01: SEP file not given') C Input file in the form of the SEP (Stanford Exploration Project) C parameter or history file must be specified. C There is no default filename. ENDIF C WRITE(*,'(A)') '+CRTPFA: Working ... ' C C Name of the model file CALL RSEP3T('MODEL',FILMOD,' ') C C Reading the data for a possibly anisotropic model: IF(FILMOD.EQ.' ') THEN C CRTPFA-02 CALL ERROR ('CRTPFA-02: Model not specified') C File MODEL must be specified. ELSE OPEN(LU0,FILE=FILMOD,STATUS='OLD') CALL MODEL1(LU0) CLOSE(LU0) C Checking whether the model is anisotropic CALL PARM4(I) IF(I.GT.0) THEN C Isotropic model C CRTPFA-03 CALL ERROR ('CRTPFA-03: Isotropic model') C The model specified by parameter MODEL must be anisotropic. END IF END IF IF(KOOR().NE.0) THEN C CRTPFA-04 CALL ERROR * ('CRTPFA-04: Anisotropic model in curvilinear coordinates') C Current version of subroutine file 'wan.for' for weakly C anisotropic models is coded only for Cartesian coordinates. END IF IF(NSRFC().NE.0) THEN C CRTPFA-05 CALL WARN('CRTPFA-05: Anisotropic model with interfaces') C This version of the program is not debugged for anisotropic C models with structural interfaces. END IF C C Prevailing frequency and other parameters for coupling ray theory: CALL RSEP3R('SINGLF',SINGLF,0.) IF (SINGLF.EQ.0.) THEN C CRTPFA-06 CALL ERROR ('CRTPFA-06: Zero prevailing frequency') C Prevailing frequency SINGLF must be nonzero. ENDIF CALL RSEP3R('ERRWAN',ERRWAN,0.001) CALL RSEP3I('QIRAY',KQIRAY,0) CALL RSEP3R('CRTANI',CRTANI, 0.) C C....................................................................... C Loop for the CRT output files (i.e. for the lines of file CRTOUT): FILOR=' ' FILOS=' ' FILOI=' ' FILORN=' ' FILOSN=' ' FILOIN=' ' CALL RSEP3T('CRTOUT',FILCRT,' ') CALL RSEP3T('CRTNEW',FILNEW,' ') DO 70 IL=1,1 C Recently only first line of CRTOUT is read. C....................................................................... C C Reading output filenames of the CRT program: IF(IL.EQ.1) THEN FILER='r01.out' FILES=' ' FILEI='r01i.out' FILERN='r01-pf.out' FILESN=' ' FILEIN='r01i-pf.out' ELSE FILER=' ' FILES=' ' FILEI=' ' FILERN=' ' FILESN=' ' FILEIN=' ' END IF LU2=0 LU2N=0 IF(FILCRT.NE.' ') THEN IF(IL.EQ.1) THEN OPEN(LU0,FILE=FILCRT,STATUS='OLD') END IF READ(LU0,*,END=90) FILER,FILES,FILEI END IF IF(FILER.EQ.' ') THEN C No new elementary wave GO TO 90 END IF LU2N=0 IF(FILNEW.NE.' ') THEN IF(IL.EQ.1) THEN OPEN(LUN,FILE=FILNEW,STATUS='OLD') END IF READ(LUN,*,END=4) FILERN,FILESN,FILEIN END IF IF((FILERN.EQ.' '.OR.FILEIN.EQ.' ').OR. * (FILES.NE.' '.AND.FILESN.EQ.' ')) THEN C No filenames given GO TO 4 END IF C C Opening output files of the CRT program for input: IF(FILER.NE.FILOR) THEN IF(IL.GT.1) THEN CLOSE(LU1) END IF OPEN(LU1,FILE=FILER,FORM='UNFORMATTED',STATUS='OLD',ERR=2) END IF IF((FILES.NE.' ').AND.(FILES.NE.FILOS)) THEN LU2=7 IF(IL.GT.1) THEN CLOSE(LU2) END IF OPEN(LU2,FILE=FILES,FORM='UNFORMATTED',STATUS='OLD',ERR=5) END IF IF(FILEI.NE.FILOI) THEN IF(IL.GT.1) THEN CLOSE(LU3) END IF OPEN(LU3,FILE=FILEI,FORM='UNFORMATTED',STATUS='OLD',ERR=3) END IF GO TO 9 C Error messages: 2 CONTINUE C CRTPFA-07 CALL ERROR('CRTPFA-07: Cannot open the input file with rays') C The default name of the file is 'r01.out'. 5 CONTINUE C CRTPFA-08 CALL ERROR * ('CRTPFA-08: Cannot open input file with points at surfaces') C The default name of the file is 's01.out'. 3 CONTINUE C CRTPFA-09 CALL ERROR * ('CRTPFA-09: Cannot open the input file with initial points') C The default name of the file is 'r01i.out'. 4 CONTINUE C CRTPFA-10 CALL ERROR * ('CRTPFA-10: No filenames of output files given') C File CRTOUT contains the names of input files CRT-R and CRT-I, C but file CRTNEW does not contain the names of the corresponding C input files. 9 CONTINUE C Opening output files: IF(FILERN.NE.FILORN) THEN IF(IL.GT.1) THEN CLOSE(LU1N) END IF OPEN(LU1N,FILE=FILERN,FORM='UNFORMATTED') END IF IF((FILESN.NE.' ').AND.(FILESN.NE.FILOSN)) THEN LU2N=8 IF(IL.GT.1) THEN CLOSE(LU2N) END IF OPEN(LU2N,FILE=FILESN,FORM='UNFORMATTED') END IF IF(FILEIN.NE.FILOIN) THEN IF(IL.GT.1) THEN CLOSE(LU3N) END IF OPEN(LU3N,FILE=FILEIN,FORM='UNFORMATTED') END IF C C....................................................................... C C For all the rays of the files LU1 and LU2, C reading the results of the complete ray tracing, C calculating new values and writing files LU1N and LU2N: CALL PFWAN(LU1,LU2,LU3,LU1N,LU2N,LU3N) C C....................................................................... C FILOR=FILER FILOS=FILES FILOI=FILEI FILORN=FILERN FILOSN=FILESN FILOIN=FILEIN 70 CONTINUE C End of loop for input/output files. 90 CONTINUE C C....................................................................... C CLOSE(LU1) CLOSE(LU2) CLOSE(LU1N) CLOSE(LU2N) IF(FILCRT.NE.' ') CLOSE(LU0) IF(FILNEW.NE.' ') CLOSE(LUN) WRITE(*,'(A)') '+CRTPFA: Done. ' STOP END C C======================================================================= C SUBROUTINE PFWAN(LU1,LU2,LU3,LU1N,LU2N,LU3N) C C---------------------------------------------------------------------- C C Subroutine PFWAN to convert a system of common S-wave rays in files C LU1,LU2 and LU3 into two systems of coupling-ray theory rays using the C prevailing frequency approximation and to write them into files C LU1N,LU2N and LU3N INTEGER LU1,LU2,LU3,LU1N,LU2N,LU3N C ------------------------------------------------------------------ C Input: C LU1 ... Number of the logical unit connected to the CRT output C file `CRT-R` with quantities along rays. C LU2 ... Number of the logical unit connected to the CRT output C file `CRT-S` with quantities at storing surface. C LU2=0 if the file is not to be read. C LU3 ... Number of the logical unit connected to the CRT output C file `CRT-I` with quantities at the initial points of the C rays. C LU1N... Number of the logical unit connected to the new file C `CRT-R` where the quantities along new coupling-ray-theory C rays will be written. C LU2N... Number of the logical unit connected to the new file C `CRT-S` where the quantities of new coupling-ray-theory C rays at storing surface will be written. C Not used if LU2=0. C LU3N... Number of the logical unit connected to the new file C `CRT-I` where the quantities at the initial points of the C new coupling-ray-theory rays will be written. C No output C C....................................................................... C Subroutines and external functions directly called: EXTERNAL AP00,APYYF,APW1,APW2,APW3,APW4,PFREA1,PFREA2,PFINTP, * PFPIGE,PFSHIF C AP00,APYYF ... File ap.for. C APW1,APW2,APW3,APW4 ... File apw.for. C PFREA1,PFREA2,PFINTP,PFPIGE,PFSHIF ... This file. C C....................................................................... C C Common block /POINTC/ to store the results of complete ray tracing: INCLUDE 'pointc.inc' C ........................... C C Common blocks /PFPTS/ and /PFSEP/ to store quantities at points C of rays and to store the input SEP parameters: INCLUDE 'crtpfa.inc' C ........................... C C Common block /RAM/ to store the points of the second coupling ray: INCLUDE 'ram.inc' INTEGER MISTOK,MSTOK PARAMETER (MISTOK=MRAM/2,MSTOK=MRAM/2) INTEGER ISTOK(MISTOK) REAL STOK(MSTOK) EQUIVALENCE (STOK,RAM),(ISTOK,RAM(MSTOK+1)) C C....................................................................... C C Auxiliary storage locations: INTEGER IRAY0 REAL PIGET1,PIGET2,PIGR1(2,2),PIGI1(2,2),PIGR2(2,2),PIGI2(2,2) REAL PI PARAMETER (PI=3.1415926) LOGICAL LAPYYF,LWLU2,LSTOR2 C C----------------------------------------------------------------------- C C NPTS=0 ISTOK(1)=0 LAPYYF=.FALSE. LWLU2=.FALSE. LSTOR2=.FALSE. C Loop for the points of rays 50 CONTINUE C Reading the results of the complete ray tracing IF (LU2.EQ.0) THEN CALL AP00(0,LU1,LU3) ELSE CALL AP00(LU1,LU2,LU3) ENDIF IF((IPT.LE.1).AND.(ISTOK(1).NE.0)) THEN C New ray - recording the duplication of the last ray CALL APW4(LU1N,ISTOK,STOK) LAPYYF=.FALSE. LWLU2=.FALSE. LSTOR2=.FALSE. END IF IF(IWAVE.LT.1) THEN C End of rays GO TO 100 END IF IF (LU2.NE.0) THEN C Reading and writing points O/S at storing surfaces: IF ((YF(1).LE.Y(1)).OR.LWLU2) THEN C The point O/F is before the point O/S, or the point O/S C was already written. Using point O/F: CALL APYYF LAPYYF=.TRUE. IF (YF(1).EQ.Y(1)) LSTOR2=.TRUE. ELSE C Point O/F is after point O/S, point O/S is to be used: LSTOR2=.TRUE. ENDIF ENDIF IF(IPT.LE.1)THEN C New ray - recording the initial point: CALL PFREA1 IF (Y(1).EQ.YI(1)) THEN C The first point on the ray has the same travel time C as the point in the source. This may happen for rays C with zero length (rays terminating in the source). C This point will be skipped: GOTO 50 ENDIF ENDIF 60 CONTINUE C Recording the new point: CALL PFREA2 C Calculating the Christoffel matrix and other values in the new C point(s), adding points by interpolation (if necessary): CALL PFINTP C Calculating the propagator matrix and travel-time corrections: CALL PFPIGE(PIGET1,PIGET2,PIGR1,PIGI1,PIGR2,PIGI2,DPIGD) IF(IPT.LE.1)THEN C Calculating the initial points of the first and C second coupling ray, recording them to the file CRT-I: IRAY0=IRAY IRAY=2*IRAY0-1 PVI(1)=RFEVEC(1,1) PVI(2)=RFEVEC(2,1) PVI(3)=RFEVEC(3,1) PVI(4)=RFEVEC(4,1) PVI(5)=RFEVEC(5,1) PVI(6)=RFEVEC(6,1) CALL APW1(LU3N) IRAY=2*IRAY0 CALL APW1(LU3N) IRAY=IRAY0 ENDIF IF ((LU2.EQ.0).OR.LAPYYF) THEN C Calculating points O/F on the two coupling rays, C recording them to the file CRT-R: IRAY0=IRAY C Calculating and recording the point of the first coupling ray: IRAY=2*IRAY0-1 PV(1)=RFEVEC(1,NPTS) PV(2)=RFEVEC(2,NPTS) PV(3)=RFEVEC(3,NPTS) PV(4)=RFEVEC(4,NPTS) PV(5)=RFEVEC(5,NPTS) PV(6)=RFEVEC(6,NPTS) Y(1)=PIGET1 Y(6)=Y(6)*(1.-DPIGD) Y(7)=Y(7)*(1.-DPIGD) Y(8)=Y(8)*(1.-DPIGD) Y(28)=PIGR1(1,1) Y(29)=PIGI1(1,1) Y(30)=PIGR1(2,1) Y(31)=PIGI1(2,1) Y(32)=PIGR1(1,2) Y(33)=PIGI1(1,2) Y(34)=PIGR1(2,2) Y(35)=PIGI1(2,2) CALL APW2(LU1N) C Calculating and recording the point of the second coupling ray: IRAY=2*IRAY0 Y(1)=PIGET2 Y(6)=Y(6)*(1.+DPIGD) Y(7)=Y(7)*(1.+DPIGD) Y(8)=Y(8)*(1.+DPIGD) Y(28)=PIGR2(1,1) Y(29)=PIGI2(1,1) Y(30)=PIGR2(2,1) Y(31)=PIGI2(2,1) Y(32)=PIGR2(1,2) Y(33)=PIGI2(1,2) Y(34)=PIGR2(2,2) Y(35)=PIGI2(2,2) CALL APW3(ISTOK,MISTOK,STOK,MSTOK) IRAY=IRAY0 IF (LAPYYF) THEN CALL APYYF LAPYYF=.FALSE. ENDIF ENDIF IF (LSTOR2) THEN C Calculating points O/S on the two coupling rays, C recording them to the file CRT-S: IF (LAPYYF) THEN C CRTPFA-11 CALL ERROR('CRTPFA-11: Missorder in points along ray') C Point O/F should be at or behind point O/S, so APYYF C should not have been called and LAPYYF should be .FALSE.. ENDIF IRAY0=IRAY C Calculating and recording the point of the first coupling ray: IRAY=2*IRAY0-1 PV(1)=RFEVEC(1,NPTS) PV(2)=RFEVEC(2,NPTS) PV(3)=RFEVEC(3,NPTS) PV(4)=RFEVEC(4,NPTS) PV(5)=RFEVEC(5,NPTS) PV(6)=RFEVEC(6,NPTS) Y(1)=PIGET1 Y(6)=Y(6)*(1.-DPIGD) Y(7)=Y(7)*(1.-DPIGD) Y(8)=Y(8)*(1.-DPIGD) Y(28)=PIGR1(1,1) Y(29)=PIGI1(1,1) Y(30)=PIGR1(2,1) Y(31)=PIGI1(2,1) Y(32)=PIGR1(1,2) Y(33)=PIGI1(1,2) Y(34)=PIGR1(2,2) Y(35)=PIGI1(2,2) CALL APW2(LU2N) C Calculating and recording the point of the second coupling ray: IRAY=2*IRAY0 Y(1)=PIGET2 Y(6)=Y(6)*(1.+DPIGD) Y(7)=Y(7)*(1.+DPIGD) Y(8)=Y(8)*(1.+DPIGD) Y(28)=PIGR2(1,1) Y(29)=PIGI2(1,1) Y(30)=PIGR2(2,1) Y(31)=PIGI2(2,1) Y(32)=PIGR2(1,2) Y(33)=PIGI2(1,2) Y(34)=PIGR2(2,2) Y(35)=PIGI2(2,2) CALL APW2(LU2N) IRAY=IRAY0 LWLU2=.TRUE. LSTOR2=.FALSE. IF (YF(1).GT.Y(1)) THEN C Point O/F is behind point O/S, point O/S was recently C processed, now continuing by processing of point O/F: CALL APYYF LAPYYF=.TRUE. CALL PFSHIF(NPTS,1) GOTO 60 ENDIF ENDIF IF (LAPYYF) THEN CALL APYYF LAPYYF=.FALSE. ENDIF C Shifting the values from "new" point to "old" point: CALL PFSHIF(NPTS,1) GO TO 50 C 100 CONTINUE RETURN END C C======================================================================= C SUBROUTINE PFREA1 C C---------------------------------------------------------------------- C Subroutine to store the unformatted output of program CRT. C Subroutine assumes that subroutine AP00 of file 'ap.for' was called C recently and the values are stored in arrays of common block POINTC, C this subroutine calculates the necessary quantities and stores them C into arrays of common block PFPTS. C C No input, the values are read from common block POINTC. C No output, calculated quantities are written to common block PFPTS. C C ........................... C Common block /POINTC/ to store the results of complete ray tracing: INCLUDE 'pointc.inc' C ........................... C Common blocks /PFPTS/ and /PFSEP/ to store quantities at points C of rays and to store the input SEP parameters: INCLUDE 'crtpfa.inc' C....................................................................... C C Auxiliary storage locations: INTEGER NPTSO,NPTSN,NPTSA,ISTART,ISTOP,ISTEP,ISHIFT,I1,I2,I3 REAL UP(10),US(10),RHO,QPA,QSA,VPA,VSA,VD(10),QL REAL HI(18),H(18),HUI(9) C Subroutines and external functions directly called: EXTERNAL ERROR,AP03,AN03,AP05,AP06,PARM2,VELOC C ERROR.. File error.for. C AP03,AP05,AP06 ... File ap.for. C AN03... File an.for. C PARM2.. File parm.for. C VELOC.. File model.for. C C----------------------------------------------------------------------- C C Recording initial point of a new ray: NPTS=1 C Index of a point: IFIND(NPTS)=NPTS C Travel time: RFTT(NPTS)=YI(1) C Storing the travel time in the initial point: RFTT(0)=YI(1) C Coordinates: RFCOOR(1,NPTS)=YI(3) RFCOOR(2,NPTS)=YI(4) RFCOOR(3,NPTS)=YI(5) C Slowness vector: RFSLO1(NPTS)=YI(6) RFSLO2(NPTS)=YI(7) RFSLO3(NPTS)=YI(8) C Polarization vector: RFPOL1(NPTS)=YI(9) RFPOL2(NPTS)=YI(10) RFPOL3(NPTS)=YI(11) C Index of complex block: IFICB(NPTS)=ICB1I C S-wave velocity and their derivatives: IF (CRTANI.EQ.0) THEN RFVS(NPTS)=YLI(2) RFVS1(NPTS)=YLI(4) RFVS2(NPTS)=YLI(5) RFVS3(NPTS)=YLI(6) ELSE CALL PARM2(IFICB(NPTS),RFCOOR(1,NPTS),UP,US,RHO,QPA,QSA) CALL VELOC(IFICB(NPTS),UP,US,QPA,QSA,VPA,VSA,VD,QL) RFVS(NPTS)=VD(1) RFVS1(NPTS)=VD(2) RFVS2(NPTS)=VD(3) RFVS3(NPTS)=VD(4) ENDIF C Matrix of geometrical spreading: RFQMAT(1,NPTS)=YI(12) RFQMAT(2,NPTS)=YI(13) RFQMAT(3,NPTS)=YI(16) RFQMAT(4,NPTS)=YI(17) C Transformation matrix P: RFPMAT(1,NPTS)=YI(14) RFPMAT(2,NPTS)=YI(15) RFPMAT(3,NPTS)=YI(18) RFPMAT(4,NPTS)=YI(19) C Index of the surface: IFISRF(NPTS)=0 C Covariant and contravariant components of the basis vectors C of the RCCS: CALL AP03(0,RFH(1,NPTS),H,HUI) IF (CRTANI.NE.0) THEN CALL AN03(ICB1I,YI,RFH(1,NPTS)) ENDIF RETURN C======================================================================= C ENTRY PFREA2 C C....................................................................... C C Entry designed to read the output of the CRT at other than initial C point of a ray. C NPTS=NPTS+1 IF (NPTS.GT.MPTS) THEN C CRTPFA-12 CALL ERROR('CRTPFA-12: NPTS greater than MPTS') C This error should not happen, NPTS should equal 2 here. ENDIF C Index of a point: IFIND(NPTS)=IFIND(NPTS-1)+1 C Travel time: RFTT(NPTS)=Y(1) C Coordinates: RFCOOR(1,NPTS)=Y(3) RFCOOR(2,NPTS)=Y(4) RFCOOR(3,NPTS)=Y(5) C Slowness vector: RFSLO1(NPTS)=Y(6) RFSLO2(NPTS)=Y(7) RFSLO3(NPTS)=Y(8) C Polarization vector RFPOL1(NPTS)=Y(9) RFPOL2(NPTS)=Y(10) RFPOL3(NPTS)=Y(11) C Index of complex block: IFICB(NPTS)=ICB1 C S-wave velocity and their derivatives: IF (CRTANI.EQ.0) THEN RFVS(NPTS)=YL(2) RFVS1(NPTS)=YL(4) RFVS2(NPTS)=YL(5) RFVS3(NPTS)=YL(6) ELSE CALL PARM2(IFICB(NPTS),RFCOOR(1,NPTS),UP,US,RHO,QPA,QSA) CALL VELOC(IFICB(NPTS),UP,US,QPA,QSA,VPA,VSA,VD,QL) RFVS(NPTS)=VD(1) RFVS1(NPTS)=VD(2) RFVS2(NPTS)=VD(3) RFVS3(NPTS)=VD(4) ENDIF C Matrix of geometrical spreading: CALL AP05(0,HUI, * RFQMAT(1,NPTS),RFQMAT(2,NPTS),RFQMAT(3,NPTS),RFQMAT(4,NPTS)) C Transformation matrix P: CALL AP06(0,HUI, * RFPMAT(1,NPTS),RFPMAT(2,NPTS),RFPMAT(3,NPTS),RFPMAT(4,NPTS)) C Index of the surface: IFISRF(NPTS)=ISRF C Covariant and contravariant components of the basis vectors C of the RCCS: CALL AP03(0,HI,RFH(1,NPTS),HUI) IF (CRTANI.NE.0) THEN CALL AN03(ICB1,Y,RFH(1,NPTS)) ENDIF RETURN C======================================================================= C ENTRY PFSHIF(NPTSO,NPTSN) C C....................................................................... C C Entry designed to shift the values in arrays of common block PFPTS. C The entry shifts all the values stored from NPTSO to NPTS so that C the values are stored from NPTSN. C Input: C NPTSO...Position in arrays of common block PFPTS from which the C values are moved. C NPTSN...Position in arrays of common block PFPTS to which the C values are moved. All the quantities stored in positions C NPTSO to NPTS are moved, and NPTS is then set to C NPTS+NPTSN-NPTSO. C Example 1: NPTSO=3, NPTSN=4, NPTS=10: C Values in positions 3 to 10 will be shifted to positions 4 to 11, C position 3 will be free, NPTS will be 11 on output. C Example 2: NPTSO=10, NPTSN=1, NPTS=10: C Value in position 10 will be shifted to position 1, C NPTS will be 1 on output, i.e. values originaly in positions 1 to C 9 will be erased on output. C IF (NPTSN.EQ.NPTSO) RETURN NPTSA=NPTS+NPTSN-NPTSO IF (NPTSA.GT.MPTS) THEN C CRTPFA-13 CALL ERROR('CRTPFA-13: Arrays in common block PFPTS small') C The dimension MPTS of the arrays in common block PFPTS C is given in the include file crtpfa.inc. ENDIF IF (NPTSN.GT.NPTSO) THEN ISTART=NPTSA ISTOP=NPTSO+1 ISTEP=-1 ISHIFT=-1 ELSE ISTART=NPTSN ISTOP=NPTSA ISTEP=1 ISHIFT=NPTSO-NPTSN ENDIF DO 30, I1=ISTART,ISTOP,ISTEP I2=I1+ISHIFT IFIND(I1)= IFIND(I2) IFICB(I1)= IFICB(I2) IFISRF(I1)= IFISRF(I2) RFTT(I1)= RFTT(I2) RFSLO1(I1)= RFSLO1(I2) RFSLO2(I1)= RFSLO2(I2) RFSLO3(I1)= RFSLO3(I2) RFPOL1(I1)= RFPOL1(I2) RFPOL2(I1)= RFPOL2(I2) RFPOL3(I1)= RFPOL3(I2) RFVS(I1)= RFVS(I2) RFVS1(I1)= RFVS1(I2) RFVS2(I1)= RFVS2(I2) RFVS3(I1)= RFVS3(I2) RFGS1(I1)= RFGS1(I2) RFGS2(I1)= RFGS2(I2) RFGP(I1)= RFGP(I2) RFANGL(I1)= RFANGL(I2) RFTINT(I1)= RFTINT(I2) DO 32, I3=1,4 RFQMAT(I3,I1)=RFQMAT(I3,I2) RFPMAT(I3,I1)=RFPMAT(I3,I2) 32 CONTINUE DO 34, I3=1,18 RFH(I3,I1)= RFH(I3,I2) 34 CONTINUE DO 35, I3=1,9 RFEVEC(I3,I1)= RFEVEC(I3,I2) 35 CONTINUE DO 36, I3=1,3 RFCOOR(I3,I1)=RFCOOR(I3,I2) 36 CONTINUE DO 38, I3=1,2 RFES1P(I3,I1)=RFES1P(I3,I2) RFES2P(I3,I1)=RFES2P(I3,I2) 38 CONTINUE 30 CONTINUE NPTS=NPTSA RETURN END C C======================================================================= C SUBROUTINE PFINTP C C----------------------------------------------------------------------- C Subroutine to check the angular change between two 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 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 C ........................... C Common blocks /PFPTS/ and /PFSEP/ to store quantities at points C of rays and to store the input SEP parameters: INCLUDE 'crtpfa.inc' C ........................... C External functions directly called: EXTERNAL ERROR,HIVD2,PARM2,PARM3,VELOC, * PFSHIF,WACHRI,WAPROJ,WAVECR,WAANGL REAL WAANGL C ERROR.. File error.for. C HIVD2.. File means.for. C PARM2,PARM3.. File parm.for. C VELOC.. File model.for. C PFSHIF..This file. C WACHRI,WAPROJ,WAVECR,WAANGL ... File wanpfa.for. C....................................................................... C C Auxiliary storage locations: INTEGER NPTSE,LPTISO REAL UP(10),US(10),RHO,QP,QS,VP,VS,VD(10),QL,AA(10,21),QQ(21) REAL EE(9),DER(9) REAL DEFI,DELTFI,PI PARAMETER (PI=3.1415926) INTEGER I1 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 SAVE ROLDP,TTOLD,XOLD,DXOLD,POLD,DPOLD,EOLD,DEOLD 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 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 IF (NPTS.NE.2) THEN C CRTPFA-14 CALL ERROR('CRTPFA-14: Strange number of points in PFINTP') C This error should not appear, PFINTP should be called C for a subinterval formed by two points along the ray. ENDIF C IF (IFIND(1).EQ.1) THEN C Initial point of the ray, calculating eigenvalues and C eigenvectors. C Reading the material parameters in the initial point: CALL PARM3(IFICB(1),RFCOOR(1,1),AA,RHO,QQ) C Computing eigenvectors and eigenvalues of C the Christoffel matrix in the initial point: CALL WACHRI(RFSLO1(1),RFSLO2(1),RFSLO3(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), * RFGP(1),RFGS1(1),RFGS2(1),EE,DER) C Storing the eigenvectors of the Christoffel matrix RFEVEC(1,0)=EE(4) RFEVEC(2,0)=EE(5) RFEVEC(3,0)=EE(6) RFEVEC(4,0)=EE(7) RFEVEC(5,0)=EE(8) RFEVEC(6,0)=EE(9) RFEVEC(7,0)=EE(1) RFEVEC(8,0)=EE(2) RFEVEC(9,0)=EE(3) RFEVEC(1,1)=EE(4) RFEVEC(2,1)=EE(5) RFEVEC(3,1)=EE(6) RFEVEC(4,1)=EE(7) RFEVEC(5,1)=EE(8) RFEVEC(6,1)=EE(9) RFEVEC(7,1)=EE(1) RFEVEC(8,1)=EE(2) RFEVEC(9,1)=EE(3) IF (ABS(RFGS2(1)-RFGS1(1)).LT.0.00001) THEN C Isotropic case, no projection of eigenvectors: RFES1P(1,1)=1. RFES1P(2,1)=0. RFES2P(1,1)=0. RFES2P(2,1)=1. LPTISO=1 ELSE C Projecting the qS eigenvectors to the plane C perpendicular to the ray: CALL WAPROJ(RFSLO1(1),RFSLO2(1),RFSLO3(1), * RFPOL1(1),RFPOL2(1),RFPOL3(1), * EE(4),EE(5),EE(6),EE(7),EE(8),EE(9), * RFES1P(1,1),RFES1P(2,1),RFES2P(1,1),RFES2P(2,1)) LPTISO=0 ENDIF RFANGL(1)=0. RFTINT(1)=0. ENDIF C C Quantities for future possible interpolation: TTOLD=RFTT(1) C Coordinates: XOLD(1)=RFCOOR(1,1) XOLD(2)=RFCOOR(2,1) XOLD(3)=RFCOOR(3,1) C dx/dt=p*v**2=p/sum(pi pi): AUX=(RFSLO1(1)**2+RFSLO2(1)**2+RFSLO3(1)**2) DXOLD(1)=RFSLO1(1)/AUX DXOLD(2)=RFSLO2(1)/AUX DXOLD(3)=RFSLO3(1)/AUX C Slowness: POLD(1)=RFSLO1(1) POLD(2)=RFSLO2(1) POLD(3)=RFSLO3(1) C dpi/dt=-1/v*dv/dxi: AUX=SQRT(AUX) DPOLD(1)=-RFVS1(1)/AUX DPOLD(2)=-RFVS2(1)/AUX DPOLD(3)=-RFVS3(1)/AUX C Polarization: EOLD(1)=RFPOL1(1) EOLD(2)=RFPOL2(1) EOLD(3)=RFPOL3(1) C dei/dt=v*sum(dv/dxk ek)*pi AUX=AUX* *(RFVS1(1)*RFPOL1(1)+RFVS2(1)*RFPOL2(1)+RFVS3(1)*RFPOL3(1)) DEOLD(1)=AUX*RFSLO1(1) DEOLD(2)=AUX*RFSLO2(1) DEOLD(3)=AUX*RFSLO3(1) C ROLDP( 1)=RFTT(1) ROLDP( 2)=RFCOOR(1,1) ROLDP( 3)=RFCOOR(2,1) ROLDP( 4)=RFCOOR(3,1) ROLDP( 5)=RFSLO1(1) ROLDP( 6)=RFSLO2(1) ROLDP( 7)=RFSLO3(1) ROLDP( 8)=RFPOL1(1) ROLDP( 9)=RFPOL2(1) ROLDP(10)=RFPOL3(1) ROLDP(11)=RFGP(1) ROLDP(12)=RFGS1(1) ROLDP(13)=RFGS2(1) ROLDP(14)=RFES1P(1,1) ROLDP(15)=RFES1P(2,1) ROLDP(16)=RFES2P(1,1) ROLDP(17)=RFES2P(2,1) ROLDP(18)=RFEVEC(1,1) ROLDP(19)=RFEVEC(2,1) ROLDP(20)=RFEVEC(3,1) ROLDP(21)=RFEVEC(4,1) ROLDP(22)=RFEVEC(5,1) ROLDP(23)=RFEVEC(6,1) ROLDP(24)=RFEVEC(7,1) ROLDP(25)=RFEVEC(8,1) ROLDP(26)=RFEVEC(9,1) ROLDP(27)=RFVS(1) ROLDP(28)=RFVS1(1) ROLDP(29)=RFVS2(1) ROLDP(30)=RFVS3(1) ROLDP(31)=RFQMAT(1,1) ROLDP(32)=RFQMAT(2,1) ROLDP(33)=RFQMAT(3,1) ROLDP(34)=RFQMAT(4,1) ROLDP(35)=RFPMAT(1,1) ROLDP(36)=RFPMAT(2,1) ROLDP(37)=RFPMAT(3,1) ROLDP(38)=RFPMAT(4,1) ROLDP(39)=RFH( 1,1) ROLDP(40)=RFH( 2,1) ROLDP(41)=RFH( 3,1) ROLDP(42)=RFH( 4,1) ROLDP(43)=RFH( 5,1) ROLDP(44)=RFH( 6,1) ROLDP(45)=RFH( 7,1) ROLDP(46)=RFH( 8,1) ROLDP(47)=RFH( 9,1) ROLDP(48)=RFH(10,1) ROLDP(49)=RFH(11,1) ROLDP(50)=RFH(12,1) ROLDP(51)=RFH(13,1) ROLDP(52)=RFH(14,1) ROLDP(53)=RFH(15,1) ROLDP(54)=RFH(16,1) ROLDP(55)=RFH(17,1) ROLDP(56)=RFH(18,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 CRTPFA-15 CALL ERROR('CRTPFA-15: 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 CRTPFA-16 CALL ERROR('CRTPFA-16: 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)=RFTT(I1) RNEWP( 2,1)=RFCOOR(1,I1) RNEWP( 3,1)=RFCOOR(2,I1) RNEWP( 4,1)=RFCOOR(3,I1) RNEWP( 5,1)=RFSLO1(I1) RNEWP( 6,1)=RFSLO2(I1) RNEWP( 7,1)=RFSLO3(I1) RNEWP( 8,1)=RFPOL1(I1) RNEWP( 9,1)=RFPOL2(I1) RNEWP(10,1)=RFPOL3(I1) RNEWP(27,1)=RFVS(I1) RNEWP(28,1)=RFVS1(I1) RNEWP(29,1)=RFVS2(I1) RNEWP(30,1)=RFVS3(I1) RNEWP(31,1)=RFQMAT(1,I1) RNEWP(32,1)=RFQMAT(2,I1) RNEWP(33,1)=RFQMAT(3,I1) RNEWP(34,1)=RFQMAT(4,I1) RNEWP(35,1)=RFPMAT(1,I1) RNEWP(36,1)=RFPMAT(2,I1) RNEWP(37,1)=RFPMAT(3,I1) RNEWP(38,1)=RFPMAT(4,I1) RNEWP(39,1)=RFH( 1,I1) RNEWP(40,1)=RFH( 2,I1) RNEWP(41,1)=RFH( 3,I1) RNEWP(42,1)=RFH( 4,I1) RNEWP(43,1)=RFH( 5,I1) RNEWP(44,1)=RFH( 6,I1) RNEWP(45,1)=RFH( 7,I1) RNEWP(46,1)=RFH( 8,I1) RNEWP(47,1)=RFH( 9,I1) RNEWP(48,1)=RFH(10,I1) RNEWP(49,1)=RFH(11,I1) RNEWP(50,1)=RFH(12,I1) RNEWP(51,1)=RFH(13,I1) RNEWP(52,1)=RFH(14,I1) RNEWP(53,1)=RFH(15,I1) RNEWP(54,1)=RFH(16,I1) RNEWP(55,1)=RFH(17,I1) RNEWP(56,1)=RFH(18,I1) NNEWP=1 CALL PARM3(IFICB(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=RFTT(I1) C Coordinates: XNEW(1)=RFCOOR(1,I1) XNEW(2)=RFCOOR(2,I1) XNEW(3)=RFCOOR(3,I1) C dx/dt=p*v**2=p/sum(pi pi): AUX=(RFSLO1(I1)**2+RFSLO2(I1)**2+RFSLO3(I1)**2) DXNEW(1)=RFSLO1(I1)/AUX DXNEW(2)=RFSLO2(I1)/AUX DXNEW(3)=RFSLO3(I1)/AUX C Slowness: PNEW(1)=RFSLO1(I1) PNEW(2)=RFSLO2(I1) PNEW(3)=RFSLO3(I1) C dpi/dt=-1/v*dv/dxi: AUX=SQRT(AUX) DPNEW(1)=-RFVS1(I1)/AUX DPNEW(2)=-RFVS2(I1)/AUX DPNEW(3)=-RFVS3(I1)/AUX C Polarization: ENEW(1)=RFPOL1(I1) ENEW(2)=RFPOL2(I1) ENEW(3)=RFPOL3(I1) C dei/dt=v*sum(dv/dxk ek)*pi AUX=AUX* * (RFVS1(I1)*RFPOL1(I1)+RFVS2(I1)*RFPOL2(I1)+RFVS3(I1)*RFPOL3(I1)) DENEW(1)=AUX*RFSLO1(I1) DENEW(2)=AUX*RFSLO2(I1) DENEW(3)=AUX*RFSLO3(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 (RFTT(NPTSE).EQ.RFTT(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)) DEFI=DEFI*PI*SINGLF*RFTT(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 (NNEWP.GE.MNEWP-2) THEN C CRTPFA-17 CALL ERROR('CRTPFA-17: 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 CRTPFA-18 CALL ERROR('CRTPFA-18: Array KNEWP too small') C This error should not appear, error17 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(IFICB(NPTSE+1),RNEWP(2,NNEWP),UP,US,RHO,QP,QS) CALL VELOC(IFICB(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(IFICB(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)) DEFI=DEFI*PI*SINGLF*RFTT(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 arrays RF* and IF*: CALL PFSHIF(NPTSE,NPTSE+1) C Recording interpolated quantities: IFIND( NPTSE)=0 RFTT( NPTSE)=RNEWP( 1,NNEWP) RFCOOR(1,NPTSE)=RNEWP( 2,NNEWP) RFCOOR(2,NPTSE)=RNEWP( 3,NNEWP) RFCOOR(3,NPTSE)=RNEWP( 4,NNEWP) RFSLO1( NPTSE)=RNEWP( 5,NNEWP) RFSLO2( NPTSE)=RNEWP( 6,NNEWP) RFSLO3( NPTSE)=RNEWP( 7,NNEWP) RFPOL1( NPTSE)=RNEWP( 8,NNEWP) RFPOL2( NPTSE)=RNEWP( 9,NNEWP) RFPOL3( NPTSE)=RNEWP(10,NNEWP) IFICB( NPTSE)=IFICB(NPTSE-1) RFVS( NPTSE)=RNEWP(27,NNEWP) RFVS1( NPTSE)=RNEWP(28,NNEWP) RFVS2( NPTSE)=RNEWP(29,NNEWP) RFVS3( NPTSE)=RNEWP(30,NNEWP) RFQMAT(1,NPTSE)=RNEWP(31,NNEWP) RFQMAT(2,NPTSE)=RNEWP(32,NNEWP) RFQMAT(3,NPTSE)=RNEWP(33,NNEWP) RFQMAT(4,NPTSE)=RNEWP(34,NNEWP) RFPMAT(1,NPTSE)=RNEWP(35,NNEWP) RFPMAT(2,NPTSE)=RNEWP(36,NNEWP) RFPMAT(3,NPTSE)=RNEWP(37,NNEWP) RFPMAT(4,NPTSE)=RNEWP(38,NNEWP) IFISRF( NPTSE)=0 RFH( 1,NPTSE)=RNEWP(39,NNEWP) RFH( 2,NPTSE)=RNEWP(40,NNEWP) RFH( 3,NPTSE)=RNEWP(41,NNEWP) RFH( 4,NPTSE)=RNEWP(42,NNEWP) RFH( 5,NPTSE)=RNEWP(43,NNEWP) RFH( 6,NPTSE)=RNEWP(44,NNEWP) RFH( 7,NPTSE)=RNEWP(45,NNEWP) RFH( 8,NPTSE)=RNEWP(46,NNEWP) RFH( 9,NPTSE)=RNEWP(47,NNEWP) RFH(10,NPTSE)=RNEWP(48,NNEWP) RFH(11,NPTSE)=RNEWP(49,NNEWP) RFH(12,NPTSE)=RNEWP(50,NNEWP) RFH(13,NPTSE)=RNEWP(51,NNEWP) RFH(14,NPTSE)=RNEWP(52,NNEWP) RFH(15,NPTSE)=RNEWP(53,NNEWP) RFH(16,NPTSE)=RNEWP(54,NNEWP) RFH(17,NPTSE)=RNEWP(55,NNEWP) RFH(18,NPTSE)=RNEWP(56,NNEWP) ENDIF C Recording quantities for computation of anisotropic corrections: RFGP(NPTSE)=RNEWP(11,NNEWP) RFGS1(NPTSE)=RNEWP(12,NNEWP) RFGS2(NPTSE)=RNEWP(13,NNEWP) RFES1P(1,NPTSE)=RNEWP(14,NNEWP) RFES1P(2,NPTSE)=RNEWP(15,NNEWP) RFES2P(1,NPTSE)=RNEWP(16,NNEWP) RFES2P(2,NPTSE)=RNEWP(17,NNEWP) RFANGL(NPTSE)=DELTFI RFTINT(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)) RFEVEC(1,NPTSE)=RNEWP(18,NNEWP) RFEVEC(2,NPTSE)=RNEWP(19,NNEWP) RFEVEC(3,NPTSE)=RNEWP(20,NNEWP) RFEVEC(4,NPTSE)=RNEWP(21,NNEWP) RFEVEC(5,NPTSE)=RNEWP(22,NNEWP) RFEVEC(6,NPTSE)=RNEWP(23,NNEWP) RFEVEC(7,NPTSE)=RNEWP(24,NNEWP) RFEVEC(8,NPTSE)=RNEWP(25,NNEWP) RFEVEC(9,NPTSE)=RNEWP(26,NNEWP) C IF (NPTSE.LE.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 IF (NPTSE.LT.NPTS) GOTO 5 ENDIF RETURN END C C======================================================================= C SUBROUTINE PFPIGE(PIGET1,PIGET2,PIGR1,PIGI1,PIGR2,PIGI2,DPIGD) C C----------------------------------------------------------------------- C Subroutine to calculate the propagator matrix C REAL PIGET1,PIGET2,PIGR1(2,2),PIGI1(2,2),PIGR2(2,2),PIGI2(2,2), * DPIGD C No input. C Output: C PIGET1,PIGET2 ... Travel times of two coupling rays. C PIGR1,PIGI1,PIGR2,PIGI2 ... Real and imaginary part of the C propagator matrix of the first and second coupling ray. C DPIGD.. Derivative of travel-time difference D with respect C to the reference travel time. C ........................... C Common block /POINTC/ to store the results of complete ray tracing: INCLUDE 'pointc.inc' C ........................... C Common blocks /PFPTS/ and /PFSEP/ to store quantities at points C of rays and to store the input SEP parameters: INCLUDE 'crtpfa.inc' C ........................... C External functions directly called: EXTERNAL ERROR,WARN,AP21,APYYF,PFPERT,WAMAT,LENGTH INTEGER LENGTH C ERROR,WARN ... File error.for. C AP21,APYYF ... File ap.for. C PFPERT..This file. C WAMAT.. File wanpfa.for. C LENGTH..File length.for. C....................................................................... C C Auxiliary storage locations: INTEGER I1,I2,I3,IPTTMP,ISRFFT,NYFTMP,ISINGL,NSINGL REAL PI PARAMETER (PI=3.1415926) REAL PIGR(2,2),PIGI(2,2),DPIGR(2,2),DPIGI(2,2), * PIGRA(2,2),PIGIA(2,2),DPIGRA(2,2),DPIGIA(2,2), * PIGRB(2,2),PIGIB(2,2),DPIG,C,S, * DPIGRO(2,2),DPIGIO(2,2),DPIGO, * DPIGRN(2,2),DPIGIN(2,2),DPIGN, * TTCOR,DTAU(7),DTLIN,STEP,FREQ,YI1TMP,YF1TMP,GAMA, * AUX,AUX0,AUX1,AUX2,AUX3,AUX4 CHARACTER*260 TXTERR LOGICAL LPWAVE,LSWAVE SAVE LPWAVE,LSWAVE,TTCOR,DTAU,DTLIN,STEP,FREQ,PIGR,PIGI, * DPIGR,DPIGI,DPIGRO,DPIGIO,DPIGO C----------------------------------------------------------------------- IF (IFIND(1).EQ.1) THEN C Initial point of the ray: LPWAVE=.TRUE. LSWAVE=.TRUE. TTCOR=0. DTAU(1)=0. DTAU(3)=0. DTLIN=0. DO 10, I1=1,NPTS IF (IFIND(I1).EQ.2) THEN STEP=RFTT(I1)-RFTT(1) ENDIF 10 CONTINUE 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. 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. ENDIF C Computing the values of travel time corrections along the ray: DO 30, I1=1,NPTS-1 I2=I1+1 IF (IFICB(I2).GT.0) THEN C CRTPFA-19 CALL ERROR('CRTPFA-19: P wave encountered') C P wave is not yet coded. CP TTCOR=TTCOR+(1./SQRT(RFGP(I2))+1./SQRT(RFGP(I1)))*0.5* CP * (RFTT(I2)-RFTT(I1)) CP LSWAVE=.FALSE. ELSE C S wave: TTCOR=TTCOR+(0.25/SQRT(RFGS1(I2)) +0.25/SQRT(RFGS2(I2)) * +0.25/SQRT(RFGS1(I1)) +0.25/SQRT(RFGS2(I1))) * *(RFTT(I2) -RFTT(I1)) LPWAVE=.FALSE. ENDIF 30 CONTINUE C CP IF (LPWAVE) THEN CP P wave along the ray from the source to the current point: CP Y(1)=TTCOR+YI(1) CP CALL AP21(0,GREEN) CP NQ=32 CP GOTO 91 (to the end of subroutine) CP ENDIF C IF (LSWAVE.AND.(KQIRAY.NE.0)) THEN C S wave along the ray from the source to the current point, C calculating second-order travel-time perturbations: TTCOR=TTCOR-0.5*(DTAU(3)+DTAU(1)) CALL APYYF IPTTMP=IPT YI1TMP=YI(1) NYFTMP=NYF ISRFFT=ISRFF YF1TMP=YF(1) C YI(1) =RFTT(0) NYF =1 DO 95, I1=1,NPTS IF (IFIND(I1).EQ.0) THEN STEP=0. GOTO 96 ENDIF 95 CONTINUE 96 CONTINUE IF (IFIND(1).EQ.1) THEN I2=1 ELSE I2=2 ENDIF DO 100, I1=I2,NPTS IF (I1.NE.1) DTLIN=DTLIN+RFTINT(I1) IPT =IFIND(I1)-1 ISRFF =IFISRF(I1) YF(1) =RFTT(I1) RFTINT(I1)=RFTINT(I1)+DTAU(1)-DTAU(3) CALL PFPERT(I1,STEP,DTAU) RFTINT(I1)=RFTINT(I1)+DTAU(3)-DTAU(1) 100 CONTINUE C IPT =IPTTMP YI(1) =YI1TMP NYF =NYFTMP ISRFF =ISRFFT YF(1) =YF1TMP CALL APYYF C TTCOR=TTCOR+0.5*(DTAU(3)+DTAU(1)) ELSEIF ((.NOT.LSWAVE).AND.(KQIRAY.NE.0)) THEN C CRTPFA-20 WRITE(TXTERR,'(2A,1I6,A)') * 'CRTPFA-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 Computing the propagator matrix PiGe along the ray: C Loop along points on the ray: DO 40, I1=2,NPTS I3=I1-1 IF ((RFTT(I1).EQ.RFTT(I3)).AND.(IFICB(I1).NE.IFICB(I3))) THEN C The ray crossed an interface: C CRTPFA-21 CALL ERROR('CRTPFA-21: Ray crossed an interface.') C Not yet coded. ENDIF C Element of the ray inside a complex block: IF (IFICB(I1).LT.0) THEN C S wave: GAMA=PI*FREQ*RFTINT(I1) AUX0=SQRT(RFANGL(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 AUX4=GAMA/(2.*PI*FREQ) C Matrix for this step along the ray: C R22, Eq.(49): PIGRA(1,1)= AUX1 PIGRA(2,1)=-RFANGL(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) C R22, Eq.(56): DPIGRA(1,1)=-AUX2*GAMA*AUX4 DPIGRA(2,1)=-AUX3*RFANGL(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) ELSE C P wave: GOTO 39 ENDIF C Matrix for all the steps along the ray to current point: 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) C R22, Eq.(47): CALL WAMAT(PIGRA,PIGIA,PIGR,PIGI) 39 CONTINUE 40 CONTINUE C C Single-frequency approximation C R22, Eq.(21): DPIG=DPIGR(1,1)**2+DPIGI(1,1)**2+DPIGR(2,1)**2+DPIGI(2,1)**2 DPIG=SQRT(DPIG) DPIGRN(1,1)=DPIGR(1,1) DPIGIN(1,1)=DPIGI(1,1) DPIGRN(2,1)=DPIGR(2,1) DPIGIN(2,1)=DPIGI(2,1) DPIGRN(1,2)=DPIGR(1,2) DPIGIN(1,2)=DPIGI(1,2) DPIGRN(2,2)=DPIGR(2,2) DPIGIN(2,2)=DPIGI(2,2) DPIGN=DPIG IF (DPIG.EQ.0.) THEN DPIGRN(1,1)= PIGI(1,1) DPIGIN(1,1)=-PIGR(1,1) DPIGRN(2,1)=-PIGI(2,1) DPIGIN(2,1)= PIGR(2,1) DPIGRN(1,2)= PIGI(1,2) DPIGIN(1,2)=-PIGR(1,2) DPIGRN(2,2)=-PIGI(2,2) DPIGIN(2,2)= PIGR(2,2) DPIGN=1. END IF IF (IFIND(1).NE.1) THEN AUX=DPIGRN(1,1)*DPIGRO(1,1)+DPIGIN(1,1)*DPIGIO(1,1) * +DPIGRN(2,1)*DPIGRO(2,1)+DPIGIN(2,1)*DPIGIO(2,1) DPIG= SIGN(DPIG,AUX*DPIGO) DPIGN=SIGN(DPIGN,AUX*DPIGO) ENDIF 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 C C Computing the propagator matrix: DO 80 ISINGL=0,NSINGL,32 C Transforming the propagator matrix from eigenvectors C to polarization vectors: IF (ISINGL.EQ.0) THEN C R22, Eq.(41): PIGET1=TTCOR+YI(1)-DPIG C R22, Eq.(22): PIGIB(1,1)=(PIGIA(1,1)+DPIGRN(1,1)/DPIGN)/2. PIGRB(1,1)=(PIGRA(1,1)-DPIGIN(1,1)/DPIGN)/2. PIGIB(2,1)=(PIGIA(2,1)+DPIGRN(2,1)/DPIGN)/2. PIGRB(2,1)=(PIGRA(2,1)-DPIGIN(2,1)/DPIGN)/2. PIGIB(1,2)=(PIGIA(1,2)+DPIGRN(1,2)/DPIGN)/2. PIGRB(1,2)=(PIGRA(1,2)-DPIGIN(1,2)/DPIGN)/2. PIGIB(2,2)=(PIGIA(2,2)+DPIGRN(2,2)/DPIGN)/2. PIGRB(2,2)=(PIGRA(2,2)-DPIGIN(2,2)/DPIGN)/2. C=COS(2.*PI*SINGLF*DPIG) S=SIN(2.*PI*SINGLF*DPIG) C R22, Eq.(43): PIGI1(1,1)=C*PIGIB(1,1)+S*PIGRB(1,1) PIGR1(1,1)=C*PIGRB(1,1)-S*PIGIB(1,1) PIGI1(2,1)=C*PIGIB(2,1)+S*PIGRB(2,1) PIGR1(2,1)=C*PIGRB(2,1)-S*PIGIB(2,1) PIGI1(1,2)=C*PIGIB(1,2)+S*PIGRB(1,2) PIGR1(1,2)=C*PIGRB(1,2)-S*PIGIB(1,2) PIGI1(2,2)=C*PIGIB(2,2)+S*PIGRB(2,2) PIGR1(2,2)=C*PIGRB(2,2)-S*PIGIB(2,2) ELSE C R22, Eq.(42): PIGET2=TTCOR+YI(1)+DPIG C R22, Eq.(23): PIGIB(1,1)=(PIGIA(1,1)-DPIGRN(1,1)/DPIGN)/2. PIGRB(1,1)=(PIGRA(1,1)+DPIGIN(1,1)/DPIGN)/2. PIGIB(2,1)=(PIGIA(2,1)-DPIGRN(2,1)/DPIGN)/2. PIGRB(2,1)=(PIGRA(2,1)+DPIGIN(2,1)/DPIGN)/2. PIGIB(1,2)=(PIGIA(1,2)-DPIGRN(1,2)/DPIGN)/2. PIGRB(1,2)=(PIGRA(1,2)+DPIGIN(1,2)/DPIGN)/2. PIGIB(2,2)=(PIGIA(2,2)-DPIGRN(2,2)/DPIGN)/2. PIGRB(2,2)=(PIGRA(2,2)+DPIGIN(2,2)/DPIGN)/2. S=-S C R22, Eq.(43): PIGI2(1,1)=C*PIGIB(1,1)+S*PIGRB(1,1) PIGR2(1,1)=C*PIGRB(1,1)-S*PIGIB(1,1) PIGI2(2,1)=C*PIGIB(2,1)+S*PIGRB(2,1) PIGR2(2,1)=C*PIGRB(2,1)-S*PIGIB(2,1) PIGI2(1,2)=C*PIGIB(1,2)+S*PIGRB(1,2) PIGR2(1,2)=C*PIGRB(1,2)-S*PIGIB(1,2) PIGI2(2,2)=C*PIGIB(2,2)+S*PIGRB(2,2) PIGR2(2,2)=C*PIGRB(2,2)-S*PIGIB(2,2) ENDIF 80 CONTINUE C C Calculating the length correction of the reference slowness C vector (R23 eq.37): DPIGD=(PIGIA(1,1)*DPIGRN(1,1)-PIGRA(1,1)*DPIGIN(1,1) * +PIGIA(1,2)*DPIGRN(1,2)-PIGRA(1,2)*DPIGIN(1,2))/DPIGN DPIGD=0.5*DPIGD*(1./SQRT(RFGS2(NPTS))-1./SQRT(RFGS1(NPTS))) C C Saving quantities for future invocation of the subroutine: DPIGRO(1,1)=DPIGRN(1,1) DPIGRO(2,1)=DPIGRN(2,1) DPIGRO(1,2)=DPIGRN(1,2) DPIGRO(2,2)=DPIGRN(2,2) DPIGIO(1,1)=DPIGIN(1,1) DPIGIO(2,1)=DPIGIN(2,1) DPIGIO(1,2)=DPIGIN(1,2) DPIGIO(2,2)=DPIGIN(2,2) DPIGO =DPIGN C RETURN END C C======================================================================= C SUBROUTINE PFPERT(I1,STEP,DTAU) C C---------------------------------------------------------------------- C Subroutine to compute second-order travel-time corrections INTEGER I1 REAL STEP,DTAU(7) C C Input: C I1 ... Position of the point in arrays of common block PFPTS. C STEP... Step of the independent variable along the ray, see C the CRT input parameter STORE. C Output: C DTAU... Second-order travel-time corrections. C C....................................................................... C Common blocks /PFPTS/ and /PFSEP/ to store quantities at points C of rays and to store the input SEP parameters: INCLUDE 'crtpfa.inc' C ........................... C Subroutines and external functions directly called: EXTERNAL ERROR,AP28,PARM3,HDER,WAIJKL,WASUM,WASUM3,WASUM4,WASUM5 REAL WASUM,WASUM3,WASUM4,WASUM5 C ERROR.. File error.for. C AP28... File ap.for. C PARM3.. File parm.for. C HDER... File hder.for. C WAIJKL,WASUM,WASUM3,WASUM4,WASUM5 ... File wanpfa.for. C....................................................................... C Auxiliary storage locations: 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,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,TT,IX,NDER,X1T,ISRF1T,NFUN11,IFUN1T,FUN1T,NFUN21, * IFUN2T,QTKT, * NSUM2,T,X1K,ISRF1K,NFUN12,IFUN1K,FUN1K,NFUN22,IFUN2K,VK 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 V0=RFVS(I1) VV0(1)=RFVS1(I1) VV0(2)=RFVS2(I1) VV0(3)=RFVS3(I1) P(1)=RFSLO1(I1) P(2)=RFSLO2(I1) P(3)=RFSLO3(I1) E1(1)=RFEVEC(1,I1) E1(2)=RFEVEC(2,I1) E1(3)=RFEVEC(3,I1) E2(1)=RFEVEC(4,I1) E2(2)=RFEVEC(5,I1) E2(3)=RFEVEC(6,I1) QQT(1)=RFQMAT(1,I1) QQT(2)=RFQMAT(2,I1) QQT(3)=RFQMAT(3,I1) QQT(4)=RFQMAT(4,I1) PPT(1)=RFPMAT(1,I1) PPT(2)=RFPMAT(2,I1) PPT(3)=RFPMAT(3,I1) PPT(4)=RFPMAT(4,I1) CALL PARM3(IFICB(I1),RFCOOR(1,I1),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)=RFH(1,I1) H1(2)=RFH(2,I1) H1(3)=RFH(3,I1) H2(1)=RFH(4,I1) H2(2)=RFH(5,I1) H2(3)=RFH(6,I1) H3(1)=RFH(7,I1) H3(2)=RFH(8,I1) H3(3)=RFH(9,I1) C IF (CRTANI.EQ.0.) THEN 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 (IFIND(I1).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 CRTPFA-22 CALL ERROR('CRTPFA-22: 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 HC1(1)=RFH(10,I1) HC1(2)=RFH(11,I1) HC1(3)=RFH(12,I1) HC2(1)=RFH(13,I1) HC2(2)=RFH(14,I1) HC2(3)=RFH(15,I1) HC3(1)=RFH(16,I1) HC3(2)=RFH(17,I1) HC3(3)=RFH(18,I1) 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) 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 (IFIND(I1).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 CRTPFA-23 CALL ERROR * ('CRTPFA-23: 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 CRTPFA-24 CALL ERROR('CRTPFA-24: 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 INCLUDE 'sep.for' C sep.for INCLUDE 'error.for' C error.for INCLUDE 'ap.for' C ap.for INCLUDE 'an.for' C an.for INCLUDE 'apw.for' C apw.for INCLUDE 'parm.for' C parm.for INCLUDE 'model.for' C model.for INCLUDE 'hder.for' C hder.for INCLUDE 'means.for' C means.for INCLUDE 'metric.for' C metric.for INCLUDE 'wanpfa.for' C wanpfa.for INCLUDE 'length.for' C length.for INCLUDE 'eigen.for' C eigen.for INCLUDE 'srfc.for' C srfc.for INCLUDE 'val.for' C val.for INCLUDE 'fit.for' C fit.for C C