C Subroutine file 'rpar.for' to control the take-off parameters of rays.
C
C By Vlastislav Cerveny, Ludek Klimes, Ivan Psencik
C
C This file consists of the following subroutines:
C     RPARB...Block data subroutine defining common block /RPARD/ to
C             store the input data for the take-off parameters of rays,
C             common block /RPARC/ to collect the quantities needed for
C             boundary-value ray tracing, common block /RPARH/
C             describing the histories of rays.  The histories of rays
C             are used by function IHIST to determine the ray-history
C             index. Different values of ray-history index are attached
C             to rays with different histories.
C     RPAR1...Subroutine designed to read the input data for the
C             take-off parameters, 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     RPAR2...Subroutine designed to specify the take-off parameters of
C             individual rays.
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     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     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     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     XFUN... Subroutine returning the value and first derivatives of
C             the X-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     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     SRPARH..Auxiliary subroutine to the subroutine RPAR31 to shift the
C             storage locations of the common block /RPARH/.
C
C Note:
C     The lines denoted by '*' in the first column contain an undebugged
C     code intended for future extensions.
C
C Input data 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 (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 crossections of the surface ISRFR and the surfaces
C             defined by the X1-function.  The receivers are the
C             crossections of the surface ISRFR and the surfaces defined
C             by the X1 and X2 functions.  ISRFR must be one of the
C             indices listed in the input data RAY-(7) (the indices for
C             storing computed quantities), including the sign, see the
C             subroutine file 'ray.for'.  For ISRFR=0 there are no
C             profiles for boundary-value ray tracing, ISRFR not listed
C             in the input data RAY-(7) has 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             For ISRFX1=-1, -2, or -3,  the X1-function coincides with
C               the -ISRFX1-th model coordinate.
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 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'.  ISRFX1 must not exceed the
C               number of the defined surfaces more than by 1.
C             Example: The most often 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 in the
C               subroutine 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     ISRFX2..Index of the function that becomes to be the X2-function.
C             ISRFX2=0: Initial-value ray traing or one-parametric
C               shooting.
C             Otherwise: Meaning of ISRFX2 is analogous to ISRFX1.
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 does not
C             exceed XERR, the ray is considered to be a boundary-value
C             ray.  The distance of the ray from the given profile is
C             measured by means of the difference of the X1-function.
C     AERR... Specifies the accuracy of the determination of the
C             boundary ray between the bundles of rays with different
C             histories: the step between two consecutive rays with
C             different histories is partitioned in order not to exceed
C             the basic step in the first direction, multiplyed by AERR.
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 alloved 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 alloved 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(I),I=3,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 positive and ISRFX2.EQ.0 (3.1),
C     For NREC positive and ISRFX2.NE.0 (3.2), for NREC negative (3.3):
C (3.1) 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.1) 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 (3.3) 'FREC'
C     'FREC'..Name of the file with profile or receiver coordinates.
C     Default: 'FREC'='rec.dat'.
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 (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     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     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 consideredred to be successful.  All
C               other rays are rejected (IMPORTANT OPTION).
C             IPOINT=large:  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             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,IPOINT:
C             Zeroes when starting the complete ray tracing program,
C             otherwise the previous values.
C     PRM1(I)=0, I=1,MPRM1.
C (5.2) /(a slash)
C
C Input file 'FREC' 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
C Storage in the memory:
C     The input data (2) are stored in common block /RPARD/ defined
C     together with common blocks /RPARC/ and /RPARH/ in the following
C     subroutine:
C     ------------------------------------------------------------------
      BLOCK DATA RPARB
      INCLUDE 'rpard.inc'
      INCLUDE 'rparc.inc'
      END
C     ------------------------------------------------------------------
C
C Date: 1996, September 30
C Coded by Ludek Klimes
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 None of the storage locations of the common block are altered.
C
C Common block /RPARD/:
      INCLUDE '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 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 RPARB,NSRFC
C     RPARB.. Block data subroutine of this file.
C     NSRFC... File 'model.for' of the package 'MODEL'.
C
C Error messages:
C     650...  NREC negative:
C             Negative input value of NREC, data line (2), is not
C             allowed.
C     651...  Too small array XREC in /RPARD/:
C             The number of given profiles (specified by the values of
C             X-functions, see input data 3) is exceeding the dimension
C             mrec of the array XREC in the common block /RPARD/.  Thus,
C             mrec should be adjusted in all declarations of /RPARD/.
C     652...  Wrong value of ISRFX1:
C             The function no. ISRFX1 specified by the routines of the
C             file 'srfc.for' is defined neither in the input data set
C             'model.dat' nor 'dcrt.dat' and cannot be defined in this
C             data set.  The index ISRFX1 is too large or is less than
C             -3.
C
C Date: 1996, September 30
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
        PRM0(2)=999999.
        READ(LUN,*) ISRFR,ISRFX(1),ISRFX(2),NREC,XERR,AERR,PRM0
c       IF(NREC.LT.0) THEN
c         PAUSE 'Error 650 in RPAR1: NREC negative'
c         STOP
c       END IF
C       (3):
        IF(NREC.GT.0) THEN
          IF(NREC.GT.MREC) THEN
            PAUSE 'Error 651 in RPAR1: Too small array XREC in /RPARD/'
            STOP
          END IF
          IF(ISRFX(2).EQ.0) THEN
C           Initial-value or 1-paremetric shooting:
            READ(LUN,*) (XREC(1,I),I=1,NREC)
          ELSE
C           2-paremetric 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
              PAUSE'Error 651 in RPAR1: Too small array XREC in /RPARD/'
              STOP
            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       (4) Auxiliary functions no. ISRFX(1) and ISRFX(2):
        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
          PAUSE 'Error 652 in RPAR1: Wrong value of ISRFX1'
          STOP
        END IF
        IF(ISRFX(2).EQ.I) THEN
          CALL SRFC1(LUN,1)
        ELSE IF(ISRFX(2).LT.-3.OR.ISRFX(2).GT.I) THEN
          PAUSE 'Error 653 in RPAR1: Wrong value of ISRFX2'
          STOP
        END IF
C
C       Default values
        PAR1L=0.
        PAR2L=0.
        PAR1A=0.
        PAR2A=0.
        PAR1B=0.
        PAR2B=0.
        ANUM=0.
        BNUM=0.
        IPOINT=0
      END IF
C
      LURPAR=LUN
      RETURN
      END
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 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 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 Error messages:
C     654...  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.
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
          PAUSE 'Error 654 in RPAR2: Zero value of ISRFX1'
          STOP
        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
      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 None of the storage locations of the common block are altered.
C
C Common block /RPARH/:
      INCLUDE '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: 1996, September 29
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
C4.10     IF(IY(6).NE.0.AND.IABS(IY(6)).LE.NSRFC()) THEN
          IF(IY(6).NE.0) THEN
C           Structural interface
            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
      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 iput 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 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 None of the storage locations of the common block are altered.
C
C Common block /RPARD/:
      INCLUDE '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 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: 1996, Septenber 30
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 X-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 X-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 X-functions in R.C.C.S.
C     XR...   Projection of the gradient of the X-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         Succesful 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 X-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 X-functions:
          DO 5 IX=1,NX
C
C           Value and gradient of the X-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 X-function along the reference surface
C           in R.C.C.S.:
            AUX1=XH1*RH2-XH2*RH1
            AUX2=1.-(RH1*RH1+RH2*RH2)/RR
            XH1=(XH1-(XR*RH1+AUX1*RH2)/RR)/AUX2
            XH2=(XH2-(XR*RH2-AUX1*RH1)/RR)/AUX2
C
C           Derivative of the X-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
      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
      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     IREC... Index of the receiver for a two-point ray,
C             otherwise zero.
C
C Common block /INITC/ (see subroutine file 'init.for'):
      INCLUDE '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 None of the storage locations of the common block are altered.
C
C Common block /RPARC/:
      INCLUDE '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: 1996, June 17
C Coded by Ludek Klimes
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 X-functions.
C
C.......................................................................
C
C     Ray-history index: positive - successful ray,
C     otherwise - unsuccessful ray
      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 directons 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) Storing the quantities 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       Determination of two-point rays (to be moved to RPMEM):
        IREC=0
        IF(JTYPE.LT.-1000) THEN
C         Two-point ray is being searched for
          IF(ISHEET.GT.0) THEN
            IREC=-JTYPE-1000
            IF(IREC.LT.1.OR.NREC.LT.IREC) THEN
              PAUSE 'Error *** in RPAR4: Wrong index of the receiver'
              STOP
            END IF
            DIST2=(X1-XREC(1,IREC))**2+(X2-XREC(2,IREC))**2
            IF(DIST2.GT.XERR**2) THEN
              JTYPE=-2
              IREC=0
            END IF
          ELSE
            JTYPE=-2
          END IF
        END IF
C       ****************************************************************
        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.
      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
      RETURN
      END
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 X-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, and 2 for the second receiver parameter
C             (X-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 X-function at the given point.
C     XDER... Array containing the derivatives of the X-function with
C             respect to the general coordinates X1, X2, X3.
C
C Common block /RPARD/:
      INCLUDE '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
      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 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
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             Succesfull 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 None of the storage locations of the common block are altered.
C
C Common block /RPARH/:
      INCLUDE '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: 1995, September 29
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:
C         KMAH index is not included in the history of an unsuccessful
C         ray.  It saves CPU time but may cause missing two-point rays,
C         especially close to the boundary of the computational volume.
          NRPARH=NRPARH+1
          CALL SRPARH
          IRPARH(NRPARH)=IEND
        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
      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 Storage locations of the common block may be altered.
C
C No subroutines and external functions required.
C
C Error messages:
C     655...  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
C             three different ray histories.  The dimension MRPARH of
C             the array IRPARH should be adjusted in all declarations of
C             /RPARH/.
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
          PAUSE 'Error 655 in SRPARH: Too small array IRPARH in /RPARH/'
          STOP
        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