C
C Program GPSTEP to calculate the discretization steps
C for (a) the decomposition of the wavefield into optimized
C 2-D Gaussian packets and (b) migration.
C
C Version: 6.20
C Date: 2008, 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     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     RYMM='string'... Name of the input file containing the
C             minimum and maximal values of the optimized and smoothed
C             initial parameters of Gaussian beams R0MIN,R0MAX,
C             Y0MIN,Y0MAX.
C             Default: RYMM='ryminmax.out'
C Input SEP parameters:
C     NX=positive integer... Number of gridpoints along the profile
C             axis X of the decomposed wavefield.
C             Default: NX=1
C     NT=positive integer... Number of gridpoints along the time axis T
C             of the decomposed wavefield.
C             Default: NT=1
C     DX=real... Grid interval along the profile axis X of
C             the decomposed wavefield.
C             Default: DX=1.
C     DT=real... Grid interval along the time axis T of the decomposed
C             wavefield.
C             Default: DT=1.
C     OT=real, OX=real ... Coordinates of the origin of the grid.
C             Default: OT=0., OX=0.
C     ISRC=positive integer...
C             Default: ISRC=1
C     VSURF=real... Velocity determined according to maximum velocity
C             in the depth of the receivers.
C             Default: VSURF=2000.
C     FMIN=real... Minimum frequency of the Gaussian packets.
C             Default: FMIN=1.
C     FMAX=real... Maximum frequency of the Gaussian packets.
C             Default: FMAX=100.
C     FREF=real... Reference frequencies of the Gaussian packets.
C             Default: FREF=50.
C     PARKSI=real... Parameter ksi with values in the interval (0,1).
C             Default: PARKSI=0.5
C     PARNY=real... Parameters ny with values in the interval (0,1).
C             Default: PARNY=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 Name of the output file:
C     GPSTEP='string'... Name of the output 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='gpstep.out'
C     GBGRD='string'... Name of the output file containing new
C             grid parameters of the Hamiltonian hypersurface (XR,TR,P)
C             computed for smoothed and optimized Im(N0) and Re(N0).
C
C=======================================================================
C
C     External functions and subroutines:
      EXTERNAL LENGTH,RSEP1,RSEP3T,RSEP3I,RSEP3R,ERROR
      INTEGER  LENGTH
C
C     Filenames and parameters:
      CHARACTER*80 FILSEP,RYMM,GPSTEP,GBGRD
      INTEGER LU1,LU2
      REAL PI
      PARAMETER (LU1=1,LU2=2,PI=3.141592654)
C
C     Other variables:
      INTEGER N,NX,NT,NXR,NTR,NP,NOMEGA,N1,N2,N3,N1NEW,N2NEW,N3NEW,ISRC
      REAL DX,DT,FMAX,FMIN,FREF,OMEGAMAX,OMEGAREF,VSURF
      REAL DXR,OXR,DTR,OTR,DP,OP,DOMEGA,OOMEGA,D1,D2,D3,O1,O2,O3
      REAL R0MIN,R0MAX,Y0MIN,Y0MAX,PARKSI,PARNY,DSRC,OXRINI
      REAL ANIMIN,ANIMAX,ANRMIN,ANRMAX,AK0,AK0T,KAPPA
      COMPLEX AN0
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 filename of the input file
      CALL RSEP3T('RYMM',RYMM,'ryminmax.out')
C
C     Reading input parameters
      CALL RSEP3I('NX',NX,116)
      CALL RSEP3R('DX',DX,25.0)
      CALL RSEP3R('OX',OX,0.)
      CALL RSEP3I('NT',NT,750)
      CALL RSEP3R('DT',DT,0.004)
      CALL RSEP3R('OT',OT,0.)
C
      CALL RSEP3I('ISRC',ISRC,1)
      CALL RSEP3R('DSRC',DSRC,1.)
      CALL RSEP3R('OXRINI',OXRINI,1.)
C     Grid parameters of the Hamiltonian hypersurface (NTR,NP,NXR)
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3R('O1',O1,1.)
      CALL RSEP3R('O2',O2,1.)
      CALL RSEP3R('O3',O3,1.)
      CALL RSEP3R('D1',D1,1.)
      CALL RSEP3R('D2',D2,1.)
      CALL RSEP3R('D3',D3,1.)
C
      CALL RSEP3R('VSURF',VSURF,2000.)
      CALL RSEP3R('FMAX',FMAX,1.)
      CALL RSEP3R('FMIN',FMIN,1.)
      CALL RSEP3R('FREF',FREF,1.)
      CALL RSEP3R('PARKSI',PARKSI,0.5)
      CALL RSEP3R('PARNY',PARNY,0.5)
      CALL RSEP3R('KAPPA',KAPPA,1.253314137)
C
C     Reading filenames of the output files
      CALL RSEP3T('GPSTEP',GPSTEP,'gpstep.out')
      CALL RSEP3T('GBGRD',GBGRD,'gbgrd.out')
C
C     Read maximum and minimum values of the optimized and smoothed
C     initial parameters of Gaussian beams R0, Y0
      OPEN(LU2,FILE=RYMM,FORM='FORMATTED',STATUS='OLD')
      READ(LU2,*) R0MIN,R0MAX,Y0MIN,Y0MAX
      CLOSE(LU2)
      IF(Y0MIN.LT.0.) THEN
C       GPSTEP-02
        CALL ERROR('GPSTEP-02: Y0 contains negative values')
      ENDIF
C
C     N=R0+iY0
      ANRMIN=R0MIN
      ANRMAX=R0MAX
      ANIMIN=Y0MIN
      ANIMAX=Y0MAX
      WRITE(*,'(A,1(G16.6,1X))') 'ANRMIN =',ANRMIN
      WRITE(*,'(A,1(G16.6,1X))') 'ANRMAX =',ANRMAX
      WRITE(*,'(A,1(G16.6,1X))') 'ANIMIN =',ANIMIN
      WRITE(*,'(A,1(G16.6,1X))') 'ANIMAX =',ANIMAX
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     Calculation of parameters N0, K0 and ~K0, equations (33), (34)
C     and (35), (see Zacek, 2003)
      AN0=PARNY*2.0*CMPLX(0.0,ANIMIN)
      AK0=PARKSI*OOMEGA*ANIMIN*VSURF**2
      AK0T=AK0*(1.0-PARNY)/PARNY
C
C     Calculation of parameters N0, K0 and ~K0 according to ????
C     AN0=CMPLX(ANRMIN,ANIMIN)
C     AK0=0.5*OOMEGA*ANIMIN*VSURF**2
C     AK0T=AK0
C
      IF(OOMEGA.LT.OMEGAMAX) THEN
C       Computation of maximum discretization step DOMEGA according to
C       equation (30), (Zacek, 2003)
        DOMEGA=KAPPA*SQRT(2.0*AK0*AK0T/(AK0+AK0T))
        N=NINT((OMEGAMAX-OOMEGA)/DOMEGA)
        NOMEGA=N+2
        DOMEGA=(OMEGAMAX-OOMEGA)/(N+1)
      ELSE
        OOMEGA=OMEGAMAX
        DOMEGA=1.0
        NOMEGA=1
      END IF
C
C     Computation of maximum discretization step DXR according to
C     equation (38), (Zacek, 2003)
      DXR=KAPPA*SQRT(-2.0*(1.0-PARNY)/OMEGAREF*AIMAG(
     *   (CMPLX(1.0,0.0)-OOMEGA/OMEGAREF*
     *    CMPLX(0.0,ANIMIN)/CMPLX(0.0,ANIMAX)*PARKSI)/
     *    CMPLX(0.0,ANIMAX)))
C
C     Computation of DXR according to ?????
C     DXR=KAPPA*SQRT(-AIMAG(0.5/(OMEGAREF*AN0)))
C
      N=NINT((NX-1)*DX/DXR)
      NXR=N+2
      DXR=(NX-1)*DX/(N+1)
C
C     Computation of maximum discretization step DP according to
C     equation (29), (Zacek, 2003)
      DP=KAPPA*SQRT(AIMAG(AN0)/OMEGAREF)
C
      N=NINT(2.0/(VSURF*DP))
      NP=N+2
      DP=2.0/(VSURF*(N+1))
C
C     Computation of maximum discretization step DTR according to
C     equation (40), (Zacek, 2003)
      DTR=KAPPA*SQRT(-2.0*PARNY*AIMAG(
     *   (CMPLX(1.0,0.0)-CMPLX(0.0,ANIMIN)/CMPLX(0.0,ANIMIN)*PARKSI)
     *   /CMPLX(0.0,AK0)))
C
C     Computation of DTR according to ?????
C     DTR=KAPPA*SQRT(0.5/AK0)
C
      N=NINT((NT-1)*DT/DTR)
      NTR=N+2
      DTR=(NT-1)*DT/(N+1)
C
      OXR=OXRINI+(ISRC-1)*DSRC
C
C     OXR=1587.5-FLOAT(NXR-1)*DXR/2.+ISRC*25.0
C     OXR=OX+FLOAT(NX-1)*DX/2.+FLOAT(ISRC-1)*DSRC
      OP=-1.0/VSURF
      OTR=0.0
C
C     Computing new grid parameters of the Hamiltonian hypersurface
      N1NEW=NTR
      N2NEW=NP
      N3NEW=NXR
      D1=FLOAT(N1-1)/FLOAT(N1NEW-1)
      D2=FLOAT(N2-1)/FLOAT(N2NEW-1)
      D3=FLOAT(N3-1)/FLOAT(N3NEW-1)
      OPEN(LU2,FILE=GBGRD,FORM='FORMATTED')
      CALL WSEP3I(2,'N1',N1NEW)
      CALL WSEP3I(2,'N2',N2NEW)
      CALL WSEP3I(2,'N3',N3NEW)
      CALL WSEP3R(2,'D1',D1)
      CALL WSEP3R(2,'D2',D2)
      CALL WSEP3R(2,'D3',D3)
      CALL WSEP3R(2,'O1',O1)
      CALL WSEP3R(2,'O2',O2)
      CALL WSEP3R(2,'O3',O3)
      CLOSE(LU2)
C
C     Writing discretization steps to the output file GPSTEP
      WRITE(*,'(A)')'+GPSTEP: Writing output file.'
      OPEN(LU2,FILE=GPSTEP,FORM='FORMATTED')
      WRITE(LU2,'(3A)') '''',GPSTEP(1:LENGTH(GPSTEP)),''''
      WRITE(LU2,*) NXR,DXR,OXR,' / (NXR,DXR,OXR)'
      WRITE(LU2,*) NTR,DTR,OTR,' / (NTR,DTR,OTR)'
      WRITE(LU2,*) NP, DP, OP ,' / (NP,DP,OP)'
      WRITE(LU2,*) NOMEGA,DOMEGA,OOMEGA,' / (NOMEGA,DOMEGA,OOMEGA)'
C
C     Writing complex-valued constants to the output file GPSTEP
      WRITE(LU2,*) REAL(AN0),AIMAG(AN0),' / (Re(N0),Im(N0))'
      WRITE(LU2,*) AK0,AK0T,ANIMIN,' / (K0,~K0,Im(N)min)'
      CLOSE(LU2)
C
C     Writing discretization steps to the screen
      WRITE(*,'(A,I6)')          ' NXR   =',NXR
      WRITE(*,'(A,1(G16.6,1X))') ' DXR   =',DXR
      WRITE(*,'(A,1(G16.6,1X))') ' OXR   =',OXR
      WRITE(*,'(A,I6)')          ' NTR   =',NTR
      WRITE(*,'(A,1(G16.6,1X))') ' DTR   =',DTR
      WRITE(*,'(A,1(G16.6,1X))') ' OTR   =',OTR
      WRITE(*,'(A,I6)')          ' NP    =',NP
      WRITE(*,'(A,1(G16.6,1X))') ' DP    =',DP
      WRITE(*,'(A,1(G16.6,1X))') ' OP    =',OP
      WRITE(*,'(A,I6)')          ' NOMEGA=',NOMEGA
      WRITE(*,'(A,1(G16.6,1X))') ' DOMEGA=',DOMEGA
      WRITE(*,'(A,1(G16.6,1X))') ' OOMEGA=',OOMEGA
C
C     Writing complex-valued constants to the screen
      WRITE(*,'(A,1(G16.6,1X))') ' Re(N0)  =',REAL(AN0)
      WRITE(*,'(A,1(G16.6,1X))') ' Im(N0)  =',AIMAG(AN0)
      WRITE(*,'(A,1(G16.6,1X))') ' K0      =',AK0
      WRITE(*,'(A,1(G16.6,1X))') ' ~K0     =',AK0T
      WRITE(*,'(A,1(G16.6,1X))') ' Im(N)min=',ANIMIN
C
      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
C
C=======================================================================
C