C 
C Subroutine file 'rp2d.for' to control 1-parametric shooting.
C
C Date: 2000, January 24
C Coded by Ludek Klimes
C
C.......................................................................
C
C This file consists of the following subroutines:
C     RP2D... Auxiliary subroutine to RPAR2, designed to control
C             one-parametric boundary-value ray tracing (as well as the
C             initial-value ray tracing).
C             RP2D
C     RP2DM...Auxiliary subroutine to RPAR4 designed to store the
C             quantities dependent on the ray history, required for
C             one-parametric boundary-value ray tracing.
C             RP2DM
C     RP2DG...Auxiliary subroutine to RPAR4 designed to determine the
C             area GG and inverse specific moment G11,G12,G22 of the
C             element of the ray-parameter surface, corresponding
C             to the ray at one-parametric boundary-value ray tracing.
C             RP2DG
C
C.......................................................................
C
C Storage in the memory:
C     Common block /RP2DC/ to store the quantities needed for
C     one-parametric boundary-value ray tracing is defined in the
C     include file 'rp2d.inc'.
C     rp2d.inc
C
C=======================================================================
C
C     
C
      SUBROUTINE RP2D(IRAY0,JRAY,G1NEW)
      INTEGER IRAY0,JRAY
      REAL G1NEW
C
C This subroutine controls one-parametric boundary-value ray tracing.
C Auxiliary subroutine to RPAR2.
C
C Input:
C     IRAY0...Difference IRAY0=IRAY-JRAY between index of a ray within
C             the elementary wave, and the index within one shooting
C             domain.
C     JRAY... Number of the already calculated rays within the shooting
C             domain.  JRAY=0 when starting to cover a new shooting
C             domain (e.g. when beginning the computation of a new
C             elementary wave.  Otherwise, the output from the previous
C             invocation of this subroutine.
C
C Output:
C     JRAY... JRAY=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     G1NEW...Shooting parameter of the ray to be calculated.
C             The domain of the shooting parameter is specified by the
C             two first executive statements of this subroutine.
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 /RP2DC/:
      INCLUDE 'rp2d.inc'
C     rp2d.inc
C All storage locations of the common block may be altered.
C
C Subroutines and external functions required:
      EXTERNAL LUWARN,WRITA
      INTEGER  LUWARN
C     WRITA... File 'writ.for'.
C     LUWARN.. File 'crt.for'.
C
C Date: 2000, January 24
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I,IAUX
      REAL G1MIN,G1MAX,G1STEP,XAUX
      CHARACTER*80 TEXT
C
C     G1MIN,G1MAX... Boundaries of the shooting domain (domain of the
C             ray parameter selected to control the shooting).
C     G1STEP... Basic step in the shooting parameter.
C
C     Specification of shooting domain and basic step:
      G1MIN=0.
      G1MAX=ANUM
      G1STEP=1.
C
C.......................................................................
C
      IF(JRAY.EQ.0) THEN
        NEW=1
        A  (NEW)=G1MIN-G1STEP
        IND(NEW)=0
        INDOLD =-1
      END IF
C
C     Control centre of the one-parametric boundary-value ray tracing
   10 CONTINUE
      IF(A(1).GT.G1MIN+G1STEP/2.) THEN
C       At least 2 rays of the elementary wave are computed - tests may
C       be performed:
        IF(IND(NEW).NE.INDOLD) THEN
C         Rays with different histories - looking for the boundary ray:
          IF(NEW.LT.NEWMAX.AND.ABS(A(NEW)-AOLD).GT.AERR) THEN
            A(NEW+1)=(A(NEW)+AOLD)/2.
            IRE(NEW+1)=0
            ITER=0
            GO TO 80
          END IF
        ELSE IF(INDOLD.GT.0) THEN
          IF(ISRFX(1).NE.0) THEN
C           Two successful rays - looking for the profiles in the
C           interval:
            XAUX=X(NEW)
            XLEFT=XOLD
            IF(IREOLD.NE.0) THEN
              IF((XREC(1,IREOLD)-X(NEW))*(XREC(1,IREOLD)-XOLD).LE.0.)
     *                                                              THEN
                XLEFT=XREC(1,IREOLD)
              END IF
            END IF
            IAUX=0
C           Loop for the given profiles:
            DO 20 I=1,NREC
              IF((XREC(1,I)-X(NEW))*(XREC(1,I)-XLEFT).LE.0.) THEN
C               There is a profile in the interval: boundary-value ray
C               tracing
C5.10           IF(ABS(XREC(1,I)-XOLD).LE.XERR) THEN
                IF(IREOLD.EQ.I) THEN
C                 Boundary-value ray successfully found - endpoint XOLD
                  ITER=0
                ELSE IF(NEW.GE.NEWMAX) THEN
C                 661
                  WRITE(TEXT,'(A,I4,A)')
     *              'Warning 661 in RPAR2: Boundary-value ray to',
     *              I,'th receiver probably inaccurate'
                  CALL WARN(TEXT)
C                 The dimension NEWMAX of the arrays in the common block
C                 /RP2DC/ is too small to allow for sufficiently fine
C                 division of the basic step in the ray parameters.
C                 Boundary-value ray is thus found with greater error
C                 than XERR.
                  ITER=0
                ELSE
                  IF((XREC(1,I)-XAUX)*(XREC(1,I)-XLEFT).LE.0.) THEN
C                   There is a boundary-value ray to be found
                    XAUX=XREC(1,I)
                    IAUX=I
                  END IF
                END IF
              END IF
   20       CONTINUE
C5.00       IF(XAUX.NE.X(NEW)) THEN
            IF(IAUX.NE.0) THEN
C             There is a boundary-value ray to be found
C5.00         IF(ABS(XAUX-X(NEW)).LE.XERR) THEN
              IF(IRE(NEW).EQ.IAUX) THEN
C               Boundary-value ray successfully found
                ITER=0
              ELSE
C               Shooting a new ray within the interval
                IRE(NEW+1)=0
                IF(MOD(ITER,2).EQ.0) THEN
C                 Regula falsi:
                  A(NEW+1)=((XAUX-XOLD)*A(NEW)+(X(NEW)-XAUX)*AOLD)
     *                     /(X(NEW)-XOLD)
                ELSE
C                 Newton-Raphson (paraxial approximation):
                  IF(ABS(XAUX-XOLD).LE.ABS(XAUX-X(NEW))) THEN
                    A(NEW+1)=AOLD+(XAUX-XOLD)/XAOLD
                  ELSE
                    A(NEW+1)=A(NEW)+(XAUX-X(NEW))/XA(NEW)
                  END IF
                END IF
                IF(A(NEW+1).LE.AOLD.OR.A(NEW).LE.A(NEW+1)) THEN
                  A(NEW+1)=(A(NEW)+AOLD)/2.
                END IF
                IF(A(NEW+1).LE.AOLD.OR.A(NEW).LE.A(NEW+1)) THEN
C                 662
                  WRITE(LUWARN(0),'(A,I4,A)')
     *              'Warning 662 in RPAR2: Boundary-value ray to',
     *              IAUX,'th receiver inaccurate'
C                 Rounding error does not allow for sufficiently fine
C                 division of the basic step in the ray parameters.
C                 Numerically, there is no ray between the two rays
C                 surrounding the profile.
C                 It is recommended to decrease the allowed error UEB of
C                 the computation of the ray (see line (3) of the input
C                 data set dcrt for the complete ray tracing described
C                 in the file 'ray.for'), or to increase the error XERR
C                 (line (2) of the input data for this file) of the
C                 boundary-value ray.
                  IF(ABS(XAUX-XOLD).LE.ABS(XAUX-X(NEW))) THEN
                    A(NEW+1)=AOLD
                  ELSE
                    A(NEW+1)=A(NEW)
                  END IF
                  IRE(NEW+1)=IAUX
                END IF
                ITER=ITER+1
                GO TO 80
              END IF
            END IF
          END IF
        END IF
        IF(IND(NEW).EQ.INDOLD) THEN
C         Leaving a homogeneous interval:
          CALL WRITA(IRAY0+IRAOLD,IRAY0+IRA(NEW),0)
        END IF
      END IF
C     Tests are performed - changing to the next interval:
C     Note: now, DAOLD is an irremovable discrepancy in the element of
C     the take-off parameter A, corresponding to the old ray.  We are
C     going to forget it.
      INDOLD=IND(NEW)
      IREOLD=IRE(NEW)
      IRAOLD=IRA(NEW)
      AOLD=A(NEW)
      DAOLD=DA(NEW)
      XOLD=X(NEW)
      XAOLD=XA(NEW)
      NEW=NEW-1
      IF(NEW.EQ.0) THEN
C       New basic interval:
        A(1)=AOLD+G1STEP
        IRE(1)=0
        ITER=0
        IF(A(1).GT.G1MAX+0.001*G1STEP) THEN
C         End of one-parametric shooting
          JRAY=0
          RETURN
        END IF
        GO TO 80
      END IF
      GO TO 10
C
C     New ray
   80 CONTINUE
      JRAY=JRAY+1
      NEW=NEW+1
      IRA(NEW)=JRAY
      G1NEW=A(NEW)
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RP2DM(IRAY,ISHEET,X1,X1A,IREC)
      INTEGER IRAY,ISHEET,IREC
      REAL X1,X1A
C
C This subroutine stores the quantities dependent on the ray history and
C required for one-parametric boundary-value ray tracing.
C Designed to be called by subroutine RPAR4.
C
C Input:
C     IRAY... Index of the ray within a single shooting domain.
C     ISHEET..Ray-history index: positive - successful ray,
C             otherwise - unsuccessful ray.
C     X1...   Value of the first X-function at the point of intersection
C             of the ray with the reference surface ISRFR.
C             Meaningful only for successful rays, i.e. ray crossing the
C             reference surface ISRFR at least IABS(IPOINT) times.
C     X1A...  Derivative of the first X-function with respect to the
C             shooting parameter A.
C             Meaningful only for successful rays.
C
C Output:
C     IREC... Index of the receiver for a two-point ray,
C             otherwise zero.
C
C Common block /RPARD/ (defined in file 'rpar.for'):
      INCLUDE 'rpard.inc'
C     rpard.inc
C None of the storage locations of the common block are altered.
C
C Common block /RP2DC/:
      INCLUDE 'rp2d.inc'
C     rp2d.inc
C Storage locations IND,X,XA of the common block may be altered.
C
C No subroutines and external functions required.
C
C Date: 1997, November 19
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I
      REAL XAUX
C
C.......................................................................
C
      IF(IRE(NEW).NE.0) THEN
        IREC=IRE(NEW)
      ELSE
C       Determination of two-point rays:
        IREC=0
        IF(ISRFX(1).NE.0) THEN
          IF(ISHEET.GT.0) THEN
            XAUX=XERR
C           Loop for the given profiles:
            DO 20 I=1,NREC
              IF((X1-XOLD)*(XREC(1,I)-XOLD).GE.0.) THEN
                IF(ABS(XREC(1,I)-X1).LE.XAUX) THEN
C                 Ray endpoint is close to the profile
                  IF(IRAY.GT.1) THEN
                    IF(INDOLD.EQ.ISHEET) THEN
                      IF(IREOLD.EQ.I) THEN
                        GO TO 10
                      END IF
                    END IF
                  END IF
                  IF(NEW.GT.1) THEN
                    IF(IND(NEW-1).EQ.ISHEET) THEN
                      IF(IRE(NEW-1).EQ.I) THEN
                        GO TO 10
                      END IF
                    END IF
                    IF((X1-X(NEW-1))*(XREC(1,I)-X(NEW-1)).LT.0.) THEN
                      GO TO 10
                    END IF
                  END IF
                  IREC=I
                  XAUX=ABS(XREC(1,I)-X1)
                END IF
              END IF
   10         CONTINUE
   20       CONTINUE
          END IF
        END IF
      END IF
C
C     Storing:
      IND(NEW)=ISHEET
      IRE(NEW)=IREC
      X(NEW)=X1
      XA(NEW)=X1A
C
      RETURN
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE RP2DG(GG,GG11,GG12,GG22)
      REAL GG,GG11,GG12,GG22
C
C This subroutine determines the area GG and inverse specific moment
C G11,G12,G22 of the element of the ray-parameter surface, corresponding
C to the ray at one-parametric boundary-value ray tracing.
C Auxiliary subroutine to RPAR4.
C
C No input.
C
C Output:
C     GG...   Area of the element of the ray-parameter surface,
C             corresponding to the ray.
C     GG11,GG12,GG22... Components of the symmetric matrix inverse to
C             the specific moment of the element of the ray-parameter
C             surface corresponding to the ray, see C.R.T.6.1.
C
C Common block /RPARD/ (defined in file 'rpar.for'):
      INCLUDE 'rpard.inc'
C     rpard.inc
C None of the storage locations of the common block are altered.
C
C Common block /RP2DC/:
      INCLUDE 'rp2d.inc'
C     rp2d.inc
C Storage locations DAOLD,DA of the common block may be altered.
C
C No subroutines and external functions required.
C
C Date: 1994, January 3
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      REAL A1,A2,B1,B2,AB
C
C     A1,A2...Step along the ray-parameter surface in the direction A.
C     B1,B2...Step along the ray-parameter surface in the direction B.
C     AB...   Area of the basic element of the ray-parameter surface for
C             a two-parametric system of rays, or the square of its
C             width for a one-parametric system of rays.  The basic
C             element of the ray-parameter surface is given by the basic
C             steps in the directions A and B.
C
C.......................................................................
C
C     Determination of the element of the parameter a corresponding to
C     the ray:
      IF(NEW.EQ.1) THEN
        DA(NEW)=1.
      ELSE
        DA(NEW)=(A(NEW-1)-AOLD)/2.
        DA(NEW-1)=DA(NEW-1)-(A(NEW)-AOLD)/2.
        DAOLD=DAOLD-(A(NEW-1)-A(NEW))/2.
        IF(IND(NEW).EQ.IND(NEW-1)) THEN
          DA(NEW)=DA(NEW)+DA(NEW-1)
          DA(NEW-1)=0.
        END IF
      END IF
      IF(IND(NEW).EQ.INDOLD) THEN
        DA(NEW)=DA(NEW)+DAOLD
        DAOLD=0.
      END IF
C
C     Determination of the area GG and inverse specific moment G11, G12,
C     G22 of the element of the ray-parameter surface, corresponding to
C     the ray:
      IF(ANUM.NE.0.) THEN
        A1=(PAR1A-PAR1L)/ANUM
        A2=(PAR2A-PAR2L)/ANUM
        IF(BNUM.NE.0.) THEN
C         Two-parametric system of rays
          B1=(PAR1B-PAR1L)/BNUM
          B2=(PAR2B-PAR2L)/BNUM
          AB=A1*B2-A2*B1
          GG  =DA(NEW)*ABS(AB)
          GG11= 12.*(A2*A2+B2*B2)/(AB*AB)
          GG12=-12.*(A1*A2+B1*B2)/(AB*AB)
          GG22= 12.*(A1*A1+B1*B1)/(AB*AB)
        ELSE
C         One-parametric system of rays
          AB=A1*A1+A2*A2
          GG  =DA(NEW)*SQRT(AB)
          GG11=12.*(A1*A1)/(AB*AB)
          GG12=12.*(A1*A2)/(AB*AB)
          GG22=12.*(A2*A2)/(AB*AB)
        END IF
      ELSE
        IF(BNUM.NE.0.) THEN
C         One-parametric system of rays
          B1=(PAR1B-PAR1L)/BNUM
          B2=(PAR2B-PAR2L)/BNUM
          AB=B1*B1+B2*B2
          GG  =SQRT(AB)
          GG11=12.*(B1*B1)/(AB*AB)
          GG12=12.*(B1*B2)/(AB*AB)
          GG22=12.*(B2*B2)/(AB*AB)
        ELSE
C         Single ray
          GG  =1.
          GG11=0.
          GG12=0.
          GG22=0.
        END IF
      END IF
      DA(NEW)=0.
      RETURN
      END
C
C=======================================================================
C