C
C Subroutine file 'raycb.for' for complete ray tracing within one C complex block C C Version: 7.10 C Date: 2014, February 16 C Coded by: Ludek Klimes C C....................................................................... C C This file consists of the following subroutines: C RAYCB...Subroutine transferring the quantities given at the C initial point of the element of a ray into the C corresponding quantities at the endpoint of the element, C see C C.R.T.5.8.1. C RAYCB C FCT... Subroutine evaluating the right-hand sides of the system C of ordinary differential equations for Y(1) to Y(27), see C C.R.T.5.8.2. C FCT C FCTA... Subroutine evaluating the right-hand sides of the system C of ordinary differential equations for Y(1) to Y(27), see C C.R.T.5.8.2. C for anisotropic ray tracing. C FCTA C OUTP... Subroutine containing various tests of the position of the C newly computed point of the ray, tests for possible C caustic points, etc., see C C.R.T.5.8.3. C OUTP C CDEF... Auxiliary subroutine to OUTP, performing steps C 5.8.3c, d, C e and f C of the algorithm. C CDEF C PSHIFT..Auxiliary subroutine to OUTP and CDEF. It corrects reduced C amplitudes with regard to phase shift due to caustics, see C 5.8.3f. C PSHIFT C PVINT...Auxiliary subroutine to OUTP. It calculates the reference C polarization vectors in anisotropic media. C It interpolates the P-wave reference polarization vector C from points X1 and X2 to point X, and then calculates the C S-wave reference polarization vectors at point X from C point X1. C C======================================================================= C C C SUBROUTINE RAYCB(YL,Y,YY,IY) REAL YL(6),Y(44),YY(5) INTEGER IY(13) C C This subroutine transfers the quantities given at the initial point of C the element of a ray into the corresponding quantities at the endpoint C of the element, see C C.R.T.5.8.1. C C Input: C YL... Array containing local quantities at the initial point of C the element of the ray. C Y... Array containing basic quantities computed along the ray C at the initial point of the element of the ray. C YY... Array containing real auxiliary quantities computed along C the ray at the initial point of the element of the ray. C IY... Array containing integer auxiliary quantities computed C along the ray at the initial point of the element of the C ray. C C Output: C YL... Array containing local quantities at the endpoint of the C element of the ray. C For the description refer to file ray.for. C Y... Array containing basic quantities computed along the ray C at the endpoint of the element of the ray. C For the description refer to file ray.for. C YY... Array containing real auxiliary quantities computed along C the ray at the endpoint of the element of the ray. C For the description refer to file ray.for. C IY... Array containing integer auxiliary quantities computed C along the ray at the endpoint of the element of the ray. C For the description refer to file ray.for. C C Common block /DCRT/ (see subroutine file 'ray.for'): INCLUDE 'dcrt.inc' C dcrt.inc C None of the storage locations of the common block are altered. C C Common block /INITC/ (see subroutine file 'init.for') is required in C subroutine OUTP. None of the storage locations of the common C block are altered there. C C Common block /RAYC/ - communication between RAYCB, FCT, and OUTP: C ------------------------------------------------------------------ INCLUDE 'raycb.inc' C raycb.inc C ------------------------------------------------------------------ C YLRC,YYRC,IYRC... Working copies of the arrays YL,YY,YI. C DELT11,DELT13,DELT33... Deviations (in absolute values of the two C computed basis vectors of the ray-centred coordinate C system from the conditions of orthonormality, see C 5.8.2d. C YLRC,YYRC,IYRC are defined in this subroutine. C C Subroutines and external functions required: EXTERNAL KOOR,METRIC,HPCG INTEGER KOOR C KOOR,METRIC... File 'metric.for' of the package 'MODEL'. C NSRFC,BLOCK,VELOC... File 'model.for' of the package 'MODEL'. C SRFC2 and subsequent routines... File 'srfc.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C PARM2 and subsequent routines... File 'parm.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C CDE,CROSS,HIVD2,SMVPRD... File 'means.for' of the package 'MODEL'. C HPCG... IBM Scientific Subroutine Package (file 'hpcg.for' of the C package 'MODEL'). C RPAR31,RPAR32... File 'rpar.for'. C WRIT31,WRIT32... File 'writ.for'. C FCT,OUTP,CDEF,PSHIFT... This file. C C Date: 2014, February 4 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations for local model parameters: FAUX(10), C G(12),GAMMA(18),GSQRD, UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL: INCLUDE 'auxmod.inc' C auxmod.inc C C....................................................................... C C External procedure names: EXTERNAL FCT,FCTA,OUTP C These subroutines are called by the subroutine HPCG. C C Auxiliary storage locations: INTEGER NDIM,IHLF,I PARAMETER (NDIM=27) REAL PRMT(6),DERY(NDIM),AUX(16,NDIM),VRED,AUX1 C C NDIM,IHLF,PRMT,DERY,AUX... See subroutine HPCG of the IBM C Scientific Subroutine Package. C I... Auxiliary loop variable. C VRED... Approximate reference velocity to estimate error weights. C AUX1... Auxiliary storage location. C C....................................................................... C C Anisotropic ray tracing: REAL CRTANI SAVE CRTANI DATA CRTANI/-999./ IF(CRTANI.EQ.-999.) THEN CALL RSEP3R('CRTANI',CRTANI,0.) IF(CRTANI.NE.0..AND.KOOR().NE.0) THEN C 5A1 CALL ERROR('5A1 in RAYCB: Anisotropy in curvilinear coord.') C Anisotropic ray tracing is now applicable only to Cartesian C coordinates. END IF END IF C C (a) Storing auxiliary input quantities in common block /RAYC/: DO 11 I=1,6 YLRC(I)=YL(I) 11 CONTINUE DO 12 I=1,5 YYRC(I)=YY(I) 12 CONTINUE DO 13 I=1,12 IYRC(I)=IY(I) 13 CONTINUE C C (b) Initial values for numerical integration: 20 PRMT(1)=YYRC(1) PRMT(2)=PRMT(1)+999999. PRMT(3)=STEP PRMT(4)=13.444*YYRC(2) CALL METRIC(Y(3),GSQRD,G,GAMMA) IF(IYRC(5).GE.0) THEN VRED=YLRC(1) ELSE VRED=YLRC(2) END IF DERY(1)=1. DERY(2)=1. DERY(3)=SQRT(G(1))/VRED DERY(4)=SQRT(G(3))/VRED DERY(5)=SQRT(G(6))/VRED AUX1=STEP*VRED**(1-NEXPS) DERY(6)=SQRT(G(7))*AUX1 DERY(7)=SQRT(G(9))*AUX1 DERY(8)=SQRT(G(12))*AUX1 DO 21 I=9,NDIM DERY(I)=0. 21 CONTINUE C C (c) Numerical integration: IF(CRTANI.EQ.0.) THEN IYRC(13)=0 CALL HPCG(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) ELSE IF(IY(5).GT.0) THEN IYRC(13)=3 ELSE IYRC(13)=6 END IF CALL HPCG(PRMT,Y,DERY,NDIM,IHLF,FCTA,OUTP,AUX) END IF C C (d) Great number of bisections: IF(PRMT(5).LT.0.) THEN YYRC(2)=2.*YYRC(2) GO TO 20 END IF C C (e) Recalling auxiliary input quantities from common block /RAYC/: DO 51 I=1,6 YL(I)=YLRC(I) 51 CONTINUE DO 52 I=1,5 YY(I)=YYRC(I) 52 CONTINUE DO 53 I=1,13 IY(I)=IYRC(I) 53 CONTINUE C RETURN END C C======================================================================= C C C SUBROUTINE FCT(X,Y,D) INTEGER NDIM PARAMETER (NDIM=27) REAL X,Y(NDIM),D(NDIM) C C This subroutine evaluates the right-hand sides of the system of C ordinary differential equations for Y(1) to Y(27). See C C.R.T.5.8.2. C C Input: C X... Value of the independent variable along the ray. C Y... Array containing basic quantities computed along the ray. C None of the input parameters except Y(3)-Y(8) is altered. C C Output: C D... Array containing derivatives of the basic quantities C computed along the ray, with respect to the independent C variable along the ray. C Y(3:8)..Renormalized orthogonal vectors. The modification should C be negligible. C C Common block /DCRT/ (see subroutine file 'ray.for'): INCLUDE 'dcrt.inc' C dcrt.inc C None of the storage locations of the common block are altered. C C Common block /RAYC/ - communication between RAYCB, FCT, and OUTP: INCLUDE 'raycb.inc' C raycb.inc C YLRC,IYRC are modified in this subroutine. DELT11,DELT13, and DELT33 C are defined in this subroutine. C C Subroutines and external functions required: EXTERNAL KOOR,METRIC,PARM2,VELOC,SMVPRD INTEGER KOOR C VELOC... File 'model.for' of the package 'MODEL'. C KOOR,METRIC... File 'metric.for' of the package 'MODEL'. C PARM2 and subsequent routines... File 'parm.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C SMVPRD... File 'means.for' of the package 'MODEL'. C C Date: 2003, May 5 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations for local model parameters: FAUX(10), C G(12),GAMMA(18),GSQRD, UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL: INCLUDE 'auxmod.inc' C auxmod.inc C C....................................................................... C C Auxiliary storage locations: C REAL H11,H21,H31,H41,H51,H61,H42,H52,H62,H13,H23,H33,H43,H53,H63 REAL D11,D13,D33,V,V1,V11,V12,V22,VNEXPS,AUX1,AUX2,AUX3,AUX4 REAL AUX11,AUX21,AUX31,AUX12,AUX22,AUX32,AUX13,AUX23,AUX33 C C H11,H21,H31,H13,H23,H33... Covariant components of the basis C vectors of ray-centred coordinate system. C H41,H51,H61,H42,H52,H62,H43,H53,H63... Contravariant components of C the basis vectors of ray-centred coordinate system. H41, C H51,H61,H43,H53,H63 are not used in Cartesian coordinates. C D11,D13,D33... Norms and scalar product of basis vectors H1, H2. C V,V1,V11,V12,V22,... Velocity and velocity derivatives in C ray-centred coordinate system. C VNEXPS..Velocity**(-NEXPS). C AUX1,AUX2,AUX3,AUX4... Auxiliary storage locations. C AUX11,AUX21,AUX31,AUX12,AUX22,AUX32,AUX13,AUX23,AUX33... Auxiliary C storage locations. Not used in Cartesian coordinates. C C....................................................................... C C (a) Number of invocations of FCT: IYRC(9)=IYRC(9)+1 C C (b) Material parameters: CALL PARM2(IABS(IYRC(5)),Y(3),UP,US,YLRC(3),QP,QS) CALL VELOC(IYRC(5),UP,US,QP,QS,YLRC(1),YLRC(2),VD,QL) YLRC(4)=VD(2) YLRC(5)=VD(3) YLRC(6)=VD(4) V=VD(1) VNEXPS=V**(-NEXPS) C C (c) Basis of the ray-centred coordinate system: C (c1) Non-unit vectors - covariant components H13=Y(6) H23=Y(7) H33=Y(8) H11=Y(9) H21=Y(10) H31=Y(11) IF(KOOR().NE.0) THEN C Curvilinear coordinates: CALL METRIC(Y(3),GSQRD,G,GAMMA) C (c2) non-unit vectors - contravariant components CALL SMVPRD(G(7),H13,H23,H33,H43,H53,H63) CALL SMVPRD(G(7),H11,H21,H31,H41,H51,H61) C (c3) Norms: D33=SQRT(H13*H43+H23*H53+H33*H63) D11=SQRT(H11*H41+H21*H51+H31*H61) D13=(H11*H43+H21*H53+H31*H63)/(D11*D33) C (c4) Orthonormal vectors: C Covariant components H13=H13/D33 H23=H23/D33 H33=H33/D33 H11=H11/D11-H13*D13 H21=H21/D11-H23*D13 H31=H31/D11-H33*D13 C Contravariant components H43=H43/D33 H53=H53/D33 H63=H63/D33 H41=H41/D11-H43*D13 H51=H51/D11-H53*D13 H61=H61/D11-H63*D13 H42=(H23*H31-H33*H21)/GSQRD H52=(H33*H11-H13*H31)/GSQRD H62=(H13*H21-H23*H11)/GSQRD C C (d) Quantities useful to test the accuracy of computations: DELT11=ABS(D11-1.) DELT33=ABS(V*D33-1.) DELT13=ABS(D13) C C (e) First and second derivatives of the velocity in ray-centred C coordinate system: C Second covariant velocity derivatives: VD(5)=VD(5)-GAMMA(1)*VD(2)-GAMMA(7)*VD(3)-GAMMA(13)*VD(4) VD(6)=VD(6)-GAMMA(2)*VD(2)-GAMMA(8)*VD(3)-GAMMA(14)*VD(4) VD(7)=VD(7)-GAMMA(3)*VD(2)-GAMMA(9)*VD(3)-GAMMA(15)*VD(4) VD(8)=VD(8)-GAMMA(4)*VD(2)-GAMMA(10)*VD(3)-GAMMA(16)*VD(4) VD(9)=VD(9)-GAMMA(5)*VD(2)-GAMMA(11)*VD(3)-GAMMA(17)*VD(4) VD(10)=VD(10)-GAMMA(6)*VD(2)-GAMMA(12)*VD(3)-GAMMA(18)*VD(4) C Ray-centred coordinate system: V1=VD(2)*H41+VD(3)*H51+VD(4)*H61 CALL SMVPRD(VD(5),H41,H51,H61,AUX1,AUX2,AUX3) V11=AUX1*H41+AUX2*H51+AUX3*H61 V12=AUX1*H42+AUX2*H52+AUX3*H62 CALL SMVPRD(VD(5),H42,H52,H62,AUX1,AUX2,AUX3) V22=AUX1*H42+AUX2*H52+AUX3*H62 C C (f) Correction of the computed quantities: Y(6)=H13/V Y(7)=H23/V Y(8)=H33/V Y(9)=H11 Y(10)=H21 Y(11)=H31 C C (g) Right-hand sides of the differential equations: AUX1=V*VNEXPS AUX2=V*AUX1 AUX3=VNEXPS/V AUX4=V1*VNEXPS D(1)=VNEXPS D(2)=0.5*QL*VNEXPS D(3)=H43*AUX1 D(4)=H53*AUX1 D(5)=H63*AUX1 CALL SMVPRD(GAMMA(1),H43,H53,H63,AUX11,AUX21,AUX31) CALL SMVPRD(GAMMA(7),H43,H53,H63,AUX12,AUX22,AUX32) CALL SMVPRD(GAMMA(13),H43,H53,H63,AUX13,AUX23,AUX33) D(6)=-VD(2)*AUX3+(AUX11*H13+AUX12*H23+AUX13*H33)*VNEXPS D(7)=-VD(3)*AUX3+(AUX21*H13+AUX22*H23+AUX23*H33)*VNEXPS D(8)=-VD(4)*AUX3+(AUX31*H13+AUX32*H23+AUX33*H33)*VNEXPS D(9) =H13*AUX4+(AUX11*H11+AUX12*H21+AUX13*H31)*AUX1 D(10)=H23*AUX4+(AUX21*H11+AUX22*H21+AUX23*H31)*AUX1 D(11)=H33*AUX4+(AUX31*H11+AUX32*H21+AUX33*H31)*AUX1 ELSE C Cartesian coordinates (20 of these statement lines are the same C as for curvilinear coordinates): C (c2) Non-unit vectors - contravariant components C H43=H13, H53=H23, H63=H33, H41=H11, H51=H21, H61=H31 C (c3) Norms: D33=SQRT(H13*H13+H23*H23+H33*H33) D11=SQRT(H11*H11+H21*H21+H31*H31) D13=(H11*H13+H21*H23+H31*H33)/(D11*D33) C (c4) Orthonormal vectors: C Covariant=contravariant components H13=H13/D33 H23=H23/D33 H33=H33/D33 H11=H11/D11-H13*D13 H21=H21/D11-H23*D13 H31=H31/D11-H33*D13 H42=(H23*H31-H33*H21) H52=(H33*H11-H13*H31) H62=(H13*H21-H23*H11) C C (d) Quantities useful to test the accuracy of computations: DELT11=ABS(D11-1.) DELT33=ABS(V*D33-1.) DELT13=ABS(D13) C C (e) First and second derivatives of the velocity in ray-centred C coordinate system: V1=VD(2)*H11+VD(3)*H21+VD(4)*H31 CALL SMVPRD(VD(5),H11,H21,H31,AUX1,AUX2,AUX3) V11=AUX1*H11+AUX2*H21+AUX3*H31 V12=AUX1*H42+AUX2*H52+AUX3*H62 CALL SMVPRD(VD(5),H42,H52,H62,AUX1,AUX2,AUX3) V22=AUX1*H42+AUX2*H52+AUX3*H62 C C (f) Correction of the computed quantities: Y(6)=H13/V Y(7)=H23/V Y(8)=H33/V Y(9)=H11 Y(10)=H21 Y(11)=H31 C C (g) Right-hand sides of the differential equations: AUX1=V*VNEXPS AUX2=V*AUX1 AUX3=VNEXPS/V AUX4=V1*VNEXPS D(1)=VNEXPS D(2)=0.5*QL*VNEXPS D(3)=H13*AUX1 D(4)=H23*AUX1 D(5)=H33*AUX1 D(6)=-VD(2)*AUX3 D(7)=-VD(3)*AUX3 D(8)=-VD(4)*AUX3 D(9)=H13*AUX4 D(10)=H23*AUX4 D(11)=H33*AUX4 END IF V11=-V11*AUX3 V12=-V12*AUX3 V22=-V22*AUX3 D(12)=AUX2*Y(14) D(13)=AUX2*Y(15) D(14)=V11*Y(12)+V12*Y(13) D(15)=V12*Y(12)+V22*Y(13) D(16)=AUX2*Y(18) D(17)=AUX2*Y(19) D(18)=V11*Y(16)+V12*Y(17) D(19)=V12*Y(16)+V22*Y(17) D(20)=AUX2*Y(22) D(21)=AUX2*Y(23) D(22)=V11*Y(20)+V12*Y(21) D(23)=V12*Y(20)+V22*Y(21) D(24)=AUX2*Y(26) D(25)=AUX2*Y(27) D(26)=V11*Y(24)+V12*Y(25) D(27)=V12*Y(24)+V22*Y(25) C RETURN END C C======================================================================= C C C SUBROUTINE FCTA(X,Y,D) INTEGER NDIM PARAMETER (NDIM=27) REAL X,Y(NDIM),D(NDIM) C C This subroutine evaluates the right-hand sides of the system of C ordinary differential equations for Y(1) to Y(27). See C C.R.T.5.8.2. C C Input: C X... Value of the independent variable along the ray. C Y... Array containing basic quantities computed along the ray. C None of the input parameters except Y(3)-Y(8) is altered. C C Output: C D... Array containing derivatives of the basic quantities C computed along the ray, with respect to the independent C variable along the ray. C Y(3:8)..Renormalized orthogonal vectors. The modification should C be negligible. C C Common block /DCRT/ (see subroutine file 'ray.for'): INCLUDE 'dcrt.inc' C dcrt.inc C None of the storage locations of the common block are altered. C C Common block /RAYC/ - communication between RAYCB, FCT(A) and OUTP: INCLUDE 'raycb.inc' C raycb.inc C YLRC,IYRC are modified in this subroutine. DELT11,DELT13, and DELT33 C are defined in this subroutine. C C Subroutines and external functions required: EXTERNAL KOOR,METRIC,PARM2,VELOC,SMVPRD INTEGER KOOR C VELOC... File 'model.for' of the package 'MODEL'. C KOOR,METRIC... File 'metric.for' of the package 'MODEL'. C PARM2 and subsequent routines... File 'parm.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C SMVPRD... File 'means.for' of the package 'MODEL'. C C Date: 2014, February 4 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: C INTEGER IERR REAL A(10,21),Q(21) REAL H11,H21,H31,H12,H22,H32 REAL H41,H51,H61,H42,H52,H62,H43,H53,H63,D11,D13,D33 REAL HHP,HHS,HH,HH1,HH2,HH3,HH4,HH5,HH6 REAL HH11,HH12,HH22,HH13,HH23,HH33,HH14,HH24,HH34,HH44 REAL HH15,HH25,HH35,HH45,HH55,HH16,HH26,HH36,HH46,HH56,HH66 REAL V1,V11,V12,V22,V14,V24,V15,V25,V44,V45,V55 REAL VNEXPS,AUX1,AUX2,AUX3 C C A,Q... Elastic moduli and their spatial derivatives. C H11,H21,H31,H12,H22,H32... Components of the contravariant basis C vectors of ray-centred coordinate system. C H41,H51,H61,H42,H52,H62,H43,H53,H63... Components of the covariant C basis vectors of ray-centred coordinate system. C D11,D13,D33... Norms and scalar product of contravariant basis C vector H1 and covariant basis vector H3. C HHP,HHS,HH,HH1,HH2,HH3,HH4,HH5,HH6,HH11,HH12,HH22,HH13,HH23,HH33, C HH14,HH24,HH34,HH44,HH15,HH25,HH35,HH45,HH55,HH16,HH26,HH36,HH46, c HH56,HH66... Hamiltonian and its phase-space derivatives. C V1,V2,V11,V12,V22,V14,V24,V15,V25,V44,V45,V55... Derivatives of C the Hamiltonian in the ray-centred coordinate system. C VNEXPS..Ray velocity**(-NEXPS). C AUX1,AUX2,AUX3... Auxiliary storage locations. C C....................................................................... C C (a) Number of invocations of FCT(A): IYRC(9)=IYRC(9)+1 C C (b) Material parameters: CALL PARM3(IABS(IYRC(5)),Y(3),A,YLRC(3),Q) CALL HDER(IYRC(5),A,Y(6),Y(7),Y(8), * HHP,HHS,HH,HH1,HH2,HH3,HH4,HH5,HH6, * HH11,HH12,HH22,HH13,HH23,HH33,HH14,HH24,HH34,HH44, * HH15,HH25,HH35,HH45,HH55,HH16,HH26,HH36,HH46,HH56,HH66, * E,IERR) IF(IERR.NE.0) THEN IYRC(6)=109 END IF C C (c) Basis of the ray-centred coordinate system: C (c1) Non-unit vectors: H43=Y(6) H53=Y(7) H63=Y(8) H11=Y(9) H21=Y(10) H31=Y(11) C (c3) Norms: D33=SQRT(H43*H43+H53*H53+H63*H63) D11=SQRT(H11*H11+H21*H21+H31*H31) D13=(H11*H43+H21*H53+H31*H63)/(D11*D33) C (c4) Orthonormal vectors: H43=H43/D33 H53=H53/D33 H63=H63/D33 H11=H11/D11-H43*D13 H21=H21/D11-H53*D13 H31=H31/D11-H63*D13 H12=H53*H31-H63*H21 H22=H63*H11-H43*H31 H32=H43*H21-H53*H11 H41=H22*HH6-H32*HH5 H51=H32*HH4-H12*HH6 H61=H12*HH5-H22*HH4 AUX1=H11*H41+H21*H51+H31*H61 H41=H41/AUX1 H51=H51/AUX1 H61=H61/AUX1 H42=HH5*H31-HH6*H21 H52=HH6*H11-HH4*H31 H62=HH4*H21-HH5*H11 AUX1=H12*H42+H22*H52+H32*H62 H42=H42/AUX1 H52=H52/AUX1 H62=H62/AUX1 C C (d) Quantities useful to test the accuracy of computations: DELT11=ABS(D11-1.) DELT33=ABS(SQRT(HH)-1.) DELT13=ABS(D13) C YLRC(1)=SQRT(HHP)/D33 YLRC(2)=SQRT(HHS)/D33 D33=D33/SQRT(HH) YLRC(4)=HH1/D33 YLRC(5)=HH2/D33 YLRC(6)=HH3/D33 C C (e) First and second derivatives of the Hamiltonian in ray-centred C coordinate system: V1=HH1*H11+HH2*H21+HH3*H31 V2=HH1*H12+HH2*H22+HH3*H32 AUX1=HH11*H11+HH12*H21+HH13*H31 AUX2=HH12*H11+HH22*H21+HH23*H31 AUX3=HH13*H11+HH23*H21+HH33*H31 V11=H11*AUX1+H21*AUX2+H31*AUX3-V1*V1 AUX1=HH11*H12+HH12*H22+HH13*H32 AUX2=HH12*H12+HH22*H22+HH23*H32 AUX3=HH13*H12+HH23*H22+HH33*H32 V12=H11*AUX1+H21*AUX2+H31*AUX3-V1*V2 V22=H12*AUX1+H22*AUX2+H32*AUX3-V2*V2 AUX1=HH14*H41+HH15*H51+HH16*H61 AUX2=HH24*H41+HH25*H51+HH26*H61 AUX3=HH34*H41+HH35*H51+HH36*H61 V14=H11*AUX1+H21*AUX2+H31*AUX3 V24=H12*AUX1+H22*AUX2+H32*AUX3 AUX1=HH14*H42+HH15*H52+HH16*H62 AUX2=HH24*H42+HH25*H52+HH26*H62 AUX3=HH34*H42+HH35*H52+HH36*H62 V15=H11*AUX1+H21*AUX2+H31*AUX3 V25=H12*AUX1+H22*AUX2+H32*AUX3 AUX1=HH44*H41+HH45*H51+HH46*H61 AUX2=HH45*H41+HH55*H51+HH56*H61 AUX3=HH46*H41+HH56*H51+HH66*H61 V44=H41*AUX1+H51*AUX2+H61*AUX3 AUX1=HH44*H42+HH45*H52+HH46*H62 AUX2=HH45*H42+HH55*H52+HH56*H62 AUX3=HH46*H42+HH56*H52+HH66*H62 V45=H41*AUX1+H51*AUX2+H61*AUX3 V55=H42*AUX1+H52*AUX2+H62*AUX3 AUX1=(H41*H43+H51*H53+H61*H63)/D33 AUX2=(H42*H43+H52*H53+H62*H63)/D33 V14=V14-V1*AUX1 V24=V24-V2*AUX1 V15=V15-V1*AUX2 V25=V25-V2*AUX2 C C (f) Correction of the computed quantities: Y(6)=Y(6)/SQRT(HH) Y(7)=Y(7)/SQRT(HH) Y(8)=Y(8)/SQRT(HH) Y(9) =H11 Y(10)=H21 Y(11)=H31 C C (g) Right-hand sides of the differential equations: IF(NEXPS.EQ.0) THEN VNEXPS=1. ELSE VNEXPS=SQRT(HH4**2+HH5**2+HH6**2)**(-NEXPS) V11=V11*VNEXPS V12=V12*VNEXPS V22=V22*VNEXPS V14=V14*VNEXPS V24=V24*VNEXPS V15=V15*VNEXPS V25=V25*VNEXPS V44=V44*VNEXPS V45=V45*VNEXPS V55=V55*VNEXPS END IF D(1)=VNEXPS D(2)=0. D(3)= HH4*VNEXPS D(4)= HH5*VNEXPS D(5)= HH6*VNEXPS D(6)=-HH1*VNEXPS D(7)=-HH2*VNEXPS D(8)=-HH3*VNEXPS AUX1=V1/D33*VNEXPS D(9)= H43*AUX1 D(10)=H53*AUX1 D(11)=H63*AUX1 D(12)= V14*Y(12)+V24*Y(13)+V44*Y(14)+V45*Y(15) D(13)= V15*Y(12)+V25*Y(13)+V45*Y(14)+V55*Y(15) D(14)=-V11*Y(12)-V12*Y(13)-V14*Y(14)-V15*Y(15) D(15)=-V12*Y(12)-V22*Y(13)-V24*Y(14)-V25*Y(15) D(16)= V14*Y(16)+V24*Y(17)+V44*Y(18)+V45*Y(19) D(17)= V15*Y(16)+V25*Y(17)+V45*Y(18)+V55*Y(19) D(18)=-V11*Y(16)-V12*Y(17)-V14*Y(18)-V15*Y(19) D(19)=-V12*Y(16)-V22*Y(17)-V24*Y(18)-V25*Y(19) D(20)= V14*Y(20)+V24*Y(21)+V44*Y(22)+V45*Y(23) D(21)= V15*Y(20)+V25*Y(21)+V45*Y(22)+V55*Y(23) D(22)=-V11*Y(20)-V12*Y(21)-V14*Y(22)-V15*Y(23) D(23)=-V12*Y(20)-V22*Y(21)-V24*Y(22)-V25*Y(23) D(24)= V14*Y(24)+V24*Y(25)+V44*Y(26)+V45*Y(27) D(25)= V15*Y(24)+V25*Y(25)+V45*Y(26)+V55*Y(27) D(26)=-V11*Y(24)-V12*Y(25)-V14*Y(26)-V15*Y(27) D(27)=-V12*Y(24)-V22*Y(25)-V24*Y(26)-V25*Y(27) C RETURN END C C======================================================================= C C C SUBROUTINE OUTP(X,Y,D,IHLF,NDIM,PRMT) INTEGER IHLF,NDIM,MY PARAMETER (MY=44) REAL X,Y(MY),D(NDIM),PRMT(5) C C This subroutine includes various tests of the position of the newly C computed point of the ray, tests for possible caustic points, etc. It C also stores the computed quantities into proper files if required. C C Input: C X... Value of the independent variable along the ray. C Y... Array containing basic quantities computed along the ray. C D... Array containing derivatives of the basic quantities C computed along the ray, with respect to the independent C variable along the ray. C IHLF... Number of bisections of the initial increment of C independent variable. C NDIM... Number of ordinary differential equations. C PRMT... Array containing parameters of the integration. C None of the input parameters except X,Y,D,PRMT(5) is altered. C C Output: C X,Y,D...Values corresponding to the point of intersection of the C ray with the boundary of the complex block or the C computational volume if the ray intersects the boundary. C Unaltered if the ray does not intersect the boundary. C PRMT(5)=1.0... The ray has intersected the boundary of the complex C block or the computational volume, or cannot be traced. C -1.0... Number ihlf of bisections of the initial increment C greater than the specified limit NHLF, C otherwise unaltered. C C Common block /DCRT/ (see subroutine file 'ray.for'): INCLUDE 'dcrt.inc' C dcrt.inc C None of the storage locations of the common block are altered. C C Common block /INITC/ (see subroutine file 'init.for'): INCLUDE 'initc.inc' C initc.inc C None of the storage locations of the common block except FSRFCA are C altered. C C Common block /RAYC/ - communication between RAYCB, FCT, and OUTP: INCLUDE 'raycb.inc' C raycb.inc C YLRC,YYRC,IYRC are modified in this subroutine. C C Subroutines and external functions required: EXTERNAL SRFC2,BLOCK,RPAR31,RPAR32,WRIT31,WRIT32,CROSS,HIVD2 EXTERNAL NSRFC,PSHIFT,CDEF INTEGER NSRFC C C NSRFC,BLOCK,VELOC... File 'model.for' of the package 'MODEL'. C SRFC2 and subsequent routines... File 'srfc.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C PARM2 and subsequent routines... File 'parm.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C CDE,CROSS,HIVD2... File 'means.for' of the package 'MODEL'. C RPAR31,RPAR32... File 'rpar.for'. C WRIT31,WRIT32... File 'writ.for'. C CDEF,PSHIFT... This file. C C Date: 2014, February 14 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations for local model parameters: FAUX(10), C G(12),GAMMA(18),GSQRD, UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL: INCLUDE 'auxmod.inc' C auxmod.inc C C....................................................................... C C Other auxiliary storage locations: C REAL TLAST,X0,X1PV,X1,X2,PV1(9),Y1(MY),Y2(MY),D1(MY),D2(MY) SAVE X2,Y2,D2 REAL XB,YB(MY),DB(MY),XAUX,YAUX(MY),DAUX(MY),AUX,ERR INTEGER I,J,K,K1,K2,IE C C TLAST...Travel time at the last but one computed point of the ray. C X0... Beginning of the ray element, for X between PRMT(1) and C X0, is not checked for crossing a structural interface. C X1PV... Beginning of the interval: the last but one computed point C of the ray. C PV1... Reference polarization vectors at point X1PV. C X1,Y1,D1... Independent variable, computed quantities and their C derivatives at the last checked point of the ray. C X2,Y2,D2... Independent variable, computed quantities and their C derivatives at the last computed point of the ray. C XB,YB,DB,XAUX,YAUX,DAUX,AUX,I,J,K,K1,K2... Auxiliary storage C locations. C ERR... Maximum error in independent variable for the C determination of the point of intersection. C C....................................................................... C C Error check during anisotropic-ray-theory S-wave ray tracing IF(IYRC(6).EQ.109) THEN C Error in FCTA during the last step of numerical itegration IF(X.NE.PRMT(1)) THEN C Returning to the previous point of numerical itegration X=X2 DO 10 I=1,NDIM Y(I)=Y2(I) D(I)=D2(I) 10 CONTINUE END IF C Terminating numerical integration PRMT(5)=1. RETURN END IF C C (a) Copying the computed quantities: TLAST=Y2(1) X1PV=X2 X1=X2 X2=X DO 11 I=1,NDIM Y1(I)=Y2(I) Y2(I)=Y(I) D1(I)=D2(I) D2(I)=D(I) 11 CONTINUE C Point X1 is going to be shifted along the ray segment towards C the a priory endpoint X2 of the segment during the processing, C whereas point X serves as the auxiliary point. C When point X1 finally reaches the actual endpoint of the segment, C it is copied to point X. C TLAST will retain the present value of Y1(1). IF(X.EQ.PRMT(1)) THEN C Beginning of the numerical integration CALL PSHIFT(27,Y,YI,I) RETURN END IF C C Reduced amplitude coefficients calculated in OUTP: DO 12 I=NDIM+1,IYRC(1) Y1(I)=Y(I) 12 CONTINUE C C Reference polarization vectors in anisotropic media: IF(IYRC(13).GE.6) THEN C S-wave reference polarization vectors DO 13 I=1,6 PV1(I)=Y(35+I) 13 CONTINUE END IF IF(IYRC(13).GE.3) THEN C P-wave reference polarization vectors for interpolation AUX=Y2(6)*E(7)+Y2(7)*E(8)+Y2(8)*E(9) IF(AUX.LT.0.) THEN E(7)=-E(7) E(8)=-E(8) E(9)=-E(9) END IF DO 14 I=7,9 PV1(I)=Y(35+I) Y2(35+I)=E(I) 14 CONTINUE END IF C C (b) Number of invocations of OUTP: IYRC(10)=IYRC(10)+1 C Initial value of the index of crossed surface bounding the complex C block IYRC(6)=0 C C (c), (d), (e) and (f) are located in the subroutine CDEF. IF(UEBMUL.LE.0.) THEN ERR=AMAX1(ABS(Y(3)),ABS(Y(4)),ABS(Y(5)))/1000000. ERR=ERR/SQRT(D(3)**2+D(4)**2+D(5)**2) ERR=AMAX1(ERR,YYRC(2)/D(1)) ELSE ERR=UEBMUL*YYRC(2)/D(1) END IF X0=PRMT(1)+ERR C C (h) Storage of the computed quantities along the ray: C XAUX... Points along ray, and the endpoint of the segment. IF(STORE.GT.0.) THEN K1=INT(X1/STORE+0.000001)+1 K2=INT(X/STORE+0.000001) ELSE K1=1 K2=0 END IF DO 84 K=K1,K2+1 IF(K.LE.K2) THEN XAUX=FLOAT(K)*STORE CALL HIVD2(NDIM,X1,Y1,D1,X2,Y2,D2,XAUX,YAUX,DAUX) ELSE XAUX=X2 DO 70 I=1,NDIM YAUX(I)=Y2(I) DAUX(I)=D2(I) 70 CONTINUE END IF C C (g) Storage of the computed quantities at given surfaces: C X... Intersections with the surfaces between X1 and XAUX. 71 CONTINUE X=XAUX DO 72 I=1,NDIM Y(I)=YAUX(I) D(I)=DAUX(I) 72 CONTINUE J=0 DO 75 I=1,NSTOR IF(KSTOR(I).GT.NSRFC().AND.KSTOR(I).LE.100) THEN DO 73 IE=1,NEND IF(KSTOR(I).EQ.KEND(IE)) THEN GO TO 74 END IF 73 CONTINUE CALL SRFC2(KSTOR(I),Y(3),FAUX) IF(FAUX(1)*FSRFCA(KSTOR(I)-NSRFC()).LT.0.) THEN J=I AUX=FAUX(1) CALL CROSS(SRFC2,KSTOR(I),3,5,NDIM,ERR,X1,Y1,D1,X2,Y2,D2, * X,Y,D,XB,YB,DB,FAUX) END IF END IF 74 CONTINUE 75 CONTINUE IF(J.NE.0) THEN CALL CDEF(NDIM,ERR,X0,X1,Y1,D1,X2,Y2,D2,X,Y,D,XB,YB,DB) IF(IYRC(6).NE.0) THEN C Boundary of the complex block is crossed C X1=X... Point of intersection with the boundary. GO TO 90 END IF C X1=X... Point of intersection with the last storing surface. CALL PVINT(IYRC(13),X1PV,PV1,X2,Y2(36),X1,Y1(36)) CALL RPAR32(J,YLRC,Y1,YYRC,IYRC) CALL WRIT32(J,YLRC,Y1,YYRC,IYRC,IYRC(1)-27,Y1) FSRFCA(KSTOR(J)-NSRFC())=AUX GO TO 71 END IF C (g-end) C CALL CDEF(NDIM,ERR,X0,X1,Y1,D1,X2,Y2,D2,XAUX,YAUX,DAUX,XB,YB,DB) IF(IYRC(6).NE.0) THEN C Boundary of the complex block is crossed C X1=XAUX... Point of intersection with the boundary. GO TO 90 END IF C X1=XAUX... The last point along the ray. IF(K.LE.K2) THEN CALL PVINT(IYRC(13),X1PV,PV1,X2,Y2(36),X1,Y1(36)) CALL RPAR31(YLRC,Y1,YYRC,IYRC) CALL WRIT31(YLRC,Y1,YYRC,IYRC) END IF 84 CONTINUE C (h-end) C 90 CONTINUE C X1... Actual endpoint of the segment. C C Copying point Y1 to array Y: X=X1 DO 92 I=1,NDIM Y(I)=Y1(I) D(I)=D1(I) 92 CONTINUE DO 93 I=NDIM+1,IYRC(1) Y(I)=Y1(I) 93 CONTINUE CALL PVINT(IYRC(13),X1PV,PV1,X2,Y2(36),X,Y(36)) C X=X1... Actual endpoint of the segment. C C (i) Accumulation of the renormalization errors for a test of C accuracy: AUX=0.5*(Y(1)-TLAST) YYRC(3)=YYRC(3)+DELT33*AUX YYRC(4)=YYRC(4)+DELT13*AUX YYRC(5)=YYRC(5)+DELT11*AUX C C (j) Great number of bisections of the initial increment: IF(IHLF.GT.NHLF) THEN PRMT(5)=-1. END IF C C Check for very small velocity: IF(IYRC(6).EQ.0) THEN IF(1.0E-12*(Y (6)**2+Y (7)**2+Y (8)**2).GT. * (YI(6)**2+YI(7)**2+YI(8)**2)) THEN IYRC(6)=108 END IF END IF C C (k) Final operations: C YYRC(1)=X is set by the subroutine CDEF IF(IYRC(6).NE.0) THEN PRMT(5)=1. END IF RETURN END C C======================================================================= C C C SUBROUTINE CDEF(NDIM,ERR,X0,X1,Y1,D1,X2,Y2,D2,X,Y,D,XB,YB,DB) INTEGER NDIM,MY PARAMETER (MY=44) REAL ERR,X0,X1,Y1(MY),D1(NDIM),X2,Y2(MY),D2(NDIM) REAL X,Y(MY),D(NDIM),XB,YB(MY),DB(NDIM) C C Auxiliary subroutine to OUTP, performing steps 5.8.3(c),(d),(e) and C (f) of the algorithm. C C Common block /DCRT/ (see subroutine file 'ray.for'): INCLUDE 'dcrt.inc' C dcrt.inc C None of the storage locations of the common block are altered. C C Common block /INITC/ (see subroutine file 'init.for'): INCLUDE 'initc.inc' C initc.inc C None of the storage locations of the common block are altered. C C Common block /RAYC/ - communication between RAYCB, FCT, and OUTP: INCLUDE 'raycb.inc' C raycb.inc C YLRC,YYRC,IYRC are modified in this subroutine. C C Subroutines and external functions required: EXTERNAL RSEP3R,CDE,PSHIFT,PARM2,VELOC C RSEP3R..File 'sep.for' of the package 'FORMS'. C NSRFC,BLOCK,VELOC... File 'model.for' of the package 'MODEL'. C SRFC2 and subsequent routines... File 'srfc.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C PARM2 and subsequent routines... File 'parm.for' and subsequent C files (like 'val.for' and 'fit.for') of the package C 'MODEL'. C CDE,CROSS,HIVD2... File 'means.for' of the package 'MODEL'. C C Date: 2014, February 16 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations for local model parameters: FAUX(10), C G(12),GAMMA(18),GSQRD, UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL: INCLUDE 'auxmod.inc' C auxmod.inc C C Auxiliary storage locations in an anisotropic model: INTEGER IERR REAL A(10,21),Q(21) REAL HHP,HHS,HH,HH1,HH2,HH3,HH4,HH5,HH6 REAL HH11,HH12,HH22,HH13,HH23,HH33,HH14,HH24,HH34,HH44 REAL HH15,HH25,HH35,HH45,HH55,HH16,HH26,HH36,HH46,HH56,HH66 C A,Q... Elastic moduli and their spatial derivatives. C HHP,HHS,HH,HH1,HH2,HH3,HH4,HH5,HH6,HH11,HH12,HH22,HH13,HH23,HH33, C HH14,HH24,HH34,HH44,HH15,HH25,HH35,HH45,HH55,HH16,HH26,HH36,HH46, C HH56,HH66... Hamiltonian and its phase-space derivatives. C C....................................................................... C C Anisotropic ray tracing: REAL CRTANI SAVE CRTANI C C Other auxiliary storage locations: REAL D33 INTEGER IBOUND(7),I DATA IBOUND/3,-3,4,-4,5,-5,-1/ C C....................................................................... C C Anisotropic ray tracing: DATA CRTANI/-999./ IF(CRTANI.EQ.-999.) THEN CALL RSEP3R('CRTANI',CRTANI,0.) END IF C C (c) Check for crossing the coordinate boundaries of the C computational volume, C (d) Check for crossing the boundary of the complex block, C (e) Check for crossing the end surfaces: CALL CDE(0,NEND,KEND,7,IBOUND,BOUNDR, * 3,5,NDIM,IYRC,ERR,X0,X1,Y1,D1,X2,Y2,D2,X,Y,D,XB,YB,DB) C C Moving point Y1 to point Y: X1=X DO 61 I=1,NDIM Y1(I)=Y(I) D1(I)=D(I) 61 CONTINUE C C (f) Phase shift due to caustics: CALL PSHIFT(IYRC(1),Y1,YI,I) IYRC(12)=IYRC(12)+I C C Quantities describing local properties of the model at point X,Y: IF(IYRC(6).NE.0) THEN IF(CRTANI.EQ.0.) THEN CALL PARM2(IABS(IYRC(5)),Y(3),UP,US,YLRC(3),QP,QS) CALL VELOC(IYRC(5),UP,US,QP,QS,YLRC(1),YLRC(2),VD,QL) YLRC(4)=VD(2) YLRC(5)=VD(3) YLRC(6)=VD(4) ELSE CALL PARM3(IABS(IYRC(5)),Y(3),A,YLRC(3),Q) CALL HDER(IYRC(5),A,Y(6),Y(7),Y(8), * HHP,HHS,HH,HH1,HH2,HH3,HH4,HH5,HH6, * HH11,HH12,HH22,HH13,HH23,HH33,HH14,HH24,HH34,HH44, * HH15,HH25,HH35,HH45,HH55,HH16,HH26,HH36,HH46,HH56, * HH66,E,IERR) IF(IERR.NE.0) THEN IYRC(6)=109 END IF D33=SQRT(Y(6)**2+Y(7)**2+Y(8)**2) YLRC(1)=SQRT(HHP)/D33 YLRC(2)=SQRT(HHS)/D33 D33=D33/SQRT(HH) YLRC(4)=HH1/D33 YLRC(5)=HH2/D33 YLRC(6)=HH3/D33 END IF END IF YYRC(1)=X C C X1=X,Y1=Y,D1=D,YLRC,YYRC,IYRC contain the quantities either at the C given point X,Y,D (for IYRC(6).EQ.0) or at the endpoint of the C computed ray element (for IYRC(6).NE.0). RETURN END C C======================================================================= C C C SUBROUTINE PSHIFT(NY,Y,YI,ISHIFT) INTEGER NY,ISHIFT REAL Y(NY),YI(21) C C This subroutine is an auxiliary routine to OUTP. It corrects reduced C amplitudes with regard to phase shift due to caustics, see C C.R.T.5.8.3f. C Hereinafter, the point 1 of the ray denotes the point corresponding C to the previous invocation of the subroutine PSHIFT, the point 2 C denotes the point corresponding to this invocation of the subroutine C PSHIFT. C C Input: C NY... Number of basic quantities describing the point of a ray. C NY=27 at the initial point of a ray element, where no C phase shift is applied to the amplitudes. C Y(1:11)... Anything. C Y(12:27)... Ray propagator matrix at the point 2 of the ray. C Y(28:NY)... Reduced amplitudes at the point 1 of the ray. C YI... Array containing quantities at the initial point of the C ray, see C C.R.T.6.1. C None of the input parameters except Y(28:NY) are altered. C C Output: C Y(28:NY)... Reduced amplitudes at the point 2 of the ray. C ISHIFT..Order of a caustic between points 1 and 2 (increment of C the KMAH index). C C Subroutines and external functions called: EXTERNAL RPAR30 C C Date: 2014, February 13 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I REAL P211,P221,P212,P222,AUX,SMALL PARAMETER (SMALL=0.000000001) REAL Q111,Q121,Q112,Q122,Q1DET REAL Q211,Q221,Q212,Q222,Q2DET SAVE Q211,Q221,Q212,Q222,Q2DET C C I,P211,P221,P212,P222,AUX...Auxiliary storage locations. C Q111,Q121,Q112,Q122,Q1DET... Matrix of ray geometrical spreading C and its determinant corresponding to the previous C invocation of PSHIFT at the point 1. C Q211,Q221,Q212,Q222,Q2DET... Matrix of ray geometrical spreading C and its determinant corresponding to the last invocation C of PSHIFT. C IF(NY.GE.28) THEN Q111=Q211 Q121=Q221 Q112=Q212 Q122=Q222 Q1DET=Q2DET END IF Q211=Y(12)*YI(12)+Y(16)*YI(13)+Y(20)*YI(14)+Y(24)*YI(15) Q221=Y(13)*YI(12)+Y(17)*YI(13)+Y(21)*YI(14)+Y(25)*YI(15) Q212=Y(12)*YI(16)+Y(16)*YI(17)+Y(20)*YI(18)+Y(24)*YI(19) Q222=Y(13)*YI(16)+Y(17)*YI(17)+Y(21)*YI(18)+Y(25)*YI(19) CALL RPAR30(YI(1),Y(1),Q211,Q221,Q212,Q222) C Regularization in order to avoid floating-point overflow: AUX=SQRT(Q211*Q211+Q221*Q221+Q212*Q212+Q222*Q222+1.) Q211=Q211/AUX Q221=Q221/AUX Q212=Q212/AUX Q222=Q222/AUX Q2DET=Q211*Q222-Q212*Q221 IF(Y(1).EQ.YI(1)) THEN IF(Q2DET.EQ.0.) THEN P211=Y(14)*YI(12)+Y(18)*YI(13)+Y(22)*YI(14)+Y(26)*YI(15) P221=Y(15)*YI(12)+Y(19)*YI(13)+Y(23)*YI(14)+Y(27)*YI(15) P212=Y(14)*YI(16)+Y(18)*YI(17)+Y(22)*YI(18)+Y(26)*YI(19) P222=Y(15)*YI(16)+Y(19)*YI(17)+Y(23)*YI(18)+Y(27)*YI(19) Q211=Q211+P211*SMALL Q221=Q221+P221*SMALL Q212=Q212+P212*SMALL Q222=Q222+P222*SMALL Q2DET=Q211*Q222-Q212*Q221 END IF END IF IF(NY.GE.28) THEN IF(Q1DET*Q2DET.LT.0.) THEN C Caustic of the first order (line caustic) between points ISHIFT=1 DO 10 I=28,NY,2 AUX=Y(I+1) Y(I+1)=-Y(I) Y(I)=AUX 10 CONTINUE ELSE AUX=(Q111*Q222-Q112*Q221+Q122*Q211-Q121*Q212)*Q2DET IF(AUX.LT.0.) THEN C Caustic of the second order (point caustic) between points ISHIFT=2 DO 20 I=28,NY Y(I)=-Y(I) 20 CONTINUE ELSE IF(AUX.EQ.0.) THEN C Stepping on a caustic ISHIFT=1 DO 30 I=28,NY,2 AUX=Y(I+1) Y(I+1)=-Y(I) Y(I)=AUX 30 CONTINUE ELSE C No caustic ISHIFT=0 END IF END IF ELSE C Initial point of the ray ISHIFT=0 END IF RETURN END C C======================================================================= C C C SUBROUTINE PVINT(NPV,X1,PV1,X2,PV2,X,PV) INTEGER NPV REAL X1,PV1(9),X2,PV2(9),X,PV(9) C C This subroutine is an auxiliary routine to OUTP calculates the C reference polarization vectors in anisotropic media. C It interpolates the P-wave reference polarization vector from points C X1 and X2 to point X, and calculates the S-wave reference polarization C vectors at point X from point X1. C C Input: C NPV... Total number of the components of the polarization vectors C to be written to the output files. C NPV=0: Isotropic medium. C NPV=3: P wave in an anisotropic medium. C NPV=6: S wave in an anisotropic medium. C X1,X2,X... Values of the independent variable at the endpoints X1 C and X2 of the given interval and at given point X. C PV1... Reference polarization vector at point X1. C PV(1:3),PV(4:6): S-wave reference polarization vectors. C Used only if NPV.GE.6. C PV1(7:9): P-wave reference polarization vector. C Used only if NPV.GE.3. C PV2... Reference polarization vector at point X2. C PV2(1:3),PV2(4:6): Not used, need not be defined. C PV2(7:9): P-wave reference polarization vector. C Used only if NPV.GE.3. C C Output: C PV... Reference polarization vector at point X. C PV(1:3),PV(4:6): S-wave reference polarization vectors. C Defined only if NPV.GE.6. C PV(7:9): P-wave reference polarization vector. C Defined only if NPV.GE.3. C C Date: 2014, February 4 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: REAL A1,A2,A3 C C....................................................................... C C Interpolation of the P-wave reference polarization vector: IF(NPV.GE.3) THEN A3=(PV2(7)-PV1(7))**2+(PV2(8)-PV1(8))**2+(PV2(9)-PV1(9))**2 IF(X.EQ.X1.OR.A3.EQ.0.) THEN PV(7)=PV1(7) PV(8)=PV1(8) PV(9)=PV1(9) ELSE IF(X.EQ.X2) THEN PV(7)=PV2(7) PV(8)=PV2(8) PV(9)=PV2(9) ELSE IF(X2-X1.LE.0.) THEN C 5A2 CALL ERROR('5A2 in PVINT: Vanishing interpolation interval') C Endpoints X1 and X2 of the given interval do not satisfy C inequality X1.LT.X2. C This error should not appear. Contact the authors. ELSE A1=0.000001*(ABS(X1)+ABS(X2)) IF(X.LT.X1-A1.OR.X2+A1.LT.X) THEN C 5A3 CALL ERROR('5A3 in PVINT: Interpolation outside interval') C Point X for interpolation of the P-wave polarization vector C is not situated between endpoints X1 and X2 of the given C interval. C This error should not appear. Contact the authors. END IF A3=2.*ASIN(0.5*SQRT(A3)) A2=A3*(X-X1)/(X2-X1) A1=A3-A2 A3=SIN(A3) A1=SIN(A1)/A3 A2=SIN(A2)/A3 PV(7)=PV1(7)*A1+PV2(7)*A2 PV(8)=PV1(8)*A1+PV2(8)*A2 PV(9)=PV1(9)*A1+PV2(9)*A2 END IF END IF C C Calculation of the S wave reference polarization vectors: IF(NPV.GE.6) THEN A1=PV1(1)*PV(7)+PV1(2)*PV(8)+PV1(3)*PV(9) A2=PV1(4)*PV(7)+PV1(5)*PV(8)+PV1(6)*PV(9) A3=PV1(7)*PV(7)+PV1(8)*PV(8)+PV1(9)*PV(9) A3=A3+1. A1=A1/A3 A2=A2/A3 PV(1)=PV1(1)-(PV1(7)+PV(7))*A1 PV(2)=PV1(2)-(PV1(8)+PV(8))*A1 PV(3)=PV1(3)-(PV1(9)+PV(9))*A1 PV(4)=PV1(4)-(PV1(7)+PV(7))*A2 PV(5)=PV1(5)-(PV1(8)+PV(8))*A2 PV(6)=PV1(6)-(PV1(9)+PV(9))*A2 A1=SQRT(PV(1)**2+PV(2)**2+PV(3)**2) A2=SQRT(PV(4)**2+PV(5)**2+PV(6)**2) PV(1)=PV(1)/A1 PV(2)=PV(2)/A1 PV(3)=PV(3)/A1 PV(4)=PV(4)/A2 PV(5)=PV(5)/A2 PV(6)=PV(6)/A2 END IF C RETURN END C C======================================================================= C