C
C Program GPSTEP to calculate the discretization steps
C for the decomposition of the wavefield into optimized
C 2-D Gaussian packets
C
C Version: 6.10
C Date: 2007, June 14
C
C Coded by: Karel Zacek
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
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 input files:
C     FGBYMIN='string'... Name of the input file containing the
C             minimum value of the optimized and smoothed
C             initial parameter of Gaussian beams Y0.
C             Default: FGBYMIN='fgbymin.out'
C     FGBYMAX='string'... Name of the input file containing the
C             maximum value of the optimized and smoothed
C             initial parameter of Gaussian beams Y0.
C             Default: FGBYMAX='fgbymax.out'
C     FGBRMIN='string'... Name of the input file containing the
C             minimum value of the optimized and smoothed
C             initial parameter of Gaussian beams R0.
C             Default: FGBRMIN='fgbrmin.out'
C     FGBRMAX='string'... Name of the input file containing the
C             maximum value of the optimized and smoothed
C             initial parameter of Gaussian beams R0.
C             Default: FGBRMAX='fgbrmax.out'
C     MODEL='string'... Input data file describing the model, it is
C             described in the subroutine file 'model.for'.
C     Description of file MODEL
C             Default:  'MODEL'='gpm-mod.dat'
C Data specifying the dimensions of the input grid:
C     NT=positive integer... Number of gridpoints along the time axis T
C             of the decomposed wavefield.
C             Default: NT=750
C     NX=positive integer... Number of gridpoints along the profile
C             axis X of the decomposed wavefield.
C             Default: NX=116
C     DT=real... Grid interval along the time axis T of the decomposed
C             wavefield.
C             Default: DT=0.004
C     DX=real... Grid interval along the profile axis X of
C             the decomposed wavefield.
C             Default: DX=25.0
C Name of input parameters:
C     SHOT=positive integer... Specifies the position of shot.
C             Default: SHOT=1
C     PARK=real... Parameter ksi with values in the interval (0,1).
C             Default: PARK=0.5
C     PARN=real... Parameters ny with values in the interval (0,1).
C             Default: PARN=0.5
C     KAPPA=real... Parameter kappa has the meaning of the ratio of
C             the step between two subsequent beams to their
C             characteristic half-width.
C             Default: KAPPA=1.253314137
C     VSURF=real... Velocity determined according to maximum velocity
C             in the depth of the receivers.
C             Default: VSURF=2000.
C     FMAX=real... Maximum frequency of the Gaussian packets.
C             Default: FMAX=1.
C     FMIN=real... Minimum frequency of the Gaussian packets.
C             Default: FMIN=1.
C     FREF=real... Reference frequencies of the Gaussian packets.
C             Default: FREF=1.
C Name of the output file:
C     GPPAR='string'... Name of the output file with parameters
C             K0, K0~, N0 - real part, N0 - imaginary part,
C             minimum of optimum initial parameter Y0.
C             Default: GPPAR='gpm-par.out'
C     GPGRID='string'... Name of the output file containing
C             discretization steps OXR,OPR,OTR,OOMEGA,DXR,DPR,DTR,
C             DOMEGA,NXR,NPR,NTR,NOMEGA
C             Default: GPGRID='gpgrid.out'
C     RPAR='string'... Name of the output file containing
C             ray parameters for back-propagation of the wavefield
C             from receivers.
C             Default: RPAR='gpm-rp.dat'
C     REC='string'... String with the name of the output data file
C             containing the names and coordinates of receivers.
C             Default: REC='gpm-rec.dat'
C     SRC='string'... String with the name of the output data file
C             containing the names and coordinates of sources.
C             Default: SRC='gpm-src.dat'
C
C-----------------------------------------------------------------------
C
C Subroutines and external functions required:
      EXTERNAL RARRAY
C     RARRAY...File 'forms.for'.
      EXTERNAL RARRAI
C     RARRAI...File 'forms.for'.
      EXTERNAL ERROR
C     ERROR...File 'error.for'.
      EXTERNAL RSEP1
C     RSEP1...File 'sep.for'.
      EXTERNAL RSEP3T
C     RSEP3T...File 'sep.for'.
      EXTERNAL RSEP3I
C     RSEP3I...File 'sep.for'.
      EXTERNAL RSEP3R
C     RSEP3R...File 'sep.for'.
      EXTERNAL PARM1
C     PARM1...File 'parm.for'.
      EXTERNAL PARM2
C     PARM2...File 'parm.for'.
C
C-----------------------------------------------------------------------
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
      INCLUDE 'sep.inc'
C     sep.inc
C
      INCLUDE 'model.inc'
C     model.inc
C
      INCLUDE 'auxmod.inc'
C     auxmod.inc
C
C-----------------------------------------------------------------------
C
C     Auxiliary storage locations:
      REAL DX,DT,FMAX,FMIN,OMEGAMAX,VSURF,FREF,OMEGAREF
C
      INTEGER LU1,LU2,LU3,LU4,LU5
      INTEGER IRAM(MRAM)
      INTEGER IXR,IPR,ISHOT
      INTEGER N,NX,NT,NXR,NTR,NPR,NOMEGA
      REAL DXR,DTR,DPR,DOMEGA,OXR,OPR,OOMEGA
      REAL XR,PR,KAPPA,PI,VEL,HUH
      REAL PRMIN,PRMAX,ANIMIN,ANIMAX,ANRMIN,ANRMAX,PARK,PARN
      REAL AK0,AK0D,THETA
      REAL COOR(3),RHO
      COMPLEX AN,ANDD,AN0,AMPL
      PARAMETER (LU1=21,LU2=22,LU3=23,LU4=24,LU5=25)
      PARAMETER (LU6=26,LU7=27,LU8=28,LU9=29,LU10=30,LU11=31)
      PARAMETER (PI=3.141592654)
      CHARACTER*80 FILSEP,FGBRMIN,FGBRMAX,FGBYMIN,FGBYMAX
      CHARACTER*80 GPGRID,GPPAR,RPAR,MODEL,REC,SRC
      EQUIVALENCE (IRAM,RAM)
C
C.......................................................................
C
C     Reading main input data:
      FILSEP=' '
      WRITE(*,'(A)')'+GPSTEP: Enter filename : '
      READ(*,*) FILSEP
      IF(FILSEP.EQ.' ') THEN
C       GPSTEP-01
        CALL ERROR('GPSTEP-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)')'+GPSTEP: 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('FGBYMIN',FGBYMIN,'fgbymin.out')
      CALL RSEP3T('FGBYMAX',FGBYMAX,'fgbymax.out')
      CALL RSEP3T('FGBRMIN',FGBRMIN,'fgbrmin.out')
      CALL RSEP3T('FGBRMAX',FGBRMAX,'fgbrmax.out')
      CALL RSEP3T('MODEL',MODEL,'gpm-mod.dat')
C
C     Reading filenames of the output files
      CALL RSEP3T('GPPAR',GPPAR,'gpm-par.out')
      CALL RSEP3T('GPGRID',GPGRID,'gpgrid.out')
      CALL RSEP3T('RPAR',RPAR,'gpm-rp.dat')
      CALL RSEP3T('REC',REC,'gpm-rec.dat')
      CALL RSEP3T('SRC',SRC,'gpm-src.dat')
C
      CALL RSEP3I('NT',NT,750)
      CALL RSEP3I('NX',NX,116)
      CALL RSEP3R('DT',DT,0.004)
      CALL RSEP3R('DX',DX,25.0)
C
      CALL RSEP3I('SHOT',ISHOT,1)
C
      CALL RSEP3R('PARK',PARK,0.5)
      CALL RSEP3R('PARN',PARN,0.5)
      CALL RSEP3R('KAPPA',KAPPA,1.253314137)
      CALL RSEP3R('VSURF',VSURF,2000.)
      CALL RSEP3R('FMAX',FMAX,1.)
      CALL RSEP3R('FMIN',FMIN,1.)
      CALL RSEP3R('FREF',FREF,1.)
C
      OPEN(LU2,FILE=FGBYMIN,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU3,FILE=FGBYMAX,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU4,FILE=FGBRMIN,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU5,FILE=FGBRMAX,FORM='FORMATTED',STATUS='OLD')
      OPEN(LU6,FILE=MODEL,FORM='FORMATTED',STATUS='OLD')
C
      OPEN(LU7,FILE=GPPAR,FORM='FORMATTED')
      OPEN(LU8,FILE=GPGRID,FORM='FORMATTED')
      OPEN(LU9,FILE=RPAR,FORM='FORMATTED')
      OPEN(LU10,FILE=REC,FORM='FORMATTED')
      OPEN(LU11,FILE=SRC,FORM='FORMATTED')
C
      READ(LU2,*) ANIMIN
      READ(LU3,*) ANIMAX
      READ(LU4,*) ANRMIN
      READ(LU5,*) ANRMAX
C
      OMEGAMAX=2.0*PI*FMAX
      OOMEGA=2.0*PI*FMIN
      OMEGAREF=2.0*PI*FREF
      IF(OOMEGA.LT.(1.0/(DT*(NT-1)))) THEN
        OOMEGA= 1.0/(DT*(NT-1))
      END IF
C
C      AK0=PARK*OOMEGA*ANIMIN*VSURF**2
C      AK0D=AK0*(1.0-PARN)/PARN
C      AN0=PARN*2.0*CMPLX(ANRMIN,ANIMIN)
      AK0=0.5*OOMEGA*ANIMIN*VSURF**2
      AK0D=AK0
      AN0=CMPLX(ANRMIN,ANIMIN)
C
      WRITE(*,'(A,1(G16.6,1X))') 'K0    =',AK0
      WRITE(*,'(A,1(G16.6,1X))') 'K0D   =',AK0D
      WRITE(*,'(A,1(G16.6,1X))') 'N0R   =',REAL(AN0)
      WRITE(*,'(A,1(G16.6,1X))') 'N0I   =',AIMAG(AN0)
      WRITE(*,'(A,1(G16.6,1X))') 'NMIN  =',ANIMIN
      WRITE(LU7,'(1(G16.6,1X))')           AK0
      WRITE(LU7,'(1(G16.6,1X))')           AK0D
      WRITE(LU7,'(1(G16.6,1X))')           REAL(AN0)
      WRITE(LU7,'(1(G16.6,1X))')           AIMAG(AN0)
      WRITE(LU7,'(1(G16.6,1X))')           ANIMIN
C
      IF(OOMEGA.LT.OMEGAMAX) THEN
        DOMEGA=KAPPA*SQRT(2.0*AK0*AK0D/(AK0+AK0D))
        N=NINT((OMEGAMAX-OOMEGA)/DOMEGA)
        NOMEGA=N+2
        DOMEGA=(OMEGAMAX-OOMEGA)/(N+1)
      ELSE
        OOMEGA=OMEGAMAX
        DOMEGA=1.0
        NOMEGA=1
      END IF
      WRITE(*,'(A,1(G16.6,1X))') 'DF    =',DOMEGA/(2.0*PI)
      WRITE(*,'(A,I6)')          'NF    =',NOMEGA
C
C      DXR=KAPPA*SQRT(-2.0*(1.0-PARN)/OMEGAREF*AIMAG(
C     *   (CMPLX(1.0,0.0)-OOMEGA/OMEGAREF*
C     *    CMPLX(0.0,ANIMIN)/CMPLX(ANRMAX,ANIMAX)*PARK)/
C     *    CMPLX(ANRMAX,ANIMAX)))
      DXR=KAPPA*SQRT(-AIMAG(0.5/(OMEGAREF*AN0)))
      N=NINT((NX-1)*DX/DXR)
      NXR=N+2
      DXR=(NX-1)*DX/(N+1)
      WRITE(*,'(A,1(G16.6,1X))') 'DXR   =',DXR
      WRITE(*,'(A,I6)')          'NXR   =',NXR
C
      OPR=-1.0/VSURF
      DPR=KAPPA*SQRT(AIMAG(AN0)/OMEGAREF)
      N=NINT(2.0/(VSURF*DPR))
      NPR=N+2
      DPR=2.0/(VSURF*(N+1))
      WRITE(*,'(A,1(G16.6,1X))') 'DP    =',DPR
      WRITE(*,'(A,I6)')          'NP    =',NPR
C
C      DTR=KAPPA*SQRT(-2.0*PARN*AIMAG(
C     *   (CMPLX(1.0,0.0)-CMPLX(0.0,ANIMIN)/CMPLX(ANRMIN,ANIMIN)*PARK)
C     *   /CMPLX(0.0,AK0)))
      DTR=KAPPA*SQRT(0.5/AK0)
      N=NINT((NT-1)*DT/DTR)
      NTR=N+2
      DTR=(NT-1)*DT/(N+1)
      WRITE(*,'(A,1(G16.6,1X))') 'DTR   =',DTR
      WRITE(*,'(A,I6)')          'NTR   =',NTR
C
      OXR=150.0+ISHOT*25.0
      OTR=0.0
C
      WRITE(LU8,'(1(G16.6,1X))')           OXR
      WRITE(LU8,'(1(G16.6,1X))')           OPR
      WRITE(LU8,'(1(G16.6,1X))')           OTR
      WRITE(LU8,'(1(G16.6,1X))')           OOMEGA
      WRITE(LU8,'(1(G16.6,1X))')           DXR
      WRITE(LU8,'(1(G16.6,1X))')           DPR
      WRITE(LU8,'(1(G16.6,1X))')           DTR
      WRITE(LU8,'(1(G16.6,1X))')           DOMEGA
      WRITE(LU8,'(I6)')                    NXR
      WRITE(LU8,'(I6)')                    NPR
      WRITE(LU8,'(I6)')                    NTR
      WRITE(LU8,'(I6)')                    NOMEGA
C
C      CALL PARM1(LU6,1)
      CALL MODEL1(LU6)
      WRITE(LU9,'(A)') '''RPAR'''
      WRITE(LU9,'(A)') '0 0 0 /'
      WRITE(LU10,'(A)') '/'
      DO 20 IXR=1,NXR
        XR=OXR+(IXR-1)*DXR
        COOR(1)=0.0
        COOR(2)=XR
        COOR(3)=12.0
        WRITE(LU10,'(A,3(G16.6,1X),A)')
     *             '''REC'' ', COOR(1),COOR(2),COOR(3),' /'
        CALL PARM2(1,COOR,UP,US,RHO,QP,QS)
        VEL=UP(1)
C        VEL=2000.0
        DO 10 IPR=1,NPR
          HUH=VEL*(OPR+(IPR-1)*DPR)
          IF(HUH.LT.-1.0) THEN
            THETA=-PI/2.0
          ELSE IF(HUH.GT.1.0) THEN
            THETA=PI/2.0
          ELSE
            THETA=ASIN(HUH)
          END IF
          WRITE(LU9,'(A,1(G16.6,1X),A)')'0 ',THETA,' 0 2 1 0 0.1 0.1 /'
   10   CONTINUE
        WRITE(LU9,'(A)') '/'
   20 CONTINUE
      WRITE(LU10,'(A)') '/'
C
      XR=3000.0+(ISHOT-1)*25.0
      COOR(1)=0.0
      COOR(2)=XR
      COOR(3)=12.0
      WRITE(LU11,'(A)') '/'
      WRITE(LU11,'(A,3(G16.6,1X),A)')
     *           '''SRC'' ', COOR(3),COOR(2),COOR(1),' /'
      WRITE(LU11,'(A)') '/'
C
      CLOSE(LU2)
      CLOSE(LU3)
      CLOSE(LU4)
      CLOSE(LU5)
      CLOSE(LU6)
      CLOSE(LU7)
      CLOSE(LU8)
      CLOSE(LU9)
      WRITE(*,'(A)') '+GPSTEP: 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