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.00
C Date: 2006, May 31
C
C Coded by: Karel Zacek
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     E-mail: zacek@karel.troja.mff.cuni.cz
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     TIMESEC='string'... Name of the input file containing the time
C             section.
C             Default: TIMESEC='timesec.out' (usually sufficient)
C     FGBY22='string'... Name of the input file with the initial
C             parameters of Gaussian beams - imaginary part.
C             Default: FGBY22='fgby22.out' (usually sufficient)
C     FGBR22='string'... Name of the input file with the initial
C             parameters of Gaussian beams - real part.
C             Default: FGBR22='fgbr22.out' (usually sufficient)
C Name of input parameters:
C     N44=real... ????????????????????.
C             Default: N44=?
C     FMIN=real... ????????????????????.
C             Default: FMIN=0.
C     FMAX=real... ????????????????????.
C             Default: FMAX=100.
C Name of the output file:
C     GPOUT='string'... Name of the output file containing
C             amplitudes of Gaussian packets.
C             Default: GPOUT='gpampl.out'
C Data specifying the dimensions of the input grid:
C     OT=real, OX=real, OP=real ... Coordinates of the origin
C             of the grid.
C             Default: OT=0., OX=0., OP=0.
C     DT=real... Grid interval along the T axis.
C             Default: DT=1.
C     DX=real... Grid interval along the X axis.
C             Default: DX=1.
C     DP=real... Grid interval along the P axis.
C             Default: DP=1.
C     NT=positive integer... Number of gridpoints along the T axis.
C             Default: NT=1
C     NX=positive integer... Number of gridpoints along the X axis.
C             Default: NX=1
C     NP=positive integer... Number of gridpoints along the P axis.
C             Default: NP=1
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 Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc

      INCLUDE 'sep.inc'
C     sep.inc

      INCLUDE 'model.inc'
C     model.inc

      INCLUDE 'auxmod.inc'
C     auxmod.inc

C-----------------------------------------------------------------------

C     Auxiliary storage locations:
      REAL DX,DT,FMAX,FMIN,OMEGAMAX,VSURF,FREF,OMEGAREF

      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,N44,N44D
      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,RAYPAR,MODEL,SRC,CSG
      EQUIVALENCE (IRAM,RAM)

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     Reading all the data from the SEP file into the memory:
      CALL RSEP1(LU1,FILSEP)

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      CALL RSEP3T('MODEL',MODEL,'model2.dat')

C     Reading filenames of the output files
      CALL RSEP3T('GPPAR',GPPAR,'gpm-par.out')
      CALL RSEP3T('GPGRID',GPGRID,'gpgrid.out')
      CALL RSEP3T('RAYPAR',RAYPAR,'gpm-rp.dat')
      CALL RSEP3T('SRC',SRC,'gpm-src.dat')
      CALL RSEP3T('CSG',CSG,'gpm-csg.dat')

      CALL RSEP3I('NT',NT,750)
      CALL RSEP3I('NX',NX,116)
      CALL RSEP3R('DT',DT,0.004)
      CALL RSEP3R('DX',DX,25.0)

      CALL RSEP3I('SHOT',ISHOT,1)

      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.)

      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')

      OPEN(LU7,FILE=GPPAR,FORM='FORMATTED')
      OPEN(LU8,FILE=GPGRID,FORM='FORMATTED')
      OPEN(LU9,FILE=RAYPAR,FORM='FORMATTED')
      OPEN(LU10,FILE=SRC,FORM='FORMATTED')
      OPEN(LU11,FILE=CSG,FORM='FORMATTED')

      READ(LU2,*) ANIMIN
      READ(LU3,*) ANIMAX
      READ(LU4,*) ANRMIN
      READ(LU5,*) ANRMAX

      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      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)

      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

      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      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

      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      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

      OXR=150.0+ISHOT*25.0
      OTR=0.0

      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      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)')
     *             '''SRC'' ', 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)') '/'

      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)') '/'

      CLOSE(LU2)
      CLOSE(LU3)
      CLOSE(LU4)
      CLOSE(LU5)
      CLOSE(LU6)
      CLOSE(LU7)
      CLOSE(LU8)
      CLOSE(LU9)
      WRITE(*,'(A)') '+GPSTEP: Done. '
      STOP
      END

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.forC
C
C=======================================================================
C