C 
C Subroutine file 'rpar.for' to control the take-off parameters of rays.
C
C Date: 1997, October 26
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     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.R.T.5.5.1).  It may be used
C             to keep the quantities useful for the determination of the
C             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.R.T.5.5.2).  It may be
C             used to keep the quantities useful for the determination
C             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).  It may be used to keep the quantities
C             useful for the determination of the take-off parameters in
C             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     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,NREC,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             Description of 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.
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 NREC.EQ.0.,
C               or IPOINT=0 at all input data lines (5.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               srfc.for
C               ISRFX1 must not exceed the number of the defined
C               surfaces more than by 1.
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 is analogous to ISRFX1:
C               For ISRFX2=-1, -2, or -3,  the X1 function coincides
C                 with the -ISRFX1-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                 srfc.for
C                 ISRFX1 must not exceed the number of the defined
C                 surfaces (including ISRFX1) more than by 1.
C     NREC... Number of the defined receivers or profiles along the
C               reference surface ISRFR for the two-point ray tracing.
C             Zero if there are no profiles.
C             Negative (e.g. -1) if the profile coordinates are to be
C               read from a separate file, see (3.2).  This is the
C               recommended option.
C             Positive values are left for the compatibility with the
C               older versions, however, they are not recommended.
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 per one step of numerical
C             integration refer to input data DCRT-(3) in 'ray.for'.
C             Description of data DCRT-(3)
C     AERR... Specifies the accuracy of the determination of the
C             boundary rays between the different ray histories.
C             The step between two consecutive rays of different
C             histories (boundary rays) is bisectioned until it does not
C             exceed the basic step in the corresponding direction
C             multiplyed by AERR.  The basic step in the direction of
C             the first normalized ray parameter (ranging from 0 to 1),
C             recalculated to the parallel great circle at the unit
C             sphere, is 1/ANUM.  The basic step in the direction of
C             the second normalized ray parameter, recalculated to the
C             parallel great circle at the unit sphere, is 1/BNUM, see
C             data (5) below.
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 3-D two-point ray tracing.
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 3-D two-point ray tracing.
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.  The boundaries
C               of the computational volume and end surfaces are not
C               taken into account.  This option limits the number of
C               thin spaces, uncovered by ray tubes, situated between
C               ray tubes corresponding to different ray histories.
C               This option is thus suitable for controlled
C               initial-value ray tracing followed be 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               (5.1) below in order to disable successful rays.
C               PRM0(3).NE.0 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.  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 be the travel-time interpolation within ray
C               tubes.
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 be 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(I),I=4,MPRM0... Reserved for future extension of the
C             numerical parameters for two-point ray tracing.
C Default values:
C     ISRFR=0, ISRFX1=0, ISRFX2=0, NREC=0, 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 (5.1).
C             The boundary-value rays are not determined if ISRFR.EQ.0
C               or NREC.EQ.0.
C               Neither are they if IPOINT=0, see input data (5.1).
C             Thus, just the initial-value ray tracing is performed if
C               (ISRFR.EQ.0 or NREC.LE.0) and AERR.GT.1.
C               Just the initial-value ray tracing is also performed if
C               IPOINT=0, see input data (5.1).
C               During the initial-value ray tracing, only the basic net
C               of rays with the given constant steps is generated.
C (3) Specification of the profiles for the one-parametric
C     boundary-value ray tracing.
C     This input is not performed if NREC=0,
C     for NREC negative (3.1),
C     for NREC positive and ISRFX2.EQ.0 (3.2),
C     for NREC positive and ISRFX2.NE.0 (3.3).
C (3.1) 'REC'
C     'REC'...Name of the file with profile or receiver coordinates.
C             Description of file REC
C     Default: 'REC'='rec.dat'.
C (3.2) XREC1(1),...,XREC1(NREC)
C     XREC1(I)... The I-th profile is the intersection of the surface
C             ISRFR with the isosurface of the X1 function,
C             corresponding to the value XREC1(I).  The values XREC1(I)
C             may be specified in any order.
C (3.3) XREC1(1),XREC2(1),...,XREC1(NREC),XREC2(NREC)
C     XREC1(I),XREC2(I)... The I-th receiver is the intersection of the
C             surface ISRFR with the isosurface of the X1 function,
C             corresponding to the value XREC1(I), and the isosurface of
C             the X2 function, corresponding to the value XREC2(I).
C (4) 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 set
C     'model.dat' or 'dcrt.dat'.  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.dat' + the number
C     ISRFCA of the auxiliary surfaces defined in the data set
C     'dcrt.dat' + 1, and so on.
C     Data for '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 (5) For each elementary wave IWAVE the following data (5.1) and (5.2):
C (5.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             in the input data INIT-(2) read by and described in
C             'init.for'.
C             Description of data INIT-(2)
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     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 basic net
C               of INT(ANUM)*INT(ANUM) rays is given by the take-off ray
C               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 basic net of rays corresponds to the initial-value
C               ray tracing.
C               Besides the basic net of rays, also additional rays are
C               computed in order to find the boundary rays between the
C               bundles of rays with different histories, and to find
C               the boundary-value rays.  All additional rays are given
C               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 basic net
C               of INT(ANUM) 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             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             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 be 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 (5.2) /(a slash)
C Example of data RPAR initial-value r.t.
C Example of data RPAR 3-D two-point r.t.
C
C                                                     
C Input file 'REC' containing profile (or receiver) coordinates:
C     This file is used only if NREC.LT.0, see the input data (2) above.
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).
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:
      INTEGER NSRFC
      EXTERNAL NSRFC
C     NSRFC... File 'model.for' of the package 'MODEL'.
C
C Date: 1997, July 19
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.1) THEN
C       (1) Reading the name of the input data:
        READ(LUN,*) TEXTP
C       (2):
        ISRFR=0
        ISRFX(1)=0
        ISRFX(2)=0
        NREC=0
        XERR=0.000
        AERR=999999.
        DO 10 I=1,MPRM0
          PRM0(I)=0.
   10   CONTINUE
        READ(LUN,*) ISRFR,ISRFX(1),ISRFX(2),NREC,XERR,AERR,PRM0
C
C       (3):
        IF(NREC.GT.0) THEN
          IF(NREC.GT.MREC) THEN
C           641
            PAUSE 'Error 641 in RPAR1: Too small array XREC in /RPARD/'
            STOP
C           The number NREC (see input data 2) of given profiles or
C           receivers is exceeding the dimension MREC of the array XREC
C           in the common block /RPARD/.  MREC should thus be adjusted
C           accordingly.
C           /RPARD/
C           Description of input data 2
          END IF
          IF(ISRFX(2).EQ.0) THEN
C           Initial-value or 1-parametric shooting:
            READ(LUN,*) (XREC(1,I),I=1,NREC)
          ELSE
C           2-parametric shooting:
            READ(LUN,*) (XREC(1,I),XREC(2,I),I=1,NREC)
          END IF
        ELSE IF(NREC.LT.0) THEN
          TEXTP='rec.dat'
          READ(LUN,*) TEXTP
          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,*) TEXTP,(AUX(J),J=1,5)
            IF(TEXTP.EQ.'$') THEN
              GO TO 12
            END IF
            IF(I.GT.MREC) THEN
C             642
              PAUSE
     *            'Error 642 in RPAR1: Too small array XREC in /RPARD/'
              STOP
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
C       (4) Auxiliary functions no. ISRFX(1) and ISRFX(2):
        IF(ISRFX(1).EQ.0.AND.ISRFX(2).NE.0) THEN
C         640
          PAUSE 'Error 640 in RPAR1: Zero ISRFX1, nonzero ISRFX2'
          STOP
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
          PAUSE 'Error 643 in RPAR1: Wrong value of ISRFX1'
          STOP
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.dat' nor 'dcrt.dat' and cannot be defined in this data
C         set.  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
          PAUSE 'Error 644 in RPAR1: Wrong value of ISRFX2'
          STOP
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.dat' nor 'dcrt.dat' and cannot be defined in this data
C         set.  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
      END IF
C
      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 RP2D
C     RP2D... One-parametric boundary-value ray tracing.  This file.
C
C Date: 1996, June 15
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
            RETURN
          END IF
        END IF
        IF(ISRFX(1).EQ.0.AND.ISRFR.NE.0.AND.NREC.NE.0.AND.IPOINT.NE.0)
     *                                                              THEN
C         645
          PAUSE 'Error 645 in RPAR2: Zero value of ISRFX1'
          STOP
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
      GO TO 10
C
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RPAR31(YL,Y,YY,IY)
      REAL YL(6),Y(35),YY(5)
      INTEGER IY(12)
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: 1997, October 15
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(12)
      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.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: 1997, July 6
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
              PAUSE
     *           'Error 648 in RPAR32: Ray tangent to reference surface'
              STOP
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(12)
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.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: 1990, January 22
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(12),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.R.T.5.4).  For a detailed description see
C             subroutine RAY (subroutine file 'ray.for').
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 Locations YI(22) to YI(25) of the common block are defined in this
C subroutine, 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
      INTEGER  IHIST
C     IHIST...This file.
C
C Date: 1997, October 17
C Coded by Ludek Klimes and Petr Bulant
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      REAL X1A
      REAL RG11,RG12,RG22,SG11,SG12,SG22,SIDE1,SIDE2,AUX1,AUX2
      REAL G1X1,G2X1,G1X2,G2X2
C
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     SIDE1,SIDE2,AUX1,AUX2... 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(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
            AUX1=YI(6)**2+YI(7)**2+YI(8)**2
            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
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
        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
        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,
     *                                        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(IEND)
      INTEGER 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.
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
C Input:
C     IEND... Reason of the termination of the computation of the ray
C             (see C.R.T.5.4).  For a detailed description see
C             subroutine RAY (subroutine file 'ray.for').
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 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
C     SRPARH... This file.
C
C Date: 1997, October 26
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I,K
C
      IF(IPOINT.EQ.0) THEN
        IHIST=0
      ELSE
        IF(JPOINT.LT.IABS(IPOINT)) THEN
C         Unsuccessful ray:
          IF(PRM0(3).GE.1.5) THEN
C           Removing the index of the surface, at which the ray
C           terminated, from the ray history:
            NRPARH=NRPARH-1
          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:
          NRPARH=NRPARH+1
          CALL SRPARH
          IRPARH(NRPARH)=KMAH
C
        ELSE
C         Successful ray:
          NRPARH=NRPARH+1
          CALL SRPARH
          IRPARH(NRPARH)=KMAH
        END IF
        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)
C
   20   CONTINUE
        IF(JPOINT.LT.IABS(IPOINT)) THEN
          IHIST=-IHIST
        END IF
        IF(IPOINT.LT.0) THEN
          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 No subroutines and external functions required.
C
C Date: 1991, May 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
      IF(NRPARH.GT.MRPARH) THEN
        IF(IRPARH(0).LE.2) THEN
C         647
          PAUSE 'Error 647 in SRPARH: Too small array IRPARH in /RPARH/'
          STOP
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 1 I=1,IRPARH(0)
          IRPARH(I)=IRPARH(I+1)-ISHIFT
    1   CONTINUE
        DO 2 I=IRPARH(0)+1,NRPARH-1
          IRPARH(I)=IRPARH(I+ISHIFT)
    2   CONTINUE
      END IF
      RETURN
      END
C
C=======================================================================
C