C
C Program CRT2D3D transforming 2-D system of rays to 3-D system of rays
C to be processed by program MTT
C
C Version: 7.00
C Date: 2013, February 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     http://sw3d.cz/staff/klimes.htm
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
C     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. The CRT output files are assumed to contain
C             a 2-D system of rays.
C             For general description of file CRTOUT refer to file
C             writ.for.
C             Only files
C             'CRT-R',
C             'CRT-I' and
C             'CRT-T'
C             are read by CRT2D3D.
C             Default: CRTOUT=' ' which means 'CRT-R'='r01.out',
C                      'CRT-I'='r01i.out', 'CRT-T'='t01.out'
C     CRTNEW='string'... File with the names of the output files
C             to contain 3-D sytem of rays composed of two 2-D systems
C             of rays shifted by vectors -(DX1,DX2,DX3) and
C             +(DX1,DX2,DX3).
C             The files 'CRT-R', 'CRT-I' and 'CRT-T' are generated.
C             The form of file CRTNEW is the same as of file CRTOUT.
C             Default: CRTNEW=' ' which means 'CRT-R'='r01-3d.out',
C                      'CRT-I'='r01i-3d.out', 'CRT-T'='t01-3d.out'
C Data describing translation vector:
C     DX1=real, DX2=real, DX3=real ... Translation vector to be used
C             to create the 3-D system of rays composed of two input 2-D
C             systems of rays shifted by vectors -(DX1,DX2,DX3) and
C             +(DX1,DX2,DX3). The vector thus should be perpendicular
C             to the input 2-D system of rays.
C             Default: DX1=0., DX2=0., DX3=0.
C
C-----------------------------------------------------------------------
C
C Memory management:
      INCLUDE 'ram.inc'
C     ram.inc
      INTEGER IRAM(MRAM)
      EQUIVALENCE (IRAM,RAM)
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-----------------------------------------------------------------------
C
      CHARACTER*80 FILSEP,FCRT,FNEW,CH
      INTEGER LU0
      PARAMETER (LU0=1)
C
C     Auxiliary storage locations:
      INTEGER LU1,LU2,LU3,LU4,MIRAM,IRAY0,I,J,I1,I2,I3,I100,N100
      PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4)
      CHARACTER*80 FILE1,FILE2,FILE3,FILE4,FILE5,FILE6
      REAL DX1,DX2,DX3,X1,X2,X3
C
C     Output format:
      INTEGER IFORM
      CHARACTER*10 FORMAT
      DATA IFORM/99999/,FORMAT/'(I6,I6,I6)'/
C
C-----------------------------------------------------------------------
C
C     Reading name of SEP file with input data:
      WRITE(*,'(A)') '+CRT2D3D: Enter input filename: '
      FILSEP=' '
      READ(*,*) FILSEP
      WRITE(*,'(A)') '+CRT2D3D: Working ...           '
C
C     Reading all data from the SEP file into the memory:
      IF (FILSEP.NE.' ') THEN
        CALL RSEP1(LU0,FILSEP)
      ELSE
C       CRT2D3D-03
        CALL ERROR('CRT2D3D-03: 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 RSEP3R('DX1',DX1,0.)
      CALL RSEP3R('DX2',DX2,0.)
      CALL RSEP3R('DX3',DX3,0.)
      CALL RSEP3T('CRTOUT',FCRT,' ')
      CALL RSEP3T('CRTNEW',FNEW,' ')
C
C     Opening input and output files:
      FILE1='r01.out'
      FILE2='r01i.out'
      FILE3='t01.out'
      FILE4='r01-3d.out'
      FILE5='r01i-3d.out'
      FILE6='t01-3d.out'
      IF (FCRT.NE.' ') THEN
        OPEN (LU1,FILE=FCRT,FORM='FORMATTED',STATUS='OLD')
        READ (LU1,*) FILE1,CH,FILE2,FILE3
        CLOSE (LU1)
      ENDIF
      IF (FNEW.NE.' ') THEN
        OPEN (LU1,FILE=FNEW,FORM='FORMATTED',STATUS='OLD')
        READ (LU1,*) FILE4,CH,FILE5,FILE6
        CLOSE (LU1)
      ENDIF
C
C     Triangles:
      WRITE(*,'(2A)') '+CRT2D3D: Writing triangles ... '
      OPEN(LU1,FILE=FILE3,STATUS='OLD')
      OPEN(LU2,FILE=FILE6)
      N100=0
C
C     Loop for the triangles
   10 CONTINUE
C       Reading the interval
        READ(LU1,*,END=20) I1,I2,I3
        N100=MAX0(I1,I2,N100)
C       Setting output format
        IF(2*I1.GT.IFORM.OR.2*I2.GT.IFORM) THEN
          IFORM=IFORM*10+9
          FORMAT(3:3)=CHAR(ICHAR(FORMAT(3:3))+1)
          FORMAT(6:6)=FORMAT(3:3)
          FORMAT(9:9)=FORMAT(3:3)
        END IF
        IF (I1.EQ.0) THEN
C         New elementary wave:
          WRITE(LU2,FORMAT) I1,I2,I3
        ELSE
          IF(I3.NE.0) THEN
C           CRT2D3D-01
C           CALL ERROR('CRT2D3D-01: Input data are not 2-D')
          END IF
C         Writing the triangles
          WRITE(LU2,FORMAT) 2*I1-1,2*I1,2*I2-1
          WRITE(LU2,FORMAT) 2*I1,2*I2-1,2*I2
        END IF
      GO TO 10
C
   20 CONTINUE
      CLOSE(LU1)
      CLOSE(LU2)
C
C     Writing rays:
      WRITE(*,'(2A)') '+CRT2D3D:  0% of rays rewritten '
      OPEN(LU1,FILE=FILE1,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU2,FILE=FILE2,FORM='UNFORMATTED',STATUS='OLD')
      OPEN(LU3,FILE=FILE4,FORM='UNFORMATTED')
      OPEN(LU4,FILE=FILE5,FORM='UNFORMATTED')
      I100=0
C
C     Loop for the points of rays
      MIRAM=MRAM/6
   50 CONTINUE
C       Reading the results of the complete ray tracing
        CALL AP00(0,LU1,LU2)
        IF(IPT.LE.1.AND.IRAM(1).GT.0) THEN
C         New ray - recording the duplication of the last ray.
C         Writing the stored ray:
          CALL APW4(LU3,IRAM(1),RAM(MIRAM+1))
C         Deleteing the stored ray from RAM:
          IRAM(1)=0
        END IF
        IF(IWAVE.LT.1) THEN
C         End of rays
          GO TO 60
        END IF
        IF(IPT.LE.1)THEN
C         New ray - recording the initial points
          IRAY0=IRAY
          X1=YI(3)
          X2=YI(4)
          X3=YI(5)
          IRAY=2*IRAY0-1
          YI(3)=X1-DX1
          YI(4)=X2-DX2
          YI(5)=X3-DX3
          CALL APW1(LU4)
          IRAY=2*IRAY0
          YI(3)=X1+DX1
          YI(4)=X2+DX2
          YI(5)=X3+DX3
          CALL APW1(LU4)
          IRAY=IRAY0
        END IF
        IRAY0=IRAY
        X1=Y(3)
        X2=Y(4)
        X3=Y(5)
        IRAY=2*IRAY0-1
        Y(3)=X1-DX1
        Y(4)=X2-DX2
        Y(5)=X3-DX3
        CALL APW2(LU3)
        IRAY=2*IRAY0
        Y(3)=X1+DX1
        Y(4)=X2+DX2
        Y(5)=X3+DX3
        CALL APW3(IRAM(1),MIRAM,RAM(MIRAM+1),MRAM-MIRAM)
        IRAY=IRAY0
        IF(100*IRAY0.GE.I100*N100) THEN
          WRITE(*,'(A,I3)') '+CRT2D3D:',I100
          I100=I100+1
        END IF
      GO TO 50
C
   60 CONTINUE
      CLOSE(LU1)
      CLOSE(LU2)
      CLOSE(LU3)
      CLOSE(LU4)
      WRITE(*,'(A)') '+CRT2D3D: Done.                 '
C
      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
      INCLUDE 'apw.for'
C     apw.for
C
C=======================================================================
C