C 
C Subroutine file 'raycb.for' for complete ray tracing within one
C complex block
C
C Date: 2004, May 19
C Coded by: Ludek Klimes
C
C.......................................................................
C
C This file consists of the following subroutines:
C     RAYCB...Subroutine transferring the quantities given at the
C             initial point of the element of a ray into the
C             corresponding quantities at the endpoint of the element,
C             see
C             C.R.T.5.8.1.
C             RAYCB
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             FCT
C     FCTA... 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             for anisotropic ray tracing.
C             FCTA
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             C.R.T.5.8.3.
C             OUTP
C     CDEF... Auxiliary subroutine to OUTP, performing steps
C             5.8.3c, d,
C             e and f
C             of the algorithm.
C             CDEF
C     PSHIFT..Auxiliary subroutine to OUTP and CDEF. It corrects reduced
C             amplitudes with regard to phase shift due to caustics, see
C             5.8.3f.
C             PSHIFT
C
C=======================================================================
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 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'):
      INCLUDE 'dcrt.inc'
C     dcrt.inc
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     ------------------------------------------------------------------
      INCLUDE 'raycb.inc'
C     raycb.inc
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             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: 2003, May 5
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations for local model parameters:  FAUX(10),
C       G(12),GAMMA(18),GSQRD,  UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL:
      INCLUDE 'auxmod.inc'
C     auxmod.inc
C
C.......................................................................
C
C     External procedure names:
      EXTERNAL FCT,FCTA,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     Anisotropic ray tracing:
      REAL CRTANI
      SAVE CRTANI
      DATA CRTANI/-999./
      IF(CRTANI.EQ.-999.) THEN
        CALL RSEP3R('CRTANI',CRTANI,0.)
        IF(CRTANI.NE.0..AND.KOOR().NE.0) THEN
C         5A1
          CALL ERROR('5A1 in RAYCB: Anisotropy in curvilinear coord.')
C         Anisotropic ray tracing is now applicable only to Cartesian
C         coordinates.
        END IF
      END IF
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:
      IF(CRTANI.EQ.0.) THEN
        CALL HPCG(PRMT,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX)
      ELSE
        CALL HPCG(PRMT,Y,DERY,NDIM,IHLF,FCTA,OUTP,AUX)
      END IF
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
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 C.R.T.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'):
      INCLUDE 'dcrt.inc'
C     dcrt.inc
C None of the storage locations of the common block are altered.
C
C Common block /RAYC/ - communication between RAYCB, FCT, and OUTP:
      INCLUDE 'raycb.inc'
C     raycb.inc
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: 2003, May 5
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations for local model parameters:  FAUX(10),
C       G(12),GAMMA(18),GSQRD,  UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL:
      INCLUDE 'auxmod.inc'
C     auxmod.inc
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.)
        DELT13=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.)
        DELT13=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
C     
C
      SUBROUTINE FCTA(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 C.R.T.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'):
      INCLUDE 'dcrt.inc'
C     dcrt.inc
C None of the storage locations of the common block are altered.
C
C Common block /RAYC/ - communication between RAYCB, FCT(A) and OUTP:
      INCLUDE 'raycb.inc'
C     raycb.inc
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: 2004, February 25
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
C
      REAL A(10,21),Q(21)
      REAL H11,H21,H31,H12,H22,H32
      REAL H41,H51,H61,H42,H52,H62,H43,H53,H63,D11,D13,D33
      REAL HHP,HHS,HH,HH1,HH2,HH3,HH4,HH5,HH6
      REAL HH11,HH12,HH22,HH13,HH23,HH33,HH14,HH24,HH34,HH44
      REAL HH15,HH25,HH35,HH45,HH55,HH16,HH26,HH36,HH46,HH56,HH66
      REAL V1,V11,V12,V22,V14,V24,V15,V25,V44,V45,V55
      REAL VNEXPS,AUX1,AUX2,AUX3
C
C     A,Q...  Elastic moduli and their spatial derivatives.
C     H11,H21,H31,H12,H22,H32... Components of the contravariant basis
C             vectors of ray-centred coordinate system.
C     H41,H51,H61,H42,H52,H62,H43,H53,H63... Components of the covariant
C             basis vectors of ray-centred coordinate system.
C     D11,D13,D33... Norms and scalar product of contravariant basis
C             vector H1 and covariant basis vector H3.
C     HHP,HHS,HH,HH1,HH2,HH3,HH4,HH5,HH6,HH11,HH12,HH22,HH13,HH23,HH33,
C     HH14,HH24,HH34,HH44,HH15,HH25,HH35,HH45,HH55,HH16,HH26,HH36,HH46,
c     HH56,HH66... Hamiltonian and its phase-space derivatives.
C     V1,V11,V12,V22,V14,V24,V15,V25,V44,V45,V55... Derivatives of the
C             Hamiltonian in ray-centred coordinate system.
C     VNEXPS..Ray velocity**(-NEXPS).
C     AUX1,AUX2,AUX3... Auxiliary storage locations.
C
C.......................................................................
C
C     (a) Number of invocations of FCT(A):
      IYRC(9)=IYRC(9)+1
C
C     (b) Material parameters:
      CALL PARM3(IABS(IYRC(5)),Y(3),A,YLRC(3),Q)
      CALL HDER(IYRC(5),A,Y(6),Y(7),Y(8),
     *          HHP,HHS,HH,HH1,HH2,HH3,HH4,HH5,HH6,
     *          HH11,HH12,HH22,HH13,HH23,HH33,HH14,HH24,HH34,HH44,
     *          HH15,HH25,HH35,HH45,HH55,HH16,HH26,HH36,HH46,HH56,HH66)
C
C     (c) Basis of the ray-centred coordinate system:
C     (c1) Non-unit vectors:
      H43=Y(6)
      H53=Y(7)
      H63=Y(8)
      H11=Y(9)
      H21=Y(10)
      H31=Y(11)
C     (c3) Norms:
      D33=SQRT(H43*H43+H53*H53+H63*H63)
      D11=SQRT(H11*H11+H21*H21+H31*H31)
      D13=(H11*H43+H21*H53+H31*H63)/(D11*D33)
C     (c4) Orthonormal vectors:
      H43=H43/D33
      H53=H53/D33
      H63=H63/D33
      H11=H11/D11-H43*D13
      H21=H21/D11-H53*D13
      H31=H31/D11-H63*D13
      H12=H53*H31-H63*H21
      H22=H63*H11-H43*H31
      H32=H43*H21-H53*H11
      H41=H22*HH6-H32*HH5
      H51=H32*HH4-H12*HH6
      H61=H12*HH5-H22*HH4
      AUX1=H11*H41+H21*H51+H31*H61
      H41=H41/AUX1
      H51=H51/AUX1
      H61=H61/AUX1
      H42=HH5*H31-HH6*H21
      H52=HH6*H11-HH4*H31
      H62=HH4*H21-HH5*H11
      AUX1=H12*H42+H22*H52+H32*H62
      H42=H42/AUX1
      H52=H52/AUX1
      H62=H62/AUX1
C
C     (d) Quantities useful to test the accuracy of computations:
      DELT11=ABS(D11-1.)
      DELT33=ABS(SQRT(HH)-1.)
      DELT13=ABS(D13)
C
      YLRC(1)=SQRT(HHP)/D33
      YLRC(2)=SQRT(HHS)/D33
      D33=D33/SQRT(HH)
      YLRC(4)=HH1/D33
      YLRC(5)=HH2/D33
      YLRC(6)=HH3/D33
C
C     (e) First and second derivatives of the Hamiltonian in ray-centred
C     coordinate system:
      V1=HH1*H11+HH2*H21+HH3*H31
      V2=HH1*H12+HH2*H22+HH3*H32
      AUX1=HH11*H11+HH12*H21+HH13*H31
      AUX2=HH12*H11+HH22*H21+HH23*H31
      AUX3=HH13*H11+HH23*H21+HH33*H31
      V11=H11*AUX1+H21*AUX2+H31*AUX3-V1*V1
      AUX1=HH11*H12+HH12*H22+HH13*H32
      AUX2=HH12*H12+HH22*H22+HH23*H32
      AUX3=HH13*H12+HH23*H22+HH33*H32
      V12=H11*AUX1+H21*AUX2+H31*AUX3-V1*V2
      V22=H12*AUX1+H22*AUX2+H32*AUX3-V2*V2
      AUX1=HH14*H41+HH15*H51+HH16*H61
      AUX2=HH24*H41+HH25*H51+HH26*H61
      AUX3=HH34*H41+HH35*H51+HH36*H61
      V14=H11*AUX1+H21*AUX2+H31*AUX3
      V24=H12*AUX1+H22*AUX2+H32*AUX3
      AUX1=HH14*H42+HH15*H52+HH16*H62
      AUX2=HH24*H42+HH25*H52+HH26*H62
      AUX3=HH34*H42+HH35*H52+HH36*H62
      V15=H11*AUX1+H21*AUX2+H31*AUX3
      V25=H12*AUX1+H22*AUX2+H32*AUX3
      AUX1=HH44*H41+HH45*H51+HH46*H61
      AUX2=HH45*H41+HH55*H51+HH56*H61
      AUX3=HH46*H41+HH56*H51+HH66*H61
      V44=H41*AUX1+H51*AUX2+H61*AUX3
      AUX1=HH44*H42+HH45*H52+HH46*H62
      AUX2=HH45*H42+HH55*H52+HH56*H62
      AUX3=HH46*H42+HH56*H52+HH66*H62
      V45=H41*AUX1+H51*AUX2+H61*AUX3
      V55=H42*AUX1+H52*AUX2+H62*AUX3
      AUX1=(H41*H43+H51*H53+H61*H63)/D33
      AUX2=(H42*H43+H52*H53+H62*H63)/D33
      V14=V14-V1*AUX1
      V24=V24-V2*AUX1
      V15=V15-V1*AUX2
      V25=V25-V2*AUX2
C
C     (f) Correction of the computed quantities:
      Y(6)=Y(6)/SQRT(HH)
      Y(7)=Y(7)/SQRT(HH)
      Y(8)=Y(8)/SQRT(HH)
      Y(9) =H11
      Y(10)=H21
      Y(11)=H31
C
C     (g) Right-hand sides of the differential equations:
      IF(NEXPS.EQ.0) THEN
        VNEXPS=1.
      ELSE
        VNEXPS=SQRT(HH4**2+HH5**2+HH6**2)**(-NEXPS)
        V11=V11*VNEXPS
        V12=V12*VNEXPS
        V22=V22*VNEXPS
        V14=V14*VNEXPS
        V24=V24*VNEXPS
        V15=V15*VNEXPS
        V25=V25*VNEXPS
        V44=V44*VNEXPS
        V45=V45*VNEXPS
        V55=V55*VNEXPS
      END IF
      D(1)=VNEXPS
      D(2)=0.
      D(3)= HH4*VNEXPS
      D(4)= HH5*VNEXPS
      D(5)= HH6*VNEXPS
      D(6)=-HH1*VNEXPS
      D(7)=-HH2*VNEXPS
      D(8)=-HH3*VNEXPS
      AUX1=V1/D33*VNEXPS
      D(9)= H43*AUX1
      D(10)=H53*AUX1
      D(11)=H63*AUX1
      D(12)= V14*Y(12)+V24*Y(13)+V44*Y(14)+V45*Y(15)
      D(13)= V15*Y(12)+V25*Y(13)+V45*Y(14)+V55*Y(15)
      D(14)=-V11*Y(12)-V12*Y(13)-V14*Y(14)-V15*Y(15)
      D(15)=-V12*Y(12)-V22*Y(13)-V24*Y(14)-V25*Y(15)
      D(16)= V14*Y(16)+V24*Y(17)+V44*Y(18)+V45*Y(19)
      D(17)= V15*Y(16)+V25*Y(17)+V45*Y(18)+V55*Y(19)
      D(18)=-V11*Y(16)-V12*Y(17)-V14*Y(18)-V15*Y(19)
      D(19)=-V12*Y(16)-V22*Y(17)-V24*Y(18)-V25*Y(19)
      D(20)= V14*Y(20)+V24*Y(21)+V44*Y(22)+V45*Y(23)
      D(21)= V15*Y(20)+V25*Y(21)+V45*Y(22)+V55*Y(23)
      D(22)=-V11*Y(20)-V12*Y(21)-V14*Y(22)-V15*Y(23)
      D(23)=-V12*Y(20)-V22*Y(21)-V24*Y(22)-V25*Y(23)
      D(24)= V14*Y(24)+V24*Y(25)+V44*Y(26)+V45*Y(27)
      D(25)= V15*Y(24)+V25*Y(25)+V45*Y(26)+V55*Y(27)
      D(26)=-V11*Y(24)-V12*Y(25)-V14*Y(26)-V15*Y(27)
      D(27)=-V12*Y(24)-V22*Y(25)-V24*Y(26)-V25*Y(27)
C
      RETURN
      END
C
C=======================================================================
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'):
      INCLUDE 'dcrt.inc'
C     dcrt.inc
C None of the storage locations of the common block are altered.
C
C Common block /INITC/ (see subroutine file 'init.for'):
      INCLUDE 'initc.inc'
C     initc.inc
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:
      INCLUDE 'raycb.inc'
C     raycb.inc
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
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: 2004, May 5
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations for local model parameters:  FAUX(10),
C       G(12),GAMMA(18),GSQRD,  UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL:
      INCLUDE 'auxmod.inc'
C     auxmod.inc
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,IE
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       Beginning 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.
      IF(UEBMUL.LE.0.) THEN
        ERR=AMAX1(ABS(Y(3)),ABS(Y(4)),ABS(Y(5)))/1000000.
        ERR=ERR/SQRT(D(3)**2+D(4)**2+D(5)**2)
        ERR=AMAX1(ERR,YYRC(2)/D(1))
      ELSE
        ERR=UEBMUL*YYRC(2)/D(1)
      END IF
      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 75 I=1,NSTOR
          IF(KSTOR(I).GT.NSRFC().AND.KSTOR(I).LE.100) THEN
            DO 73 IE=1,NEND
              IF(KSTOR(I).EQ.KEND(IE)) THEN
                GO TO 74
              END IF
   73       CONTINUE
            CALL SRFC2(KSTOR(I),Y(3),FAUX)
            IF(FAUX(1)*FSRFCA(KSTOR(I)-NSRFC()).LT.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
   74     CONTINUE
   75   CONTINUE
        IF(J.NE.0) THEN
          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)
          FSRFCA(KSTOR(J)-NSRFC())=AUX
          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.
C>    DO 91 K=1,NSTOR
C>      IF(KSTOR(K).GT.NSRFC().AND.KSTOR(K).LE.100) THEN
C>        IF(KSTOR(K).EQ.IYRC(6)) THEN
C>          CALL RPAR32(K,YLRC,Y1,YYRC,IYRC)
C>          CALL WRIT32(K,YLRC,Y1,YYRC,IYRC,IYRC(1)-27,Y1)
C>        END IF
C>      END IF
C> 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     Check for very small velocity:
      IF(IYRC(6).EQ.0) THEN
        IF(1.0E-12*(Y (6)**2+Y (7)**2+Y (8)**2).GT.
     *             (YI(6)**2+YI(7)**2+YI(8)**2)) THEN
          IYRC(6)=108
        END IF
      END IF
C
C     (k) Final operations:
C     YYRC(1)=X  is set by the subroutine CDEF
      IF(IYRC(6).NE.0) THEN
        PRMT(5)=1.
      END IF
      RETURN
      END
C
C=======================================================================
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'):
      INCLUDE 'dcrt.inc'
C     dcrt.inc
C None of the storage locations of the common block are altered.
C
C Common block /INITC/ (see subroutine file 'init.for'):
      INCLUDE 'initc.inc'
C     initc.inc
C None of the storage locations of the common block are altered.
C
C Common block /RAYC/ - communication between RAYCB, FCT, and OUTP:
      INCLUDE 'raycb.inc'
C     raycb.inc
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:  FAUX(10),
C       G(12),GAMMA(18),GSQRD,  UP(10),US(10),RO,QP,QS, VP,VS,VD(10),QL:
      INCLUDE 'auxmod.inc'
C     auxmod.inc
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     Quantities 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
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, see
C C.R.T.5.8.3f.
C Hereinafter, the point 1 of the ray denotes the point corresponding
C to the previous invocation of the subroutine PSHIFT, the point 2
C denotes the point corresponding to this invocation of the subroutine
C PSHIFT.
C
C Input:
C     NY...   Number of basic quantities describing the point of a ray.
C             NY=27 at the initial point of a ray element, where no
C             phase shift is applied to the amplitudes.
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             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 Subroutines and external functions called:
      EXTERNAL RPAR30
C
C Date: 1999, May 25
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER I
      REAL P211,P221,P212,P222,AUX,SMALL
      PARAMETER (SMALL=0.000000001)
      REAL Q111,Q121,Q112,Q122,Q1DET
      REAL Q211,Q221,Q212,Q222,Q2DET
      SAVE Q211,Q221,Q212,Q222,Q2DET
C
C     I,P211,P221,P212,P222,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
      IF(NY.GE.28) THEN
        Q111=Q211
        Q121=Q221
        Q112=Q212
        Q122=Q222
        Q1DET=Q2DET
      END IF
      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)
      CALL RPAR30(YI(1),Y(1),Q211,Q221,Q212,Q222)
C     Regularization in order to avoid floating-point overflow:
      AUX=SQRT(Q211*Q211+Q221*Q221+Q212*Q212+Q222*Q222+1.)
      Q211=Q211/AUX
      Q221=Q221/AUX
      Q212=Q212/AUX
      Q222=Q222/AUX
      Q2DET=Q211*Q222-Q212*Q221
      IF(Y(1).EQ.YI(1)) THEN
        IF(Q2DET.EQ.0.) THEN
          P211=Y(14)*YI(12)+Y(18)*YI(13)+Y(22)*YI(14)+Y(26)*YI(15)
          P221=Y(15)*YI(12)+Y(19)*YI(13)+Y(23)*YI(14)+Y(27)*YI(15)
          P212=Y(14)*YI(16)+Y(18)*YI(17)+Y(22)*YI(18)+Y(26)*YI(19)
          P222=Y(15)*YI(16)+Y(19)*YI(17)+Y(23)*YI(18)+Y(27)*YI(19)
          Q211=Q211+P211*SMALL
          Q221=Q221+P221*SMALL
          Q212=Q212+P212*SMALL
          Q222=Q222+P222*SMALL
          Q2DET=Q211*Q222-Q212*Q221
        END IF
      END IF
      IF(NY.GE.28) THEN
        IF(Q1DET*Q2DET.LT.0.) THEN
C         Caustic of the first order (line caustic) between points
          ISHIFT=1
          DO 10 I=28,NY,2
            AUX=Y(I+1)
            Y(I+1)=-Y(I)
            Y(I)=AUX
   10     CONTINUE
        ELSE
          AUX=(Q111*Q222-Q112*Q221+Q122*Q211-Q121*Q212)*Q2DET
          IF(AUX.LT.0.) THEN
C           Caustic of the second order (point caustic) between points
            ISHIFT=2
            DO 20 I=28,NY
              Y(I)=-Y(I)
   20       CONTINUE
          ELSE IF(AUX.EQ.0.) THEN
C           Stepping on a caustic
            ISHIFT=1
            DO 30 I=28,NY,2
              AUX=Y(I+1)
              Y(I+1)=-Y(I)
              Y(I)=AUX
   30       CONTINUE
          ELSE
C           No caustic
            ISHIFT=0
          END IF
        END IF
      END IF
      RETURN
      END
C
C=======================================================================
C