C SUBROUTINE FILE 'RAYCB.FOR' FOR COMPLETE RAY TRACING WITHIN ONE COMPLEX
C BLOCK
C
C BY VLASTISLAV CERVENY, LUDEK KLIMES, IVAN PSENCIK
C
C THIS FILE CONSISTS OF THE FOLLOWING SUBROUTINES:
C     RAYCB...SUBROUTINE TRANSFERING THE QUANTITIES GIVEN AT THE INITIAL
C             POINT OF THE ELEMENT OF A RAY INTO THE CORRESPONDING
C             QUANTITIES AT THE ENDPOINT OF THE ELEMENT (SEE
C             C.R.T.5.8.1).
C     FCT...  SUBROUTINE EVALUATING THE RIGHT-HAND SIDES OF THE SYSTEM
C             OF ORDINARY DIFFERENTIAL EQUATIONS FOR Y(1) TO Y(27) (SEE
C             C.R.T.5.8.2).
C     OUTP... SUBROUTINE CONTAINING VARIOUS TESTS OF THE POSITION OF THE
C             NEWLY COMPUTED POINT OF THE RAY, TESTS FOR POSSIBLE
C             CAUSTIC POINTS, ETC. (SEE C.R.T.5.8.3).
C     CDEF... AUXILIARY SUBROUTINE TO OUTP, PERFORMING STEPS 5.8.3(C),
C             (D), (E) AND (F) OF THE ALGORITHM.
C     PSHIFT..AUXILIARY SUBROUTINE TO OUTP AND CDEF. IT CORRECTS REDUCED
C             AMPLITUDES WITH REGARD TO PHASE SHIFT DUE TO CAUSTICS (SEE
C             C.R.T.5.8.3F).
C
C DATE: 1992, NOVEMBER 2
C
C=======================================================================
C
      SUBROUTINE RAYCB(YL,Y,YY,IY)
      REAL YL(6),Y(35),YY(5)
      INTEGER IY(12)
C
C THIS SUBROUTINE TRANSFERS THE QUANTITIES GIVEN AT THE INITIAL POINT OF
C THE ELEMENT OF A RAY INTO THE CORRESPONDING QUANTITIES AT THE ENDPOINT
C OF THE ELEMENT (SEE C.R.T.5.8.1).
C
C INPUT:
C     YL...   ARRAY CONTAINING LOCAL QUANTITIES AT THE INITIAL POINT OF
C             THE ELEMENT OF THE RAY.
C     Y...    ARRAY CONTAINING BASIC QUANTITIES COMPUTED ALONG THE RAY
C             AT THE INITIAL POINT OF THE ELEMENT OF THE RAY.
C     YY...   ARRAY CONTAINING REAL AUXILIARY QUANTITIES COMPUTED ALONG
C             THE RAY AT THE INITIAL POINT OF THE ELEMENT OF THE RAY.
C     IY...   ARRAY CONTAINING INTEGER AUXILIARY QUANTITIES COMPUTED
C             ALONG THE RAY AT THE INITIAL POINT OF THE ELEMENT OF THE
C             RAY.
C
C OUTPUT:
C     YL...   ARRAY CONTAINING LOCAL QUANTITIES AT THE ENDPOINT OF THE
C             ELEMENT OF THE RAY.
C     Y...    ARRAY CONTAINING BASIC QUANTITIES COMPUTED ALONG THE RAY
C             AT THE ENDPOINT OF THE ELEMENT OF THE RAY.
C     YY...   ARRAY CONTAINING REAL AUXILIARY QUANTITIES COMPUTED ALONG
C             THE RAY AT THE ENDPOINT OF THE ELEMENT OF THE RAY.
C     IY...   ARRAY CONTAINING INTEGER AUXILIARY QUANTITIES COMPUTED
C             ALONG THE RAY AT THE ENDPOINT OF THE ELEMENT OF THE RAY.
C
C COMMON BLOCK /DCRT/ (SEE 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/ (SEE SUBROUTINE FILE 'INIT.FOR') IS REQUIRED IN
C     SUBROUTINE OUTP.  NONE OF THE STORAGE LOCATIONS OF THE COMMON
C     BLOCK ARE ALTERED THERE.
C
C COMMON BLOCK /RAYC/ -  COMMUNICATION BETWEEN RAYCB, FCT, AND OUTP:
C     ------------------------------------------------------------------
      REAL YLRC(6),YYRC(5),DELT11,DELT13,DELT33
      INTEGER IYRC(12)
      COMMON/RAYC/YLRC,YYRC,IYRC,DELT11,DELT13,DELT33
C     ------------------------------------------------------------------
C     YLRC,YYRC,IYRC... WORKING COPIES OF THE ARRAYS YL,YY,YI.
C     DELT11,DELT13,DELT33... DEVIATIONS (IN ABSOLUTE VALUES OF THE TWO
C             COMPUTED BASIS VECTORS OF THE RAY-CENTRED COORDINATE
C             SYSTEM FROM THE CONDITIONS OF ORTHONORMALITY, SEE
C             C.R.T.5.8.2D.
C YLRC,YYRC,IYRC ARE DEFINED IN THIS SUBROUTINE.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL METRIC,HPCG
C     NSRFC,BLOCK,VELOC... FILE 'MODEL.FOR' OF THE PACKAGE 'MODEL'.
C     KOOR,METRIC... FILE 'METRIC.FOR' OF THE PACKAGE 'MODEL'.
C     SRFC2 AND SUBSEQUENT ROUTINES... FILE 'SRFC.FOR' AND SUBSEQUENT
C             FILES (LIKE 'VAL.FOR' AND 'FIT.FOR') OF THE PACKAGE
C             'MODEL'.
C     PARM2 AND SUBSEQUENT ROUTINES... FILE 'PARM.FOR' AND SUBSEQUENT
C             FILES (LIKE 'VAL.FOR' AND 'FIT.FOR') OF THE PACKAGE
C             'MODEL'.
C     CDE,CROSS,HIVD2,SMVPRD... FILE 'MEANS.FOR' OF THE PACKAGE 'MODEL'.
C     HPCG... IBM SCIENTIFIC SUBROUTINE PACKAGE (FILE 'HPCG.FOR' OF THE
C             PACKAGE 'MODEL').
C     RPAR31,RPAR32... FILE 'RPAR.FOR'.
C     WRIT31,WRIT32... FILE 'WRIT.FOR'.
C     FCT,OUTP,CDEF,PSHIFT... THIS FILE.
C
C DATE: 1993, JANUARY 15
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS FOR LOCAL MODEL PARAMETERS:
      REAL G(12),GAMMA(18),GSQRD,FAUX(10)
      REAL UP(10),US(10),RO,QP,QS,VP,VS,VD(10),QL
C     THESE AUXILIARY VARIABLES AND ARRAYS NEED NOT BE LOCATED IN A
C     COMMON BLOCK.  THERE IS NO REASON TO LOCATE THEM IN THE AUXILIARY
C     COMMON BLOCK /AUXMOD/ BUT TO SHARE THE MEMORY.
      COMMON/AUXMOD/ G,GAMMA,GSQRD,FAUX,UP,US,RO,QP,QS,VP,VS,VD,QL
C
C     G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE FILE 'METRIC.FOR'.
C     FAUX... AUXILIARY ARRAY TO STORE A FUNCTIONAL VALUE AND ITS
C             DERIVATIVES.
C     UP,US,RO,QP,QS... SEE SUBROUTINE PARM2 OF THE FILE 'PARM.FOR'.
C     UP,US,QP,QS,VP,VS,VD,QL... SEE SUBROUTINE VELOC OF THE FILE
C             'MODEL.FOR'.
C
C.......................................................................
C
C     EXTERNAL PROCEDURE NAMES:
      EXTERNAL FCT,OUTP
C     THESE SUBROUTINES ARE CALLED BY THE SUBROUTINE HPCG.
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER NDIM,IHLF,I
      PARAMETER (NDIM=27)
      REAL PRMT(6),DERY(NDIM),AUX(16,NDIM),VRED,AUX1
C
C     NDIM,IHLF,PRMT,DERY,AUX... SEE SUBROUTINE HPCG OF THE IBM
C             SCIENTIFIC SUBROUTINE PACKAGE.
C     I...    AUXILIARY LOOP VARIABLE.
C     VRED... APPROXIMATE REFERENCE VELOCITY TO ESTIMATE ERROR WEIGHTS.
C     AUX1... AUXILIARY STORAGE LOCATION.
C
C.......................................................................
C
C     (A) STORING AUXILIARY INPUT QUANTITIES IN COMMON BLOCK /RAYC/:
      DO 11 I=1,6
        YLRC(I)=YL(I)
   11 CONTINUE
      DO 12 I=1,5
        YYRC(I)=YY(I)
   12 CONTINUE
      DO 13 I=1,12
        IYRC(I)=IY(I)
   13 CONTINUE
C
C     (B) INITIAL VALUES FOR NUMERICAL INTEGRATION:
   20 PRMT(1)=YYRC(1)
      PRMT(2)=PRMT(1)+999999.
      PRMT(3)=STEP
      PRMT(4)=13.444*YYRC(2)
      CALL METRIC(Y(3),GSQRD,G,GAMMA)
      IF(IYRC(5).GE.0) THEN
        VRED=YLRC(1)
      ELSE
        VRED=YLRC(2)
      END IF
      DERY(1)=1.
      DERY(2)=1.
      DERY(3)=SQRT(G(1))/VRED
      DERY(4)=SQRT(G(3))/VRED
      DERY(5)=SQRT(G(6))/VRED
      AUX1=STEP*VRED**(1-NEXPS)
      DERY(6)=SQRT(G(7))*AUX1
      DERY(7)=SQRT(G(9))*AUX1
      DERY(8)=SQRT(G(12))*AUX1
      DO 21 I=9,NDIM
        DERY(I)=0.
   21 CONTINUE
C
C     (C) NUMERICAL INTEGRATION:
      CALL HPCG(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)
C
C     (D) GREAT NUMBER OF BISECTIONS:
      IF(PRMT(5).LT.0.) THEN
        YYRC(2)=2.*YYRC(2)
        GO TO 20
      END IF
C
C     (E) RECALLING AUXILIARY INPUT QUANTITIES FROM COMMON BLOCK /RAYC/:
      DO 51 I=1,6
        YL(I)=YLRC(I)
   51 CONTINUE
      DO 52 I=1,5
        YY(I)=YYRC(I)
   52 CONTINUE
      DO 53 I=1,12
        IY(I)=IYRC(I)
   53 CONTINUE
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE FCT(X,Y,D)
      INTEGER NDIM
      PARAMETER (NDIM=27)
      REAL X,Y(NDIM),D(NDIM)
C
C THIS SUBROUTINE EVALUATES THE RIGHT-HAND SIDES OF THE SYSTEM OF
C ORDINARY DIFFERENTIAL EQUATIONS FOR Y(1) TO Y(27). SEE C.R.T., SECTION
C 5.8.2.
C
C INPUT:
C     X...    VALUE OF THE INDEPENDENT VARIABLE ALONG THE RAY.
C     Y...    ARRAY CONTAINING BASIC QUANTITIES COMPUTED ALONG THE RAY.
C NONE OF THE INPUT PARAMETERS EXCEPT Y(3)-Y(8) IS ALTERED.
C
C OUTPUT:
C     D...    ARRAY CONTAINING DERIVATIVES OF THE BASIC QUANTITIES
C             COMPUTED ALONG THE RAY, WITH RESPECT TO THE INDEPENDENT
C             VARIABLE ALONG THE RAY.
C     Y(3:8)..RENORMALIZED ORTHOGONAL VECTORS.  THE MODIFICATION SHOULD
C             BE NEGLIGIBLE.
C
C COMMON BLOCK /DCRT/ (SEE 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 /RAYC/ -  COMMUNICATION BETWEEN RAYCB, FCT, AND OUTP:
      REAL YLRC(6),YYRC(5),DELT11,DELT13,DELT33
      INTEGER IYRC(12)
      COMMON/RAYC/YLRC,YYRC,IYRC,DELT11,DELT13,DELT33
C YLRC,IYRC ARE MODIFIED IN THIS SUBROUTINE.  DELT11,DELT13, AND DELT33
C ARE DEFINED IN THIS SUBROUTINE.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL KOOR,METRIC,PARM2,VELOC,SMVPRD
      INTEGER KOOR
C     VELOC... FILE 'MODEL.FOR' OF THE PACKAGE 'MODEL'.
C     KOOR,METRIC... FILE 'METRIC.FOR' OF THE PACKAGE 'MODEL'.
C     PARM2 AND SUBSEQUENT ROUTINES... FILE 'PARM.FOR' AND SUBSEQUENT
C             FILES (LIKE 'VAL.FOR' AND 'FIT.FOR') OF THE PACKAGE
C             'MODEL'.
C     SMVPRD... FILE 'MEANS.FOR' OF THE PACKAGE 'MODEL'.
C
C DATE: 1993, JANUARY 23
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS FOR LOCAL MODEL PARAMETERS:
      REAL G(12),GAMMA(18),GSQRD,FAUX(10)
      REAL UP(10),US(10),RO,QP,QS,VP,VS,VD(10),QL
C     THESE AUXILIARY VARIABLES AND ARRAYS NEED NOT BE LOCATED IN A
C     COMMON BLOCK.  THERE IS NO REASON TO LOCATE THEM IN THE AUXILIARY
C     COMMON BLOCK /AUXMOD/ BUT TO SHARE THE MEMORY.
      COMMON/AUXMOD/ G,GAMMA,GSQRD,FAUX,UP,US,RO,QP,QS,VP,VS,VD,QL
C
C     G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE FILE 'METRIC.FOR'.
C     FAUX... AUXILIARY ARRAY TO STORE A FUNCTIONAL VALUE AND ITS
C             DERIVATIVES.
C     UP,US,RO,QP,QS... SEE SUBROUTINE PARM2 OF THE FILE 'PARM.FOR'.
C     UP,US,QP,QS,VP,VS,VD,QL... SEE SUBROUTINE VELOC OF THE FILE
C             'MODEL.FOR'.
C
C.......................................................................
C
C     AUXILIARY STORAGE LOCATIONS:
C
      REAL H11,H21,H31,H41,H51,H61,H42,H52,H62,H13,H23,H33,H43,H53,H63
      REAL D11,D13,D33,V,V1,V11,V12,V22,VNEXPS,AUX1,AUX2,AUX3,AUX4
      REAL AUX11,AUX21,AUX31,AUX12,AUX22,AUX32,AUX13,AUX23,AUX33
C
C     H11,H21,H31,H13,H23,H33... COVARIANT COMPONENTS OF THE BASIS
C             VECTORS OF RAY-CENTRED COORDINATE SYSTEM.
C     H41,H51,H61,H42,H52,H62,H43,H53,H63... CONTRAVARIANT COMPONENTS OF
C             THE BASIS VECTORS OF RAY-CENTRED COORDINATE SYSTEM.  H41,
C             H51,H61,H43,H53,H63 ARE NOT USED IN CARTESIAN COORDINATES.
C     D11,D13,D33... NORMS AND SCALAR PRODUCT OF BASIS VECTORS H1, H2.
C     V,V1,V11,V12,V22,... VELOCITY AND VELOCITY DERIVATIVES IN
C             RAY-CENTRED COORDINATE SYSTEM.
C     VNEXPS..VELOCITY**(-NEXPS).
C     AUX1,AUX2,AUX3,AUX4... AUXILIARY STORAGE LOCATIONS.
C     AUX11,AUX21,AUX31,AUX12,AUX22,AUX32,AUX13,AUX23,AUX33... AUXILIARY
C             STORAGE LOCATIONS.  NOT USED IN CARTESIAN COORDINATES.
C
C.......................................................................
C
C     (A) NUMBER OF INVOCATIONS OF FCT:
      IYRC(9)=IYRC(9)+1
C
C     (B) MATERIAL PARAMETERS:
      CALL PARM2(IABS(IYRC(5)),Y(3),UP,US,YLRC(3),QP,QS)
      CALL VELOC(IYRC(5),UP,US,QP,QS,YLRC(1),YLRC(2),VD,QL)
      YLRC(4)=VD(2)
      YLRC(5)=VD(3)
      YLRC(6)=VD(4)
      V=VD(1)
      VNEXPS=V**(-NEXPS)
C
C     (C) BASIS OF THE RAY-CENTRED COORDINATE SYSTEM:
C     (C1) NON-UNIT VECTORS - COVARIANT COMPONENTS
      H13=Y(6)
      H23=Y(7)
      H33=Y(8)
      H11=Y(9)
      H21=Y(10)
      H31=Y(11)
      IF(KOOR().NE.0) THEN
C       CURVILINEAR COORDINATES:
        CALL METRIC(Y(3),GSQRD,G,GAMMA)
C       (C2) NON-UNIT VECTORS - CONTRAVARIANT COMPONENTS
        CALL SMVPRD(G(7),H13,H23,H33,H43,H53,H63)
        CALL SMVPRD(G(7),H11,H21,H31,H41,H51,H61)
C       (C3) NORMS:
        D33=SQRT(H13*H43+H23*H53+H33*H63)
        D11=SQRT(H11*H41+H21*H51+H31*H61)
        D13=(H11*H43+H21*H53+H31*H63)/(D11*D33)
C       (C4) ORTHONORMAL VECTORS:
C       COVARIANT COMPONENTS
        H13=H13/D33
        H23=H23/D33
        H33=H33/D33
        H11=H11/D11-H13*D13
        H21=H21/D11-H23*D13
        H31=H31/D11-H33*D13
C       CONTRAVARIANT COMPONENTS
        H43=H43/D33
        H53=H53/D33
        H63=H63/D33
        H41=H41/D11-H43*D13
        H51=H51/D11-H53*D13
        H61=H61/D11-H63*D13
        H42=(H23*H31-H33*H21)/GSQRD
        H52=(H33*H11-H13*H31)/GSQRD
        H62=(H13*H21-H23*H11)/GSQRD
C
C       (D) QUANTITIES USEFUL TO TEST THE ACCURACY OF COMPUTATIONS:
        DELT11=ABS(D11-1.)
        DELT33=ABS(V*D33-1.)
        DELT11=ABS(D13)
C
C       (E) FIRST AND SECOND DERIVATIVES OF THE VELOCITY IN RAY-CENTRED
C       COORDINATE SYSTEM:
C       SECOND COVARIANT VELOCITY DERIVATIVES:
        VD(5)=VD(5)-GAMMA(1)*VD(2)-GAMMA(7)*VD(3)-GAMMA(13)*VD(4)
        VD(6)=VD(6)-GAMMA(2)*VD(2)-GAMMA(8)*VD(3)-GAMMA(14)*VD(4)
        VD(7)=VD(7)-GAMMA(3)*VD(2)-GAMMA(9)*VD(3)-GAMMA(15)*VD(4)
        VD(8)=VD(8)-GAMMA(4)*VD(2)-GAMMA(10)*VD(3)-GAMMA(16)*VD(4)
        VD(9)=VD(9)-GAMMA(5)*VD(2)-GAMMA(11)*VD(3)-GAMMA(17)*VD(4)
        VD(10)=VD(10)-GAMMA(6)*VD(2)-GAMMA(12)*VD(3)-GAMMA(18)*VD(4)
C       RAY-CENTRED COORDINATE SYSTEM:
        V1=VD(2)*H41+VD(3)*H51+VD(4)*H61
        CALL SMVPRD(VD(5),H41,H51,H61,AUX1,AUX2,AUX3)
        V11=AUX1*H41+AUX2*H51+AUX3*H61
        V12=AUX1*H42+AUX2*H52+AUX3*H62
        CALL SMVPRD(VD(5),H42,H52,H62,AUX1,AUX2,AUX3)
        V22=AUX1*H42+AUX2*H52+AUX3*H62
C
C       (F) CORRECTION OF THE COMPUTED QUANTITIES:
        Y(6)=H13/V
        Y(7)=H23/V
        Y(8)=H33/V
        Y(9)=H11
        Y(10)=H21
        Y(11)=H31
C
C       (G) RIGHT-HAND SIDES OF THE DIFFERENTIAL EQUATIONS:
        AUX1=V*VNEXPS
        AUX2=V*AUX1
        AUX3=VNEXPS/V
        AUX4=V1*VNEXPS
        D(1)=VNEXPS
        D(2)=0.5*QL*VNEXPS
        D(3)=H43*AUX1
        D(4)=H53*AUX1
        D(5)=H63*AUX1
        CALL SMVPRD(GAMMA(1),H43,H53,H63,AUX11,AUX21,AUX31)
        CALL SMVPRD(GAMMA(7),H43,H53,H63,AUX12,AUX22,AUX32)
        CALL SMVPRD(GAMMA(13),H43,H53,H63,AUX13,AUX23,AUX33)
        D(6)=-VD(2)*AUX3+(AUX11*H13+AUX12*H23+AUX13*H33)*VNEXPS
        D(7)=-VD(3)*AUX3+(AUX21*H13+AUX22*H23+AUX23*H33)*VNEXPS
        D(8)=-VD(4)*AUX3+(AUX31*H13+AUX32*H23+AUX33*H33)*VNEXPS
        D(9) =H13*AUX4+(AUX11*H11+AUX12*H21+AUX13*H31)*AUX1
        D(10)=H23*AUX4+(AUX21*H11+AUX22*H21+AUX23*H31)*AUX1
        D(11)=H33*AUX4+(AUX31*H11+AUX32*H21+AUX33*H31)*AUX1
      ELSE
C       CARTESIAN COORDINATES (20 OF THESE STATEMENT LINES ARE THE SAME
C       AS FOR CURVILINEAR COORDINATES):
C       (C2) NON-UNIT VECTORS - CONTRAVARIANT COMPONENTS
C       H43=H13, H53=H23, H63=H33, H41=H11, H51=H21, H61=H31
C       (C3) NORMS:
        D33=SQRT(H13*H13+H23*H23+H33*H33)
        D11=SQRT(H11*H11+H21*H21+H31*H31)
        D13=(H11*H13+H21*H23+H31*H33)/(D11*D33)
C       (C4) ORTHONORMAL VECTORS:
C       COVARIANT=CONTRAVARIANT COMPONENTS
        H13=H13/D33
        H23=H23/D33
        H33=H33/D33
        H11=H11/D11-H13*D13
        H21=H21/D11-H23*D13
        H31=H31/D11-H33*D13
        H42=(H23*H31-H33*H21)
        H52=(H33*H11-H13*H31)
        H62=(H13*H21-H23*H11)
C
C       (D) QUANTITIES USEFUL TO TEST THE ACCURACY OF COMPUTATIONS:
        DELT11=ABS(D11-1.)
        DELT33=ABS(V*D33-1.)
        DELT11=ABS(D13)
C
C       (E) FIRST AND SECOND DERIVATIVES OF THE VELOCITY IN RAY-CENTRED
C       COORDINATE SYSTEM:
        V1=VD(2)*H11+VD(3)*H21+VD(4)*H31
        CALL SMVPRD(VD(5),H11,H21,H31,AUX1,AUX2,AUX3)
        V11=AUX1*H11+AUX2*H21+AUX3*H31
        V12=AUX1*H42+AUX2*H52+AUX3*H62
        CALL SMVPRD(VD(5),H42,H52,H62,AUX1,AUX2,AUX3)
        V22=AUX1*H42+AUX2*H52+AUX3*H62
C
C       (F) CORRECTION OF THE COMPUTED QUANTITIES:
        Y(6)=H13/V
        Y(7)=H23/V
        Y(8)=H33/V
        Y(9)=H11
        Y(10)=H21
        Y(11)=H31
C
C       (G) RIGHT-HAND SIDES OF THE DIFFERENTIAL EQUATIONS:
        AUX1=V*VNEXPS
        AUX2=V*AUX1
        AUX3=VNEXPS/V
        AUX4=V1*VNEXPS
        D(1)=VNEXPS
        D(2)=0.5*QL*VNEXPS
        D(3)=H13*AUX1
        D(4)=H23*AUX1
        D(5)=H33*AUX1
        D(6)=-VD(2)*AUX3
        D(7)=-VD(3)*AUX3
        D(8)=-VD(4)*AUX3
        D(9)=H13*AUX4
        D(10)=H23*AUX4
        D(11)=H33*AUX4
      END IF
      V11=-V11*AUX3
      V12=-V12*AUX3
      V22=-V22*AUX3
      D(12)=AUX2*Y(14)
      D(13)=AUX2*Y(15)
      D(14)=V11*Y(12)+V12*Y(13)
      D(15)=V12*Y(12)+V22*Y(13)
      D(16)=AUX2*Y(18)
      D(17)=AUX2*Y(19)
      D(18)=V11*Y(16)+V12*Y(17)
      D(19)=V12*Y(16)+V22*Y(17)
      D(20)=AUX2*Y(22)
      D(21)=AUX2*Y(23)
      D(22)=V11*Y(20)+V12*Y(21)
      D(23)=V12*Y(20)+V22*Y(21)
      D(24)=AUX2*Y(26)
      D(25)=AUX2*Y(27)
      D(26)=V11*Y(24)+V12*Y(25)
      D(27)=V12*Y(24)+V22*Y(25)
C
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE OUTP(X,Y,D,IHLF,NDIM,PRMT)
      INTEGER IHLF,NDIM,MY
      PARAMETER (MY=35)
      REAL X,Y(MY),D(NDIM),PRMT(5)
C
C THIS SUBROUTINE INCLUDES VARIOUS TESTS OF THE POSITION OF THE NEWLY
C COMPUTED POINT OF THE RAY, TESTS FOR POSSIBLE CAUSTIC POINTS, ETC.  IT
C ALSO STORES THE COMPUTED QUANTITIES INTO PROPER FILES IF REQUIRED.
C
C INPUT:
C     X...    VALUE OF THE INDEPENDENT VARIABLE ALONG THE RAY.
C     Y...    ARRAY CONTAINING BASIC QUANTITIES COMPUTED ALONG THE RAY.
C     D...    ARRAY CONTAINING DERIVATIVES OF THE BASIC QUANTITIES
C             COMPUTED ALONG THE RAY, WITH RESPECT TO THE INDEPENDENT
C             VARIABLE ALONG THE RAY.
C     IHLF... NUMBER OF BISECTIONS OF THE INITIAL INCREMENT OF
C             INDEPENDENT VARIABLE.
C     NDIM... NUMBER OF ORDINARY DIFFERENTIAL EQUATIONS.
C     PRMT... ARRAY CONTAINING PARAMETERS OF THE INTEGRATION.
C NONE OF THE INPUT PARAMETERS EXCEPT X,Y,D,PRMT(5) IS ALTERED.
C
C OUTPUT:
C     X,Y,D...VALUES CORRESPONDING TO THE POINT OF INTERSECTION OF THE
C             RAY WITH THE BOUNDARY OF THE COMPLEX BLOCK OR THE
C             COMPUTATIONAL VOLUME IF THE RAY INTERSECTS THE BOUNDARY.
C             UNALTERED IF THE RAY DOES NOT INTERSECT THE BOUNDARY.
C     PRMT(5)=1.0... THE RAY HAS INTERSECTED THE BOUNDARY OF THE COMPLEX
C               BLOCK OR THE COMPUTATIONAL VOLUME,
C             -1.0... NUMBER IHLF OF BISECTIONS OF THE INITIAL INCREMENT
C               GREATER THAN THE SPECIFIED LIMIT NHLF,
C             OTHERWISE UNALTERED.
C
C COMMON BLOCK /DCRT/ (SEE 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/ (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 NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK EXCEPT FSRFCA ARE
C ALTERED.
C
C COMMON BLOCK /RAYC/ -  COMMUNICATION BETWEEN RAYCB, FCT, AND OUTP:
      REAL YLRC(6),YYRC(5),DELT11,DELT13,DELT33
      INTEGER IYRC(12)
      COMMON/RAYC/YLRC,YYRC,IYRC,DELT11,DELT13,DELT33
C YLRC,YYRC,IYRC ARE MODIFIED IN THIS SUBROUTINE.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL SRFC2,BLOCK,RPAR31,RPAR32,WRIT31,WRIT32,CROSS,HIVD2
      EXTERNAL NSRFC,PSHIFT,CDEF
      INTEGER NSRFC

C     NSRFC,BLOCK,VELOC... FILE 'MODEL.FOR' OF THE PACKAGE 'MODEL'.
C     SRFC2 AND SUBSEQUENT ROUTINES... FILE 'SRFC.FOR' AND SUBSEQUENT
C             FILES (LIKE 'VAL.FOR' AND 'FIT.FOR') OF THE PACKAGE
C             'MODEL'.
C     PARM2 AND SUBSEQUENT ROUTINES... FILE 'PARM.FOR' AND SUBSEQUENT
C             FILES (LIKE 'VAL.FOR' AND 'FIT.FOR') OF THE PACKAGE
C             'MODEL'.
C     CDE,CROSS,HIVD2... FILE 'MEANS.FOR' OF THE PACKAGE 'MODEL'.
C     RPAR31,RPAR32... FILE 'RPAR.FOR'.
C     WRIT31,WRIT32... FILE 'WRIT.FOR'.
C     CDEF,PSHIFT... THIS FILE.
C
C DATE: 1993, JANUARY 15
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS FOR LOCAL MODEL PARAMETERS:
      REAL G(12),GAMMA(18),GSQRD,FAUX(10)
      REAL UP(10),US(10),RO,QP,QS,VP,VS,VD(10),QL
C     THESE AUXILIARY VARIABLES AND ARRAYS NEED NOT BE LOCATED IN A
C     COMMON BLOCK.  THERE IS NO REASON TO LOCATE THEM IN THE AUXILIARY
C     COMMON BLOCK /AUXMOD/ BUT TO SHARE THE MEMORY.
      COMMON/AUXMOD/ G,GAMMA,GSQRD,FAUX,UP,US,RO,QP,QS,VP,VS,VD,QL
C
C     G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE FILE 'METRIC.FOR'.
C     FAUX... AUXILIARY ARRAY TO STORE A FUNCTIONAL VALUE AND ITS
C             DERIVATIVES.
C     UP,US,RO,QP,QS... SEE SUBROUTINE PARM2 OF THE FILE 'PARM.FOR'.
C     UP,US,QP,QS,VP,VS,VD,QL... SEE SUBROUTINE VELOC OF THE FILE
C             'MODEL.FOR'.
C
C.......................................................................
C
C     OTHER AUXILIARY STORAGE LOCATIONS:
C
      REAL TLAST,X0,X1,X2,Y1(MY),Y2(MY),D1(MY),D2(MY)
      SAVE X2,Y2,D2
      REAL XB,YB(MY),DB(MY),XAUX,YAUX(MY),DAUX(MY),AUX,ERR
      INTEGER I,J,K,K1,K2
C
C     TLAST...TRAVEL TIME AT THE LAST BUT ONE COMPUTED POINT OF THE RAY.
C     X0...   BEGINNING OF THE RAY ELEMENT, FOR X BETWEEN PRMT(1) AND
C             X0, IS NOT CHECKED FOR CROSSING A STRUCTURAL INTERFACE.
C     X1,Y1,D1... INDEPENDENT VARIABLE, COMPUTED QUANTITIES AND THEIR
C             DERIVATIVES AT THE LAST CHECKED POINT OF THE RAY.
C     X2,Y2,D2... INDEPENDENT VARIABLE, COMPUTED QUANTITIES AND THEIR
C             DERIVATIVES AT THE LAST COMPUTED POINT OF THE RAY.
C     XB,YB,DB,XAUX,YAUX,DAUX,AUX,I,J,K,K1,K2... AUXILIARY STORAGE
C             LOCATIONS.
C     ERR...  MAXIMUM ERROR IN INDEPENDENT VARIABLE FOR THE
C             DETERMINATION OF THE POINT OF INTERSECTION.
C
C.......................................................................
C
C     (A) COPYING THE COMPUTED QUANTITIES:
      TLAST=Y2(1)
      X1=X2
      X2=X
      DO 11 I=1,NDIM
        Y1(I)=Y2(I)
        Y2(I)=Y(I)
        D1(I)=D2(I)
        D2(I)=D(I)
   11 CONTINUE
      DO 12 I=NDIM+1,IYRC(1)
        Y1(I)=Y(I)
   12 CONTINUE
C     X1 IS GOING TO BE SHIFTED ALONG THE RAY SEGMENT TOWARDS TO THE
C     ENDPOINT OF THE SEGMENT DURING THE PROCESSING, WHEREAS TLAST WILL
C     RETAIN THE PRESENT VALUE OF Y1(1).
      IF(X.EQ.PRMT(1)) THEN
C       BEGINNIG OF THE NUMERICAL INTEGRATION
        CALL PSHIFT(27,Y,YI,I)
        RETURN
      END IF
C
C     (B) NUMBER OF INVOCATIONS OF OUTP:
      IYRC(10)=IYRC(10)+1
C     INITIAL VALUE OF THE INDEX OF CROSSED SURFACE BOUNDING THE COMPLEX
C     BLOCK
      IYRC(6)=0
C
C     (C), (D), (E) AND (F) ARE LOCATED IN THE SUBROUTINE CDEF.
      ERR=YYRC(2)*D(1)
      X0=PRMT(1)+ERR
C
C     (H) STORAGE OF THE COMPUTED QUANTITIES ALONG THE RAY:
C     XAUX... POINTS ALONG RAY, AND THE ENDPOINT OF THE SEGMENT.
      IF(STORE.GT.0.) THEN
        K1=INT(X1/STORE+0.000001)+1
        K2=INT(X/STORE+0.000001)
      ELSE
        K1=1
        K2=0
      END IF
      DO 84 K=K1,K2+1
        IF(K.LE.K2) THEN
          XAUX=FLOAT(K)*STORE
          CALL HIVD2(NDIM,X1,Y1,D1,X2,Y2,D2,XAUX,YAUX,DAUX)
        ELSE
          XAUX=X2
          DO 70 I=1,NDIM
            YAUX(I)=Y2(I)
            DAUX(I)=D2(I)
   70     CONTINUE
        END IF
C
C       (G) STORAGE OF THE COMPUTED QUANTITIES AT GIVEN SURFACES:
C       X... INTERSECTIONS WITH THE SURFACES BETWEEN X1 AND XAUX.
   71   CONTINUE
        X=XAUX
        DO 72 I=1,NDIM
          Y(I)=YAUX(I)
          D(I)=DAUX(I)
   72   CONTINUE
        J=0
        DO 73 I=1,NSTOR
          IF(KSTOR(I).GT.NSRFC()) THEN
            CALL SRFC2(KSTOR(K),Y(3),FAUX)
            IF(FAUX(1)*FSRFCA(KSTOR(I)-NSRFC()).LE.0.)  THEN
              J=I
              AUX=FAUX(1)
              CALL CROSS(SRFC2,KSTOR(I),3,5,NDIM,ERR,X1,Y1,D1,X2,Y2,D2,
     *                                              X,Y,D,XB,YB,DB,FAUX)
            END IF
          END IF
   73   CONTINUE
        IF(J.NE.0) THEN
          FSRFCA(KSTOR(K)-NSRFC())=AUX
          CALL CDEF(NDIM,ERR,X0,X1,Y1,D1,X2,Y2,D2,X,Y,D,XB,YB,DB)
          IF(IYRC(6).NE.0) THEN
C           BOUNDARY OF THE COMPLEX BLOCK IS CROSSED
C           X1=X... POINT OF INTERSECTION WITH THE BOUNDARY.
            GO TO 90
          END IF
C         X1=X... POINT OF INTERSECTION WITH THE LAST STORING SURFACE.
          CALL RPAR32(J,YLRC,Y1,YYRC,IYRC)
          CALL WRIT32(J,YLRC,Y1,YYRC,IYRC,IYRC(1)-27,Y1)
          GO TO 71
        END IF
C       (G-END)
C
        CALL CDEF(NDIM,ERR,X0,X1,Y1,D1,X2,Y2,D2,XAUX,YAUX,DAUX,XB,YB,DB)
        IF(IYRC(6).NE.0) THEN
C         BOUNDARY OF THE COMPLEX BLOCK IS CROSSED
C         X1=XAUX... POINT OF INTERSECTION WITH THE BOUNDARY.
          GO TO 90
        END IF
C       X1=XAUX... THE LAST POINT ALONG THE RAY.
        IF(K.LE.K2) THEN
          CALL RPAR31(YLRC,Y1,YYRC,IYRC)
          CALL WRIT31(YLRC,Y1,YYRC,IYRC)
        END IF
   84 CONTINUE
C     (H-END)
C
   90 CONTINUE
C     X1... ENDPOINT OF THE SEGMENT.
      DO 91 K=1,NSTOR
        IF(KSTOR(K).GT.NSRFC()) THEN
          IF(KSTOR(K).EQ.IYRC(6)) THEN
            CALL RPAR32(K,YLRC,Y1,YYRC,IYRC)
            CALL WRIT32(K,YLRC,Y1,YYRC,IYRC,IYRC(1)-27,Y1)
          END IF
        END IF
   91 CONTINUE
      X=X1
      DO 92 I=1,NDIM
        Y(I)=Y1(I)
        D(I)=D1(I)
   92 CONTINUE
      DO 93 I=NDIM+1,IYRC(1)
        Y(I)=Y1(I)
   93 CONTINUE
C     X=X1... ENDPOINT OF THE SEGMENT.
C
C     (I) ACCUMULATION OF THE RENORMALIZATION ERRORS FOR A TEST OF
C     ACCURACY:
      AUX=0.5*(Y(1)-TLAST)
      YYRC(3)=YYRC(3)+DELT33*AUX
      YYRC(4)=YYRC(4)+DELT13*AUX
      YYRC(5)=YYRC(5)+DELT11*AUX
C
C     (J) GREAT NUMBER OF BISECTIONS OF THE INITIAL INCREMENT:
      IF(IHLF.GT.NHLF) PRMT(5)=-1.
C
C     (K) FINAL OPERATIONS:
C     YYRC(1)=X  IS SET BY THE SUBROUTINE CDEF
      IF(IYRC(6).NE.0) PRMT(5)=1.
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE CDEF(NDIM,ERR,X0,X1,Y1,D1,X2,Y2,D2,X,Y,D,XB,YB,DB)
      INTEGER NDIM,MY
      PARAMETER (MY=35)
      REAL ERR,X0,X1,Y1(MY),D1(NDIM),X2,Y2(MY),D2(NDIM)
      REAL X,Y(MY),D(NDIM),XB,YB(MY),DB(NDIM)
C
C AUXILIARY SUBROUTINE TO OUTP, PERFORMING STEPS 5.8.3(C),(D),(E) AND
C (F) OF THE ALGORITHM.
C
C COMMON BLOCK /DCRT/ (SEE 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/ (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 NONE OF THE STORAGE LOCATIONS OF THE COMMON BLOCK EXCEPT FSRFCA ARE
C ALTERED.
C
C COMMON BLOCK /RAYC/ -  COMMUNICATION BETWEEN RAYCB, FCT, AND OUTP:
      REAL YLRC(6),YYRC(5),DELT11,DELT13,DELT33
      INTEGER IYRC(12)
      COMMON/RAYC/YLRC,YYRC,IYRC,DELT11,DELT13,DELT33
C YLRC,YYRC,IYRC ARE MODIFIED IN THIS SUBROUTINE.
C
C SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED:
      EXTERNAL CDE,PSHIFT,PARM2,VELOC
C     NSRFC,BLOCK,VELOC... FILE 'MODEL.FOR' OF THE PACKAGE 'MODEL'.
C     SRFC2 AND SUBSEQUENT ROUTINES... FILE 'SRFC.FOR' AND SUBSEQUENT
C             FILES (LIKE 'VAL.FOR' AND 'FIT.FOR') OF THE PACKAGE
C             'MODEL'.
C     PARM2 AND SUBSEQUENT ROUTINES... FILE 'PARM.FOR' AND SUBSEQUENT
C             FILES (LIKE 'VAL.FOR' AND 'FIT.FOR') OF THE PACKAGE
C             'MODEL'.
C     CDE,CROSS,HIVD2... FILE 'MEANS.FOR' OF THE PACKAGE 'MODEL'.
C
C DATE: 1993, JANUARY 15
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS FOR LOCAL MODEL PARAMETERS:
      REAL G(12),GAMMA(18),GSQRD,FAUX(10)
      REAL UP(10),US(10),RO,QP,QS,VP,VS,VD(10),QL
C     THESE AUXILIARY VARIABLES AND ARRAYS NEED NOT BE LOCATED IN A
C     COMMON BLOCK.  THERE IS NO REASON TO LOCATE THEM IN THE AUXILIARY
C     COMMON BLOCK /AUXMOD/ BUT TO SHARE THE MEMORY.
      COMMON/AUXMOD/ G,GAMMA,GSQRD,FAUX,UP,US,RO,QP,QS,VP,VS,VD,QL
C
C     G,GAMMA,GSQRD... SEE SUBROUTINE METRIC OF THE FILE 'METRIC.FOR'.
C     FAUX... AUXILIARY ARRAY TO STORE A FUNCTIONAL VALUE AND ITS
C             DERIVATIVES.
C     UP,US,RO,QP,QS... SEE SUBROUTINE PARM2 OF THE FILE 'PARM.FOR'.
C     UP,US,QP,QS,VP,VS,VD,QL... SEE SUBROUTINE VELOC OF THE FILE
C             'MODEL.FOR'.
C
C.......................................................................
C
C     OTHER AUXILIARY STORAGE LOCATIONS:
      INTEGER IBOUND(7),I
      DATA IBOUND/3,-3,4,-4,5,-5,-1/
C
C.......................................................................
C
C     (C) CHECK FOR CROSSING THE COORDINATE BOUNDARIES OF THE
C     (D) CHECK FOR CROSSING THE BOUNDARY OF THE COMPLEX BLOCK
C     (E) CHECK FOR CROSSING THE END SURFACES
      CALL CDE(0,NEND,KEND,7,IBOUND,BOUNDR,
     *            3,5,NDIM,IYRC,ERR,X0,X1,Y1,D1,X2,Y2,D2,X,Y,D,XB,YB,DB)
C
C     (F) PHASE SHIFT DUE TO CAUSTICS:
      X1=X
      DO 61 I=1,NDIM
        Y1(I)=Y(I)
        D1(I)=D(I)
   61 CONTINUE
      CALL PSHIFT(IYRC(1),Y1,YI,I)
      IYRC(12)=IYRC(12)+I
C
C     QUANTITIIES DESCRIBING LOCAL PROPERTIES OF THE MODEL AT POINT X,Y:
      IF(IYRC(6).NE.0) THEN
        CALL PARM2(IABS(IYRC(5)),Y(3),UP,US,YLRC(3),QP,QS)
        CALL VELOC(IYRC(5),UP,US,QP,QS,YLRC(1),YLRC(2),VD,QL)
        YLRC(4)=VD(2)
        YLRC(5)=VD(3)
        YLRC(6)=VD(4)
      END IF
      YYRC(1)=X
C
C     X1=X,Y1=Y,D1=D,YLRC,YYRC,IYRC CONTAIN THE QUANTITIES EITHER AT THE
C     GIVEN POINT X,Y,D (FOR IYRC(6).EQ.0) OR AT THE ENDPOINT OF THE
C     COMPUTED RAY ELEMENT (FOR IYRC(6).NE.0).
      RETURN
      END
C
C=======================================================================
C
      SUBROUTINE PSHIFT(NY,Y,YI,ISHIFT)
      INTEGER NY,ISHIFT
      REAL Y(NY),YI(21)
C
C THIS SUBROUTINE IS AN AUXILIARY ROUTINE TO OUTP. IT CORRECTS REDUCED
C AMPLITUDES WITH REGARD TO PHASE SHIFT DUE TO CAUSTICS.
C (SEE C.R.T.5.8.3F).  HEREINAFTER, THE POINT 1 OF THE RAY DENOTES THE
C POINT CORRESPONDING TO THE PREVIOUS INVOCATION OF THE SUBROUTINE
C PSHIFT, THE POINT 2 DENOTES THE POINT CORRESPONDING TO THIS INVOCATION
C OF THE SUBROUTINE PSHIFT.
C
C INPUT:
C     NY...   NUMBER OF BASIC QUANTITIES DESCRIBING THE POINT OF A RAY.
C     Y(1:11)... ANYTHING.
C     Y(12:27)... RAY PROPAGATOR MATRIX AT THE POINT 2 OF THE RAY.
C     Y(28:NY)... REDUCED AMPLITUDES AT THE POINT 1 OF THE RAY.
C     YI...   ARRAY CONTAINING QUANTITIES AT THE INITIAL POINT OF THE
C             RAY (SEE C.R.T.6.1)
C NONE OF THE INPUT PARAMETERS EXCEPT Y(28:NY) ARE ALTERED.
C
C OUTPUT:
C     Y(28:NY)... REDUCED AMPLITUDES AT THE POINT 2 OF THE RAY.
C     ISHIFT..ORDER OF A CAUSTIC BETWEEN POINTS 1 AND 2 (INCREMENT OF
C             THE KMAH INDEX).
C
C NO SUBROUTINES AND EXTERNAL FUNCTIONS REQUIRED.
C
C DATE: 1990, NOVEMBER 15
C CODED BY LUDEK KLIMES
C
C-----------------------------------------------------------------------
C
C     AUXILIARY STORAGE LOCATIONS:
      INTEGER I
      REAL AUX
      REAL Q111,Q121,Q112,Q122,Q1DET
      REAL Q211,Q221,Q212,Q222,Q2DET
      SAVE Q211,Q221,Q212,Q222,Q2DET
C
C     I,AUX...AUXILIARY STORAGE LOCATIONS.
C     Q111,Q121,Q112,Q122,Q1DET... MATRIX OF RAY GEOMETRICAL SPREADING
C             AND ITS DETERMINANT CORRESPONDING TO THE PREVIOUS
C             INVOCATION OF PSHIFT AT THE POINT 1.
C     Q211,Q221,Q212,Q222,Q2DET... MATRIX OF RAY GEOMETRICAL SPREADING
C             AND ITS DETERMINANT CORRESPONDING TO THE LAST INVOCATION
C             OF PSHIFT.
C
      Q111=Q211
      Q121=Q221
      Q112=Q212
      Q122=Q222
      Q1DET=Q2DET
      Q211=Y(12)*YI(12)+Y(16)*YI(13)+Y(20)*YI(14)+Y(24)*YI(15)
      Q221=Y(13)*YI(12)+Y(17)*YI(13)+Y(21)*YI(14)+Y(25)*YI(15)
      Q212=Y(12)*YI(16)+Y(16)*YI(17)+Y(20)*YI(18)+Y(24)*YI(19)
      Q222=Y(13)*YI(16)+Y(17)*YI(17)+Y(21)*YI(18)+Y(25)*YI(19)
      Q2DET=Q211*Q222-Q212*Q221
      IF(Q1DET*Q2DET.LT.0.) THEN
C       CAUSTIC OF THE FIRST ORDER (LINE CAUSTIC)
        ISHIFT=1
        DO 10 I=28,NY,2
          AUX=Y(I+1)
          Y(I+1)=-Y(I)
          Y(I)=AUX
   10   CONTINUE
      ELSE IF(((Q111*Q222-Q112*Q221)+(Q211*Q122-Q212*Q121))*Q2DET.LT.0.)
     *     THEN
C       CAUSTIC OF THE SECOND ORDER (POINT CAUSTIC)
        ISHIFT=2
        DO 20 I=28,NY
          Y(I)=-Y(I)
   20   CONTINUE
      ELSE
C       NO CAUSTIC
        ISHIFT=0
      END IF
      RETURN
      END
C
C=======================================================================
C