C
C Subroutine file 'raycb.for' for complete ray tracing within one C complex block C C Date: 1999, May 25 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.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 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.R.T.5.8.3). C OUTP C CDEF... Auxiliary subroutine to OUTP, performing steps 5.8.3(c), C (d), (e) and (f) 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 C.R.T.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.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 C.R.T.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: 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 External procedure names: EXTERNAL FCT,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 (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: CALL HPCG(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX) 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.R.T., Section C 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: 1993, January 23 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.) DELT11=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.) DELT11=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 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: 1997, November 13 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. ERR=YYRC(2)*D(1) 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 (k) Final operations: C YYRC(1)=X is set by the subroutine CDEF IF(IYRC(6).NE.0) PRMT(5)=1. 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. C (see C.R.T.5.8.3f). Hereinafter, the point 1 of the ray denotes the C point corresponding to the previous invocation of the subroutine C PSHIFT, the point 2 denotes the point corresponding to this invocation C of the subroutine 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.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