C 
C Program CRTRAY converting the unformatted output of program CRT into
C formatted files with rays suitable for plotting.
C
C Version: 5.10
C Date: 1996, September 26
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 'Lines' described in 'formsdoc.htm' and designed to
C store 3-D curves, e.g., for the purposes of plotting.
C
C.......................................................................
C
C                                                    
C Description of data files:
C
C Main input data read from external interactive device (*):
C     The data consist of character strings and integers, read by list
C     directed (free format) input.  The strings have thus to be
C     enclosed in apostrophes.  The interactive * external unit may be
C     redirected to the file containing the data.
C (1) 'REC','SRC','RAYPOINTS','INITIALPOINTS','RAYS',NQ,KALL
C     'REC'...If non-blank, the name of the file with the names of the
C             receiver points.  The names are then used within the
C             strings describing the two-point rays.  Otherwise, the
C             two-point rays are denoted by the receiver index.
C             Description of file REC
C     'SRC'...If non-blank, the name of the file with the name of the
C             source point.  The name is then used within the strings
C             describing the rays.
C             Description of file SRC
C     'RAYPOINTS'... File with the quantities stored along rays (see
C             C.R.T.5.5.1).
C     'INITIALPOINTS'... File with the quantities at the initial points
C             of rays, corresponding to file 'RAYPOINTS' (see
C             C.R.T.6.1).
C     'RAYS'..Name of the output formatted file with rays.  This file is
C             designed, for instance, as input file for plotting rays.
C             Description of file RAYS
C     NQ...   Number of reals in each output line.
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 of ray amplitude, normalized to 1 at an
C                initial surface or along on a unit sphere around a
C                point source, corresponding to P- or S1-polarization at
C                the initial point of the ray.
C             9. Imaginary part of ray amplitude corresponding to P- or
C                S1-polarization at the initial point of the ray.
C                Printed only if greater than 0.000001 in abs value.
C            10. Real part of ray amplitude corresponding to
C                S2-polarization at the initial point of the ray.
C                Printed only if greater than 0.000001 in abs value.
C            11. Imaginary part of ray amplitude corresponding
C                S2-polarization at the initial point of the ray.
C                Printed only if greater than 0.000001 in abs value.
C     KALL... 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: 'REC'=' ', 'SRC'=' ', 'RAYPOINTS'='r01.out',
C     'INITIALPOINTS'='r01i.out', 'RAYS'='rays.out', NQ=4, KALL=0.
C
C Input unformatted file RAYPOINTS:
C     See the description within source code file 'writ.for'.
C     Description of file RAYPOINTS
C
C Input unformatted file INITIALPOINTS:
C     See the description within source code file 'writ.for'.
C     Description of file INITIALPOINTS
C
C                                                    
C Output formatted file RAYS (format LIN):
C (1) / (a slash).
C (2) For each ray (2.1), (2.2), and (2.3):
C (2.1) 'RAYTXT',/
C     'RAYTXT'... String in apostrophes describing the ray.  See the
C             description of input parameter KALL.
C     /...    An obligatory slash after the string, in place of the
C             coordinates of the reference point.
C (2.2) For each point of the ray (2.1.1):
C (2.2.1) (OUT(I),I=1,NQ),/
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=4:
C     (2.2.1) X1,X2,X3,TT,/
C     X1,X2,X3... Coordinates of the point of the ray.
C     TT...   Arrival time at the point.
C (2.3) / (a slash).
C (3) / (a slash).
C Description of format LIN
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,FORM2
C     CRTOUT..File 'crtout.for'.
C     AP00... File 'ap.for' (called by CRTOUT).
C     FORM2...File 'forms.for'.
C     FORM1...File 'forms.for' (called by FORM2).
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      INTEGER LU1,LU2,LU3,MQ,MPTS,MOUT
      PARAMETER (LU1=1,LU2=2,LU3=3)
      PARAMETER (MQ=11,MPTS=2,MOUT=MQ*MPTS)
C     1:X1, 2:X2, 3:X3, 4:TT
C     Optional:           5:P1, 6:P2, 7:P3, 8:AR1, 9:AI1, 10:AR2, 11:AI2
      CHARACTER*80 FILE1,FILE2,FILE3,FILSRC,FILREC
      CHARACTER*(2+8*MQ) FORMAT
      CHARACTER*80 RAYTXT
      INTEGER NQ,KALL,NPTS,LENTXT,IQ
      REAL OUT(MOUT),OUTMIN(MQ),OUTMAX(MQ)
C
C.......................................................................
C
C     Opening input and output files:
      FILREC=' '
      FILSRC=' '
      FILE1='r01.out'
      FILE2='r01i.out'
      FILE3='rays.out'
C     Number of output quantities:
      NQ  =4
C     Switch between all rays and only two-point rays:
      KALL=0
      WRITE(*,'(2A)')
     *  ' Enter 5 filenames (REC,SRC,R01,R01I,RAYS),',
     *  ' and 2 integers (NQ,KALL): '
      READ(*,*) FILREC,FILSRC,FILE1,FILE2,FILE3,NQ,KALL
      NQ=MIN0(NQ,MQ)
      CALL TXT1(LU3,FILSRC,FILREC)
      OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU3,FILE=FILE3)
C
C.......................................................................
C
C     Loop for the points of rays
   10 CONTINUE
C       Reading the results of the complete ray tracing
        NPTS=0
        CALL CRTOUT(LU1,LU2,KALL,0,1,MQ,MPTS,NPTS,OUT,OUTMIN,OUTMAX)
        IF(IWAVE.LT.1)THEN
C         End of rays
          GO TO 80
        END IF
C       Writing the results of the complete ray tracing
        FORMAT(1:1)='('
        CALL FORM2(NQ,OUTMIN,OUTMAX,FORMAT(2:2+8*NQ))
        IF(IPT.LE.1)THEN
C         New ray - recording the initial point
          CALL TXT2(KALL,0,IWAVE,IRAY,IREC,LENTXT,RAYTXT)
          WRITE(LU3,'(A)') '/'
          WRITE(LU3,'(2A)') RAYTXT(1:LENTXT),' /'
          WRITE(LU3,FORMAT) OUT(1),(' ',OUT(IQ),IQ=2,NQ),' /'
        END IF
        WRITE(LU3,FORMAT) OUT(1+MQ*(NPTS-1)),
     *        (' ',OUT(IQ),IQ=2+MQ*(NPTS-1),MQ*(NPTS-1)+NQ),' /'
      GO TO 10
C
   80 CONTINUE
      WRITE(LU3,'(A)') '/'
      WRITE(LU3,'(A)') '/'
      CLOSE(LU3)
      CLOSE(LU2)
      CLOSE(LU1)
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'ap.for'
C     ap.for
      INCLUDE 'crtout.for'
C     crtout.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'length.for'
C     length.for
C
C=======================================================================
C