C Subroutine file 'rpar.for' to control the take-off parameters of rays. C C Version: 7.20 C Date: 2015, June 1 C Coded by Ludek Klimes and Petr Bulant C C....................................................................... C C This file consists of the following subroutines: C RPAR1...Subroutine designed to read the input data for the C take-off parameters of rays and to store them in common C block /RPARD/, to prepare the parameters for the C subroutine RPAR2 and to store them in the common block C /RPARC/. It is called when starting the computation of a C new elementary wave. C RPAR1 C RPAR2...Subroutine designed to specify the take-off parameters of C individual rays. C RPAR2 C RPAR30..Subroutine performing the numerical quadrature of the C squared matrix of geometrical spreading, useful for the C estimination of the L2 norm of the distance between C neighbouring rays in order to control the fatness of ray C tubes. This subroutine is designed to be called by C subroutine PSHIFT of file 'raycb.for'. Subroutine PSHIFT C is called both with steps STEP and STORE along the ray. C RPAR30 C RPAR31..Subroutine called with constant step STORE of the C independent variable along the ray, and at the points of C intersection with interfaces either before and after the C transformation, i.e. called simultaneously with the C subroutine WRIT31, see also C C.R.T.5.5.1. C Subroutine RPAR31 accumulates the ray histories and C collects other quantities useful for the determination of C the take-off parameters in the memory. C RPAR31 C RPAR32..Subroutine called at the points of intersection of the ray C with specified interfaces, i.e. called simultaneously with C the subroutine WRIT32, see also C C.R.T.5.5.2. C It may be used to keep the quantities useful for the C determination of the take-off parameters in the memory. C RPAR32 C RPAR33..Subroutine called at the endpoints of the elements of C rays, situated at structural interfaces, i.e. called C simultaneously with the subroutine WRIT33, see also C C.R.T.5.5.3. C It may be used to keep the quantities useful for the C determination of the take-off parameters in the memory. C RPAR33 C RPAR4...Subroutine called after finishing the computation of each C ray. It defines four storage locations of the common C block /INITC/ introduced in the file 'init.for'. C RPAR4 C XFUN... Subroutine returning the value and first derivatives of C the Xi functions describing the profiles. It is called by C the subroutine RPAR32, and may be user-defined. In this C version, it calls the subroutine SRFC2 from the subroutine C file 'srfc.for', or directly returns one of the model C coordinates. C XFUN C IHIST...Auxiliary integer function to the subroutine RPAR4 to C determine the ray-history index. The ray-history index is C used to distinguish the rays with different histories for C the purposes of the subroutine RPAR2. C IHIST C SRPARH..Auxiliary subroutine to the subroutine RPAR31 to shift the C storage locations of the common block /RPARH/. C SRPARH C C Note: C The lines denoted by '*' in the first column contain an undebugged C code intended for future extensions. C C....................................................................... C C Description of data files: C C Input data RPAR for the take-off parameters of rays: C Data (1)-(3) common to all elementary waves should always be C located in the file specified by input SEP parameter C RPAR. C If optional input SEP parameter RPAREW=' ' (default), C data (4), repeated for all source points, should also be located C in file RPAR. C If optional input SEP parameter RPAREW is specified and is not C blank, data (4) for a single source should be located in file C RPAREW. C If there are several source points, data (4) are repeatedly used C for all source points. 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 TEXTP), the input parameter is of the C type REAL. C (1) TEXTP C String describing the data. Only the first 80 characters of the C string are significant. C C (2) ISRFR,ISRFX1,ISRFX2,XERR,AERR,PRM0(I),I=1,MPRM0,/ C ISRFR...Index of the reference surface at which the computed C quantities are stored and the given profiles or receivers C for two-point ray tracing are situated. The profiles are C the crosssections of the surface ISRFR and the isosurfaces C of the X1 function corresponding to given values (profile C coordinates). The receivers are the crosssections of the C surface ISRFR and the isosurfaces of the X1 and X2 C functions corresponding to given pairs of values (receiver C coordinates). C ISRFR should be one of the indices listed in the input C data DCRT-(7) (the indices for storing computed C quantities), including the sign, see the subroutine file C 'ray.for': C data DCRT-(7) C For ISRFR=0 there are no profiles for boundary-value ray C tracing, ISRFR not listed in the input data DCRT-(7) has C the same effect. C ISRFR=101,102,103,104,105, or 106 specify the boundary of C the computational volume. C Example: if the profile for the boundary-value ray tracing C is situated at the model (earth's) surface, ISRFR is the C index of the model surface. Its sign must correspond to C the side of the surface ISRFR facing the material block. C ISRFX1..Index of the function that becomes to be the X1 function. C The X1 function is defined here by the subroutine XFUN. C XFUN C ISRFX1=0: No boundary-value (two-point) ray tracing and C no determination of the boundaries between the ray C histories. C The X1 function is zero and has no meaning. C This option may be used only if the boundary-value rays C are not searched for, i.e. if ISRFR.EQ.0 or REC.EQ.' ' C or IPOINT=0 at all input data lines (4.1). C For ISRFX1=-1, -2, or -3, the X1 function coincides with C the -ISRFX1-th model coordinate. C For ISRFX1 positive, the subroutine XFUN is an interface C subroutine to the subroutine SRFC2, called with value C ISRFX1 assigned to the dummy argument ISRFC, see the C subroutine file 'srfc.for'. C C srfc.for C ISRFX1 must not exceed the number of the defined C surfaces more than by 1. C Hint: Even if there is no two-point ray tracing, specify, C e.g., ISRFX1=-1 (default). C ISRFX2..Index of the function that becomes to be the X2 function. C ISRFX2=0: Initial-value ray tracing (for ISRFX1=0) or C one-parametric shooting to profiles corresponding to C isosurfaces of the X1 function. The one-parametric C shooting to profiles reduces to the 2-D two-point ray C tracing in 2-D models in which the receivers at the C reference surface are determined by a single coordinate C (X1 function). C Otherwise: Two-parametric 3-D two-point ray tracing to C the receivers defined as the crosssections of the C surface ISRFR and the isosurfaces of the X1 and X2 C functions corresponding to given pairs of values C (receiver coordinates). C Meaning of ISRFX2 is then analogous to ISRFX1: C For ISRFX2=-1, -2, or -3, the X2 function coincides C with the -ISRFX2-th model coordinate. C For ISRFX2 positive, the subroutine XFUN is an interface C subroutine to the subroutine SRFC2, called with value C ISRFX2 assigned to the dummy argument ISRFC, see the C subroutine file 'srfc.for'. C C srfc.for C ISRFX2 must not exceed the number of the defined C surfaces (including ISRFX1) more than by 1. C Hint: Even if there is no two-point ray tracing, specify, C e.g., ISRFX2=-2 (default) for two-parametric shooting. C XERR... Specifies the accuracy of the boundary-value ray tracing. C If the distance of the ray from the given profile or C receiver does not exceed XERR, the ray is considered to be C a boundary-value ray. The distance of the ray from the C given profile is measured by means of the difference of C the X1 function. Similarly, the distance of the ray from C the given receiver is measured in units of the X1 and X2 C functions (e.g., length units for negative both ISRFX1 and C ISRFX2). C XERR should be reasonably greater than the maximum error C along any ray, converted to the units of XERR. C For the maximum error UEB per one step of numerical C integration refer to input data DCRT-(3) in 'ray.for'. C data DCRT-(3) C For two-parametric control of the take-off ray parameters C also rays distant more than XERR may be considered, see C parameter PRM0(5) below. C AERR... Specifies the accuracy of the determination of the C boundary rays between the different ray histories. C The distance between rays in the ray-parameter domain C is measured with respect to the metric tensor determined C by the input data (4) described below. The distance C between two consecutive rays of different histories C (boundary rays) is bisectioned until it does not exceed C AERR. Roughly speaking, the width of a demarcation belt C between different ray histories should not exceed AERR C times the size of the basic triangles in the ray-parameter C domain. C No boundary rays are determined for AERR.GT.1 (e.g. C AERR=999999.). Note that IPOINT=0 has the same effect as C AERR.GT.1. C PRM0(1)... Maximum allowed distance of the boundary ray from the C shadow zone (measured on the reference surface). C The distance is not measured if PRM0(1)=0 (default). C Used only by two-parametric control of the take-off ray C parameters. C PRM0(2)... Maximum allowed length of sides of homogeneous C triangles (measured on the reference surface). C The triangles are not measured if PRM0(2)=0 (default). C Used only by two-parametric control of the take-off ray C parameters. C PRM0(3)... Controls the determination of the ray histories. C PRM0(3).EQ.0: C Rays terminating at different boundaries of the C computational volume have different histories. C This option is suitable for two-point ray tracing. C PRM0(3).EQ.0 corresponds to CRT version 5.00 and is the C default in version 5.10. C PRM0(3).NE.0 (e.g., PRM0(3)=1): C Only structural interfaces and the reference surface C (if specified) influence ray histories of unsuccessful C rays. The boundaries of the computational volume and C end surfaces are not taken into account. This option C limits the number of thin spaces, uncovered by ray C tubes, situated between ray tubes corresponding to C different ray histories. C This option is thus suitable for controlled C initial-value ray tracing followed by the travel-time C interpolation within ray tubes. In such a case, it is C also recommended to set IPOINT=999999 (default) in data C (4.1) below in order to disable successful rays. C PRM0(3).EQ.1 corresponds to CRT version 4.10, C PRM0(3).GE.1.5 (e.g., PRM0(3)=2): C Index of the surface, at which the ray terminates, is C removed from ray histories of unsuccessful rays. C This option further limits the number of uncovered thin C spaces between ray tubes corresponding to different ray C histories and is thus suitable for controlled C initial-value ray tracing followed by the travel-time C interpolation within ray tubes. C It is recommended to set IPOINT=999999 (default) in data C (4.1) below in order to disable successful rays. C PRM0(3).GE.2.5 (e.g., PRM0(3)=4): C Reason of the termination of the computation of a ray is C removed from ray histories. This option further limits C the number of uncovered thin spaces between ray tubes C corresponding to different ray histories and is thus C suitable for controlled initial-value ray tracing C followed by the travel-time interpolation within ray C tubes. C There are 2 suboptions: C PRM0(3).LT.3.5 (e.g., PRM0(3)=3): C Distinguishing overcritical incidence (IEND=25) from C other reasons of the termination of the computation of C a ray (recommended for interpolations within ray C cells). C PRM0(3).GE.3.5 (e.g., PRM0(3)=4): C Reason of the termination of the computation of a ray C is removed from ray histories at all. C PRM0(4)... Controls the maximum length of the sides of homogeneous C triangles in the ray-tube metric. The ray-tube C metric contains an information about the geometrical C spreading corresponding to the ray. Maximum length C of the homogeneous triangle's sides in the C ray-tube metric thus controls a width of the C corresponding ray tube. C PRM0(4) is thus intended to be used for controlled C initial-value ray tracing followed by the interpolation C within ray tubes. C There are 3 options: C PRM0(4)=0: Reccommended for two-point ray tracing. C Measuring of triangles in the ray-tube metric is C disabled. C PRM0(4).GT.0: Reccommended for controlled initial-value C ray tracing. Maximum length of the sides of triangles C in the ray-tube metric is PRM0(4), the metric C describes an average squared thickness of ray tubes. C This means, that in a homogeneous model, e.g., the C maximum thickness of the ray tubes will be about C PRM0(4)*SQRT(3). C PRM0(4).LT.0: Maximum length of the sides of triangles C in the ray-tube metric is ABS(PRM0(4)), the metric C describes a thickness of ray tubes at their ends. This C may be useful, e.g., for controlled initial-value C ray tracing of a weakly refracted wave. C Used only by two-parametric control of the take-off ray C parameters. C PRM0(5)... Controls the maximum distance of the boundary-value C rays from the given receivers. If the program fails C to find a ray distant less then XERR, the nearest ray C is found. If the nearest ray is distant less than C PRM0(5), it is considered to be a two-point ray and C a warning is generated. C PRM0(5), if specified, should be positive and greater C than XERR. The default PRM0(5)=0. is understood as C infinite PRM0(5), and the nearest ray is always C considered to be a two-point ray. C Used only by two-parametric two-point ray tracing. C PRM0(I),I=6,MPRM0... Reserved for future extension of the C numerical parameters for two-point ray tracing. C Default values: C ISRFR=0, ISRFX1=-1, ISRFX2=-2, XERR=0.000, AERR=999999, C PRM0(I)=0,I=1,MPRM0. C The boundary rays are not determined if AERR.GT.1. C Neither are they if IPOINT=0, see input data (4.1). C The boundary-value rays are not determined if ISRFR.EQ.0, C if ISRFR is not listed in the input data DCRT-(7), C or if REC.EQ.' '. C Neither are they if IPOINT=0, see input data (4.1). C Thus, just the initial-value ray tracing is performed if C (ISRFR.EQ.0 or ISRFR not in DCRT-(7) or REC.EQ.' ') C and AERR.GT.1. C Just the initial-value ray tracing is also performed if C IPOINT=0, see input data (4.1). C During the initial-value ray tracing, only the system of C basic rays with the given constant steps is generated. C (3) The data specifying the functions No. ISRFX1 and ISRFX2. This C auxiliary function is specified by the subroutine SRFC2 and the C data are read by subroutine SRFC1. For the description of the C data refer to the subroutine SRFC1 of the file 'srfc.for'. This C input is performed only if ISRFX1 or ISRFX2 is positive and the C function ISRFX1 or ISRFX2 is not defined in the data sets C MODEL or DCRT. C In such a case, the index of the first C undefined function of ISRFX1 or ISRFX2 must equal the number NSRFC C of the surfaces defined in the data set MODEL + the number C ISRFCA of the auxiliary surfaces defined in the data set C DCRT + 1, and so on. C 'srfc.for' C Example: Imagine that the profile will be constructed as the C intersection of the surface ISRFR with an auxiliary C vertical plane. Using the basic version of the file C 'srfc.for', the auxiliary vertical plane passing through C the points (X1,X2,X3)=(X1A,X2A,X3A) and C (X1,X2,X3)=(X1B,X2B,X3B) may be defined using subroutine C file 'srfc.for' by the following input data: C For X1A.NE.X1B: C 1 -2 0 0 C 2 C X1A X1B C X2A X2B C The interpolated function then returns the oriented C distance from the vertical plane, measured in the C direction of the X2-axis but with opposite sign. C For X2A.NE.X2B: C 2 -1 0 0 C 2 C X2A X2B C X1A X1B C The interpolated function then returns the oriented C distance from the vertical plane, measured in the C direction of the X1-axis but with opposite sign. C C (4) For each elementary wave IWAVE the following data (4.1) and (4.2): C (4.1) Any times: C PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM,IPOINT, C PRM1(I),I=1,MPRM1,/ C PAR1L,PAR2L... Ray parameters corresponding to the origin (0,0) of C the domain of the normalized ray parameters. C The kind of ray parameters used is given by index INIPAR C read by init.for C and described in C crtin.for. C Imagine the normalized ray parameters as two coordinates C ranging from 0 to 1. The normalized ray-parameter domain C is then a unit square. Normalized ray parameters (0,0) C are mapped to ray parameters (PAR1L,PAR2L), normalized ray C parameters (1,0) are mapped to ray parameters C (PAR1A,PAR2A) and normalized ray parameters (0,1) are C mapped to ray parameters (PAR1B,PAR2B). The mapping is C linear. Any couple of normalized ray parameters can thus C be simply mapped to the ray parameters (e.g., to the C colatitude and longitude). C PAR1A,PAR2A... Ray parameters corresponding to the right bottom C corner (1,0) of the domain of the normalized ray C parameters. C PAR1B,PAR2B... Ray parameters corresponding to the left top corner C (0,1) of the domain of the normalized ray parameters. C ANUM,BNUM... Specify the 2*2 metric tensor G in the normalized ray C parameter domain. C The unit square of the normalized ray parameters is C divided into ANUM times BNUM rectangles. C For initial-value ray tracing or one-parametric shooting, C the normalized ray parameters corresponding to the C vertices of the rectangles are mapped to the ray C parameters and the corresponding basic rays are traced. C For two-parametric shooting, the size of the rectangles C define the metric for the coverage of the ray-parameter C domain by a triangularized system of basic rays. C ANUM,BNUM in RP2D (Initial-value ray tracing or one-parametric C shooting): C The components of the metric tensor are C G11=ANUM*ANUM, G12=0., G22=BNUM*BNUM. The system of basic C rays forms a regular rectangular grid in the normalized C ray parameter domain. The grid cells are unit squares C with respect to the above metric: C For ANUM.GT.0 and BNUM.GT.0: the two-parametric system of C INT(ANUM)*INT(ANUM) basic rays is given by the take-off C ray parameters C PAR1=PAR1L+A*(PAR1A-PAR1L)/ANUM+B*(PAR1B-PAR1L)/BNUM, C PAR2=PAR2L+A*(PAR2A-PAR2L)/ANUM+B*(PAR2B-PAR2L)/BNUM, C where C A=0,1,2,...,INT(ANUM) and B=0,1,2,...,INT(BNUM). C The system of basic rays corresponds to the C initial-value ray tracing. C Besides the system of basic rays, also additional rays C are computed in order to find the boundary rays between C the bundles of rays with different histories, and to C find the boundary-value rays. All additional rays are C given by C B=0,1,2,...,INT(BNUM) C and real A from the interval 0 to INT(ANUM), forming C the cuts of the ray-parameter surface given by fixed B C and parametrized by A. The boundary rays and the C boundary-value rays are searched for along the mentioned C cuts. Under the boundary-value rays we understand here C the rays from the above mentioned cuts passing through C the XERR-vicinity of the given profiles. C For ANUM.GT.0 and BNUM.EQ.0: The one-parametric system of C INT(ANUM) basic rays is given by the take-off ray C parameters C PAR1=PAR1L+A*(PAR1A-PAR1L)/ANUM, C PAR2=PAR2L+A*(PAR2A-PAR2L)/ANUM, C where C A=0,1,2,...,INT(ANUM). C Additional rays are computed as for ANUM.GT.0 and C BNUM.GT.0. C For ANUM.EQ.0 and BNUM.GT.0: Similarly for ANUM.GT.0 and C BNUM.EQ.0, but no additional rays are computed. C For ANUM.EQ.0 and BNUM.EQ.0: Just a single ray given by C PAR1=PAR1L and PAR2=PAR2L is computed. C ANUM,BNUM in RP3D (two-parametric shooting): C ANUM,BNUM must be positive for two-parametric shooting. C The components of the metric tensor are C G11=ANUM*ANUM*E11, G12=ANUM*BNUM*E12, G22=BNUM*BNUM*E22, C where E is some basic metric tensor, e.g. the unit matrix C or the Cartesian metric tensor on the unit sphere around C the source. The system of basic rays forms basic C triangles covering the normalized ray parameter domain. C The basic triangles are constructed in a way to be roughly C equilateral with unit sides with respect to metric G. C Recommendations for the choice of ANUM and BNUM: C For strong lateral heterogeneities, the ratio of ANUM C and BNUM should roughly correspond to the ratio of the C lengths of sides of the ray-parameter domain. C For moderate lateral heterogeneities, the number C corresponding to the longitude in spherical coordinates C may be smaller. C The initial value of ANUM*BNUM may be chosen, e.g., C between 20*20 to 40*40, or roughly corresponding to C ANUM*BNUM=(1/AERR)**2 for reasonable values of AERR. C ANUM and BNUM may then be adjusted according to the C complexity of the ray tracing. C If ANUM*BNUM is great, especially if it is approaching C 30 or even 50 per cent of the number of all rays traced C for the elementary wave, ANUM and BNUM should be C decreased. C If ANUM*BNUM is small with respect to the number of the C rays traced for the elementary wave, ANUM and BNUM may C be increased. C IPOINT..IPOINT=0: No ray history is recorded, all rays have the C same index 0 and no boundary-value ray tracing is C performed. C Otherwise: The IABS(IPOINT)-th point of intersection of C a ray with the surface ISRFR is considered for the C boundary-value ray tracing. C IPOINT negative: Only those rays satisfying the whole C code (IEND=10) are considered to be successful. All C other rays are rejected (IMPORTANT OPTION). C Another effect of this option: C If a ray terminates at a coordinate plane limiting the C computational volume (IEND=50), subroutine CODE is C called once again to check whether the ray satisfies C the whole code. C If yes, IEND=50 is replaced with IEND=10. C IPOINT=large (e.g. IPOINT=999999): C Ray history is recorded, boundary rays may C be determined but no boundary-value ray tracing is C performed because all rays are unsuccessful. C This is the recommended option for controlled C initial-value ray tracing followed by the travel-time C interpolation within ray tubes. C In most cases, user will wish the value of IPOINT=1. C For the structural interface (not auxiliary surface), the C positive or negative side of the surface ISRFR is C specified, then: C No point of intersection is accounted if the reflected C ray is situated at the other side of the surface ISRFR C than the specified one, and C two points of intersection are accounted if the C reflected ray is situated at the specified side of the C surface ISRFR. C Just the initial-value ray tracing is performed if C IPOINT=0. C PRM1(I),I=1,MPRM1... Reserved for future extension of the C numerical parameters for two-point ray tracing, specific C to the calculated elementary wave. C Default values: C PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM C Zeroes when starting the complete ray tracing program, C otherwise the previous values. C IPOINT=999999 when starting the complete ray tracing program, C otherwise the previous value. C PRM1(I)=0, I=1,MPRM1. C (4.2) /(a slash) C C Example of data RPAR for initial-value ray tracing. C C Example of data RPAR for 3-D two-point ray tracing. C C C Input file 'REC' containing profile (or receiver) coordinates: C (1) Several strings terminated by / (a slash). C The simplest way is to submit just the /. C (2) Several times (2.1): C (2.1) 'NAMER(IR)',XR(1,IR),XR(2,IR),XR(3,IR),XR(4,IR),XR(5,IR),/ C 'NAMER(IR)'... String reserved for the name of the receiver. C No meaning here, anything in apostrophes may be submitted. C XR(1,IR),XR(2,IR),XR(3,IR)... Model coordinates of the receiver. C For ISRFX1.LT.0: XREC1(IR)=XR(-ISRFX1,IR) is the first C coordinate of the point along the reference surface. C For ISRFX2.LT.0: XREC2(IR)=XR(-ISRFX2,IR) is the second C coordinate of the point along the reference surface. C Values of other model coordinates have no meaning here. C XR(4,IR),XR(5,IR)... Coordinates of the receiver along the C reference surface other than the model coordinates. C For ISRFX1.GT.0: XREC1(IR)=XR(4,IR) is the first C coordinate of the point along the reference surface. C For ISRFX2.GT.0: XREC2(IR)=XR(5,IR) is the second C coordinate of the point along the reference surface. C Other values have no meaning here and need not be C specified. C Default: X1R(IR)=0, X2R(IR)=0, X3R(IR)=0, X4R(IR)=0, X5R(IR)=0. C (3) / (a slash) or end of file. C C Example of data file REC C C....................................................................... C C Storage in the memory: C The input data (2) are stored in common block /RPARD/ defined C in the include file 'rpard.inc'. C rpard.inc C Common block /RPARC/ to collect the quantities needed for C boundary-value ray tracing and common block /RPARH/ describing the C histories of rays are defined in the include file 'rparc.inc'. C rparc.inc C The histories of rays are used by function IHIST to determine the C ray-history index. Different values of ray-history index are C attached to rays with different histories. C C======================================================================= C C C SUBROUTINE RPAR1(LUN,IWAVE) INTEGER LUN,IWAVE C C This subroutine is called when starting the computation of a new C elementary wave. It prepares the parameters for the subroutine RPAR2 C to specify the take-off parameters of the individual rays. C C Input: C LUN... Logical unit number of the external input device C containing the input data. C IWAVE...Index of the elementary wave that is going to be computed, C i.e. the output from the last invocation of the subroutine C CODE1. C C No output. C C Common block /DCRT/: INCLUDE 'dcrt.inc' C dcrt.inc C None of the storage locations of the common block are altered. C C Common block /RPARD/: INCLUDE 'rpard.inc' C rpard.inc C Nearly all the storage locations of the common block are defined in C this subroutine. C C Common block /RPARC/: INCLUDE 'rparc.inc' C rparc.inc C None of the storage locations of the common block, except LURPAR, are C altered. LURPAR is set to LUN. C C Subroutines and external functions required: EXTERNAL NSRFC,RSEP3T INTEGER NSRFC C NSRFC... File 'model.for' of the package 'MODEL'. C C Date: 1999, September 23 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER LUREC PARAMETER (LUREC=9) CHARACTER*80 TEXTP INTEGER I,J REAL AUX(5) C C LUREC...Logical unit number of the receiver file. The file is C opened and closed here. C TEXTP...The name of the data. String of 80 characters. C I,J... Auxiliary loop variables, or temporary variables. C AUX... Temporary storage locations to read the receiver C coordinates. C C....................................................................... C IF(IWAVE.EQ.0) THEN C C (1) Reading the name of the input data: READ(LUN,*) TEXTP C (2): ISRFR=0 ISRFX(1)=-1 ISRFX(2)=-2 XERR=0.000 AERR=999999. DO 10 I=1,MPRM0 PRM0(I)=0. 10 CONTINUE READ(LUN,*) ISRFR,ISRFX(1),ISRFX(2),XERR,AERR,PRM0 C C (3) Auxiliary functions no. ISRFX(1) and ISRFX(2): IF(ISRFX(1).EQ.0.AND.ISRFX(2).NE.0) THEN C 642 CALL ERROR('642 in RPAR1: Zero ISRFX1, nonzero ISRFX2') C For initial-value ray tracing given by ISRFX1=0, there must be C ISRFX2=0 in the input data set RPAR. END IF I=NSRFC()+NSRFCA+1 IF(ISRFX(1).EQ.I) THEN CALL SRFC1(LUN,1) I=I+1 ELSE IF(ISRFX(1).LT.-3.OR.ISRFX(1).GT.I) THEN C 643 CALL ERROR('643 in RPAR1: Wrong value of ISRFX1') C The function no. ISRFX1 specified by the routines of the file C 'srfc.for' is defined neither in the input data set C MODEL nor DCRT C and cannot be defined in this data set. C The index ISRFX1 is too large or is less than -3. END IF IF(ISRFX(2).EQ.I) THEN CALL SRFC1(LUN,1) ELSE IF(ISRFX(2).LT.-3.OR.ISRFX(2).GT.I) THEN C 644 CALL ERROR('644 in RPAR1: Wrong value of ISRFX2') C The function no. ISRFX2 specified by the routines of the file C 'srfc.for' is defined neither in the input data set C MODEL nor DCRT C and cannot be defined in this data set. C The index ISRFX2 is too large or is less than -3. END IF C C Default values of data (5): PAR1L=0. PAR2L=0. PAR1A=0. PAR2A=0. PAR1B=0. PAR2B=0. ANUM=0. BNUM=0. IPOINT=999999 C C Receivers: CALL RSEP3T('REC',TEXTP,' ') IF(TEXTP.EQ.' ') THEN NREC=0 ELSE OPEN(LUREC,FILE=TEXTP,STATUS='OLD') READ(LUREC,*) (TEXTP,I=1,20) DO 11 I=1,MREC+1 TEXTP='$' AUX(1)=0. AUX(2)=0. AUX(3)=0. AUX(4)=0. AUX(5)=0. READ(LUREC,*,END=12) TEXTP,(AUX(J),J=1,5) IF(TEXTP.EQ.'$') THEN GO TO 12 END IF IF(I.GT.MREC) THEN C 641 CALL ERROR * ('641 in RPAR1: Too small array XREC in /RPARD/') C The number of given receivers or profiles is exceeding the C dimension MREC of the array XREC in the common block C /RPARD/. MREC should thus be adjusted accordingly. C /RPARD/ END IF IF(ISRFX(1).LT.0) THEN XREC(1,I)=AUX(-ISRFX(1)) ELSE IF(ISRFX(1).GT.0) THEN XREC(1,I)=AUX(4) END IF IF(ISRFX(2).LT.0) THEN XREC(2,I)=AUX(-ISRFX(2)) ELSE IF(ISRFX(2).GT.0) THEN XREC(2,I)=AUX(5) END IF 11 CONTINUE 12 CONTINUE NREC=I-1 CLOSE(LUREC) END IF C END IF LURPAR=LUN RETURN END C C======================================================================= C C C SUBROUTINE RPAR2(IRAY,PAR1,PAR2) REAL PAR1,PAR2 INTEGER IRAY C C This subroutine determines the take-off parameters of the ray. C C Input: C IRAY... Number of the already computed rays. IRAY=0 when C beginning the computation of a new elementary wave. C Otherwise, the output from the previous invocation of C RPAR2. C C Output: C IRAY... IRAY=0 if all rays are computed and the computation of the C elementary wave will be terminated. C Otherwise, input value increased by 1. C PAR1,PAR2... Take-off parameters of the ray. C C Common block /RPARD/: INCLUDE 'rpard.inc' C rpard.inc C Storage locations PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM, C IPOINT,PRM1 of the common block may be redefined in this subroutine. C C Common blocks /RPARC/ and /RPARH/: INCLUDE 'rparc.inc' C rparc.inc C Storage locations IRAY0,JRAY,JTYPE,G1,G2 of common block /RPARC/ are C altered. C Storage locations JPOINT and KMAH of common block /RPARH/ are C annulled. C Storage locations NRPARH and IRPARH(0) of common block /RPARH/ may be C altered. C C Subroutines and external functions required: EXTERNAL LUWARN,RP2D INTEGER LUWARN C LUWARN..File error.for. C RP2D... One-parametric boundary-value ray tracing. This file. C C Date: 2013, February 19 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I,IAUX REAL AUX1,AUX2,AUX3,AUX4,AUX5,AUX6,AUX7,AUX8 C C....................................................................... C JPOINT=0 KMAH=0 C IF(IRAY.EQ.0) THEN C New elementary wave IRAY0=0 IF(ISRFX(2).EQ.0) THEN G2=0. END IF GO TO 90 END IF C C New ray - determination of shooting parameters G1,G2: 10 CONTINUE JRAY=IRAY-IRAY0 IF(ISRFX(2).EQ.0) THEN C One-parametric shooting (or initial-value ray tracing): CALL RP2D(IRAY0,JRAY,G1) IF(JRAY.EQ.0) THEN C End of shooting - new one-parametric slice or new data G2=G2+1. IF(G2.GT.BNUM+0.001) THEN C New data or new elementary wave GO TO 90 ELSE C New one-parametric slice IRAY0=IRAY GO TO 10 END IF END IF ELSE C Two-parametric shooting: CALL RP3D(JRAY,JTYPE,G1,G2) IF(JRAY.EQ.0) THEN C End of shooting - new data or new elementary wave GO TO 90 END IF END IF IRAY=IRAY0+JRAY C C New ray - determination of take-off parameters PAR1,PAR2: NRPARH=IRPARH(IRPARH(0)) IF(ISRFX(2).EQ.0) THEN C One-parametric shooting: C Projection of the shooting domain (0,ANUM)*(0,BNUM) PAR1=PAR1L PAR2=PAR2L IF(ANUM.NE.0.) THEN PAR1=PAR1+(PAR1A-PAR1L)*G1/ANUM PAR2=PAR2+(PAR2A-PAR2L)*G1/ANUM END IF IF(BNUM.NE.0.) THEN PAR1=PAR1+(PAR1B-PAR1L)*G2/BNUM PAR2=PAR2+(PAR2B-PAR2L)*G2/BNUM END IF ELSE C Two-parametric shooting: C Projection of the shooting domain (0,1)*(0,1) PAR1=PAR1L+(PAR1A-PAR1L)*G1+(PAR1B-PAR1L)*G2 PAR2=PAR2L+(PAR2A-PAR2L)*G1+(PAR2B-PAR2L)*G2 END IF RETURN C C Reading input data 90 CONTINUE AUX1=PAR1L AUX2=PAR2L AUX3=PAR1A AUX4=PAR2A AUX5=PAR1B AUX6=PAR2B AUX7=ANUM AUX8=BNUM IAUX=IPOINT DO 91 I=1,MPRM1 PRM1(I)=0. 91 CONTINUE READ(LURPAR,*) * PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM,IPOINT,PRM1 IF(IRAY.NE.0) THEN IF(PAR1L.EQ.AUX1.AND.PAR2L.EQ.AUX2.AND.PAR1A.EQ.AUX3.AND. * PAR2A.EQ.AUX4.AND.PAR1B.EQ.AUX5.AND.PAR2B.EQ.AUX6.AND. * ANUM.EQ.AUX7.AND.BNUM.EQ.AUX8.AND.IAUX.EQ.IPOINT) THEN C New elementary wave IRAY=0 C Writing the table of ray histories: WRITE(LUWARN(0),'(A)') 'Table of ray histories:' DO 92 I=1,IRPARH(0) WRITE(LUWARN(0),'(I4,A,I3,19I4,(4X,19I4))') * INDMIN+I,':',(IRPARH(IAUX),IAUX=IRPARH(I-1)+1,IRPARH(I)) 92 CONTINUE RETURN END IF END IF IF(ISRFX(1).NE.0) THEN IF(PAR1L.EQ.PAR1A.AND.PAR2L.EQ.PAR2A) THEN C 645 CALL ERROR('645 in RPAR2: Zero ray-parameter domain') C (PAR1A,PAR2A) cannot equal (PAR1L,PAR2L) for one or two C parametric boundary-value ray tracing indicated by non-zero C ISRFX1. END IF END IF IF(ISRFX(2).NE.0) THEN IF(PAR1L.EQ.PAR1B.AND.PAR2L.EQ.PAR2B) THEN C 646 CALL ERROR('646 in RPAR2: Zero ray-parameter domain') C (PAR1B,PAR2B) cannot equal (PAR1L,PAR2L) for two-parametric C boundary-value ray tracing indicated by non-zero ISRFX2. END IF IF(ANUM.LE.0..OR.BNUM.LE.0.) THEN C 653 CALL ERROR('653 in RPAR2: Non-positive ANUM or BNUM') C Both ANUM and BNUM must be positive for two-parametric C boundary-value ray tracing indicated by non-zero ISRFX2. END IF END IF IF(ISRFX(1).EQ.0.AND.ISRFR.NE.0.AND.NREC.NE.0.AND.IPOINT.NE.0) * THEN C 647 CALL ERROR('647 in RPAR2: Zero value of ISRFX1') C ISRFX1 must not be zero if the boundary-value ray tracing C is to be performed, i.e. if ISRFR.NE.0 and NREC.NE.0 and C IPOINT.NE.0 in the input data. END IF C Initializing ray histories INDMIN=0 IRPARH(0)=0 C New initial ray IRAY0=IRAY G2=0. GO TO 10 C END C C======================================================================= C C C SUBROUTINE RPAR30(TTINI,TT,Q11,Q21,Q12,Q22) REAL TTINI,TT,Q11,Q21,Q12,Q22 C C This subroutine performs the numerical quadrature of the squared C matrix of geometrical spreading, useful for the estimination of the C L2 norm of the distance between neighbouring rays in order to control C the fatness of ray tubes. This subroutine is designed to be called C by subroutine PSHIFT of file 'raycb.for'. Subroutine PSHIFT is called C both with steps STEP and STORE along the ray. C C Input: C TTINI...Travel time at the initial point of the ray, used to C identify the inititial point of a ray. C TT... Travel time (independent variable for the numerical C quadrature along the ray). C Q11,Q21,Q12,Q22... Matrix of geometrical spreading. C None of the input parameters are altered. C C No output. C C Common block /RPARD/: INCLUDE 'rpard.inc' C rpard.inc C None of the storage locations of the common block are altered. C C Common block /RPARC/: INCLUDE 'rparc.inc' C rparc.inc C Storage locations of the common block may be altered. C C No subroutines and external functions required. C C Date: 1999, May 25 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C REAL TT1,GG111,GG112,GG122,AUX REAL TT2,GG211,GG212,GG222 SAVE TT2,GG211,GG212,GG222 C IF(TT.EQ.TTINI) THEN GG11=0. GG12=0. GG22=0. ELSE TT1 =TT2 GG111=GG211 GG112=GG212 GG122=GG222 END IF TT2 =TT GG211=Q11*Q11+Q21*Q21 GG212=Q11*Q12+Q21*Q22 GG222=Q12*Q12+Q22*Q22 IF(TT.NE.TTINI) THEN IF(PRM0(4).LT.0) THEN C Metric at the last point GG11=GG211 GG12=GG212 GG22=GG222 ELSE C Average squared metric AUX=(TT2-TT1)/2. GG11=GG11+(GG111+GG211)*AUX GG12=GG12+(GG112+GG212)*AUX GG22=GG22+(GG122+GG222)*AUX END IF END IF RETURN END C C======================================================================= C C C SUBROUTINE RPAR31(YL,Y,YY,IY) REAL YL(6),Y(35),YY(5) INTEGER IY(13) C C This subroutine may be used to store the quantities useful for the C determination of the take-off parameters in the memory. It is called C with constant step STORE (see the input data in the subroutine file C 'ray.for') of the independent variable along the ray (if STORE.NE.0), C and at the points of intersection with interfaces either before and C after the transformation (in any case). C C Input: C YL... Array containing local quantities at the point of the ray. C Y... Array containing basic quantities computed along the ray. C YY... Array containing real auxiliary quantities computed along C the ray. C IY... Array containing integer auxiliary quantities computed C along the ray. C None of the input parameters are altered. C C No output. C C Common block /RPARD/: INCLUDE 'rpard.inc' C rpard.inc C None of the storage locations of the common block are altered. C C Common block /RPARH/: INCLUDE 'rparc.inc' C rparc.inc C Storage locations of the common block may be altered. C C Subroutines and external functions required: EXTERNAL NSRFC INTEGER NSRFC C NSRFC... File 'model.for' of the package 'MODEL'. C C Date: 2014, February 13 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C No auxiliary storage locations. C IF(IPOINT.NE.0) THEN IF(NRPARH.EQ.IRPARH(IRPARH(0))) THEN C Initial element of the ray NRPARH=NRPARH+1 CALL SRPARH IRPARH(NRPARH)=IY(5) END IF IF(JPOINT.LT.IABS(IPOINT)) THEN IF(IY(6).NE.0.AND.(IABS(IY(6)).LE.NSRFC().OR. * PRM0(3).EQ.0.0.OR. * PRM0(3).GE.1.5 )) THEN C Structural interface (or model boundary) NRPARH=NRPARH+1 CALL SRPARH IF(MOD(NRPARH-IRPARH(IRPARH(0)),2).EQ.0) THEN C Incidence at a structural interface IRPARH(NRPARH)=IY(6) ELSE C Leaving a structural interface IRPARH(NRPARH)=IY(5) END IF END IF END IF END IF RETURN END C C======================================================================= C C C SUBROUTINE RPAR32(ISTOR,YL,Y,YY,IY) INTEGER ISTOR,IY(13) REAL YL(6),Y(27),YY(5) C C This subroutine may be used to store the quantities useful for the C determination of the take-off parameters in the memory. It is called C at the points of intersection of the ray with interfaces specified in C the input data in the subroutine file 'ray.for' for storing the C computed quantities, see also C C.R.T.5.5.2. C C Input: C ISTOR...The sequential number in the input data of the specified C surface. C YL... Array containing local quantities at the point of the ray. C Y... Array containing basic quantities computed along the ray. C YY... Array containing real auxiliary quantities computed along C the ray. C IY... Array containing integer auxiliary quantities computed C along the ray. C None of the input parameters are altered. C C No output. C C Common block /DCRT/ (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 /RPARD/: INCLUDE 'rpard.inc' C rpard.inc C None of the storage locations of the common block are altered. C C Common block /RPARC/ and /RPARH/: INCLUDE 'rparc.inc' C rparc.inc C Storage locations X1,X2,X1G1,X2G1,X1G2,X2G2 of common block /RPARC/ C may be altered. C Storage locations JPOINT and KMAH of common block /RPARH/ may be C altered. C C Subroutines and external functions required: EXTERNAL METRIC,SMVPRD,XFUN C METRIC... File 'metric.for' of the package 'MODEL'. C SMVPRD... File 'means.for' of the package 'MODEL'. C XFUN... This file. C C Date: 2014, February 13 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Substitution in the common block /RPARC/: REAL X(2),XG(4) EQUIVALENCE (X(1),X1),(XG(1),X1G1) C C Auxiliary storage locations: INTEGER NX,IX REAL GSQRD,G(12),AUX(18),AUX1,AUX2,AUX3 EQUIVALENCE (AUX(11),AUX1),(AUX(12),AUX2),(AUX(13),AUX3) REAL H42,H52,H62,RH1,RH2,RR,XH1,XH2,XR,SIDE1,SIDE2 C C NX... Number of defined Xi functions. C IX... Auxiliary loop variable C GSQRD,G... See subroutine METRIC of the file 'metric.for'. C AUX,AUX1,AUX2,AUX3... Auxiliary storage locations: C AUX(1:18)... Christoffel symbols, C AUX(1:4)... Value and gradient of the function specifying C the reference surface ISRFR, C AUX(5:10)... Second derivatives of the function specifying C the reference surface ISRFR, C AUX(5:7)... Gradient of the Xi function. C AUX(11:13)... AUX1,AUX2,AUX3, C AUX1,AUX2,AUX3... Contravariant derivatives the functions C etc. C H42,H52,H62... Contravariant components of the 2-nd basis vector C of R.C.C.S. C RH1,RH2... Derivatives of the function describing reference the C surface ISRFR in R.C.C.S. C RR... Norm of the gradient of the function describing reference C the surface ISRFR in R.C.C.S. C XH1,XH2... Derivatives of the Xi functions in R.C.C.S. C XR... Projection of the gradient of the Xi functions onto the C normal to the reference surface. C SIDE1,SIDE2... Vectorial side of the shooting domain, in the ray C parameters PAR1,PAR2. C C....................................................................... C IF(KSTOR(ISTOR).EQ.ISRFR) THEN DO 1 IX=1,ISTOR-1 IF(KSTOR(IX).EQ.ISRFR) THEN C The surface is already encountered GO TO 9 END IF 1 CONTINUE JPOINT=JPOINT+1 IF(JPOINT.EQ.IABS(IPOINT)) THEN C Successful ray: C C KMAH index: KMAH=IY(12) C C Metric tenzor CALL METRIC(Y(3),GSQRD,G,AUX) C C Contravariant components of the 2-nd basis vector of R.C.C.S.: CALL SMVPRD(G(7),Y(6),Y(7),Y(8),AUX1,AUX2,AUX3) AUX1=GSQRD*SQRT(AUX1*Y(6)+AUX2*Y(7)+AUX3*Y(8)) H42=(Y(7)*Y(11)-Y(8)*Y(10))/AUX1 H52=(Y(8)*Y( 9)-Y(6)*Y(11))/AUX1 H62=(Y(6)*Y(10)-Y(7)*Y( 9))/AUX1 C C Gradient of the function describing reference surface ISRFR: C Covariant components IF(IABS(ISRFR).LE.100) THEN CALL SRFC2(IABS(ISRFR),Y(3),AUX) ELSE C Reference surface limits the computational volume AUX(2)=0. AUX(3)=0. AUX(4)=0. AUX((IABS(ISRFR)-97)/2)=1. END IF C Contravariant components CALL SMVPRD(G(7),AUX(2),AUX(3),AUX(4),AUX1,AUX2,AUX3) C Gradient in R.C.C.S. RH1=AUX1*Y(9)+AUX2*Y(10)+AUX3*Y(11) RH2=AUX(2)*H42+AUX(3)*H52+AUX(4)*H62 C Projection onto the normal to the reference surface RR =AUX1*AUX(2)+AUX2*AUX(3)+AUX3*AUX(4) C C Number of Xi functions: IF(ISRFX(2).EQ.0) THEN IF(ISRFX(1).EQ.0) THEN NX=0 ELSE NX=1 END IF ELSE NX=2 END IF C C Loop for Xi functions: DO 5 IX=1,NX C C Value and gradient of the Xi functions: C Covariant components CALL XFUN(IX,Y(3),X(IX),AUX(5)) C Contravariant components CALL SMVPRD(G(7),AUX(5),AUX(6),AUX(7),AUX1,AUX2,AUX3) C Gradient in R.C.C.S. XH1=AUX1*Y(9)+AUX2*Y(10)+AUX3*Y(11) XH2=AUX(5)*H42+AUX(6)*H52+AUX(7)*H62 C Projection onto the normal to the reference surface XR =AUX1*AUX(2)+AUX2*AUX(3)+AUX3*AUX(4) C C Derivatives of the Xi function along the reference surface C in R.C.C.S.: AUX1=XH1*RH2-XH2*RH1 AUX2=1.-(RH1*RH1+RH2*RH2)/RR IF(AUX2.LE.0.) THEN C AUX2 is the square of cosine of the incidence angle at the C reference surface. C 648 CALL ERROR * ('648 in RPAR32: Ray tangent to reference surface') C If the source is situated at the reference surface, the C rays tangent to the reference surface should be avoided C by adjusting the ray-parameter domain properly. C The angular deviation of rays from the reference surface C should be at least 0.001 radians. C Data for the ray-parameter domain C If this is not the case, please, contact the author. END IF XH1=(XH1-(XR*RH1+AUX1*RH2)/RR)/AUX2 XH2=(XH2-(XR*RH2-AUX1*RH1)/RR)/AUX2 C C Derivative of the Xi function with respect to the shooting C parameters: SIDE1=PAR1A-PAR1L SIDE2=PAR2A-PAR2L XG(IX)= ((XH1*Y(12)+XH2*Y(13))*(YI(12)*SIDE1+YI(16)*SIDE2) * +(XH1*Y(16)+XH2*Y(17))*(YI(13)*SIDE1+YI(17)*SIDE2) * +(XH1*Y(20)+XH2*Y(21))*(YI(14)*SIDE1+YI(18)*SIDE2) * +(XH1*Y(24)+XH2*Y(25))*(YI(15)*SIDE1+YI(19)*SIDE2)) SIDE1=PAR1B-PAR1L SIDE2=PAR2B-PAR2L XG(IX+2)=((XH1*Y(12)+XH2*Y(13))*(YI(12)*SIDE1+YI(16)*SIDE2) * +(XH1*Y(16)+XH2*Y(17))*(YI(13)*SIDE1+YI(17)*SIDE2) * +(XH1*Y(20)+XH2*Y(21))*(YI(14)*SIDE1+YI(18)*SIDE2) * +(XH1*Y(24)+XH2*Y(25))*(YI(15)*SIDE1+YI(19)*SIDE2)) C 5 CONTINUE C C Undefined functional values: DO 6 IX=NX+1,2 X(IX)=0. XG(IX)=0. XG(IX+2)=0. 6 CONTINUE C END IF END IF 9 CONTINUE RETURN END C C======================================================================= C C C SUBROUTINE RPAR33(YL,Y,YY,IY) REAL YL(6),Y(35),YY(5) INTEGER IY(13) C C This subroutine may be used to store the quantities useful for the C determination of the take-off parameters in the memory. It is called C at the endpoints of the elements of rays, situated at structural C interfaces, see also C C.R.T.5.5.3. C C Input: C YL... Array containing local quantities at the point of the ray. C Y... Array containing basic quantities computed along the ray. C YY... Array containing real auxiliary quantities computed along C the ray. C IY... Array containing integer auxiliary quantities computed C along the ray. For the definition of arrays YL, Y, YY and C IY see the file 'ray.for'. C None of the input parameters are altered. C C No output. C C No subroutines and external functions required. C C Date: 2014, February 13 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C RETURN END C C======================================================================= C C C SUBROUTINE RPAR4(IRAY,PAR1,PAR2,YL,Y,YY,IY,IEND,ISHEET,IREC) C REAL PAR1,PAR2,YL(6),Y(35),YY(5) INTEGER IRAY,IY(13),IEND,ISHEET,IREC C C This subroutine is called after finishing the computation of a ray. C It defines the storage locations YI(22) to YI(25) of the common block C /INITC/. C C Input: C IRAY... Index of the ray, i.e. the number of the already computed C rays, i.e. the output from the last invocation of RPAR2. C PAR1,PAR2... Take-off parameters of the ray. C YL... Array containing local quantities at the point of the ray. C Y... Array containing basic quantities computed along the ray. C YY... Array containing real auxiliary quantities computed along C the ray. C IY... Array containing integer auxiliary quantities computed C along the ray. For the definition of arrays YL, Y, YY and C IY see the file 'ray.for'. C IEND... Reason of the termination of the computation of the ray, C see C C.R.T.5.4. C See subroutine RAY in the subroutine file 'ray.for' for C detailed description of C IEND. C C Output: C ISHEET..Ray-history index. The different ray histories are c consecutively indexed by positive integers 1,2,3,... C according to their appearance during ray tracing. C The ray histories are indexed independently within each C elementary wave. C The ray-history indices are complemented with sign: C Positive - successful ray (crossing reference surface), C negative - unsuccessful ray (terminating before crossing C reference surface). C C IREC... Index of the receiver for a two-point ray, C 0 for a basic ray (and other rays in 2-D), C -1 for other rays. C C Common block /INITC/ (see subroutine file 'init.for'): INCLUDE 'initc.inc' C initc.inc C ------------------------------------------------------------------ C YI(22)..Area of the element of the ray-parameter surface, C corresponding to the ray, C YI(23),YI(24),YI(25)... Components 11, 12, 22 of the symmetric C matrix inverse to the specific moment of the element of C the ray-parameter surface corresponding to the ray, see C C.R.T.6.1. C For one-parametric shooting, locations YI(22) to YI(25) of the common C block are defined in this subroutine. C For two-parametric shooting, locations YI(22) to YI(25) of the common C block are set to zeros in this subroutine. C Other storage locations remain unchanged. C C Common block /RPARD/: INCLUDE 'rpard.inc' C rpard.inc C None of the storage locations of the common block are altered. C C Common block /RPARC/: INCLUDE 'rparc.inc' C rparc.inc C Storage locations of the common block are not altered. C C Subroutines and external functions required: EXTERNAL IHIST,KOOR,METRIC,SMVPRD,RPMEM INTEGER IHIST,KOOR C IHIST...This file. C KOOR,METRIC... File 'metric.for' of the package 'MODEL'. C SMVPRD..File 'means.for' of the package 'MODEL'. C RPMEM...File 'rp3d.for' of the package 'CRT'. C C Date: 2014, February 13 C Coded by Ludek Klimes and Petr Bulant C C----------------------------------------------------------------------- C C Auxiliary storage locations: REAL GSQRD,G(12),AUX(18) REAL X1A REAL RG11,RG12,RG22,SG11,SG12,SG22,SIDE1,SIDE2,AUX1,AUX2,AUX3,AUX4 REAL XG11,XG12,XG22,G1X1,G2X1,G1X2,G2X2 C C GSQRD,G,AUX... For subroutine METRIC. C X1A... Derivative of the X1 function with respect to the C shooting parameter A. C RG11,RG12,RG22... Ray take-off parameter metric tensor. C SG11,SG12,SG22... Shooting parameter metric tensor. C SCALE,SIDE1,SIDE2,AUX1,AUX2,AUX3... Temporary variables. C G1X1,G1X2,G2X1,G2X2... Derivatives of shooting parameters with C respect to the Xi functions. C C....................................................................... C C Ray-history index: positive - successful ray, C otherwise - unsuccessful ray IF(JPOINT.LT.IABS(IPOINT)) THEN C KMAH index of an unsuccessful ray: KMAH=IY(12) END IF ISHEET=IHIST(IY,IEND) C C Calculating shooting-parameter metric tensor, storing the results: IF(ISRFX(2).EQ.0) THEN C Storing the quantities for one-parametric shooting: IF(ANUM.GT.0.) THEN X1A=X1G1/ANUM ELSE X1A=0. END IF CALL RP2DM(JRAY,ISHEET,X1,X1A,IREC) C X1,X1G1... Determined in RPAR32 (if defined), C ISHEET... Determined in this subroutine. C IREC... Determined in RP2DM. ELSE C (a) metric tensor per ray take-off parameters: C Check for non-existing rays: IF(IEND.LT.70) THEN C Cartesian metric tensor (zero for a point source) RG11= YI(12)*YI(12)+YI(13)*YI(13) RG12= YI(12)*YI(16)+YI(13)*YI(17) RG22= YI(16)*YI(16)+YI(17)*YI(17) IF(ABS(YI(12)*YI(17)-YI(13)*YI(14)).LT.0.000001) THEN C Plus unit-sphere metric tensor for point or line source IF(KOOR().EQ.0) THEN AUX1=YI(6)**2+YI(7)**2+YI(8)**2 ELSE CALL METRIC(YI(3),GSQRD,G,AUX) CALL SMVPRD(G(7),YI(6),YI(7),YI(8),AUX1,AUX2,AUX3) AUX1=AUX1*YI(6)+AUX2*YI(7)+AUX3*YI(8) END IF RG11=RG11+(YI(14)*YI(14)+YI(15)*YI(15))/AUX1 RG12=RG12+(YI(14)*YI(18)+YI(15)*YI(19))/AUX1 RG22=RG22+(YI(18)*YI(18)+YI(19)*YI(19))/AUX1 END IF C The above two metric tensors (equal along an exploding unit C sphere source) should be combined in a more sophisticated way. C ************************************************************** C The metric tensor should be determined in 'init.for' ! C ************************************************************** ELSE RG11=1. RG12=0. RG22=1. END IF IF(IEND.GE.70) THEN C The ray has not been traced, GGij have not been initialized. GG11=0. GG12=0. GG22=0. END IF IF(PRM0(4).GE.0.) THEN IF(Y(1)-YI(1).GT.0.) THEN GG11=GG11/(Y(1)-YI(1)) GG12=GG12/(Y(1)-YI(1)) GG22=GG22/(Y(1)-YI(1)) ELSE GG11=0. GG12=0. GG22=0. END IF END IF C (b) Metric tensor per shooting parameters: C SCALE renormalizes the unit-sphere metric in radians to the C metric in some scaled units, proportional to radians by the C factors of scale in the directions of both the shooting C parameters. SIDE1=PAR1A-PAR1L SIDE2=PAR2A-PAR2L AUX3=GG11*SIDE1+GG12*SIDE2 AUX4=GG12*SIDE1+GG22*SIDE2 XG11=SIDE1*AUX3+SIDE2*AUX4 SCALE=ANUM/SQRT(SIDE1*SIDE1+SIDE2*SIDE2) SIDE1=SIDE1*SCALE SIDE2=SIDE2*SCALE AUX1=RG11*SIDE1+RG12*SIDE2 AUX2=RG12*SIDE1+RG22*SIDE2 SG11=SIDE1*AUX1+SIDE2*AUX2 SIDE1=PAR1B-PAR1L SIDE2=PAR2B-PAR2L XG12=SIDE1*AUX3+SIDE2*AUX4 AUX3=GG11*SIDE1+GG12*SIDE2 AUX4=GG12*SIDE1+GG22*SIDE2 XG22=SIDE1*AUX3+SIDE2*AUX4 SCALE=BNUM/SQRT(SIDE1*SIDE1+SIDE2*SIDE2) SIDE1=SIDE1*SCALE SIDE2=SIDE2*SCALE SG12=SIDE1*AUX1+SIDE2*AUX2 AUX1=RG11*SIDE1+RG12*SIDE2 AUX2=RG12*SIDE1+RG22*SIDE2 SG22=SIDE1*AUX1+SIDE2*AUX2 C Another extremely simple metric, conformal with the metric used C at one-parametric shooting, and coinciding with the above scaled C metric for two-parametric shooting in the case of a small C rectangular shooting domain situated at the equator: C SG11=ANUM*ANUM C SG12=0. C SG22=BNUM*BNUM C Regularization near the poles: SG11=SG11+ANUM*ANUM/100. SG22=SG22+BNUM*BNUM/100. C (c) Calculating the derivatives for two-parametric shooting: IF(ISHEET.GT.0) THEN AUX1=X1G1*X2G2-X2G1*X1G2 IF(AUX1.NE.0.) THEN G1X1= X2G2/AUX1 G2X1=-X2G1/AUX1 G1X2=-X1G2/AUX1 G2X2= X1G1/AUX1 ELSE C Caustic: AUX1=X1G1*X1G1+X2G1*X2G1+X1G2*X1G2+X2G2*X2G2 IF(AUX1.NE.0.) THEN C Line caustic - generalized inverse matrix: G1X1=X1G1/AUX1 G2X1=X1G2/AUX1 G1X2=X2G1/AUX1 G2X2=X2G2/AUX1 ELSE C Point caustic: G1X1=0. G2X1=0. G1X2=0. G2X2=0. END IF END IF ELSE X1=0. X2=0. G1X1=0. G2X1=0. G1X2=0. G2X2=0. END IF C (c) Storing the quantities for two-parametric shooting: CALL RPMEM(JRAY,JTYPE,ISHEET,G1,G2,SG11,SG12,SG22, * XG11,XG12,XG22,X1,X2,G1X1,G2X1,G1X2,G2X2) C JRAY,JTYPE,G1,G2... Determined in RPAR2 (output of RP3D), C X1,X2,X1G1,X2G1,X1G2,X2G2... Determined in RPAR32 (if defined), C ISHEET,SG11,SG12,SG22,G1X1,G2X1,G1X2,G2X2... Determined in this C subroutine. C (d) Identifying basic and two-point rays: IF (JTYPE.LT.-1000) THEN C Two-point ray: IREC=-JTYPE-1000 ELSEIF (JTYPE.EQ.0) THEN C Basic ray: IREC=0 ELSE C Auxiliary or boundary ray: IREC=-1 END IF END IF C C Determination of the area YI(22) and inverse specific moment C YI(23:25) of element of the ray-parameter surface, corresponding C to the ray: IF(ISRFX(2).EQ.0) THEN CALL RP2DG(YI(22),YI(23),YI(24),YI(25)) ELSE YI(22)=0. YI(23)=0. YI(24)=0. YI(25)=0. C CALL RP3DG(YI(22),YI(23),YI(24),YI(25)) *** for future extension END IF C C Storing additional quantities related to the shooting algorithm: YI(26)=G1 YI(27)=G2 YI(28)=X1 YI(29)=X2 C RETURN END C C======================================================================= C C C SUBROUTINE XFUN(IFUN,COOR,X,XDER) INTEGER IFUN REAL COOR(3),X,XDER(3) C C Interface subroutine designed to return the value and first C derivatives of the Xi functions describing the profiles. It is called C by the subroutine RPAR32, and may be user-defined. C C Input: C IFUN... 1 for the first (X1 function), and 2 for the second C receiver parameter (X2 function). C COOR... Array containing coordinates X1, X2, X3 of the given C point. C None of the input parameters are altered. C C Output: C X... Value of the corresponding Xi function at the given point. C XDER... Array containing the derivatives of the Xi function with C respect to the general coordinates X1, X2, X3. C C Common block /RPARD/: INCLUDE 'rpard.inc' C rpard.inc C None of the storage locations of the common block are altered. C C Subroutines and external functions required: EXTERNAL SRFC2 C SRFC2... File 'srfc.for' of the package 'MODEL'. C C Date: 1993, December 20 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: REAL FAUX(10) C IF(ISRFX(IFUN).GT.0) THEN CALL SRFC2(ISRFX(IFUN),COOR,FAUX) X=FAUX(1) XDER(1)=FAUX(2) XDER(2)=FAUX(3) XDER(3)=FAUX(4) ELSE X=0. XDER(1)=0. XDER(2)=0. XDER(3)=0. IF(ISRFX(IFUN).LT.0) THEN X=COOR(-ISRFX(IFUN)) XDER(-ISRFX(IFUN))=1. END IF END IF RETURN END C C======================================================================= C C C INTEGER FUNCTION IHIST(IY,IEND) INTEGER IY(13),IEND C C Auxiliary integer function to the subroutine RPAR4 designed to C determine the ray-history index. The ray-history index is used to C distinguish between the rays with different histories for the purposes C of the subroutine RPAR2. C C The ray history is recorded as a sequence of integers. In this C version, the history consists of: C (A) Indices of complex blocks along the ray (including the sign), C (B) Indices of structural interfaces intersected by the ray (including C the sign), and index of the end surface limiting the computational C volume, if intersected. C (C) Reason IEND (see subroutine RAY2 of 'ray.for' and INIT2 of C 'init.for') if the ray is not successful (does not reach the C reference surface). C (D) KMAH index. It may be excluded by input SEP parameteter KMAHI. C Note - differences in versions 4.10 and older: C (B) Index of the end surface not considered in the history. C (C) IEND not considered if IEND.GE.40. C Note - differences in versions 5.00 and older: C (D) KMAH index is not included in the history of an unsuccessful ray. C It saves CPU time but may cause missing two-point rays, especially C close to the boundary of the computational volume. This option C have been discarded in ver. 5.10 because of introducing controlled C initial-value ray tracing followed by interpolation within ray C cells. C Note - differences in versions 5.10-7.00: C (D) KMAH index is always included. C C Input: C IY... Array containing integer auxiliary quantities computed C along the ray. For the definition of array IY see file C 'ray.for'. C IEND... Reason of the termination of the computation of the ray, C see C C.R.T.5.4. C See subroutine RAY in the subroutine file 'ray.for' for C detailed description of C IEND. C None of the input parameters are altered. C C Output: C IHIST...Positive (1,2,3,...) for successful ray, C otherwise (0 or -1,-2,-3,...) for unsuccessful ray. C Successful ray is understood here to be the ray having at C least IABS(IPOINT) points of intersection with the C specified surface ISRFR, see the input data (2). C Different values of IHIST correspond to rays with C different histories. C C Input SEP parameter: C KMAHI='integer'... Optionally excluding the KMAH index from the C ray-history index. C KMAHI=0: The KMAH index is included in the ray-history. C KMAHI=-1: The KMAH index is excluded from the ray-history. C Default: KMAHI=0 C C Common block /RPARD/: INCLUDE 'rpard.inc' C rpard.inc C None of the storage locations of the common block are altered. C C Common block /RPARH/: INCLUDE 'rparc.inc' C rparc.inc C Storage locations of the common block may be altered. C C Subroutines and external functions required: EXTERNAL SRPARH,CODE,RSEP3I C SRPARH... This file. C CODE... File 'code.for' of package 'CRT'. C RSEP3I... File 'sep.for'. C C Date: 2014, May 5 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER I,K C C Optionally excluding the KMAH index from the ray-history index: INTEGER KMAHI SAVE KMAHI DATA KMAHI/-999999/ C IF(KMAHI.EQ.-999999) THEN CALL RSEP3I('KMAHI',KMAHI,0) IF(KMAHI.LT.-1.OR.KMAHI.GT.0) THEN C 657 CALL ERROR('657 in RPAR: Incorrect SEP parameter KMAHI') END IF END IF C IF(IPOINT.EQ.0) THEN IHIST=0 ELSE IF(JPOINT.LT.IABS(IPOINT)) THEN C Unsuccessful ray: IF(PRM0(3).GE.1.5.OR.IRPARH(NRPARH).GE.108) THEN C Removing the index of the surface, at which the ray C terminated, from the ray history: IF(IEND.LT.70) THEN IF(NRPARH-IRPARH(IRPARH(0)).EQ.0) THEN C 651 CALL ERROR('651 in IHIST: Empty ray history') C The history of a ray, which has defined initial C conditions (i.e., IEND.LT.70), should contain at C least two items. C This error should not appear. Contact the authors. END IF IF(MOD(NRPARH-IRPARH(IRPARH(0)),2).EQ.1) THEN C 652 CALL ERROR('652 in IHIST: Odd length of ray history') C At this moment, the history of the ray should contain C an even number of items. C This error should not appear. Contact the authors. END IF NRPARH=NRPARH-1 END IF END IF IF(PRM0(3).LT.2.5) THEN C Including reason IEND of the termination of the computation C of a ray into the ray history: NRPARH=NRPARH+1 CALL SRPARH IRPARH(NRPARH)=IEND ELSE IF(PRM0(3).LT.3.5) THEN C Distinguishing overcritical incidence (IEND=25) from other C reasons of the termination of the computation of a ray: IF(IEND.EQ.25) THEN NRPARH=NRPARH+1 CALL SRPARH IRPARH(NRPARH)=IEND END IF END IF C Following 3 lines include KMAH index also in the history of an C unsuccessful ray: IF(IEND.LT.70) THEN IF(KMAHI.EQ.0) THEN NRPARH=NRPARH+1 CALL SRPARH IRPARH(NRPARH)=KMAH END IF IF(IY(12).NE.KMAH) THEN C 654 CALL ERROR('654 in RPAR2: Incorrect KMAH index') C This error should not appear. Contact the authors. END IF END IF ELSE C Successful ray: IF(KMAHI.EQ.0) THEN NRPARH=NRPARH+1 CALL SRPARH IRPARH(NRPARH)=KMAH END IF IF(IY(12).NE.KMAH) THEN C 655 CALL ERROR('655 in RPAR2: Incorrect KMAH index') C This error should not appear. Contact the authors. END IF END IF C DO 13 K=1,IRPARH(0) IF(IRPARH(K)-IRPARH(K-1).EQ.NRPARH-IRPARH(IRPARH(0))) THEN DO 11 I=IRPARH(K-1)+1,IRPARH(K) IF(IRPARH(I).NE.IRPARH(I-IRPARH(K)+NRPARH)) THEN GO TO 12 END IF 11 CONTINUE C The history of the last ray is equal to the K-th ray history IHIST=INDMIN+K GO TO 20 END IF 12 CONTINUE 13 CONTINUE C New ray history: NRPARH=NRPARH+1 CALL SRPARH DO 18 I=NRPARH,IRPARH(0)+2,-1 IRPARH(I)=IRPARH(I-1) 18 CONTINUE DO 19 I=1,IRPARH(0) IRPARH(I)=IRPARH(I)+1 19 CONTINUE IRPARH(0)=IRPARH(0)+1 IRPARH(IRPARH(0))=NRPARH IHIST=INDMIN+IRPARH(0) CH WRITE(LUWARN(0),'(A/I4,A,I3,19I4,(4X,19I4))') CH * 'New ray history:', CH * IHIST,':',(IRPARH(I),I=IRPARH(IRPARH(0)-1)+1,NRPARH) C 20 CONTINUE IF(JPOINT.LT.IABS(IPOINT)) THEN IHIST=-IHIST END IF IF(IPOINT.LT.0) THEN IF(IEND.EQ.50) THEN CALL CODE(IY,K,I,IEND) IF(IEND.NE.10) THEN IEND=50 END IF END IF IF(IEND.NE.10) THEN IHIST=-IABS(IHIST) END IF END IF END IF RETURN END C C======================================================================= C C C SUBROUTINE SRPARH() C C Auxiliary subroutine to the subroutine RPAR31 designed to shift the C storage locations of the common block /RPARH/. C C No input. C C No output. C C Common block /RPARH/: INCLUDE 'rparc.inc' C rparc.inc C Storage locations of the common block may be altered. C C Subroutines and external functions required: EXTERNAL LUWARN INTEGER LUWARN C LUWARN..File error.for. C C Date: 2013, February 19 C Coded by Ludek Klimes C C----------------------------------------------------------------------- C C Auxiliary storage locations: INTEGER ISHIFT,I C C ISHIFT..Shift in the array IRPARH. C I... Auxiliary loop variable. C C....................................................................... C 1 CONTINUE IF(NRPARH.GT.MRPARH) THEN IF(IRPARH(0).LE.2) THEN C 649 CALL ERROR('649 in SRPARH: Too small array IRPARH in /RPARH/') C Array IRPARH of the common block /RPARH/ designed to store C different ray histories is too small to contain the last three C different ray histories. The dimension MRPARH of array C IRPARH should be adjusted in all declarations of /RPARH/. END IF INDMIN=INDMIN+1 IRPARH(0)=IRPARH(0)-1 ISHIFT=IRPARH(1)-IRPARH(0) NRPARH=NRPARH-ISHIFT DO 2 I=1,IRPARH(0) IRPARH(I)=IRPARH(I+1)-ISHIFT 2 CONTINUE DO 3 I=IRPARH(0)+1,NRPARH-1 IRPARH(I)=IRPARH(I+ISHIFT) 3 CONTINUE C 650 WRITE(LUWARN(0),'(A/I4,A,I3,19I4,(4X,19I4))') * 'Warning 650 in SRPARH: Discarded ray history:', * INDMIN,':',(IRPARH(I),I=IRPARH(0)+1,IRPARH(1)) C Array IRPARH of the common block /RPARH/ designed to store C different ray histories is too small for all ray histories. C Some ray histories have thus been discarded. Dimension MRPARH C of array IRPARH may be adjusted to prevent problems due to C discarded ray histories. GO TO 1 END IF RETURN END C C======================================================================= C