C FILE 'RAY.FOR' FOR THE COMPLETE TRACING OF A RAY FROM A GIVEN POINT C C BY VLASTISLAV CERVENY, LUDEK KLIMES, IVAN PSENCIK C C THIS FILE CONSISTS OF: C RAYB... BLOCK DATA SUBROUTINE DEFINING COMMON BLOCKS /RAYT/ AND C /DCRT/ TO STORE THE INPUT DATA FOR COMPLETE RAY TRACING. C DESCRIPTION OF THE QUANTITIES COMPUTED ALONG A RAY (SEE C C.R.T.5.2). C RAY1... SUBROUTINE DESIGNED TO READ THE INPUT DATA FOR COMPLETE C RAY TRACING AND TO STORE THEM IN THE COMMON BLOCKS /RAYT/ C AND /DCRT/ (SEE C.R.T.5.6). C RAY2... SUBROUTINE DESIGNED TO CONTINUE IN THE COMPLETE RAY C TRACING OF A RAY FROM THE GIVEN POINT (SEE C.R.T.5.7). C C INPUT DATA FOR COMPLETE RAY TRACING: C THE DATA ARE READ IN BY THE LIST DIRECTED INPUT (FREE FORMAT). IN C THE LIST OF INPUT DATA BELOW, EACH NUMBERED PARAGRAPH INDICATES C THE BEGINNING OF A NEW INPUT OPERATION (NEW READ STATEMENT). IF C THE FIRST LETTER OF THE SYMBOLIC NAME OF THE INPUT VARIABLE IS C I-N, THE CORRESPONDING VALUE IN INPUT DATA MUST BE OF THE TYPE C INTEGER. OTHERWISE (EXCEPT TEXTR), THE INPUT PARAMETER IS OF THE C TYPE REAL. C (1) TEXTR C STRING DESCRIBING THE DATA. ONLY THE FIRST 80 CHARACTERS OF THE C STRING ARE SIGNIFICANT. C (2) KSTORE,NEXPS,NHLF,MODCRT C QUANTITIES CONTROLING NUMERICAL INTEGRATION. C KSTORE..SPECIFIES WHETHER THE CONVERSION COEFFICIENTS ARE TO BE C CONSIDERED (SEE C.R.T.5.5.4 AND 5.6E): C KSTORE.LE.1... NO AMPLITUDE CONVERSION COEFFICIENTS ARE C APPLIED AT THE POINT OF INTERSECTION WITH AN INTERFACE. C KSTORE.GE.2... AMPLITUDE CONVERSION COEFFICIENTS ARE C APPLIED AT THE POINT OF INTERSECTION WITH AN INTERFACE. C WHENEVER KSTORE IS CHANGED, DATA (8) AND (9) SHOULD BE C ADJUSTED PROPERLY. C NEXP... SPECIFIES INDEPENDENT VARIABLE ALONG RAYS, SEE C C.R.T.(5.1). C NHLF... MAXIMUM ALLOWED NUMBER OF HALVINGS (BISECTIONS) OF INITIAL C INCREMENT OF INDEPENDENT VARIABLE DURING NUMERICAL C INTEGRATION. C MODCRT..SPECIFIES THE MODE OF COMPLETE RAY TRACING, AND THE RAY C ELEMENTS ALONG WHICH THE COMPUTED QUANTITIES ARE STORED. C 0...INITIAL MODE, C STORING THE COMPUTED QUANTITIES ALONG THE WHOLE RAYS, C NO STORING AT THE ENDPOINTS OF THE RAY ELEMENTS. C 1...INITIAL MODE, C STORING THE QUANTITIES JUST ALONG NEW RAY ELEMENTS, C NO STORING AT THE ENDPOINTS OF THE RAY ELEMENTS. C 2...INITIAL MODE, C STORING THE QUANTITIES JUST ALONG NEW RAY ELEMENTS, C STORING AT THE ENDPOINTS OF THE NEW RAY ELEMENTS. C 3...INTERPOLATION MODE, C STORING THE QUANTITIES JUST ALONG NEW RAY ELEMENTS, C STORING AT THE ENDPOINTS OF THE NEW RAY ELEMENTS. C *** MODCRT=3 IS NOT ENABLED BY THE CURRENT VERSION OF C 'INIT.FOR' ***. C (3) STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT C QUANTITIES CONTROLING NUMERICAL INTEGRATION. C STORE...STEP OF INDEPENDENT VARIABLE FOR STORING THE COMPUTED C QUANTITIES ALONG A RAY (C.R.T.5.5.1). FOR STORE=0 THE C QUANTITIES ARE NOT STORED ALONG RAYS. C STEP... INITIAL INCREMENT OF INDEPENDENT VARIABLE FOR NUMERICAL C INTEGRATION. C UEB... UPPER ERROR BOUND OF TRAVEL TIME PER ONE STEP OF NUMERICAL C INTEGRATION. ERRORS IN THE COORDINATES OF POINTS ALONG C THE RAY ARE APPROXIMATELY TRANSFORMED TO UNITS OF TRAVEL C TIME AND ARE ALSO BOUNDED BY UEB. THE ERROR PER STEP OF C NUMERICAL INTEGRATION IS AUTOMATICALLY KEPT WITHIN THE C LIMIT UEB IF THIS DOES NOT REQUIRE MORE THAN NHLF C BISECTIONS OF THE INITIAL INCREMENT STEP. IN THE OPPOSITE C CASE THE UPPER ERROR BOUND IS 2,4,8... TIMES GREATER FOR C THE COMPUTATION OF THE REST OF THE RAY. THUS THE C COMPUTATION OF EACH RAY IS COMPLETED. C UEBPP,UEBPH,UEBHH... MAXIMUM ALLOWED ACCUMULATED DEVIATIONS OF THE C TWO POLARIZATION VECTORS FROM THE CONDITIONS OF C ORTHONORMALITY (C.R.T.5.8.2D AND 5.8.3I). THE ACCUMULATED C DEVIATIONS ARE EXPRESSED IN TIME UNITS. IF ANY OF THE C ACCUMULATED DEVIATIONS OVERRFLOWS ITS UPPER BOUND, A C WARNING IS GENERATED. C UEBDRT..MAXIMUM ALLOWED DEVIATION OF L.H.S. AND R.H.S. IN C C.R.T.,EQ.(5.11). IF THE DEVIATION OF ANY COMPONENT OF C THE 4*4 MATRIX EXCEEDS UEBDRT, A WARNING IS GENERATED. C (4) X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,TMAX C THE BOUNDARIES OF THE COMPUTATIONAL VOLUME. C X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX... COORDINATES SPECIFYING THE C COORDINATE PLANES BOUNDING THE COMPUTATIONAL VOLUME. THE C COORDINATE PLANES ARE INDEXED 101, 102, 103, 104, 105 AND C 106. DEFAULT VALUES ARE THE BOUNDARIES OF THE MODEL, SEE C SUBROUTINE MODEL1 AND THE INPUT DATA FOR THE MODEL. C TMAX... MAXIMUM TRAVEL TIME FOR THE COMPUTATION OF RAYS. THE C CORRESPONDING ISOCHRONE IS INDEXED BY 107. DEFAULT VALUE C IS 999999. C (5) NSRFCA C NUMBER OF AUXILIARY SURFACES, SEE C.R.T.5.3. THE SURFACES ARE C INDEXED SEQUENTIALLY BY POSITIVE INTEGERS FROM NSRFC+1 TO C NSRFC+NSRFCA. DEFAULT VALUE IS 0. C (6) THE INDICES OF END SURFACES. C EACH INDEX MUST HAVE THE SIGN CORRESPONDING TO THE SIDE OF THE C SURFACE AT WHICH THE COMPUTATIONAL VOLUME IS SITUATED. C THE INDICES MAY BE SPECIFIED IN AN ARBITRARY ORDER AND MUST BE C TERMINATED BY A SLASH. SEE C.R.T.5.4F. C (7) THE INDICES OF SURFACES FOR STORING COMPUTED QUANTITIES. C THE INDICES OF SURFACES 1 TO NSRFC COVERING STRUCTURAL INTERFACES C MUST INCLUDE THE PROPER SIGN. THE INDICES MAY BE SPECIFIED IN AN C ARBITRARY ORDER AND MUST BE TERMINATED BY A SLASH. NOTE THAT THE C SAME SURFACE MAY BE SPECIFIED SEVERAL TIMES, ESPECIALLY IF THE C QUANTITIES AT DIFFERENT POINTS OF INTERSECTIONS WITH THE SURFACE C ARE STORED IN DIFFERENT FILES, SEE (8) AND (9) BELOW. C FOR DETAILS SEE C.R.T.5.5.2. C (8) FOR EACH ABOVE SURFACE, THE SEQUENTIAL NUMBER OF THE FIRST POINT C OF INTERSECTION OF A RAY WITH THE SURFACE IN WHICH THE COMPUTED C QUANTITIES ARE STORED INTO THE CORRESPONDING FILE. A RAY MAY HAVE C SEVERAL POINTS OF INTERSECTION WITH A SPECIFIED SURFACE. THE C QUANTITIES NEED NOT BE STORED AT ALL POINTS OF INTERSECTION. THEY C MAY BE STORED JUST FROM THE POINT OF INTERSECTION SPECIFIED IN C THIS INPUT (8) TO THE ONE SPECIFIED IN THE INPUT (9) BELOW. C USUALLY, USER WILL WISH TO STORE THE QUANTITIES AT ALL POINTS OF C INTERSECTION, THAT CORRESPONDS TO THE VALUES OF 1 IN THIS INPUT C AND IS THE DEFAULT (INDICATED BY A SLASH IN THE INPUT DATA). C FOR A STRUCTURAL INTERFACE (NOT AUXILIARY SURFACE), THE POSITIVE C OR NEGATIVE SIDE OF THE SURFACE ISRFR IS SPECIFIED, THEN: C NO POINT OF INTERSECTION IS ACCOUNTED IF THE REFLECTED RAY IS C SITUATED AT THE OTHER SIDE OF THE SURFACE ISRFR THAN THE C SPECIFIED ONE, AND C ONE POINT OF INTERSECTION IS ACCOUNTED IF THE REFLECTED RAY IS C SITUATED AT THE SPECIFIED SIDE OF THE SURFACE ISRFR AND C KSTORE.GE.2 IN THE DATA (2). C TWO POINTS OF INTERSECTION ARE ACCOUNTED IF THE REFLECTED RAY IS C SITUATED AT THE SPECIFIED SIDE OF THE SURFACE ISRFR AND C KSTORE.LE.1 IN THE DATA (2). C (9) FOR EACH ABOVE SURFACE, THE SEQUENTIAL NUMBER OF THE LAST POINT C OF INTERSECTION OF A RAY WITH THE SURFACE IN WHICH THE COMPUTED C QUANTITIES ARE STORED INTO THE CORRESPONDING FILE. DEFAULT VALUES C OF THIS INPUT ARE 999999. C (10) THE DATA SPECIFYING NSRFCA FUNCTIONS DESCRIBING THE AUXIIARY C SURFACES. THE DATA ARE READ BY SUBROUTINE SRFC1. FOR THEIR C DESCRIPTION REFER TO SUBROUTINE SRFC1. C THIS INPUT IS PERFORMED ONLY IF NSRFCA.GT.0, SEE (5). C C STORAGE IN THE MEMORY: C THE INPUT DATA (1) TO (9) ARE STORED IN COMMON BLOCKS /RAYT/ AND C /DCRT/ DEFINED IN THE FOLLOWING SUBROUTINE: C ------------------------------------------------------------------ BLOCK DATA RAYB CHARACTER*80 TEXTR COMMON/RAYT/TEXTR SAVE /RAYT/ INTEGER MEND,MSTOR PARAMETER (MEND=128) PARAMETER (MSTOR=128) INTEGER KSTORE,NEXPS,NHLF,MODCRT REAL STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT,BOUNDR(7) INTEGER NSRFCA,NEND,KEND(MEND),NSTOR,KSTOR(MSTOR) COMMON/DCRT/ KSTORE,NEXPS,NHLF,MODCRT,STORE,STEP,UEB,UEBPP,UEBPH, * UEBHH,UEBDRT,BOUNDR,NSRFCA,NEND,KEND,NSTOR,KSTOR SAVE /DCRT/ END C ------------------------------------------------------------------ C TEXTR...THE NAME OF THE DATA. STRING OF 80 CHARACTERS. C KSTORE,NEXPS,NHLF,MODCRT... INPUT DATA (2). C STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT... INPUT DATA (3). C BOUNDR..BOUNDARIES X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,TMAX OF THE C COMPUTATIONAL VOLUME, SEE INPUT DATA (4). C NSRFCA..NUMBER OF AUXILIARY SURFACES, SEE INPUT DATA (5). C NEND... NUMBER OF END SURFACES, SEE INPUT DATA (6). C KEND... CONTAINS THE INDICES OF END SURFACES, SEE INPUT DATA (6). C NSTOR...NUMBER OF SURFACES FOR STORING COMPUTED QUANTITIES, SEE C INPUT DATA (7). C KSTOR...CONTAINS THE INDICES OF SURFACES FOR STORING COMPUTED C QUANTITIES, SEE INPUT DATA (7), THE INPUT DATA (8) AND C FINALY, THE INPUT DATA (9). THE TOTAL OF 3*NSTOR C INTEGERS. C COMMON BLOCK /DCRT/ IS ALSO INCLUDED IN SUBROUTINE FILES C 'RAYCB.FOR', 'INIT.FOR', 'WRIT.FOR' AND 'SCROPC.FOR'. C ALL THE INPUT DATA ARE STORED SEQUENTIALLY IN THE SAME ORDER AS C THEY WERE READ. THE ONLY EXCEPTION ARE LOCATIONS NEND AND NSTOR C WHICH ARE INSERTED WHEN READING THE INPUT DATA. THE INDEX OF THE C LAST ALLOCATED NUMERIC UNIT OF ARRAY KEND IS NAMED MEND. THE C INDEX OF THE LAST ALLOCATED NUMERIC UNIT OF ARRAY KSTOR IS NAMED C MSTOR. IF MEND OR MSTOR IS CHANGED, IT MUST BE ADJUSTED IN ALL C SUBROUTINES WHICH INCLUDE COMMON BLOCK /DCRT/. C C DATE: 1994, JANUARY 15 C CODED BY LUDEK KLIMES C C======================================================================= C C DESCRIPTION OF THE QUANTITIES COMPUTED ALONG A RAY (SEE C.R.T.5.2) C C LOCAL QUANTITIES (C.R.T.5.5.4): C YL(1)=VP... VELOCITY OF P-WAVES AT THE POINT. C YL(2)=VS... VELOCITY OF S-WAVES AT THE POINT. C YL(3)=RO... DENSITY AT THE POINT. C YL(4)=V1,YL(5)=V2,YL(6)=V3... VELOCITY DERIVATIVES IN GENERAL C COORDINATES. C C QUANTITIES COMPUTED ALONG A RAY (C.R.T.5.2): C BASIC QUANTITIES COMPUTED ALONG A RAY (C.R.T.5.2.1): C Y(1)... TRAVEL TIME. C Y(2)... IMAGINARY PART OF THE COMPLEX-VALUED TRAVEL TIME. C Y(3),Y(4),Y(5)... COORDINATES OF POINTS ALONG THE RAY. C Y(6),Y(7),Y(8)... COVARIANT COMPONENTS OF THE SLOWNESS VECTOR. C Y(9),Y(10),Y(11)... COVARIANT COMPONENTS OF THE POLARIZATION C VECTOR E1 PERPENDICALAR TO THE RAY. C ( Y(12),Y(16),Y(20),Y(24) ) RAY PROPAGATOR MATRIX (I.E. THE C ( Y(13),Y(17),Y(21),Y(25) ) ...MATRIX OF FUNDAMENTAL SOLUTIONS C ( Y(14),Y(18),Y(22),Y(26) ) OF THE DYNAMIC RAY TRACING SYSTEM. C ( Y(15),Y(19),Y(23),Y(27) ) C Y(28) TO Y(NY), WHERE NY=27+NAMPL... NAMPL REAL QUANTITIES C REPRESENTING COMPLEX-VALUED VECTORIAL REDUCED AMPLITUDES. C THE VECTORIAL REDUCED AMPLITUDES ARE SPECIFIED IN THE C RAY-CENTRED COORDINATE SYSTEM. C P WAVE AT THE INITIAL POINT OF THE RAY, C P WAVE AT THE POINT UNDER CONSIDERATION: C NAMPL=2, C Y(28)=REAL(A33), Y(29)=AIMAG(A33). C P WAVE AT THE INITIAL POINT OF THE RAY, C S WAVE AT THE POINT UNDER CONSIDERATION: C NAMPL=4, C Y(28)=REAL(A13), Y(29)=AIMAG(A13), C Y(30)=REAL(A23), Y(31)=AIMAG(A23). C S WAVE AT THE INITIAL POINT OF THE RAY, C P WAVE AT THE POINT UNDER CONSIDERATION: C NAMPL=4, C Y(28)=REAL(A31), Y(29)=AIMAG(A31), C Y(30)=REAL(A32), Y(31)=AIMAG(A32). C S WAVE AT THE INITIAL POINT OF THE RAY, C S WAVE AT THE POINT UNDER CONSIDERATION: C NAMPL=8, C Y(28)=REAL(A11), Y(29)=AIMAG(A11), C Y(30)=REAL(A21), Y(31)=AIMAG(A21), C Y(32)=REAL(A12), Y(33)=AIMAG(A12), C Y(34)=REAL(A22), Y(35)=AIMAG(A22). C Y(NY+1) TO Y(35), WHERE NY=27+NAMPL... UNDEFINED. C AUXILIARY QUANTITIES COMPUTED ALONG A RAY (C.R.T.5.2.2): C YY(1)...INDEPENDENT VARIABLE ALONG A RAY. C YY(2)=UEBRAY... UPPER ERROR BOUND FOR RAY TRACING, WHICH IS EQUAL C TO THE INPUT VALUE UEB AT THE INITIAL POINT OF THE RAY C (C.R.T.5.6J). IT IS ALWAYS DOUBLED WHEN THE NUMERICAL C INTEGRATION REQUIRES MORE THAN NHLF BISECTIONS OF THE C INITIAL INCREMENT STEP (C.R.T.5.6G). UEBRAY.GT.UEB AT C THE ENDPOINT OF THE RAY INDICATES A DECREASED ACCURACY C OF COMPUTATION. C YY(3)=ERRPP,YY(4)=ERRPH,YY(5)=ERRHH... DEVIATIONS (IN ABSOLUTE C VALUES OF THE TWO COMPUTED BASIS VECTORS OF THE C RAY-CENTRED COORDINATE SYSTEM FROM THE CONDITIONS OF C ORTHONORMALITY (C.R.T.5.8.2D AND 5.8.3I), ACCUMULATED C ALONG THE RAY. ANY OF THEM MAY BE COMPARED WITH THE C CORRESPONDING SPECIFIED LIMIT UEBPP, UEBPH, UEBHH AT THE C ENDPOINT OF THE RAY. C IY(1)=NY=27+NAMPL... NUMBER OF THE BASIC QUANTITIES C DESCRIBING THE POINT OF A RAY, SEE ABOVE. C IY(2)=KODIND... POSITION IN THE CODE (INDEX IN ARRAY KODE) C CORRESPONDING TO THE CONSIDERED ELEMENT OF A RAY. ITS C VALUE IS DETERMINED BY SUBROUTINE CODE (C.R.T.4.3). C IY(3)=ICB0... INDEX OF THE COMPLEX BLOCK, FROM WHICH THE RAY C ENTERED THE COMPLEX BLOCK IN WHICH THE COMPUTED C ELEMENT OF THE RAY IS SITUATED. C IY(3)=0 BEFORE LEAVING THE COMPLEX BLOCK, IN WHICH THE C INITIAL POINT OF THE RAY IS SITUATED. C IY(4)=ISB1... INDEX OF THE SIMPLE BLOCK CONTAINING THE C COMPUTED ELEMENT OF THE RAY. C IY(5)=ICB1... INDEX OF THE COMPLEX BLOCK CONTAINING THE C COMPUTED ELEMENT OF THE RAY, SUPPLEMENTED BY A SIGN '+' C FOR P WAVE AND SIGN '-' FOR S WAVE. C IY(6)=ISRF... INDEX OF THE SURFACE AT WHICH THE ENDPOINT OF THE C COMPUTED ELEMENT OF THE RAY IS SITUATED, SUPPLEMENTED BY C A SIGN '+' OR '-' FOR THE ENDPOINT SITUATED AT THE C POSITIVE OR NEGATIVE SIDE OF THE SURFACE, RESPECTIVELY. C ZERO INSIDE THE COMPLEX BLOCK. NOTE: THE SIGN OF THIS C QUANTITY IS THE EXEPTION TO THE DEFINITION IN THE ORIGINAL C PAPER ON C.R.T. C IY(7)=ISB2... INDEX OF THE SIMPLE BLOCK TOUCHING THE COMPLEX C BLOCK ICB1 FROM THE OTHER SIDE OF THE SURFACE ISRF AT C THE ENDPOINT OF THE COMPUTED ELEMENT OF THE RAY. C ISB2=0 FOR A FREE SPACE ON THE OTHER SIDE OF ISRF. C UNDEFINED INSIDE THE COMPLEX BLOCK, DEFINED ONLY AT C THE ENDPOINT OF THE ELEMENT OF THE RAY. C IY(8)=ICB2... INDEX OF THE COMPLEX BLOCK TOUCHING THE COMPLEX C BLOCK ICB1 FROM THE OTHER SIDE OF THE SURFACE ISRF AT C THE ENDPOINT OF THE COMPUTED ELEMENT OF THE RAY. C ICB2=0 FOR A FREE SPACE ON THE OTHER SIDE OF ISRF. C UNDEFINED INSIDE THE COMPLEX BLOCK, DEFINED ONLY AT C THE ENDPOINT OF THE ELEMENT OF THE RAY. C IY(9)=IFCT... NUMBER OF INVOCATIONS OF SUBROUTINE FCT EVALUATING C THE RIGHT-HAND SIDES OF THE ORDINARY DIFFERENTIAL C EQUATIONS, ALONG THE COMPUTED PART OF THE RAY. C IY(10)=IOUTP... NUMBER OF SUCCESSFUL STEPS OF THE NUMERICAL C INTEGRATION, ALONG THE RAY (I.E. THE NUMBER OF INVOCATIONS C OF SUBROUTINE OUTP DECREASED BY THE NUMBER OF INVOCATIONS C OF SUBROUTINE HPCG). C IY(11)=ITRANS... NUMBER OF TRANSFORMATIONS AT AN INTERFACE (I.E. C THE NUMBER OF INVOCATIONS OF SUBROUTINE TRANS). C IY(12)=KMAH... NUMBER OF CAUSTIC POINTS ALONG THE RAY (THE INDEX C OF THE RAY TRAJECTORY). C C DATE: 1990, JANUARY 23 C WRITTEN BY VLASTISLAV CERVENY, LUDEK KLIMES, IVAN PSENCIK C C======================================================================= C SUBROUTINE RAY1(LUN) INTEGER LUN C C SUBROUTINE RAY1 READS THE INPUT DATA FOR COMPLETE RAY TRACING (SEE C C.R.T.5.6) AND STORES THEM IN COMMON BLOCKS /RAYT/ AND /DCRT/. IT IS C CALLED ONCE, BEFORE THE COMPLETE TRACING OF INDIVIDUAL RAYS BEGINS. C SUBROUTINE RAY1 CALLS SUBROUTINE SRFC1 IN ORDER TO READ THE INPUT DATA C FOR THE AUXILIARY SURFACES, AND TO STORE THEM IN THE MEMORY. C C INPUT: C LUN... LOGICAL UNIT NUMBER OF THE EXTERNAL INPUT DEVICE C CONTAINING THE INPUT DATA. C THE INPUT PARAMETER IS NOT ALTERED. C C NO OUTPUT. C C COMMON BLOCK /MODELC/ (SUBROUTINE FILE 'MODEL.FOR'): INTEGER MSB,MCB PARAMETER (MSB=128) PARAMETER (MCB=128) INTEGER NEXPV,NEXPQ REAL BOUNDM(6) INTEGER NSRFCS,NSB,KSB(0:MSB),NCB,KCB(0:MSB) EQUIVALENCE (KSB(0),NSB),(KCB(0),NCB) COMMON/MODELC/NEXPV,NEXPQ,BOUNDM,NSRFCS,KSB,KCB C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED. C C COMMON BLOCK /RAYT/: CHARACTER*80 TEXTR COMMON/RAYT/TEXTR C COMMON BLOCK /DCRT/: INTEGER MEND,MSTOR PARAMETER (MEND=128) PARAMETER (MSTOR=128) INTEGER KSTORE,NEXPS,NHLF,MODCRT REAL STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT,BOUNDR(7) INTEGER NSRFCA,NEND,KEND(MEND),NSTOR,KSTOR(MSTOR) COMMON/DCRT/ KSTORE,NEXPS,NHLF,MODCRT,STORE,STEP,UEB,UEBPP,UEBPH, * UEBHH,UEBDRT,BOUNDR,NSRFCA,NEND,KEND,NSTOR,KSTOR C ALL THE STORAGE LOCATIONS OF THE COMMON BLOCKS ARE DEFINED IN THIS C SUBROUTINE. C C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED: INTEGER NSRFC EXTERNAL RAYB,NSRFC,SRFC1 C RAYB... BLOCK DATA SUBROUTINE OF THIS FILE. C NSRFC...FILE 'MODEL.FOR' OF THE PACKAGE 'MODEL'. C SRFC1 AND SUBSEQUENT ROUTINES... FILE 'SRFC.FOR' AND SUBSEQUENT C FILES (LIKE 'VAL.FOR' AND 'FIT.FOR') OF THE PACKAGE C 'MODEL'. C C ERROR MESSAGES: C 561... INSUFFICIENT MEMORY IN /DCRT/: C INSUFFICIENT MEMORY FOR THE INPUT DATA IN COMMON BLOCK C /DCRT/. THE DIMENSION MEND OF ARRAY KEND MUST BE C ENLARGED. SEE THE BLOCK DATA SUBROUTINE RAYB. C 562... INSUFFICIENT MEMORY IN /DCRT/: C INSUFFICIENT MEMORY FOR THE INPUT DATA IN COMMON BLOCK C /DCRT/. THE DIMENSION MSTOR OF ARRAY KSTOR MUST BE C ENLARGED. SEE THE BLOCK DATA SUBROUTINE RAYB. C 563... INADMISSIBLE END SURFACE: C THE INDEX OF AN END SURFACE IN INPUT DATA IS GREATER THAN C THE NUMBER OF SPECIFIED SURFACES. C 564... INADMISSIBLE STORE SURFACE: C THE INDEX OF A SURFACE FOR STORING THE COMPUTED QUANTITIES C IS GREATER THAN THE NUMBER OF SPECIFIED SURFACES AND WHILE C NOT SPECIFYING THE COMPUTATIONAL VOLUME BOUNDARY, OR THE C INDEX OF AN AUXILIARY SURFACE FOR STORING IS NEGATIVE. C C DATE: 1993, DECEMBER 18 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS: INTEGER I C C (1) STRING DESCRIBING THE DATA READ(LUN,*) TEXTR C (2,3) QUANTITIES CONTROLING NUMERICAL INTEGRATION MODCRT=0 READ(LUN,*) KSTORE,NEXPS,NHLF,MODCRT READ(LUN,*) STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT C (4) COORDINATE PLANES AND ISOCHRONE BOUNDING THE COMPUTATIONAL C VOLUME DO 10 I=1,6 BOUNDR(I)=BOUNDM(I) 10 CONTINUE BOUNDR(7)=999999. READ(LUN,*) BOUNDR C (5) NUMBER OF AUXILIARY SURFACES NSRFCA=0 READ(LUN,*) NSRFCA C C (6) END SURFACES: C INITIALIZING MEMORY FOR INDICES OF END SURFACES DO 11 I=1,MEND KEND(I)=0 11 CONTINUE C READING INDICES OF END SURFACES: READ(LUN,*) (KEND(I),I=1,MEND) DO 12 I=1,MEND IF(IABS(KEND(I)).GT.NSRFC()+NSRFCA) THEN PAUSE 'ERROR 563 IN RAY1: INADMISSIBLE END SURFACE' STOP ELSE IF(KEND(I).EQ.0) THEN NEND=I-1 GO TO 13 END IF 12 CONTINUE PAUSE 'ERROR 561 IN RAY1: INSUFFICIENT MEMORY IN /DCRT/' STOP 13 CONTINUE C C (7) SURFACES FOR STORING COMPUTED QUANTITIES: C INITIALIZING MEMORY FOR INDICES OF THE SURFACES DO 21 I=1,MSTOR KSTOR(I)=0 21 CONTINUE C READING INDICES OF THE SURFACES: READ(LUN,*) (KSTOR(I),I=1,MSTOR) DO 22 I=1,MSTOR/3+1 IF(KSTOR(I).LT.-NSRFC().OR.KSTOR(I).GT.NSRFC()+NSRFCA) THEN IF(KSTOR(I).LT.101.OR.KSTOR(I).GT.107) THEN PAUSE 'ERROR 564 IN RAY1: INADMISSIBLE STORE SURFACE' STOP END IF ELSE IF(KSTOR(I).EQ.0) THEN NSTOR=I-1 GO TO 23 END IF 22 CONTINUE PAUSE 'ERROR 562 IN RAY1: INSUFFICIENT MEMORY IN /DCRT/' STOP 23 CONTINUE C (8) SEQUENTIAL NUMBERS OF THE FIRST POINTS OF INTERSECTION DO 28 I=NSTOR+1,2*NSTOR KSTOR(I)=1 28 CONTINUE READ(LUN,*) (KSTOR(I),I=NSTOR+1,2*NSTOR) C (9) SEQUENTIAL NUMBERS OF THE LAST POINTS OF INTERSECTION DO 29 I=2*NSTOR+1,3*NSTOR KSTOR(I)=999999 29 CONTINUE READ(LUN,*) (KSTOR(I),I=2*NSTOR+1,3*NSTOR) C C (10) AUXILIARY SURFACES: IF(NSRFCA.GT.0) THEN CALL SRFC1(LUN,NSRFCA) END IF RETURN END C C======================================================================= C SUBROUTINE RAY2(YL,Y,YY,IY,IEND) REAL YL(6),Y(35),YY(5) INTEGER IY(12),IEND C C THIS SUBROUTINE CONTINUES IN THE COMPLETE TRACING OF A RAY FROM THE C GIVEN POINT (SEE C.R.T.5.7). FOR MODCRT.LE.2 (SEE THE INPUT DATA C (2)), THE WHOLE RAY IS COMPUTED. FOR MODCRT.GE.3, JUST ONE RAY C ELEMENT IS COMPUTED. C C INPUT: C YL... ARRAY CONTAINING LOCAL QUANTITIES AT THE GIVEN POINT. C Y... ARRAY CONTAINING BASIC QUANTITIES COMPUTED ALONG THE RAY C AT THE GIVEN POINT. C YY... ARRAY CONTAINING REAL AUXILIARY QUANTITIES COMPUTED ALONG C THE RAY AT THE GIVEN POINT. C IY... ARRAY CONTAINING INTEGER AUXILIARY QUANTITIES COMPUTED C ALONG THE RAY AT THE GIVEN POINT. C C OUTPUT: C YL... ARRAY CONTAINING LOCAL QUANTITIES AT THE ENDPOINT OF THE C RAY. C Y... ARRAY CONTAINING BASIC QUANTITIES COMPUTED ALONG THE RAY C AT THE ENDPOINT OF THE RAY. C YY... ARRAY CONTAINING REAL AUXILIARY QUANTITIES COMPUTED ALONG C THE RAY AT THE ENDPOINT OF THE RAY. C IY... ARRAY CONTAINING INTEGER AUXILIARY QUANTITIES COMPUTED C ALONG THE RAY AT THE ENDPOINT OF THE RAY. C IEND... REASON OF THE TERMINATION OF THE COMPUTATION OF A RAY (SEE C C.R.T.5.4). IEND=10,21,22,23,30 ARE INDICATED IN C SUBROUTINE CODE (FILE 'CODE.FOR'). IEND=24,25,26,32,33 C ARE INDICATED IN SUBROUTINE TRANS (FILE 'TRANS.FOR'). C THE REASONS IEND=40,50,60 ARE CHECKED AND INDICATED BY C IY(6) IN SUBROUTINE RAYCB (FILE 'RAYCB.FOR'). IEND=61 C IS INDICATED IN THIS SUBROUTINE RAY2 USING IY(6). C 0... THE COMPUTATION OF THE RAY MAY CONTINUE. C 10... RAY SATISFIES THE WHOLE CODE. C 21... THE POINT OF INCIDENCE IS SITUATED AT A DIFFERENT C SURFACE THAN THAT SPECIFIED BY THE CODE. C 22... THE NEXT ELEMENT OF THE RAY IS REQUIRED BY THE CODE C TO BE SITUATED IN A COMPLEX BLOCK THAT DOES NOT TOUCH C THE POINT OF INCIDENCE. C 23,24... TRANSMISSION IS REQUIRED BY THE CODE AT A FREE C SURFACE. C 25... RAY OF THE REQUIRED REFLECTED OR TRANSMITTED WAVE C IS NOT REAL-VALUED (OVERCRITICAL REFLECTION OR C TRANSMISSION). C 26... S WAVE IN A LIQUID BLOCK IS REQUIRED BY THE CODE. C 30,32... REFLECTION OR TYPE CONVERSION AT THE FICTITIOUS C PART OF THE INTERFACE. C 33... ZERO REFLECTION/TRANSMISSION COEFFICIENT. C 40... TRAVEL TIME GREATER THAN THE SPECIFIED LIMIT. C 50... RAY HAS INTERSECTED A COORDINATE PLANE LIMITING THE C COMPUTATIONAL VOLUME FOR COMPLETE RAY TRACING. C 60,61... RAY HAS INTERSECTED A SMOOTH SURFACE LIMITING THE C COMPUTATIONAL VOLUME FOR COMPLETE RAY TRACING. C OTHER VALUES OF IEND MAY INDICATE THE PROGRAM FAILURE. C C COMMON BLOCK /DCRT/: INTEGER MEND,MSTOR PARAMETER (MEND=128) PARAMETER (MSTOR=128) INTEGER KSTORE,NEXPS,NHLF,MODCRT REAL STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT,BOUNDR(7) INTEGER NSRFCA,NEND,KEND(MEND),NSTOR,KSTOR(MSTOR) COMMON/DCRT/ KSTORE,NEXPS,NHLF,MODCRT,STORE,STEP,UEB,UEBPP,UEBPH, * UEBHH,UEBDRT,BOUNDR,NSRFCA,NEND,KEND,NSTOR,KSTOR C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED. C C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED: EXTERNAL NSRFC,KOOR,CODE,RAYCB,TRANS,CONV,RPAR31,RPAR32,RPAR33 EXTERNAL WRIT31,WRIT32,WRIT33 INTEGER NSRFC,KOOR 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 CODE... FILE 'CODE.FOR'. C RAYCB AND SUBSEQUENT ROUTINES... FILE 'RAYCB.FOR'. C TRANS,CONV AND SUBSEQUENT ROUTINES... FILE 'TRANS.FOR'. C RPAR31,RPAR32,RPAR33... FILE 'RPAR.FOR'. C WRIT31,WRIT32,WRIT33 AND SUBSEQUENT ROUTINES... FILE 'WRIT.FOR'. C FCT,OUTP,CDEF,PSHIFT... THIS FILE. C C ERROR MESSAGES: C 570... WRONG FUNCTION OF SUBROUTINE CODE: C REASON 37 OR 38 OF THE TERMINATION OF THE RAY COMPUTATION C IS REPORTED BY THE SUBROUTINE TRANS, WHILE IT SHOULD BE C REPORTED PREVIOUSLY BY THE SUBROUTINE CODE AS THE REASON C 22 OR 23, RESPECTIVELY. C C DATE: 1990, NOVEMBER 18 C CODED BY LUDEK KLIMES C C----------------------------------------------------------------------- C C AUXILIARY STORAGE LOCATIONS: INTEGER KODIND,ICBNEW,I,NAMPL REAL YC(12) C C KODIND..POSITION IN THE CODE (INDEX IN THE ARRAY KODE) CORRES- C PONDING TO THE NEXT ELEMENT OF THE RAY. C ICBNEW..THE INDEX OF THE COMPLEX BLOCK IN WHICH THE NEXT C ELEMENT OF THE RAY IS TO BE SITUATED, SUPPLEMENTED C BY THE SIGN "+" FOR P WAVE OR "-" FOR S WAVE. C I... AUXILIARY LOOP VARIABLE. C NAMPL...NUMBER OF REAL QUANTITIES REPRESENTING COMPLEX-VALUED C VECTORIAL REDUCED AMPLITUDES. SEE THE SUBROUTINE CONV OF C THE SUBROUTINE FILE 'TRANS.FOR'. C YC... ARRAY CONTAINING REAL QUANTITIES REPRESENTING COMPLEX- C -VALUED VECTORIAL REDUCED AMPLITUDES. SEE C.R.T.5.5.4 AND C THE SUBROUTINE CONV OF THE SUBROUTINE FILE 'TRANS.FOR'. C C....................................................................... C 10 CONTINUE IF(IY(6).EQ.0) THEN C C (A1) ELEMENT OF THE RAY INSIDE A COMPLEX BLOCK CALL RAYCB(YL,Y,YY,IY) C C (B1) STORAGE OF THE COMPUTED QUANTITIES AT THE INTERFACE CALL RPAR31(YL,Y,YY,IY) IF(STORE.NE.0) THEN CALL WRIT31(YL,Y,YY,IY) END IF DO 21 I=1,NSTOR IF(KSTOR(I).EQ.IY(6)) THEN CALL RPAR32(I,YL,Y,YY,IY) CALL CONV(KSTORE,YL,Y,IY,NAMPL,YC) CALL WRIT32(I,YL,Y,YY,IY,NAMPL,YC) END IF IF(KSTOR(I).EQ.-IY(6).AND.KSTORE.GE.2) THEN IY(6)=-IY(6) CALL CONV(KSTORE,YL,Y,IY,NAMPL,YC) CALL WRIT32(I,YL,Y,YY,IY,NAMPL,YC) IY(6)=-IY(6) END IF 21 CONTINUE IF(IABS(IY(6)).LE.NSRFC()) THEN CALL RPAR33(YL,Y,YY,IY) CALL WRIT33(YL,Y,YY,IY) END IF C C (A2) CHECK FOR THE BOUNDARY OF THE COMPUTATIONAL VOLUME IF(IABS(IY(6)).GE.107) THEN IEND=40 GO TO 99 ELSE IF(IABS(IY(6)).GE.101) THEN IEND=50 GO TO 99 ELSE IF(IABS(IY(6)).GT.NSRFC()) THEN IEND=60 GO TO 99 ELSE IF(MODCRT.GE.3) THEN IEND=0 GO TO 99 END IF C ELSE C C (B2) INTERPRETATION OF THE ELEMENTARY WAVES CODE CALL CODE(IY,KODIND,ICBNEW,IEND) IF(IEND.NE.0) THEN GO TO 99 END IF C C (B3) CHECK FOR THE END SURFACES C NOTE: REFLECTION FROM AN END SURFACE IS ALLOWED. IF(IABS(IY(5)).NE.IABS(ICBNEW)) THEN DO 23 I=1,NEND IF(IABS(KEND(I)).EQ.IABS(IY(6))) THEN IEND=61 GO TO 99 END IF 23 CONTINUE END IF C C (B4) TRANSFORMATION ACROSS THE INTERFACE CALL TRANS(YL,Y,YY,IY,KODIND,ICBNEW,IEND) IF(IEND.GE.37) THEN PAUSE 'ERROR 570 IN RAY2: WRONG FUNCTION OF SUBROUTINE CODE' STOP ELSE IF(IEND.NE.0) THEN GO TO 99 END IF C C (B5) STORAGE OF THE COMPUTED QUANTITIES AT THE INTERFACE CALL RPAR31(YL,Y,YY,IY) IF(STORE.NE.0) THEN CALL WRIT31(YL,Y,YY,IY) END IF DO 25 I=1,NSTOR IF(KSTOR(I).EQ.IY(6)) THEN CALL RPAR32(I,YL,Y,YY,IY) IF(KSTORE.LE.1) THEN CALL CONV(KSTORE,YL,Y,IY,NAMPL,YC) CALL WRIT32(I,YL,Y,YY,IY,NAMPL,YC) END IF END IF 25 CONTINUE C C (B6) THE RAY MAY CONTINUE IN THE NEXT COMPLEX BLOCK IY(6)=0 C END IF GO TO 10 C 99 CONTINUE RETURN END C C======================================================================= C