C Subroutine file 'rp2d.for' to control 1-parametric shooting.
C
C By Vlastislav Cerveny, Ludek Klimes, Ivan Psencik
C
C This file consists of the following subroutines:
C     RP2DB...Block data subroutine defining common block /RP2DC/
C             to store the quantities needed for one-parametric
C             boundary-value ray tracing.
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     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     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
C Storage in the memory:
C     Common block /RP2DC/ is declared in the following subroutine:
C     ------------------------------------------------------------------
      BLOCK DATA RP2DB
      INCLUDE 'rp2d.inc'
      END
C     ------------------------------------------------------------------
C
C Date: 1997, July 8
C Coded by Ludek Klimes
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/ (defined in file 'rpar.for'):
      INCLUDE 'rpard.inc'
C None of the storage locations of the common block are altered.
C
C Common block /RP2DC/:
      INCLUDE 'rp2d.inc'
C All storage locations of the common block may be altered.
C
C Subroutines and external functions required:
      EXTERNAL WRITA
C     WRITA... File 'writ.for'.
C
C Warning messages:
C     656...  Boundary-value ray cannot be accurately found:
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     657...  Boundary-value ray cannot be accurately determined:
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 the
C             computation of the ray (see line (3) of the input data set
C             dcrt for the complete ray tracing described in the file
C             'ray.for'), or to increase the error XERR (line (2) of the
C             input data for this file) of the boundary-value ray.
C
C Date: 1996, September 29
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I
      REAL G1MIN,G1MAX,G1STEP,XAUX
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.
            ITER=0
            GO TO 80
          END IF
        ELSE IF(INDOLD.GT.0) THEN
          IF(ISRFX(1).NE.0) THEN
C           Two succesful rays - looking for the profiles in the
C           interval:
            XAUX=X(NEW)
C           Loop for the given profiles:
            DO 20 I=1,NREC
              IF((XREC(1,I)-X(NEW))*(XREC(1,I)-XOLD).LE.0.) THEN
C               There is a profile in the interval: boundary-value ray
C               tracing
                IF(ABS(XREC(1,I)-XOLD).LE.XERR) THEN
C                 Boundary-value ray succesfully found - endpoint XOLD
                  ITER=0
                ELSE IF(NEW.GE.NEWMAX) THEN
C                 Boundary-value ray inaccurately found
                  WRITE(*,'('' Warning 656 in RPAR2: Boundary-value'',
     *                      '' ray cannot be accurately found.'')')
                  ITER=0
                ELSE
                  IF((XREC(1,I)-XAUX)*(XREC(1,I)-XOLD).LE.0.) THEN
C                   There is a boundary-value ray to be found
                    XAUX=XREC(1,I)
                  END IF
                END IF
              END IF
   20       CONTINUE
            IF(XAUX.NE.X(NEW)) THEN
C             There is a boundary-value ray to be found
              IF(ABS(XAUX-X(NEW)).LE.XERR) THEN
C               Boundary-value ray succesfully found - endpoint X(NEW)
                ITER=0
              ELSE
C               Shooting a new ray within the interval
                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                 Boundary-value ray inaccurately found
                  WRITE(*,'('' Warning 657 in RPAR2: Boundary-value'',
     *                      '' ray cannot be accurately determined.'')')
                  ITER=0
                  GO TO 70
                END IF
                ITER=ITER+1
                GO TO 80
              END IF
            END IF
          END IF
        END IF
      END IF
   70 CONTINUE
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.
      IF(IND(NEW).EQ.INDOLD) THEN
C       Leaving a homogeneous interval:
        CALL WRITA(IRAY0+IRAOLD,IRAY0+IRA(NEW),0)
      END IF
      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
        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
      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 succesful 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 succesful 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 None of the storage locations of the common block are altered.
C
C Common block /RP2DC/:
      INCLUDE '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: 1996, September 29
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I
      REAL XAUX
C
C.......................................................................
C
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(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
              END IF
              IREC=I
              XAUX=ABS(XREC(1,I)-X1)
            END IF
   10       CONTINUE
   20     CONTINUE
        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
      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 None of the storage locations of the common block are altered.
C
C Common block /RP2DC/:
      INCLUDE '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
