C 
C Program CRTLEW to calculate the directions of the projections of the
C ray segments into the given 2-D section and to collect them according
C to the given angular intervals.
C
C Version: 5.70
C Date: 2003, May 18
C
C Coded by: Petr Bulant
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: bulant@seis.karlov.mff.cuni.cz
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-R', and 'CRT-I' are
C             read by CRTLEW.
C             For general description of file CRTOUT refer to
C             file writ.for.
C             Default: CRTOUT=' ' which means 'CRT-R'='r01.out',
C                      'CRT-I'='r01i.out'
C     MODLEW='string' ... Name of the output formatted file with the
C             distribution of the ray segments directions according
C             to the given angular intervals.
C             Description of file MODLEW
C             Default: MODLEW='modlew.out'
C Switch between all rays and two-point rays only:
C     KALL=integer:
C             KALL.LE.0: only two-point rays are considered,
C             KALL.GE.1: all rays are considered.
C             Default: KALL=0
C Specification of the 2-D section for the projection of the rays:
C     KOOR1=integer... Index of the first coordinate determining the
C             2-D section for the calculation of the ray directions.
C             KOOR1 must be 1, 2 or 3.
C             Default: KOOR1=1
C     KOOR2=integer... Index of the second coordinate determining the
C             2-D section for the calculation of the ray directions.
C             KOOR2 must be 1, 2 or 3 and must differ from KOOR1.
C             Default: KOOR2=2
C Data to specify the angular intervals, according which the ray
C directions should be sorted:
C     NA=integer... Number of angular directions.
C             Default: NA=90
C     DA...   Angular step.
C             The default value is recommended.
C             Default: DA=3.141592/NA
C     OA...   The first angle. The angles are defined in radians,
C             -pi/2 for positive half-axis KOOR1, 0 for positive
C             half-axis KOOR2, pi/2 for negative half-axis KOOR1.
C             The default value is usually sufficient.
C             Default: OA=-1.570796
C
C                                                  
C Output formatted file MODLEW:
C     The program calculates lenghts of the ray segments projected
C     to the plane defined by COOR1 and COOR2.  (By the ray segment
C     we understand the line connecting two successive points in the
C     file CRT-R.)  For each projected ray segment, the slowness vector
C     at its endpoints is averaged, and the average is considered to be
C     the direction of the segment.  For each angular interval given by
C     OA, DA and NA the lengths of the ray segments of the direction
C     corresponding to the interval are accumulated.  Finally, the
C     accumulated lengths are divided by the length af all the ray
C     segments, multiplied by NA, and the resulting number is written to
C     the file MODLEW.  The file has form
C     LINes:
C (1) String identifying the file, terminated by / (a slash).
C (2) The header 2.1 and for each angular interval (NA times)
C     the values 2.2:
C (2.1) 'W' / ... Reference text, terminated by / (a slash).
C (2.2) A W / ... Mean angle A of the interval, and the relative length
C         of the ray segments corresponding to the interval.  The data
C         2.2 repeat NA times.
C (3) / Checksum CHESUM ... CHESUM is the sum of all the lenghts W.
C (4) / Average angular change AVDFI ... AVDFI is the average difference
C         of the directions of the ray along ray segments.
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:
C     AP00... File 'ap.for'.
C
C-----------------------------------------------------------------------
C
      CHARACTER*80 FILSEP,FCRT,CH
      INTEGER LU0
      PARAMETER (LU0=1)
C
C     Auxiliary storage locations:
      INTEGER LU1,LU2,LU3,MA
      PARAMETER (LU1=1,LU2=2,LU3=3,MA=360)
      CHARACTER*80 FILE1,FILE2,FILE3
      INTEGER KALL,KOOR1,KOOR2,NA,NFI,I1
      REAL RAYLEN(MA),FI,CFI,DFI,SUMLEN,SUMDFI,SUMSEG,SEGLEN,AVDFI,
     *     FIOLD,YOLD1,YOLD2,OA,DA,A,CHESUM
      DATA RAYLEN/MA*0./
C
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+CRTLEW: Enter input filename: '
      FILSEP=' '
      READ(*,*) FILSEP
      WRITE(*,'(A)') '+CRTLEW: Working ...           '
C
C     Reading all data from the SEP file into the memory:
      IF (FILSEP.NE.' ') THEN
        CALL RSEP1(LU0,FILSEP)
      ELSE
C       CRTLEW-01
        CALL ERROR('CRTLEW-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='r01.out'
      FILE2='r01i.out'
      IF (FCRT.NE.' ') THEN
        OPEN (LU1,FILE=FCRT,FORM='FORMATTED',STATUS='OLD')
        READ (LU1,*) FILE1,CH,FILE2,CH
        CLOSE (LU1)
      ENDIF
      CALL RSEP3T('MODLEW',FILE3,'modlew.out')
C
C     Switch between all rays and only two-point rays:
      CALL RSEP3I('KALL',KALL,0)
C     Coordinates of the section:
      CALL RSEP3I('KOOR1',KOOR1,1)
      CALL RSEP3I('KOOR2',KOOR2,2)
      IF (KOOR1.EQ.KOOR2) THEN
C       CRTLEW-02
        CALL ERROR('CRTLEW-02: KOOR1 equal to KOOR2')
C       KOOR2 must differ from KOOR1.
      ENDIF
C     Angular intervals:
      CALL RSEP3I('NA',NA,90)
      IF(NA.GT.MA) THEN
C       CRTLEW-03
        CALL ERROR('CRTLEW-03: Too many directions')
C       The dimension MA of the array RAYLEN should be increased.
      END IF
      CALL RSEP3R('DA',DA,3.141592/FLOAT(NA))
      CALL RSEP3R('OA',OA,-1.570796+0.5*DA)
      IF(OA.LT.-1.570796.OR.OA+DA*FLOAT(NA-1).GT.1.570796) THEN
C       CRTLEW-04
        CALL ERROR('CRTLEW-04: Wrong direction')
      END IF
C
      OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU3,FILE=FILE3)
      WRITE(LU3,'(9(A,F8.5))') '''Directions and relative weights'' /'
C
C.......................................................................
C
C     Loop for the points of rays
      SUMDFI=0.
      SUMSEG=0.
   10 CONTINUE
        CALL AP00(0,LU1,LU2)
        IF (IWAVE.LT.1) THEN
C         End of rays:
          GOTO 80
        ELSEIF (KALL.GT.0.OR.IREC.GT.0) THEN
C         This ray is to be processed:
C         Calculating the direction
C         at the point of the ray:
          IF (Y(5+KOOR2).EQ.0.) THEN
            IF (Y(5+KOOR1).LT.0.) THEN
              FI=1.570796
            ELSE
              FI=-1.570796
            ENDIF
          ELSE
            FI=ATAN(-Y(5+KOOR1)/Y(5+KOOR2))
          ENDIF
          IF (IPT.GT.1) THEN
C           Length of the ray segment:
            SEGLEN=SQRT((YOLD1-Y(2+KOOR1))**2 + (YOLD2-Y(2+KOOR2))**2)
C           Angular change along the ray segment and the direction
C           in the center of the segment:
            DFI=ABS(FI-FIOLD)
            CFI=(FI+FIOLD)/2.
            IF (DFI.GT.1.570796) THEN
              DFI=3.141592-DFI
              IF (CFI.LT.0.) THEN
                CFI=1.570796+CFI
              ELSE
                CFI=-1.570796+CFI
              ENDIF
            ENDIF
            SUMDFI=SUMDFI+DFI
            SUMSEG=SUMSEG+1.
C           Number of the interval for the center of the ray segment:
            NFI=NINT((CFI-OA)/DA)+1
            IF (NFI.EQ.NA+1) NFI=NA
            IF (NFI.EQ.0)    NFI=1
            RAYLEN(NFI)=RAYLEN(NFI)+SEGLEN
            SUMLEN=SUMLEN+SEGLEN
          ENDIF
C         Storing the quantities at the point of the ray:
          FIOLD=FI
          YOLD1=Y(2+KOOR1)
          YOLD2=Y(2+KOOR2)
        ENDIF
      GOTO 10
C
  80  CONTINUE
      IF ((SUMSEG.EQ.0.).OR.(SUMLEN.EQ.0.)) GOTO 100
      AVDFI=SUMDFI/SUMSEG
      WRITE(LU3,'(9(A,F8.5))') '''W'' /'
      CHESUM=0.
      DO 90, I1=1,NA
        RAYLEN(I1)=RAYLEN(I1)/SUMLEN*FLOAT(NA)
        A=OA+DA*FLOAT(I1-1)
        WRITE(LU3,'(9(A,F8.5))') ' ',A,' ',RAYLEN(I1),' /'
        CHESUM=CHESUM+RAYLEN(I1)
  90  CONTINUE
C
  100 CONTINUE
      WRITE(LU3,'(9(A,F8.5))') '/ Checksum ',CHESUM
      WRITE(LU3,'(9(A,F8.5))') '/ Average angular change ',AVDFI,' rad,'
      WRITE(LU3,'(9(A,F8.3))') '/ i.e. ',180*AVDFI/3.141592,' deg.'
      CLOSE(LU3)
      CLOSE(LU2)
      CLOSE(LU1)
      WRITE(*,'(A)') '+CRTLEW: Done.                 '
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'ap.for'
C     ap.for
C
C=======================================================================
C