C
C Program CRTPTS converting the unformatted output of program CRT into a
C formatted file containing coordinates, travel times, slowness vectors,
C and amplitudes at the endpoints of (usually two-point) rays.
C
C Version: 5.80
C Date: 2004, June 10
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: klimes@seis.karlov.mff.cuni.cz
C
C.......................................................................
C
C This simple conversion program may serve as an example how to read and
C process output files of the complete ray tracing program 'CRT'.
C Reading the output files is completed by a simple invocation of
C subroutine AP00 of file 'ap.for', called by means of subroutine CRTOUT
C of file 'crtout.for'.
C
C The structure of the output file with rays is an extension of the
C general file form POINTS
C designed to store 3-D points.
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Names of the input and output files:
C     CRTOUT='string'... File with the names of the output files of
C             program CRT. A single set of CRT output files is read
C             from CRTOUT.  Only the files 'CRT-S', and 'CRT-I' are
C             read by CRTRAY.
C             For general description of file CRTOUT refer to
C             file writ.for.
C             Default: CRTOUT=' ' which means 'CRT-S'='s01.out',
C                      'CRT-I'='s01i.out'
C     REC='string' ... If non-blank, the name of the file with the
C             names of the receiver points.  The names are then used
C             within the strings describing the points of two-point
C             rays.  Otherwise, the two-point rays are denoted by the
C             receiver index.
C             Description of file REC
C             Default: REC=' '
C     SRC='string' ... If non-blank, the name of the file with the
C             name and coordinates of the source point.  The source name
C             may then be used within the strings describing the points
C             of the rays (see parameter KALL), and the coordinates may
C             be used to modify the initial point of a ray (see
C             parameter KSRC) and two-point travel time (see parameter
C             KTWO).   In this version, the first point of file SRC
C             is considered to be the source point.
C             Description of file SRC
C             Default: SRC=' '
C     PTS='string' .. Name of the output formatted file with ray points.
C             Description of file PTS
C             Default: PTS='pts.out'
C Parameters describing the form of output file PTS:
C     NQ=integer ... Number of reals in each output line of file PTS.
C             If NQ exceeds parameter MQ (see the code below), it is set
C             to MQ.
C             The output reals represent:
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 or absolute value of complex-valued amplitude
C                A11 or A13 or A31 or A33, depending on a P or S wave.
C                See parameter KPHASE below.
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             9. Imaginary part or phase of complex-valued amplitude
C                A11 or A31 or A13 or A33, depending on a P or S wave.
C                See parameter KPHASE below.
C            10. Real part or absolute value of complex-valued amplitude
C                A21 or A32 or A23, or 0, depending on a P or S wave.
C                See parameter KPHASE below.
C            11. Imaginary part or phase of complex-valued amplitude
C                A21 or A32 or A23, or 0, depending on a P or S wave.
C                See parameter KPHASE below.
C            12. Real part or absolute value of complex-valued amplitude
C                A12, or 0, depending on a P or S wave.
C                See parameter KPHASE below.
C            13. Imaginary part or phase of complex-valued amplitude
C                A12, or 0, depending on a P or S wave.
C                See parameter KPHASE below.
C            14. Real part or absolute value of complex-valued amplitude
C                A22, or 0, depending on a P or S wave.
C                See parameter KPHASE below.
C            15. Imaginary part or phase of complex-valued amplitude
C                A22, or 0, depending on a P or S wave.
C                See parameter KPHASE below.
C             Default: NQ=11
C     KPHASE=integer:
C             0: Real and imaginary parts of the complex-valued
C                amplitude are written to the output file PTS,
C                see the description of NQ.
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                Imaginary part of ray amplitude is printed
C                only if it is greater (in absolute value) than
C                0.000001 times the absolute value of the complex-valued
C                amplitude.
C             1: Absolute value and phase of the complex-valued
C                amplitude are written, instead of the real and
C                imaginary parts.
C             Default: KPHASE=0
C     KALL=integer:
C             KALL.LE.0: only two-point rays are considered.
C             KALL.GE.1: all rays are considered.
C             Absolute value specifies the form of the strings
C             describing individual points.  Here are some examples:
C             ABS(KALL)=0: 'REC   13'
C                          'recnam'
C                          'srcnam TO recnam'
C             ABS(KALL)=1: 'RAY  112'
C                          'RAY  112, REC   13'
C                          'RAY  112 TO recnam'
C                          'RAY  112 FROM srcnam'
C                          'RAY  112 FROM srcnam TO recnam'
C             ABS(KALL)=2: 'WAVE   1, RAY  112'
C                          'WAVE   1, RAY  112, REC   13'
C                          'WAVE   1, RAY  112 TO recnam'
C                          'WAVE   1, RAY  112 FROM srcnam'
C                          'WAVE   1, RAY  112 FROM srcnam TO recnam'
C             Values KALL=0 and KALL=1 specify the briefest strings.
C             Default: KALL=0
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 sourcer 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                extrapolated from the initial point to the source
C                point.
C                In this version, the first point of file SRC is
C                considered to be the source point.
C             Default: KSRC=0
C     KREC=integer... Modification of the surface point of the ray.
C            -1: Initial points of rays are written to the output file
C                instead of ray points situated on the storing surface.
C             0: Ray points situated on the storing surface are written
C                to the output file.
C             1: If the receiver file is specified, the coordinates
C                of ray endpoints are replaced by receiver coordinates
C                and the travel time is linearly extrapolated to the
C                receivers.
C             2: If the receiver file is specified, the coordinates
C                of ray endpoints are replaced by receiver coordinates
C                and the travel time is quadratically extrapolated to
C                the receivers.
C             Default: KREC=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                The travel time at the initial point of the ray depends
C                on parameter KSRC.
C             Default: KTWO=0
C     KTT=integer:
C             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.  This option yields
C                a meaningful output only for two-point rays, with
C                KALL=0 and KREC.GE.0.
C             Default: KTT=0
C
C                                                     
C Output formatted file PTS:
C (1) / (a slash).
C (2) For each ray endpoint (or initial point) (2.1):
C (2.1) 'RAYTXT',(OUT(I),I=1,NQ),/
C     'RAYTXT'... One or two strings in apostrophes describing the ray.
C             See the description of input parameters KALL and KTT.
C             One string:  output format PTS.
C             Two strings: output format FTT.
C     (OUT(I),I=1,NQ)... Output quantities at the ray point, see the
C             description of input parameter NQ.
C     /...    An obligatory slash after at the end of line, in place
C             where the slowness vector components could be written.
C     For default NQ=11 and P-wave at the source one of:
C     (2.1) 'RAYTXT',X1,X2,X3,TT,P1,P2,P3,AR,AI,/
C     (2.1) 'RAYTXT',X1,X2,X3,TT,P1,P2,P3,AR,/
C     X1,X2,X3... Coordinates of the point of the ray.
C     TT...   Arrival time at the point.
C     P1,P2,P3... Slowness vector.
C     AR...   Real part or absolute value of the complex-valued
C             amplitude, see parameter KPHASE.
C     AI...   Imaginary part or phase of the complex-valued amplitude,
C             see parameter KPHASE.
C     /...    An obligatory slash after at the end of line.
C (3) / (a slash).
C Description of format PTS
C Description of format FTT
C
C-----------------------------------------------------------------------
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 CRTOUT,TXT1,TXT2,TTSORT,FORM2
C     CRTOUT,TXT1,TXT2... File 'crtout.for'.
C     AP00... File 'ap.for' (called by CRTOUT).
C     LENGTH..File 'length.for' (called by CRTOUT).
C     TTSORT..File 'ttsort.for'.
C     INDEXX..File 'indexx.for' (called by TTSORT).
C     FORM2...File 'forms.for'.
C     FORM1...File 'forms.for' (called by FORM2).
C
C-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C     Allocation of working arrays:
      INTEGER MPTS,MOUT
      PARAMETER (MPTS=MRAM/9,MOUT=MRAM-5*MPTS)
      INTEGER IRECS(MPTS),IWAVES(MPTS),IRAYS(MPTS),INDX(MPTS)
      INTEGER ISRCS(MPTS)
      REAL OUT(MOUT)
      EQUIVALENCE (IRECS ,RAM          )
      EQUIVALENCE (IWAVES,RAM(  MPTS+1))
      EQUIVALENCE (IRAYS ,RAM(2*MPTS+1))
      EQUIVALENCE (INDX  ,RAM(3*MPTS+1))
      EQUIVALENCE (ISRCS ,RAM(4*MPTS+1))
      EQUIVALENCE (OUT   ,RAM(5*MPTS+1))
      REAL RECS(MPTS)
      EQUIVALENCE (IRECS,RECS)
C
C-----------------------------------------------------------------------
C
      CHARACTER*80 FILSEP,FCRT,CH
      INTEGER LU0
      PARAMETER (LU0=1)
C
C     Auxiliary storage locations:
      INTEGER LU1,LU2,LU3,MQ
      PARAMETER (LU1=1,LU2=2,LU3=3)
      PARAMETER (MQ=15)
C     1:X1, 2:X2, 3:X3, 4:TT, 5:P1, 6:P2, 7:P3,
C     8-15: 1, 2 or 4 complex-valued amplitudes
      CHARACTER*80 FILREC,FILSRC,FILE1,FILE2,FILE3
      CHARACTER*(4+8*MQ) FORMAT
      CHARACTER*80 RAYTXT
      INTEGER NQ,KALL,KSRC,KREC,INI,KTT
      INTEGER MPTS1,NPTS,LENTXT,IQ,II,I,J
      REAL OUTMIN(MQ),OUTMAX(MQ)
C
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+CRTPTS: Enter input filename: '
      FILSEP=' '
      READ(*,*) FILSEP
      WRITE(*,'(A)') '+CRTPTS: Working ...           '
C
C     Reading all data from the SEP file into the memory:
      IF (FILSEP.NE.' ') THEN
        CALL RSEP1(LU0,FILSEP)
      ELSE
C       CRTPTS-01
        CALL ERROR('CRTPTS-01: SEP file not given')
C       Input file in the form of the SEP (Stanford Exploration Project)
C       parameter or history file must be specified.
C       There is no default filename.
      ENDIF
C
C     Reading input parameters from the SEP file:
      CALL RSEP3T('CRTOUT',FCRT,' ')
      FILE1='s01.out'
      FILE2='s01i.out'
      IF (FCRT.NE.' ') THEN
        OPEN (LU1,FILE=FCRT,FORM='FORMATTED',STATUS='OLD')
        READ (LU1,*) CH,FILE1,FILE2,CH
        CLOSE (LU1)
      ENDIF
      CALL RSEP3T('REC',FILREC,' ')
      CALL RSEP3T('SRC',FILSRC,' ')
      CALL RSEP3T('PTS',FILE3,'pts.out')
C
C     Number of output quantities:
      CALL RSEP3I('NQ',NQ,11)
C     Switch between all rays and only two-point rays:
      CALL RSEP3I('KALL',KALL,0)
      CALL RSEP3I('KSRC',KSRC,0)
      CALL RSEP3I('KREC',KREC,0)
      CALL RSEP3I('KTT',KTT,0)
C
      NQ=MIN0(NQ,MQ)
      MPTS1=MIN0(MPTS,MOUT/NQ)
      INI =MAX0(0,MIN0(-KREC,1))
      CALL TXT1(LU3,FILSRC,FILREC)
      FORMAT(1:1)='('
      OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU3,FILE=FILE3)
      WRITE(LU3,'(A)') '/'
C
C.......................................................................
C
C     Loop for the points of rays:
      NPTS=0
   10 CONTINUE
C       Reading the results of the complete ray tracing
        CALL CRTOUT
     *           (LU1,LU2,KALL,KREC,INI,NQ,MPTS1,NPTS,OUT,OUTMIN,OUTMAX)
        IF(IWAVE.LT.1)THEN
C         End of rays
          GO TO 60
        END IF
        NPTS=NPTS-INI
        IRECS(NPTS)=IREC
        ISRCS(NPTS)=ISRC
        IWAVES(NPTS)=IWAVE
        IRAYS(NPTS)=IRAY
      GO TO 10
   60 CONTINUE
C
C.......................................................................
C
C     Sorting two-point rays:
      IF(KALL.LE.0) THEN
        CALL TTSORT(NQ,NPTS,4,OUT,IRECS,RECS,INDX)
      ELSE
        DO 71 I=1,NPTS
          INDX(I)=I
   71   CONTINUE
      END IF
C
C.......................................................................
C
C     Writing ray points:
C
C     Text describing the point:
C5.70 IF(KREC.LT.0) THEN
C5.70   LENTXT=9
C5.70   RAYTXT(1:LENTXT)='''000-000'''
C5.70   WRITE(RAYTXT(2:4),'(I3.3)') -KREC
C5.70 END IF
C
C     Writing:
      FORMAT(1:4)='(2A,'
      CALL FORM2(NQ,OUTMIN,OUTMAX,FORMAT(5:4+8*NQ))
      DO 89 I=1,NPTS
        J=INDX(I)
C
C       Text describing the point:
C5.70   IF(KREC.LT.0) THEN
C5.70     WRITE(RAYTXT(LENTXT-3:LENTXT-1),'(I3.3)') IRECS(J)
C5.70   ELSE
          CALL TXT2(KALL,KTT,ISRCS(J),IWAVES(J),IRAYS(J),IRECS(J),
     *                                                    LENTXT,RAYTXT)
C5.70   END IF
C
        J=NQ*(J-1)
        DO 81 IQ=NQ,1,-1
C         Not writing zeros:
          IF(ABS(OUT(IQ+J)).GT.0.) THEN
            GO TO 82
          END IF
   81   CONTINUE
   82   CONTINUE
        WRITE(LU3,FORMAT)
     *                   RAYTXT(1:LENTXT),(' ',OUT(II),II=1+J,IQ+J),' /'
   89 CONTINUE
C
      WRITE(LU3,'(A)') '/'
      CLOSE(LU3)
      CLOSE(LU2)
      CLOSE(LU1)
      WRITE(*,'(A)') '+CRTPTS: Done.                 '
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'ap.for'
C     ap.for
      INCLUDE 'crtout.for'
C     crtout.for
      INCLUDE 'ttsort.for'
C     ttsort.for
      INCLUDE 'indexx.for'
C     indexx.for of Numerical Recipes
C
C=======================================================================
C