C
C Program GPRPAR to calculate the take-off parameters of rays
C for complete ray tracing from initial point.
C
C Version: 6.30
C Date: 2009, June 5
C
C Coded by: Karel Zacek
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz
C
C The code of this program was a part of the original gpstep.for
C program (SW3D-CD-10).
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, read by list directed (free
C     format) input.  The strings have thus to be enclosed in
C     apostrophes.  The interactive * external unit may be redirected to
C     the file containing the data.
C (1) 'SEP'/
C     'SEP'...Name of the file with input parameters.
C             Description of file SEP
C             If blank, default values of the corresponding data are
C             considered.
C Default: 'SEP'=' '.
C
C                                                     
C Input file 'SEP' in the SEP format:
C     The file has the form of a SEP parameter file.
C     For the description of the SEP format refer to file
C     'sep.for'.
C Names of the input files:
C     MODEL='string'... Name of the input file with the data specifying
C             the model.
C             Description of 
C             MODEL.
C             No default, 'MODEL' must be specified and cannot be blank.
C     GPSTEP='string'... Name of the input file containing
C             discretization steps NXR,DXR,OXR,NTR,DTR,OTR,NP,DP,OP,
C             NOMEGA,DOMEGA,OOMEGA, complex-valued constants
C             N0, K0, ~K0, Im(N)min (same as Y0min).
C             Default: GPSTEP=' '
C Names of input parameters:
C     X1SRC=real, X2SRC=real, X3SRC=real... Initial coordinates of
C             the source.
C             Default: X1SRC=0., X2SRC=0., X3SRC=0.
C     X1REC=real, X2REC=real, X3REC=real... Initial coordinates of
C             the receivers.
C             Default: X1REC=0., X2REC=0., X3REC=0.
C     ISRC=positive integer... Index of the source (shot).
C             Default: ISRC=1
C     DSRC=real... Distance between sources (shots).
C             Default: DSRC=1.
C     OSRC=real... Position of the first source (shot).
C             Default: OSRC=1.
C     OXRINI=real... Initial value of OXR.
C             Default: OXRINI=1.
C     KOLPAR=positive integer... Column where the take-off parameters
C             will be written in the file RPAR. Allowable values
C             are 1 and 2.
C             Default: KOLPAR=3.
C Names of the output files:
C     RPAR='string'... Name of the output file containing
C             the take-off parameters of rays for complete ray tracing
C             from initial point.
C             Default: RPAR=' '
C     SRC='string'... String with the name of the output data file
C             with the names and coordinates of the sources.
C             Default: SRC=' '
C     REC='string'... String with the name of the output data file
C             with the names and coordinates of the receivers.
C             Default: REC=' '
C
C=======================================================================
C
C Common block /AUXMOD/:
      INCLUDE 'auxmod.inc'
C     auxmod.inc
C
C.......................................................................
C
C     External functions and subroutines:
      EXTERNAL RSEP1,RSEP3T,RSEP3I,RSEP3R,ERROR,PARM2
C
C     Filenames and parameters:
      CHARACTER*80 FILSEP,RPAR,MODEL,SRC,REC,GPSTEP
      INTEGER LU1,LU2,LU3
      PARAMETER (LU1=1,LU2=2,LU3=3)
C
C     Other variables:
      INTEGER IXR,NXR,NTR,IP,NP,ISRC,KOLPAR
      REAL DXR,OXR,DTR,OTR,DP,OP,DSRC,OSRC,OXRINI,XR,VXR,RATP,PARL
      REAL X1SRC,X2SRC,X3SRC,X1REC,X2REC,X3REC
      REAL COOR(3)
      CHARACTER*80 TEXT
C
C.......................................................................
C
C     Reading main input data:
      FILSEP=' '
      WRITE(*,'(A)')'+GPRPAR: Enter filename : '
      READ(*,*) FILSEP
      IF(FILSEP.EQ.' ') THEN
C       GPRPAR-01
        CALL ERROR('GPRPAR-01: No input file specified')
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
      WRITE(*,'(A)')'+GPRPAR: Reading, calculating... '
C
C     Reading all the data from the SEP file into the memory:
      CALL RSEP1(LU1,FILSEP)
C
C     Reading filenames of the input files
      CALL RSEP3T('MODEL',MODEL,' ')
      CALL RSEP3T('GPSTEP',GPSTEP,' ')
C
C     Reading input parameters
      CALL RSEP3R('X1SRC',X1SRC,1.)
      CALL RSEP3R('X2SRC',X2SRC,1.)
      CALL RSEP3R('X3SRC',X3SRC,1.)
      CALL RSEP3R('X1REC',X1REC,1.)
      CALL RSEP3R('X2REC',X2REC,1.)
      CALL RSEP3R('X3REC',X3REC,1.)
      CALL RSEP3I('ISRC',ISRC,1)
      CALL RSEP3R('DSRC',DSRC,1.)
      CALL RSEP3R('OSRC',OSRC,1.)
      CALL RSEP3R('OXRINI',OXRINI,1.)
      CALL RSEP3I('KOLPAR',KOLPAR,3)
C
C     Reading filenames of the output files
      CALL RSEP3T('RPAR',RPAR,' ')
      CALL RSEP3T('SRC',SRC,' ')
      CALL RSEP3T('REC',REC,' ')
      CLOSE(LU1)
C
C     Reading discretization steps
      OPEN(LU1,FILE=GPSTEP,FORM='FORMATTED',STATUS='OLD')
      READ(LU1,*) TEXT
      READ(LU1,*) NXR,DXR,OXR
      READ(LU1,*) NTR,DTR,OTR
      READ(LU1,*) NP,DP,OP
      CLOSE(LU1)
C
      WRITE(*,'(A,I6)')          ' ISRC=',ISRC
      WRITE(*,'(A,1(G16.6,1X))') ' DSRC=',DSRC
      WRITE(*,'(A,1(G16.6,1X))') ' OSRC=',OSRC
      WRITE(*,'(A,I6)')          ' NP   =',NP
      WRITE(*,'(A,1(G16.6,1X))') ' DP   =',DP
      WRITE(*,'(A,1(G16.6,1X))') ' OP   =',OP
      WRITE(*,'(A,I6)')          ' NXR  =',NXR
      WRITE(*,'(A,1(G16.6,1X))') ' DXR  =',DXR
C
      OXR=OXRINI+(ISRC-1)*DSRC
      WRITE(*,'(A,1(G16.6,1X))') ' OXR  =',OXR
C
      OPEN(LU1,FILE=MODEL,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU2,FILE=REC,FORM='FORMATTED')
      OPEN(LU3,FILE=RPAR,FORM='FORMATTED')
      CALL MODEL1(LU1)
      WRITE(LU3,'(A)') '''RPAR'''
      WRITE(LU3,'(A)') '0 0 0 /'
      WRITE(LU2,'(A)') '/'
C     Loop over intersection of central ray of Gaussian paket
C     with profile XR
      DO 20 IXR=1,NXR
        XR=OXR+(IXR-1)*DXR
        COOR(1)=X1REC
        COOR(2)=XR
        COOR(3)=X3REC
        WRITE(LU2,'(A,3(G16.6,1X),A)')
     *             '''REC'' ', COOR(1),COOR(2),COOR(3),' /'
        CALL PARM2(1,COOR,UP,US,RHO,QP,QS)
        VXR=UP(1)
        WRITE(*,'(A,1(G16.6,1X))') 'VXR  =',VXR
C       Loop over slowness vector along profile P
        DO 10 IP=1,NP
          RATP=VXR*(OP+(IP-1)*DP)
          IF((RATP.LT.-1.).OR.(RATP.GT.1.)) THEN
C           GPRPAR-02
            CALL ERROR('GPRPAR-02: Decrease value of OP')
          ELSE
            PARL=ASIN(RATP)
          END IF
          IF (KOLPAR.EQ.1) THEN
            WRITE(LU3,'(1(G16.6,1X),A)')PARL,' 0 0 0 0 0 0 0 0 /'
          ELSE IF (KOLPAR.EQ.2) THEN
            WRITE(LU3,'(A,1(G16.6,1X),A)')'0 ',PARL,' 0 0 0 0 0 0 0 /'
          ELSE
C           GPRPAR-03
            CALL ERROR('GPRPAR-03: KOLPAR should be equal to 1 or 2')
          END IF
   10   CONTINUE
      WRITE(LU3,'(A)') '/'
   20 CONTINUE
      WRITE(LU2,'(A)') '/'
      CLOSE(LU1)
      CLOSE(LU2)
      CLOSE(LU3)
C
C     Source coordinates
      XR=OSRC+(ISRC-1)*DSRC
      COOR(1)=X1SRC
      COOR(2)=XR
      COOR(3)=X3SRC
      OPEN(LU1,FILE=SRC,FORM='FORMATTED')
      WRITE(LU1,'(A)') '/'
      WRITE(LU1,'(A,3(G16.6,1X),A)')
     *           '''SRC'' ', COOR(1),COOR(2),COOR(3),' /'
      WRITE(LU1,'(A)') '/'
      CLOSE(LU1)
C
      WRITE(*,'(A)') '+GPRPAR: Done. '
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'model.for'
C     model.for
      INCLUDE 'srfc.for'
C     srfc.for
      INCLUDE 'metric.for'
C     metric.for
      INCLUDE 'parm.for'
C     parm.for
      INCLUDE 'val.for'
C     val.for
      INCLUDE 'fit.for'
C     fit.for
C
C=======================================================================
C