C
C Subroutine file 'raycb.for' for complete ray tracing within one C complex block C C Date: 2004, May 19 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 C======================================================================= C C C SUBROUTINE RAYCB(YL,Y,YY,IY) REAL YL(6),Y(35),YY(5) INTEGER IY(12) 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 Y... Array containing basic quantities computed along the ray C at the endpoint of the element of the ray. C YY... Array containing real auxiliary quantities computed along C the ray at the endpoint of the element of the ray. C IY... Array containing integer auxiliary quantities computed C along the ray at the endpoint of the element of the ray. 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 METRIC,HPCG C NSRFC,BLOCK,VELOC... File 'model.for' of the package 'MODEL'. C KOOR,METRIC... File 'metric.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: 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 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 CALL HPCG(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) ELSE 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,12 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: 2004, February 25 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: C 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,V11,V12,V22,V14,V24,V15,V25,V44,V45,V55... Derivatives of the C Hamiltonian in 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) 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=35) 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, 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: 2004, 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 Other auxiliary storage locations: C REAL TLAST,X0,X1,X2,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 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 (a) Copying the computed quantities: TLAST=Y2(1) 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 DO 12 I=NDIM+1,IYRC(1) Y1(I)=Y(I) 12 CONTINUE C X1 is going to be shifted along the ray segment towards to the C endpoint of the segment during the processing, whereas TLAST will C 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 (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 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 RPAR31(YLRC,Y1,YYRC,IYRC) CALL WRIT31(YLRC,Y1,YYRC,IYRC) END IF 84 CONTINUE C (h-end) C 90 CONTINUE C X1... Endpoint of the segment. C> DO 91 K=1,NSTOR C> IF(KSTOR(K).GT.NSRFC().AND.KSTOR(K).LE.100) THEN C> IF(KSTOR(K).EQ.IYRC(6)) THEN C> CALL RPAR32(K,YLRC,Y1,YYRC,IYRC) C> CALL WRIT32(K,YLRC,Y1,YYRC,IYRC,IYRC(1)-27,Y1) C> END IF C> END IF C> 91 CONTINUE 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 C X=X1... 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) PRMT(5)=-1. 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=35) 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 CDE,PSHIFT,PARM2,VELOC 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: 1993, January 15 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: INTEGER IBOUND(7),I DATA IBOUND/3,-3,4,-4,5,-5,-1/ C C....................................................................... C C (c) Check for crossing the coordinate boundaries of the 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 (f) Phase shift due to caustics: X1=X DO 61 I=1,NDIM Y1(I)=Y(I) D1(I)=D(I) 61 CONTINUE 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 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) 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: 1999, May 25 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 END IF RETURN END C C======================================================================= C