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, AND COMMON BLOCK /RP2DC/
C             TO STORE THE QUANTITIES NEEDED FOR ONE-PARAMETRIC
C             BOUNDARY-VALUE RAY TRACING.  THE HISTORIES OF RAYS ARE
C             USED BY THE 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-FUNCTION 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     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 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
C     ISRFR...INDEX OF THE REFERENCE SURFACE AT WHICH THE COMPUTED
C             QUANTITIES ARE STORED AND THE GIVEN PROFILES FOR
C             BOUNDARY-VALUE RAY TRACING ARE SITUATED.  THE PROFILES ARE
C             THE CROSSECTIONS OF THE SURFACE ISRFR AND THE SURFACES
C             DEFINED BY THE X-FUNCTION.  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             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 X-FUNCTION.
C             THE X-FUNCTION IS DEFINED HERE BY THE SUBROUTINE XFUN.
C             FOR ISRFX1=-1, -2, OR -3,  THE X-FUNCTION COINCIDES WITH
C               THE -ISRFX1-TH MODEL COORDINATE.
C             FOR ISRFX1=0,  THE X-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.LE.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 THE 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..ZERO.
C     NREC... NUMBER OF THE DEFINED PROFILES ALONG THE SURFACE ISRFR FOR
C               THE ONE-PARAMETRIC BOUNDARY-VALUE 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.
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 X-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 (OR NEGATIVE) HAS THE
C             SAME EFFECT AS AERR.GT.1.
C DEFAULT VALUES:
C     ISRFR=0, ISRFX1=0, ISRFX2=0, NREC=0, XERR=0.000, AERR=999999.
C             THE BOUNDARY RAYS ARE NOT DETERMINED IF AERR.GT.1.
C               NEITHER ARE THEY IF IPOINT.LE.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.LE.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.LE.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 (3.1), FOR NREC NEGATIVE (3.2):
C (3.1) XREC(1),...,XREC(NREC)
C     XREC(I)... THE I-TH PROFILE IS THE INTERSECTION OF THE SURFACE
C             ISRFR WITH THE ISOSURFACE OF THE X-FUNCTION, CORRESPONDING
C             TO THE VALUE XREC(I).  THE VALUES XREC(I) MAY BE SPECIFIED
C             IN ANY ORDER.
C (3.2) 'FREC'
C     'FREC'..NAME OF THE FILE WITH PROFILE COORDINATES XREC(1),...,
C             XREC(NREC).
C     DEFAULT: 'FREC'='REC.DAT'.
C (4) THE DATA SPECIFYING THE FUNCTION NO. ISRFX1.  THIS AUXILIARY
C     FUNCTION IS SPECIFIED BY THE SUBROUTINE SRFC2 AND THE DATA ARE
C     READ BY SUBROUTINE SRFC1.  FOR THE DESCRIPTION OF THE DATA REFER
C     TO THE SUBROUTINE SRFC1 OF THE FILE 'SRFC.FOR'.  THIS INPUT IS
C     PERFORMED ONLY IF ISRFX IS POSITIVE AND THE FUNCTION ISRFX IS NOT
C     DEFINED IN THE DATA SET 'MODEL.DAT' OR 'DCRT.DAT'.  IN SUCH A
C     CASE, ISRFX MUST EQUAL TO THE NUMBER NSRFC OF THE SURFACES DEFINED
C     IN THE DATA SET 'MODEL.DAT' + THE NUMBER ISRFCA OF THE AUXILIARY
C     SURFACES DEFINED IN THE DATA SET 'DCRT.DAT' + 1.
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             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=PARL1+A*(PAR1A-PAR1L)/ANUM+B*(PAR1B-PAR1L)/BNUM,
C                 PAR2=PARL2+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=PARL1+A*(PAR1A-PAR1L)/ANUM,
C                 PAR2=PARL2+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     IPOINT..THE IPOINT-TH POINT OF INTERSECTION OF A RAY WITH THE
C             SURFACE ISRFR IS CONSIDERED FOR THE BOUNDARY-VALUE RAY
C             TRACING.  IN MOST CASES, USER WILL WISH THE VALUE
C             IPOINT=1.  FOR IPOINT.LE.0, NO RAY HISTORY IS RECORDED,
C             ALL RAYS HAVE THE SAME INDEX 0 AND NO BOUNDARY-VALUE RAY
C             RAY TRACING IS PERFORMED.  FOR THE STRUCTURAL INTERFACE
C             (NOT AUXILIARY SURFACE), THE POSITIVE OR NEGATIVE SIDE OF
C             THE SURFACE ISRFR IS 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.LE.0.
C     DEFAULT VALUES: ZEROES WHEN STARTING THE COMPLETE RAY TRACING
C     PROGRAM, OTHERWISE THE PREVIOUS VALUES.
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 SIMLEST WAY IS TO SUBMIT JUST THE /.
C (2) SEVERAL TIMES (2.1):
C (2.1) 'NAMER(IR)',XREC(IR),X2R(IR),X3R(IR),/
C     'NAMER(IR)'... STRING RESERVED FOR THE NAME OF THE RECEIVER.
C             NO MEANING HERE, ANYTHING IN APOSTROPHES MAY BE SUBMITTED.
C     XREC(IR)... X-FUNCTION (COORDINATE) OF THE IR-TH PROFILE
C             (RECEIVER).
C     X2R(IR),X3R(IR)... NOT READ IN THIS VERSION.  RESERVED FOR OTHER
C             COORDINATES (X-FUNCTIONS) OF THE IR-TH RECEIVER.
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
      INTEGER MREC
      PARAMETER (MREC=128)
      INTEGER ISRFR,IPOINT,ISRFX(2),NREC
      REAL XERR,AERR,XREC(2,MREC)
      REAL PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
      COMMON/RPARD/ISRFR,IPOINT,ISRFX,NREC,AERR,XERR,XREC,
     *             PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
      SAVE /RPARD/
      INTEGER LURPAR,JRAY,JTYPE
      REAL G1,G2,X1,X2,X1G1,X2G1,X1G2,X2G2
      COMMON/RPARC/LURPAR,JRAY,JTYPE,G1,G2,X1,X2,X1G1,X2G1,X1G2,X2G2
      SAVE /RPARC/
      INTEGER MRPARH
      PARAMETER (MRPARH=1024)
      INTEGER JPOINT,INDMIN,NRPARH,IRPARH(0:MRPARH)
      COMMON/RPARH/JPOINT,INDMIN,NRPARH,IRPARH
      SAVE /RPARH/
      INTEGER NEWMAX
      PARAMETER (NEWMAX=20)
      INTEGER ITER,NEW,INDOLD,IND(NEWMAX)
      REAL AOLD,A(NEWMAX),DAOLD,DA(NEWMAX)
      REAL XOLD,X(NEWMAX),XAOLD,XA(NEWMAX)
      COMMON/RP2DC/ITER,NEW,INDOLD,IND,AOLD,A,DAOLD,DA,XOLD,X,XAOLD,XA
      SAVE /RP2DC/
      END
C     ------------------------------------------------------------------
C     ISRFR,IPOINT,ISRFX,NREC,XERR,AERR,XREC... INPUT DATA (2) AND (3).
C     PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM... LAST VALUES OF
C             THE INPUT DATA (5.1).
C     LURPAR..LOGICAL UNIT NUMBER OF THE EXTERNAL INPUT DEVICE
C             CONTAINING THE INPUT DATA.
C     JRAY... REFERENCE INDEX OF THE RAY WITHIN ONE SHOOTING DOMAIN
C             (NOT WITHIN ONE ELEMENTARY WAVE).  TO BE USED IN FUTURE
C             VERSIONS.
C     JTYPE...RESERVED FOR FUTURE USE.
C     G1,G2...NORMALIZED SHOOTING PARAMETERS.  HERE, NORMALIZED MEANS
C             THAT THE SHOOTING DOMAIN IS (0,1)*(0,1) IN THE NORMALIZED
C             SHOOTING PARAMETERS.
C             FOR ONE-PARAMETRIC SHOOTING:
C               G1,G2: MULTIPLIED BY ANUM AND BNUM, RESPECTIVELY.
C               G1: THE FIRST, SHOOTING PARAMETER.
C               G2: THE SECOND, FIXED PARAMETER.
C     X1,X2 ..VALUES OF THE X-FUNCTIONS AT THE POINT OF INTERSECTION
C             WITH THE REFERENCE SURFACE (E.G., TWO SELECTED COORDINATES
C             OF THE POINT OF INTERSECTION).
C     X1G1,X2G1,X1G2,X2G2... DERIVATIVES OF THE X-FUNCTIONS WITH RESPECT
C             TO THE NORMALIZED SHOOTING PARAMETERS.
C     JPOINT..THE NUMBER OF CROSSECTIONS OF THE COMPUTED RAY WITH THE
C             SURFACE ISRFR.  BECAUSE, FOR THE STRUCTURAL INTERFACES,
C             THE POSITIVE OR NEGATIVE SIDE OF THE SURFACE IS SPECIFIED,
C             EITHER ZERO OR TWO CROSSECTIONS CORRESPOND TO THE
C             REFLECTION FROM THE SURFACE ISRFR.
C     INDMIN..ABSOLUTE VALUE OF THE RAY-HISTORY INDEX MINUS THE
C             LOCATION OF THE CORRESPONDING RAY HISTORY IN THE ARRAY
C             IRPARH OF THE COMMON BLOCK /RPARH/.
C     NRPARH..INDEX IN THE ARRAY IRPARH OF THE LAST ELEMENT OF THE
C             HISTORY OF THE COMPUTED RAY, SEE BELOW.
C     IRPARH..ARRAY CONTAINING DIFFERENT HISTORIES OF RAYS.  UNDER THE
C             RAY HISTORY WE UNDERSTAND HERE THE BLOCK-SURFACE CODE (SEE
C             THE FILE 'CODE.FOR') DESCRIBING THE BEHAVIOUR OF
C             INDIVIDUAL RAYS, COMPLEMENTED BY THE INDICATION IEND
C             SPECIFYING THE REASON OF THE TERMINATION OF THE
C             COMPUTATION OF THE RAY (SEE THE SUBROUTINE RAY2 OF THE
C             FILE 'RAY.FOR').
C             IRPARH(0)... NUMBER OF DIFFERENT RAY HISTORIES ENCOUNTERED
C               YET.
C             IRPARH(I),I=1,...,IRPARH(0)... INDICES IN THE ARRAY IRPARH
C               OF THE LAST ELEMENTS OF INDIVIDUAL RAY HISTORIES.
C             IRPARH(I),I=IRPARH(0)+1,...,IRPARH(IRPARH(0))...
C               INDIVIDUAL ENCOUNTERED RAY HISTORIES.
C             IRPARH(I),I=IRPARH(IRPARH(0)),...,NRPARH... HISTORY OF THE
C               COMPUTED RAY.  AFTER FINISHING COMPUTATION OF THE RAY,
C               IT IS COMPARED WITH THE RECORDED HISTORIES.  IF IT
C               APPEARS TO BE A NEW HISTORY, IT IS RECORDED, OTHERWISE
C               IT IS FORGOTTEN.
C             IF THE ARRAY IRPARH IS TOO SMALL TO KEEP ALL DIFFERENT RAY
C             HISTORIES CORRESPONDING TO THE COMPUTED ELEMENTARY WAVE,
C             THE OLDEST HISTORIES ARE FORGOTTEN.
C     ITER... NUMBER OF ITERATIONS WHEN SEARCHING FOR THE LAST
C             BOUNDARY-VALUE RAY.
C     NEW...  NUMBER OF COMPUTED RAYS IN THE CURRENT CUT (I.E.
C             CORRESPONDING TO THE LAST VALUE OF B) TO THE RIGHT (I.E.
C             FOR GREATER VALUES OF A) FROM THE INTERVAL UNDER STUDY.
C             THE INTERVAL UNDER STUDY IS A = AOLD TO A(NEW).  NEW=1
C             WHEN BEGINNING TO STUDY A NEW BASIC INTERVAL.  NEW=NEW+1
C             IF THE INTERVAL UNDER STUDY IS DIVIDED BY A NEW RAY.
C             AOLD=A(NEW) AND THEN NEW=NEW-1 IF WE CHANGE TO THE NEXT
C             INTERVAL.
C     INDOLD..INDEX OF THE LEFT-END RAY OF THE INTERVAL UNDER STUDY.
C     IND...  INDICES OF THE RAYS TO THE RIGHT FROM THE INTERVAL UNDER
C             STUDY.
C     AOLD... VALUE OF A CORRESPONDING TO THE LEFT-END RAY OF THE
C             INTERVAL UNDER STUDY.
C     A...    VALUES OF A CORRESPONDING TO THE COMPUTED RAYS SITUATED TO
C             THE RIGHT OF THE INTERVAL UNDER STUDY.
C     DAOLD...ELEMENT OF A CORRESPONDING TO THE RAY A=AOLD.
C     DA...   DA(I) IS THE ELEMENT OF A CORRESPONDING TO THE RAY A=A(I).
C     XOLD... VALUE OF X-FUNCTION CORRESPONDING TO THE RAY A=AOLD.
C             DEFINED ONLY FOR SUCCESFUL RAYS, I.E. RAY CROSSING THE
C             SURFACE ISRFR AT LEAST IPOINT TIMES.
C     X...    X(I) IS THE VALUE OF X-FUNCTION CORRESPONDING TO THE RAY
C             A=A(I).  DEFINED ONLY FOR SUCCESFUL RAYS.
C     XAOLD...DERIVATIVE OF X-FUNCTION WITH RESPECT TO A CORRESPONDING
C             TO THE RAY A=AOLD.  DEFINED ONLY FOR SUCCESFUL RAYS.
C     XA...   XA(I) IS THE DERIVATIVE OF X-FUNCTION WITH RESPECT TO A
C             CORRESPONDING TO THE RAY A=A(I).  DEFINED ONLY FOR
C             SUCCESFUL RAYS.
C     COMMON BLOCK /RPARD/ IS INCLUDED IN EXTERNAL PROCEDURES RPAR1,
C     RPAR2, RPAR32, RPAR4, XFUN, AND RP2DG.
C     COMMON BLOCK /RPARC/ IS INCLUDED IN EXTERNAL PROCEDURES RPAR1,
C     RPAR2, RPAR32, AND RPAR4.
C     COMMON BLOCK /RPARH/ IS INCLUDED IN EXTERNAL PROCEDURES RPAR2,
C     RPAR31, RPAR32, IHIST, AND SRPARH.
C     COMMON BLOCK /RP2DC/ IS INCLUDED IN EXTERNAL PROCEDURES RP2D,
C     RP2DM, AND RP2DG.
C
C DATE: 1994, JANUARY 15
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 SUBROTINE 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/:
      INTEGER MEND,MSTOR
      PARAMETER (MEND=128)
      PARAMETER (MSTOR=128)
      INTEGER KSTORE,NEXPS,NHLF,MODCRT
      REAL STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT,BOUNDR(7)
      INTEGER NSRFCA,NEND,KEND(MEND),NSTOR,KSTOR(MSTOR)
      COMMON/DCRT/ KSTORE,NEXPS,NHLF,MODCRT,STORE,STEP,UEB,UEBPP,UEBPH,
     *             UEBHH,UEBDRT,BOUNDR,NSRFCA,NEND,KEND,NSTOR,KSTOR
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C COMMON BLOCK /RPARD/:
      INTEGER MREC
      PARAMETER (MREC=128)
      INTEGER ISRFR,IPOINT,ISRFX(2),NREC
      REAL XERR,AERR,XREC(2,MREC)
      REAL PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
      COMMON/RPARD/ISRFR,IPOINT,ISRFX,NREC,AERR,XERR,XREC,
     *             PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
C ALL THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE DEFINED IN THIS
C SUBROUTINE.
C
C COMMON BLOCK /RPARC/:
      INTEGER LURPAR,JRAY,JTYPE
      REAL G1,G2,X1,X2,X1G1,X2G1,X1G2,X2G2
      COMMON/RPARC/LURPAR,JRAY,JTYPE,G1,G2,X1,X2,X1G1,X2G1,X1G2,X2G2
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-FUNCTION, 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     653...  WRONG VALUE OF ISRFX2:
C             ISRFX2 IN THE INPUT DATA, LINE (2), HAS TO EQUAL ZERO.
C
C DATE: 1994, JANUARY 15
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER LUREC
      PARAMETER (LUREC=9)
      CHARACTER*80 TEXTP
      INTEGER I
      REAL AUX1,AUX2
C
C     TEXTP...THE NAME OF THE DATA. STRING OF 80 CHARACTERS.
C     I...    AUXILIARY LOOP VARIABLE, OR TEMPORARY VARIABLE.
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.
        READ(LUN,*) ISRFR,ISRFX(1),ISRFX(2),NREC,XERR,AERR
        IF(NREC.LT.0) THEN
          PAUSE 'ERROR 650 IN RPAR1: NREC NEGATIVE'
          STOP
        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
          READ(LUN,*) (XREC(1,I),I=1,NREC)
        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 12 I=1,MREC+1
            AUX1=-999999.
            AUX2=0.
            READ(LUREC,*) TEXTP,AUX1,AUX2
            IF(AUX1.EQ.-999999.) 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
            XREC(1,I)=AUX1
            XREC(2,I)=AUX2
   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
*       GO TO 19
        IF(ISRFX(2).NE.0) THEN
          PAUSE 'ERROR 653 IN RPAR1: WRONG VALUE OF ISRFX2'
          STOP
        END IF
*  19   CONTINUE
*       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/:
      INTEGER MREC
      PARAMETER (MREC=128)
      INTEGER ISRFR,IPOINT,ISRFX(2),NREC
      REAL XERR,AERR,XREC(2,MREC)
      REAL PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
      COMMON/RPARD/ISRFR,IPOINT,ISRFX,NREC,AERR,XERR,XREC,
     *             PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
C STORAGE LOCATIONS PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM OF THE
C COMMON BLOCK MAY BE REDEFINED IN THIS SUBROUTINE.
C
C COMMON BLOCK /RPARC/:
      INTEGER LURPAR,JRAY,JTYPE
      REAL G1,G2,X1,X2,X1G1,X2G1,X1G2,X2G2
      COMMON/RPARC/LURPAR,JRAY,JTYPE,G1,G2,X1,X2,X1G1,X2G1,X1G2,X2G2
C STORAGE LOCATIONS JRAY,JTYPE,G1,G2 OF THE COMMON BLOCK ARE ALTERED.
C
C COMMON BLOCK /RPARH/:
      INTEGER MRPARH
      PARAMETER (MRPARH=1024)
      INTEGER JPOINT,INDMIN,NRPARH,IRPARH(0:MRPARH)
      COMMON/RPARH/JPOINT,INDMIN,NRPARH,IRPARH
C STORAGE LOCATIONS NRPARH AND IRPARH(0) OF THE COMMON BLOCK 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.GT.0 IN THE INPUT DATA.
C
C DATE: 1994, JANUARY 3
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER IRAY0,I
      REAL AUX1,AUX2,AUX3,AUX4,AUX5,AUX6,AUX7,AUX8
      SAVE IRAY0
C
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
C.......................................................................
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:
        CALL RP2D(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
      I=IPOINT
      READ(LURPAR,*)PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM,IPOINT
      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.I.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.GT.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 /RPARH/:
      INTEGER MRPARH
      PARAMETER (MRPARH=1024)
      INTEGER JPOINT,INDMIN,NRPARH,IRPARH(0:MRPARH)
      COMMON/RPARH/JPOINT,INDMIN,NRPARH,IRPARH
C ALL STORAGE LOCATIONS OF THE COMMON BLOCK MAY BE ALTERED.
C
C COMMON BLOCK /RPARD/:
      INTEGER MREC
      PARAMETER (MREC=128)
      INTEGER ISRFR,IPOINT,ISRFX(2),NREC
      REAL XERR,AERR,XREC(2,MREC)
      REAL PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
      COMMON/RPARD/ISRFR,IPOINT,ISRFX,NREC,AERR,XERR,XREC,
     *             PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL NSRFC
      INTEGER NSRFC
C     NSRFC... FILE 'MODEL.FOR' OF THE PACKAGE 'MODEL'.
C
C DATE: 1993, DECEMBER 20
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     NO AUXILIARY STORAGE LOCATIONS.
C
      IF(NRPARH.EQ.IRPARH(IRPARH(0))) THEN
C       INITIAL ELEMENT OF THE RAY
        JPOINT=0
        NRPARH=NRPARH+1
        CALL SRPARH
        IRPARH(NRPARH)=IY(5)
      END IF
      IF(JPOINT.LT.IPOINT) THEN
        IF(IY(6).NE.0.AND.IABS(IY(6)).LE.NSRFC()) 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
      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'):
      INTEGER MEND,MSTOR
      PARAMETER (MEND=128)
      PARAMETER (MSTOR=128)
      INTEGER KSTORE,NEXPS,NHLF,MODCRT
      REAL STORE,STEP,UEB,UEBPP,UEBPH,UEBHH,UEBDRT,BOUNDR(7)
      INTEGER NSRFCA,NEND,KEND(MEND),NSTOR,KSTOR(MSTOR)
      COMMON/DCRT/ KSTORE,NEXPS,NHLF,MODCRT,STORE,STEP,UEB,UEBPP,UEBPH,
     *             UEBHH,UEBDRT,BOUNDR,NSRFCA,NEND,KEND,NSTOR,KSTOR
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C COMMON BLOCK /INITC/ (SUBROUTINE FILE 'INIT.FOR'):
      INTEGER MSRFCA
      PARAMETER (MSRFCA=128)
      INTEGER ISB1I,ICB1I
      REAL YLI(6),YI(25),FSRFCA(MSRFCA)
      COMMON/INITC/ISB1I,ICB1I,YLI,YI,FSRFCA
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C COMMON BLOCK /RPARD/:
      INTEGER MREC
      PARAMETER (MREC=128)
      INTEGER ISRFR,IPOINT,ISRFX(2),NREC
      REAL XERR,AERR,XREC(2,MREC)
      REAL PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
      COMMON/RPARD/ISRFR,IPOINT,ISRFX,NREC,AERR,XERR,XREC,
     *             PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C COMMON BLOCK /RPARC/:
      INTEGER LURPAR,JRAY,JTYPE
      REAL G1,G2,X1,X2,X1G1,X2G1,X1G2,X2G2
      COMMON/RPARC/LURPAR,JRAY,JTYPE,G1,G2,X1,X2,X1G1,X2G1,X1G2,X2G2
C STORAGE LOCATIONS X1,X2,X1G1,X2G1,X1G2,X2G2 OF THE COMMON BLOCK MAY BE
C ALTERED.
C
C COMMON BLOCK /RPARH/:
      INTEGER MRPARH
      PARAMETER (MRPARH=1024)
      INTEGER JPOINT,INDMIN,NRPARH,IRPARH(0:MRPARH)
      COMMON/RPARH/JPOINT,INDMIN,NRPARH,IRPARH
C JUST THE STORAGE LOCATION JPOINT OF THE COMMON BLOCK IS 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: 1994, JANUARY 3
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 IX
      REAL GSQRD,G(12),AUX(18),AUX1,AUX2,AUX3
      REAL H42,H52,H62,RH1,RH2,RR,XH1,XH2,XR,SIDE1,SIDE2
      EQUIVALENCE (AUX(11),AUX1),(AUX(12),AUX2),(AUX(13),AUX3)
C
C     I...    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 FUNCTION X.
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 FUNCTION X IN R.C.C.S.
C     XR...   PROJECTION OF THE GRADIENT OF THE X-FUNCTION 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.IPOINT) THEN
C         SUCCESFUL RAY:
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
          CALL SRFC2(ISRFR,Y(3),AUX)
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
              IX=0
            ELSE
              IX=1
            END IF
*         ELSE
*           IX=2
*         END IF
C
C         LOOP FOR X-FUNCTIONS:
          DO 5 IX=1,IX
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=IX,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)
C
      REAL PAR1,PAR2,YL(6),Y(35),YY(5)
      INTEGER IRAY,IY(12),IEND,ISHEET
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 COMMON BLOCK /INITC/ (SEE SUBROUTINE FILE 'INIT.FOR'):
      INTEGER MSRFCA
      PARAMETER (MSRFCA=128)
      INTEGER ISB1I,ICB1I
      REAL YLI(6),YI(25),FSRFCA(MSRFCA)
      COMMON/INITC/ISB1I,ICB1I,YLI,YI,FSRFCA
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/:
      INTEGER MREC
      PARAMETER (MREC=128)
      INTEGER ISRFR,IPOINT,ISRFX(2),NREC
      REAL XERR,AERR,XREC(2,MREC)
      REAL PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
      COMMON/RPARD/ISRFR,IPOINT,ISRFX,NREC,AERR,XERR,XREC,
     *             PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C COMMON BLOCK /RPARC/:
      INTEGER LURPAR,JRAY,JTYPE
      REAL G1,G2,X1,X2,X1G1,X2G1,X1G2,X2G2
      COMMON/RPARC/LURPAR,JRAY,JTYPE,G1,G2,X1,X2,X1G1,X2G1,X1G2,X2G2
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: 1994, JANUARY 8
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 FIRST X-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(ISHEET,X1,X1A)
C       X1,X1G1... DETERMINED IN RPAR32 (IF DEFINED),
C       ISHEET... DETERMINED IN THIS SUBROUTINE.
*     ELSE
C       (A) UNIT-SPHERE METRIC TENSOR PER RAY TAKE-OFF PARAMETERS:
*       IF(YI(12).EQ.0..AND.YI(13).EQ.0..AND.
*    *     YI(16).EQ.0..AND.YI(17).EQ.0.) THEN
C         POINT SOURCE
*         AUX1=YI(6)**2+YI(7)**2+YI(8)**2
*         RG11=(YI(14)*YI(14)+YI(15)*YI(15))/AUX1
*         RG12=(YI(14)*YI(18)+YI(15)*YI(19))/AUX1
*         RG22=(YI(18)*YI(18)+YI(19)*YI(19))/AUX1
*       ELSE
*         PAUSE 'ERROR *** IN RPAR4: ONLY POINT SOURCE CONSIDERED NOW'
*       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=ANUM/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
*         G1X1= X2G2/AUX1
*         G2X1=-X2G1/AUX1
*         G1X2=-X1G2/AUX1
*         G2X2= X1G1/AUX1
*       ELSE
*         X1=0.
*         X2=0.
*         G1X1=0.
*         G2X1=0.
*         G1X2=0.
*         G2X2=0.
*       END IF
*       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
*       if(iray.eq.1) open(40,file='rp.out')
*       write(40,'(3i6,2f10.6,3f10.3,2f10.5,4f10.6)')
*    *  jray,jtype,isheet,g1,g2,sg11,sg12,sg22,x1,x2,g1x1,g2x1,g1x2,g2x2
*     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/:
      INTEGER MREC
      PARAMETER (MREC=128)
      INTEGER ISRFR,IPOINT,ISRFX(2),NREC
      REAL XERR,AERR,XREC(2,MREC)
      REAL PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
      COMMON/RPARD/ISRFR,IPOINT,ISRFX,NREC,AERR,XERR,XREC,
     *             PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
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 THE RAYS WITH DIFFERENT HISTORIES FOR THE PURPOSES OF THE
C SUBROUTINE RPAR2.
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 IPOINT POINTS OF INTERSECTION WITH THE SPECIFIED
C             SURFACE ISRFR, SEE THE INPUT DATA (2).  DIFFERENT VALUES
C             OF IHIST CORRESPOND TO RAYS WITH DIFFERENT HISTORIES.
C
C COMMON BLOCK /RPARD/:
      INTEGER MREC
      PARAMETER (MREC=128)
      INTEGER ISRFR,IPOINT,ISRFX(2),NREC
      REAL XERR,AERR,XREC(2,MREC)
      REAL PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
      COMMON/RPARD/ISRFR,IPOINT,ISRFX,NREC,AERR,XERR,XREC,
     *             PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C COMMON BLOCK /RPARH/:
      INTEGER MRPARH
      PARAMETER (MRPARH=1024)
      INTEGER JPOINT,INDMIN,NRPARH,IRPARH(0:MRPARH)
      COMMON/RPARH/JPOINT,INDMIN,NRPARH,IRPARH
C ALL 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: 1994, JANUARY 8
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER I,K
C
      IF(JPOINT.LT.IPOINT.AND.IEND.LT.40) THEN
        NRPARH=NRPARH+1
        CALL SRPARH
        IRPARH(NRPARH)=IEND
      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.IPOINT) THEN
        IHIST=-IHIST
      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/:
      INTEGER MRPARH
      PARAMETER (MRPARH=1024)
      INTEGER JPOINT,INDMIN,NRPARH,IRPARH(0:MRPARH)
      COMMON/RPARH/JPOINT,INDMIN,NRPARH,IRPARH
C ALL 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
      SUBROUTINE RP2D(IRAY,G1NEW)
      INTEGER IRAY
      REAL G1NEW
C
C THIS SUBROUTINE CONTROLS ONE-PARAMETRIC BOUNDARY-VALUE RAY TRACING.
C AUXILIARY SUBROUTINE TO RPAR2.
C
C INPUT:
C     IRAY... INDEX OF THE ALREADY CALCULATED RAYS WITHIN THE SHOOTING
C             DOMAIN.  IRAY=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     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     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/:
      INTEGER MREC
      PARAMETER (MREC=128)
      INTEGER ISRFR,IPOINT,ISRFX(2),NREC
      REAL XERR,AERR,XREC(2,MREC)
      REAL PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
      COMMON/RPARD/ISRFR,IPOINT,ISRFX,NREC,AERR,XERR,XREC,
     *             PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C COMMON BLOCK /RP2DC/:
      INTEGER NEWMAX
      PARAMETER (NEWMAX=20)
      INTEGER ITER,NEW,INDOLD,IND(NEWMAX)
      REAL AOLD,A(NEWMAX),DAOLD,DA(NEWMAX)
      REAL XOLD,X(NEWMAX),XAOLD,XA(NEWMAX)
      COMMON/RP2DC/ITER,NEW,INDOLD,IND,AOLD,A,DAOLD,DA,XOLD,X,XAOLD,XA
C ALL STORAGE LOCATIONS OF THE COMMON BLOCK MAY BE ALTERED.
C
C NO SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED.
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: 1994, JANUARY 3
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(IRAY.EQ.0) THEN
        NEW=1
        A(1)=-G1STEP
      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
C         TWO SUCCESFUL RAYS - LOOKING FOR THE PROFILES IN THE 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
   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.
      INDOLD=IND(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
          IRAY=0
          RETURN
        END IF
        GO TO 80
      END IF
      GO TO 10
C
C     NEW RAY
   80 CONTINUE
      IRAY=IRAY+1
      NEW=NEW+1
      G1NEW=A(NEW)
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE RP2DM(ISHEET,X1,X1A)
      INTEGER ISHEET
      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     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 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 NO OUTPUT.
C
C COMMON BLOCK /RP2DC/:
      INTEGER NEWMAX
      PARAMETER (NEWMAX=20)
      INTEGER ITER,NEW,INDOLD,IND(NEWMAX)
      REAL AOLD,A(NEWMAX),DAOLD,DA(NEWMAX)
      REAL XOLD,X(NEWMAX),XAOLD,XA(NEWMAX)
      COMMON/RP2DC/ITER,NEW,INDOLD,IND,AOLD,A,DAOLD,DA,XOLD,X,XAOLD,XA
C STORAGE LOCATIONS IND,X,XA 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     NO AUXILIARY STORAGE LOCATIONS.
C
      IND(NEW)=ISHEET
      X(NEW)=X1
      XA(NEW)=X1A
      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/:
      INTEGER MREC
      PARAMETER (MREC=128)
      INTEGER ISRFR,IPOINT,ISRFX(2),NREC
      REAL XERR,AERR,XREC(2,MREC)
      REAL PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
      COMMON/RPARD/ISRFR,IPOINT,ISRFX,NREC,AERR,XERR,XREC,
     *             PAR1L,PAR2L,PAR1A,PAR2A,PAR1B,PAR2B,ANUM,BNUM
C NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK ARE ALTERED.
C
C COMMON BLOCK /RP2DC/:
      INTEGER NEWMAX
      PARAMETER (NEWMAX=20)
      INTEGER ITER,NEW,INDOLD,IND(NEWMAX)
      REAL AOLD,A(NEWMAX),DAOLD,DA(NEWMAX)
      REAL XOLD,X(NEWMAX),XAOLD,XA(NEWMAX)
      COMMON/RP2DC/ITER,NEW,INDOLD,IND,AOLD,A,DAOLD,DA,XOLD,X,XAOLD,XA
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