C 
      SUBROUTINE CRTOUT
     *            (LU1,LU2,KALL,KREC,INI,NQ,MPTS,NPTS,OUT,OUTMIN,OUTMAX)
C
      INTEGER LU1,LU2,KALL,KREC,INI,NQ,MPTS,NPTS
      REAL OUT(NQ,MPTS),OUTMIN(NQ),OUTMAX(NQ)
C
C Subroutine designed to prepare some output quantities of the CRT
C program for printing.
C
C Input:
C     LU1,LU2... Logical unit numbers corresponding to files with ray
C             points and ray initial points.
C     KALL... KALL.LE.0: only two-point rays are considered,
C             KALL.GE.1: all rays are considered.
C     KREC... 0: No Taylor expansion of travel time.
C             1: Linear Taylor expansion of travel time to the receiver.
C             2: Quadratic Taylor expansion of travel time and linear
C                Taylor expansion of slowness vector to the receiver.
C     INI...  INI.LE.0: initial points of rays are not considered,
C             INI.GE.1: initial points of rays are considered.
C     NQ...   Number of reals in each output line.
C     MPTS... Maximum total number of ray points.
C     NPTS... Number of ray points already stored in array OUT.
C     OUT...  Quantities at ray points already stored in the memory
C             during previous invocations.
C     OUTMIN,OUTMAX... Minimum and maximum values of corresponding
C             quantities stored in array OUT.
C
C Output:
C     NPTS... Number of ray points stored in array OUT.
C             Input value increased by 1 or 2 (if also the initial
C             point of a ray has been stored).
C     OUT...  For each ray point, the first NQ quantities of the
C             following ones:
C             1. X1-coordinate.
C             2. X2-coordinate.
C             3. X3-coordinate.
C             4. Travel time.
C             5. P1 slowness-vector component.
C             6. P2 slowness-vector component.
C             7. P3 slowness-vector component.
C             8. Real part of complex-valued amplitude
C                A11 or A13 or A31 or A33, depending on a P or S wave.
C                The first index of amplitude Aij corresponds to the
C                endpoint of the ray, the second index corresponds to
C                initial point of the ray.  Indices i,j=1,2 correspond
C                to two S-wave polarizations, indices i,j=3 correspond
C                to a P-wave.
C                Ray amplitudes A11, A22 or A33 are normalized to 1 at
C                an initial surface or on a unit sphere around a point
C                source.
C             9. Imaginary part of complex-valued amplitude
C                A11 or A31 or A13 or A33, depending on a P or S wave.
C            10. Real part of complex-valued amplitude
C                A21 or A32 or A23, or 0, depending on a P or S wave.
C            11. Imaginary part of complex-valued amplitude
C                A21 or A32 or A23, or 0, depending on a P or S wave.
C            12. Real part of complex-valued amplitude A12, or 0,
C                depending on a P or S wave.
C            13. Imaginary part of complex-valued amplitude A12, or 0,
C                depending on a P or S wave.
C            14. Real part of complex-valued amplitude A22, or 0,
C                depending on a P or S wave.
C            15. Imaginary part of complex-valued amplitude A22, or 0,
C                depending on a P or S wave.
C     OUTMIN,OUTMAX... Minimum and maximum values of corresponding
C             quantities stored in array OUT.
C
C SEP parameter:
C     KSRC=integer... Modification of the initial point of the ray.
C             0: No modification of the initial point of the ray.
C             1: If the source file is specified, the coordinates of
C                the initial point of a ray are replaced by the source
C                coordinates and the travel time is linearly
C                interpolated from the initial point to the source
C                point.
C             Default: KSRC=0
C     KTWO=integer... Converting initial-value to two-point travel time.
C             0: No modification of the initial-value travel time.
C             1: The initial-value travel time is replaced by the
C                two-point travel time.
C                The two-point travel time is the initial-value travel
C                time minus the travel time at the initial point of the
C                ray.
C             Default: KTWO=0
C
C Common block /POINTC/ to store the results of complete ray tracing:
      INCLUDE 'pointc.inc'
C     pointc.inc
C None of the storage locations of the common block are altered.
C
C Subroutines and external functions required:
      EXTERNAL AP00
C     AP00... File 'ap.for'.
C
C Date: 2004, June 10
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
C     Local storage locations:
      REAL GIANT
      PARAMETER (GIANT=1.E20)
      INTEGER IQ
      REAL QDETI,QDET,VI,V,RHOI,RHO,AUX
      REAL HI(18),H(18),HUI(9),RM(6),RN(6),R(3),P(3)
      REAL TTINI
      SAVE TTINI
      INTEGER KSRC,KTWO,KPHASE
      SAVE    KSRC,KTWO,KPHASE
      DATA    KSRC/-100/
C
C.......................................................................
C
      IF(KSRC.EQ.-100) THEN
        CALL RSEP3I('KSRC',KSRC,0)
        CALL RSEP3I('KTWO',KTWO,0)
        CALL RSEP3I('KPHASE',KPHASE,0)
      END IF
C
      IF(NPTS.EQ.0) THEN
        DO 10 IQ=1,NQ
          OUTMIN(IQ)= GIANT
          OUTMAX(IQ)=-GIANT
   10   CONTINUE
      END IF
C
   20 CONTINUE
C       Reading the results of the complete ray tracing
        CALL AP00(0,LU1,LU2)
        IF(IWAVE.LT.1)THEN
C         End of rays
          RETURN
        ELSE IF (KALL.GT.0.OR.IREC.GT.0) THEN
C         Initial travel time
          IF(IPT.LE.1)THEN
            TTINI=YI(1)
            IF(KSRC.GE.1) THEN
              R(1)=YI(3)
              R(2)=YI(4)
              R(3)=YI(5)
C             The first point in the source file is selected:
              CALL SRC(1,R(1),R(2),R(3))
C             Linear Taylor expansion of travel time:
              R(1)=R(1)-YI(3)
              R(2)=R(2)-YI(4)
              R(3)=R(3)-YI(5)
              TTINI=TTINI+YI(6)*R(1)+YI(7)*R(2)+YI(8)*R(3)
            END IF
          END IF
C         ..............................................................
C         Initial point of the ray:
          IF(INI.GT.0.AND.IPT.LE.1)THEN
C           New ray - recording the initial point
            IF(NPTS.GE.MPTS) THEN
C             CRTOUT-01
              CALL ERROR
     *               ('CRTOUT-01: Too many ray points to fit in memory')
            END IF
            NPTS=NPTS+1
C
C           Coordinates:
            R(1)=YI(3)
            R(2)=YI(4)
            R(3)=YI(5)
            IF(KSRC.GE.1) THEN
C             The first point in the source file is selected:
              CALL SRC(1,R(1),R(2),R(3))
            END IF
            DO 31 IQ=1,MIN0(NQ,3)
              OUT(IQ,NPTS)=R(IQ)
   31       CONTINUE
C           Travel time:
            DO 32 IQ=4,MIN0(NQ,4)
              OUT(IQ,NPTS)=TTINI
              IF(KTWO.GE.1) THEN
C               Converting initial-value to two-point travel time
                OUT(IQ,NPTS)=OUT(IQ,NPTS)-TTINI
              END IF
   32       CONTINUE
C           Slowness vector:
            DO 33 IQ=5,MIN0(NQ,7)
              OUT(IQ,NPTS)=YI(1+IQ)
   33       CONTINUE
C           Amplitudes:
            DO 34 IQ=8,NQ
              OUT(IQ,NPTS)=0.
   34       CONTINUE
            DO 35 IQ=8,MIN0(NQ,NY-20),6
              OUT(IQ,NPTS)=1.
   35       CONTINUE
C
            DO 39 IQ=1,NQ
              OUTMIN(IQ)=AMIN1(OUTMIN(IQ),OUT(IQ,NPTS))
              OUTMAX(IQ)=AMAX1(OUTMAX(IQ),OUT(IQ,NPTS))
   39       CONTINUE
          END IF
C         ..............................................................
C         Other than initial ray point:
          IF(NPTS.GE.MPTS) THEN
C           CRTOUT-02
            CALL ERROR
     *               ('CRTOUT-02: Too many ray points to fit in memory')
          END IF
          NPTS=NPTS+1
C
C         Coordinates:
          R(1)=Y(3)
          R(2)=Y(4)
          R(3)=Y(5)
          IF(KREC.GE.1) THEN
            CALL REC(IREC,R(1),R(2),R(3))
          END IF
          DO 41 IQ=1,MIN0(NQ,3)
            OUT(IQ,NPTS)=R(IQ)
   41     CONTINUE
C         Travel time:
          DO 42 IQ=4,MIN0(NQ,4)
            OUT(IQ,NPTS)=Y(1)
            IF(KREC.GE.1) THEN
C             Linear Taylor expansion of travel time:
              R(1)=R(1)-Y(3)
              R(2)=R(2)-Y(4)
              R(3)=R(3)-Y(5)
              OUT(IQ,NPTS)=OUT(IQ,NPTS)+Y(6)*R(1)+Y(7)*R(2)+Y(8)*R(3)
              IF(KREC.GE.2) THEN
C               Quadratic Taylor expansion of travel time:
                CALL AP03(0,HI,H,HUI)
                CALL AP08(0,H,HUI,RM,RN)
                P(1)=RN(1)*R(1)+RN(2)*R(2)+RN(4)*R(3)
                P(2)=RN(2)*R(1)+RN(3)*R(2)+RN(5)*R(3)
                P(3)=RN(4)*R(1)+RN(5)*R(2)+RN(6)*R(3)
                OUT(IQ,NPTS)=OUT(IQ,NPTS)
     *                              +0.5*(P(1)*R(1)+P(2)*R(2)+P(3)*R(3))
              END IF
            END IF
            IF(KTWO.GE.1) THEN
C             Converting initial-value to two-point travel time
              OUT(IQ,NPTS)=OUT(IQ,NPTS)-TTINI
            END IF
   42     CONTINUE
C         Slowness vector:
          DO 43 IQ=5,MIN0(NQ,7)
            OUT(IQ,NPTS)=Y(1+IQ)
            IF(KREC.GE.2) THEN
C             Linear Taylor expansion of slowness vector:
              IF(KREC.LT.2) THEN
                OUT(IQ,NPTS)=Y(1+IQ)
              ELSE
                OUT(IQ,NPTS)=Y(1+IQ)+P(IQ-4)
              END IF
            END IF
   43     CONTINUE
C         Amplitudes:
          IF(NQ.GT.7) THEN
C wrong!    AUX=SQRT( YLI(1)**3*YLI(3)*ABS(YI(14)*YI(19)-YI(15)*YI(18))/
C old       AUX=SQRT( YLI(1)**3*YLI(3)/
C old*                (YL(1)   *YL(3) *ABS( Y(20)*Y(25) - Y(21)*Y(24))))
            CALL AP07(QDETI,QDET,VI,V,RHOI,RHO,IQ)
            IF(QDET.NE.0.) THEN
              AUX=SQRT(QDETI*VI*RHOI/(QDET*V*RHO))
            ELSE
              AUX=0.
            END IF
            DO 44 IQ=8,MIN0(NQ,NY-20)
              OUT(IQ,NPTS)=Y(20+IQ)*AUX
   44       CONTINUE
            DO 45 IQ=NY-19,NQ
              OUT(IQ,NPTS)=0.
   45       CONTINUE
            DO 46 IQ=9,NQ,2
C             Modifying amplitudes:
              AUX=SQRT(OUT(IQ-1,NPTS)**2+OUT(IQ,NPTS)**2)
              IF(KPHASE.EQ.0) THEN
C               Real and imaginary parts of the complex-valued amplitude
                IF(ABS(OUT(IQ-1,NPTS)).LT.0.000001*AUX) THEN
                  OUT(IQ-1,NPTS)=0.
                END IF
                IF(ABS(OUT(IQ,NPTS)).LT.0.000001*AUX) THEN
                  OUT(IQ,NPTS)=0.
                END IF
              ELSE
C               Absolute value and phase of the complex-valued amplitude
                IF(AUX.GT.0.) THEN
                  OUT(IQ,NPTS)=ATAN2(OUT(IQ,NPTS),OUT(IQ-1,NPTS))
                  OUT(IQ-1,NPTS)=AUX
                END IF
                IF(ABS(OUT(IQ,NPTS)).LT.0.000001) THEN
                  OUT(IQ,NPTS)=0.
                END IF
              END IF
   46       CONTINUE
          END IF
C
          DO 49 IQ=1,NQ
            OUTMIN(IQ)=AMIN1(OUTMIN(IQ),OUT(IQ,NPTS))
            OUTMAX(IQ)=AMAX1(OUTMAX(IQ),OUT(IQ,NPTS))
   49     CONTINUE
C         ..............................................................
          RETURN
        END IF
      GO TO 20
C
      END
C
C=======================================================================
C
C     
C
      SUBROUTINE TXT1(LU,SRCFIL,RECFIL)
C
C Subroutine designed to read the source and receiver names and prepare
C them for entry TXT2.
C
C     ENTRY      TXT2(KALL,KTT,ISRC,IWAVE,IRAY,IREC,LENTXT,TXT)
C
C Entry designed to generate the string describing the ray or its
C endpoint.
C
C     ENTRY REC(IREC,R1,R2,R3)
C
C     ENTRY SRC(IREC,R1,R2,R3)
C
      INTEGER LU,KALL,KTT,ISRC,IWAVE,IRAY,IREC,LENTXT
      CHARACTER*(*) SRCFIL,RECFIL,TXT
      REAL R1,R2,R3
C
C Input:
C     KALL... IABS(KALL)=0: only source (if the source file has been
C               submitted) and receiver (for two-point rays) information
C               is contained within the output string.
C             IABS(KALL)=1: in addition, the string is prefixed with
C               the index of the ray.
C             IABS(KALL)=2: in addition, the string is prefixed also
C               with the index of the elementary wave.
C     KTT...  0: creates a single string.
C             1: separates the string into 2 strings: the source and
C                receiver parts.  This option is intended to generate
C                files with synthetic travel times.
C     ISRC... Index of the source.
C     IWAVE...Index of the elementary wave.
C     IRAY... Index of the ray within the elementary wave.
C     IREC... Index of the receiver for a two-point ray, determined in
C             subroutine RPAR4.
C
C Output:
C     LENTXT..Length of the string, including also apostrophes.
C     TXT...  The string, beginning and terminating with apostrophes.
C             If KTT=1, the string is separated by
C             apostrophe,blank,apostrophe into the source and receiver
C             parts.
C             Examples:                                  KTT  IREC   REC
C             ---------                                    KALL   SRC
C             ' '                                         0  0  0  n
C             'REC   13'                                  0  0  +  n  n
C             'recnam'                                    0  0  +  n  y
C             'srcnam'                                    0  0  0  y
C             'srcnam, REC   13'                          0  0  +  y  n
C             'srcnam TO recnam'                          0  0  +  y  y
C             'RAY  112'                                  0  1  0  n
C             'RAY  112, REC   13'                        0  1  +  n  n
C             'RAY  112 TO recnam'                        0  1  +  n  y
C             'RAY  112 FROM srcnam'                      0  1  0  y
C             'RAY  112 FROM srcnam, REC   13'            0  1  +  y  n
C             'RAY  112 FROM srcnam to recnam'            0  1  +  y  y
C             'WAVE   1, RAY  112'                        0  2  0  n
C             'WAVE   1, RAY  112, REC   13'              0  2  +  n  n
C             'WAVE   1, RAY  112 TO recnam'              0  2  +  n  y
C             'WAVE   1, RAY  112 FROM srcnam'            0  2  0  y
C             'WAVE   1, RAY  112 FROM srcnam, rec   13'  0  2  +  y  n
C             'WAVE   1, RAY  112 FROM srcnam TO recnam'  0  2  +  y  y
C             ' ' ' '                                     1  0  0  n
C             ' ' 'REC   13'                              1  0  +  n  n
C             ' ' 'recnam'                                1  0  +  n  y
C             'srcnam' ' '                                1  0  0  y
C             'srcnam' 'REC   13'                         1  0  +  y  n
C             'srcnam' 'recnam'                           1  0  +  y  y
C             'RAY  112' ' '                              1  1  0  n
C             'RAY  112' 'REC   13'                       1  1  +  n  n
C             'RAY  112' 'recnam'                         1  1  +  n  y
C             'RAY  112 FROM srcnam' ' '                  1  1  0  y
C             'RAY  112 FROM srcnam' 'REC   13'           1  1  +  y  n
C             'RAY  112 FROM srcnam' 'recnam'             1  1  +  y  y
C             'WAVE   1, RAY  112' ' '                    1  2  0  n
C             'WAVE   1, RAY  112' 'REC   13'             1  2  +  n  n
C             'WAVE   1, RAY  112' 'recnam'               1  2  +  n  y
C             'WAVE   1, RAY  112 FROM srcnam' ' '        1  2  0  y
C             'WAVE   1, RAY  112 FROM srcnam' 'REC   13' 1  2  +  y  n
C             'WAVE   1, RAY  112 FROM srcnam' 'recnam'   1  2  +  y  y
C
C Subroutines and external functions required:
      INTEGER  LENGTH
      EXTERNAL LENGTH
C     LENGTH..File 'length.for'.
C
C Date: 2003, August 14
C Coded by Ludek Klimes
C
C-----------------------------------------------------------------------
C
      INTEGER MPTS,NSRC,NREC,LENSRC,LENREC
      PARAMETER (MPTS=1000)
      CHARACTER*8 POINTS(MPTS)
      REAL UNDEF,COOR1(MPTS),COOR2(MPTS),COOR3(MPTS)
      PARAMETER (UNDEF=-999999.)
      SAVE NSRC,NREC,LENSRC,LENREC,POINTS,COOR1,COOR2,COOR3
C
      INTEGER I
C
C-----------------------------------------------------------------------
C
      NSRC=0
      LENSRC=0
      IF(SRCFIL.NE.' ') THEN
        OPEN(LU,FILE=SRCFIL,STATUS='OLD')
        READ(LU,*) (POINTS(1),I=1,20)
   11   CONTINUE
          IF(NSRC.GE.MPTS) THEN
C           CRTOUT-03
            CALL ERROR
     *             ('CRTOUT-03: Too many source names to fit in memory')
          END IF
          NSRC=NSRC+1
          POINTS(NSRC)='$'
          COOR1(NSRC)=UNDEF
          COOR2(NSRC)=UNDEF
          COOR3(NSRC)=UNDEF
          READ(LU,*,END=12)
     *                  POINTS(NSRC),COOR1(NSRC),COOR2(NSRC),COOR3(NSRC)
          IF(POINTS(NSRC).EQ.'$') THEN
            GO TO 12
          END IF
          LENSRC=MAX0(LENGTH(POINTS(NSRC)),LENSRC)
        GO TO 11
   12   CONTINUE
        NSRC=NSRC-1
        CLOSE(LU)
      END IF
C
      NREC=NSRC
      LENREC=0
      IF(RECFIL.NE.' ') THEN
        OPEN(LU,FILE=RECFIL,STATUS='OLD')
        READ(LU,*) (POINTS(NREC+1),I=1,20)
   21   CONTINUE
          IF(NREC.GE.MPTS) THEN
C           CRTOUT-04
            CALL ERROR
     *           ('CRTOUT-04: Too many receiver names to fit in memory')
          END IF
          NREC=NREC+1
          POINTS(NREC)='$'
          COOR1(NREC)=UNDEF
          COOR2(NREC)=UNDEF
          COOR3(NREC)=UNDEF
          READ(LU,*) POINTS(NREC),COOR1(NREC),COOR2(NREC),COOR3(NREC)
          IF(POINTS(NREC).EQ.'$') THEN
            NREC=NREC-1
            GO TO 22
          END IF
          LENREC=MAX0(LENGTH(POINTS(NREC)),LENREC)
        GO TO 21
   22   CONTINUE
        CLOSE(LU)
      END IF
C
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY REC(IREC,R1,R2,R3)
C
      IF(LENREC.GT.0.AND.IREC.GT.0.AND.IREC.LE.NREC-NSRC) THEN
        IF(COOR1(NSRC+IREC).NE.UNDEF) R1=COOR1(NSRC+IREC)
        IF(COOR2(NSRC+IREC).NE.UNDEF) R2=COOR2(NSRC+IREC)
        IF(COOR3(NSRC+IREC).NE.UNDEF) R3=COOR3(NSRC+IREC)
      END IF
C
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY SRC(IREC,R1,R2,R3)
C
      IF(LENSRC.GT.0.AND.IREC.GT.0.AND.IREC.LE.NSRC) THEN
        IF(COOR1(IREC).NE.UNDEF) R1=COOR1(IREC)
        IF(COOR2(IREC).NE.UNDEF) R2=COOR2(IREC)
        IF(COOR3(IREC).NE.UNDEF) R3=COOR3(IREC)
      END IF
C
      RETURN
C
C-----------------------------------------------------------------------
C
      ENTRY TXT2(KALL,KTT,ISRC,IWAVE,IRAY,IREC,LENTXT,TXT)
C
C     Initial apostrophe:
      LENTXT=1
      TXT=''''
C
C     Index of the elementary wave:
      IF(IABS(KALL).GE.2) THEN
        TXT(LENTXT+1:LENTXT+10)='WAVE0000, '
        WRITE(TXT(LENTXT+5:LENTXT+8),'(I4)') IWAVE
        LENTXT=LENTXT+10
      END IF
C
C     Index of the ray:
      IF(IABS(KALL).GE.1) THEN
        TXT(LENTXT+1:LENTXT+8)='RAY00000'
        WRITE(TXT(LENTXT+4:LENTXT+8),'(I5)') IRAY
        LENTXT=LENTXT+8
      END IF
C
C     Name of the source:
      IF(LENSRC.GT.0) THEN
        IF(ISRC.GT.0) THEN
          IF(ISRC.GT.NSRC) THEN
C           CRTOUT-05
            CALL ERROR
     *           ('CRTOUT-05: Source index exceeding number of sources')
          END IF
C         Separator:
          IF(KALL.NE.0) THEN
            TXT(LENTXT+1:LENTXT+6)=' FROM '
            LENTXT=LENTXT+6
          END IF
C         Name:
          TXT(LENTXT+1:LENTXT+LENSRC)=POINTS(ISRC)(1:LENSRC)
          LENTXT=LENTXT+LENSRC
        END IF
      END IF
C
C     Separation of source and receiver strings:
      IF(KTT.NE.0) THEN
        IF(LENTXT.LE.1) THEN
          LENTXT=2
        END IF
        TXT(LENTXT+1:LENTXT+3)=''' '''
        LENTXT=LENTXT+3
      END IF
C
C     Name of the receiver:
      IF(LENREC.GT.0) THEN
C       Separator:
        IF(KTT.EQ.0.AND.LENTXT.GT.1) THEN
          IF(IREC.GT.0) THEN
            TXT(LENTXT+1:LENTXT+4)=' TO '
          END IF
          LENTXT=LENTXT+4
        END IF
C       Name:
        IF(IREC.GT.0) THEN
          IF(IREC.GT.NREC-NSRC) THEN
C           CRTOUT-06
            CALL ERROR
     *       ('CRTOUT-06: Receiver index exceeding number of receivers')
          END IF
          TXT(LENTXT+1:LENTXT+LENREC)=POINTS(NSRC+IREC)(1:LENREC)
        END IF
        LENTXT=LENTXT+LENREC
      END IF
C
C     Index of the receiver:
      IF(LENREC.EQ.0) THEN
C       Separator:
        IF(KTT.EQ.0.AND.LENTXT.GT.1) THEN
          IF(IREC.GT.0) THEN
            TXT(LENTXT+1:LENTXT+2)=', '
          END IF
          LENTXT=LENTXT+2
        END IF
C       Index:
        IF(IREC.GT.0) THEN
          TXT(LENTXT+1:LENTXT+8)='REC00000'
          WRITE(TXT(LENTXT+4:LENTXT+8),'(I5)') IREC
        END IF
        LENTXT=LENTXT+8
      END IF
C
C     Terminating apostrophe:
      LENTXT=LENTXT+1
      TXT(LENTXT:LENTXT)=''''
C
      RETURN
      END
C
C=======================================================================
C