C 
C Program GPANAL to calculate the amplitudes of 2-D Gaussian packets
C
C Version: 8.00
C Date: 2022, April 24
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.......................................................................
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     SS='string'... Name of the input file with synthetic
C             seismograms to be decomposed.
C             Default: SS='ss.out'
C     FGBR22='string'... Name of the input file with the resampled
C             smoothed optimized parameters of Gaussian beams Re(N0).
C             Default: FGBR22=' '
C     FGBY22='string'... Name of the input file with the resampled
C             smoothed optimized parameters of Gaussian beams Im(N0).
C             Default: FGBY22=' '
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 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     PKSI=real... Parameter ksi with values in the interval (0,1).
C             Default: PKSI=0.5
C     PNY=real... Parameters ny with values in the interval (0,1).
C             Default: PNY=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     GWM=real.. Gaussian beam decay at the halfwidth is EXP(-GWM*GWM/2)
C             Default: GWM=3.
C     RELAMP=real... Relative decay of the Gaussian envelope at which
C             the loop over the points of the seismograms is terminated.
C             The relative error due to this economizing roughly
C             corresponds to the value of RELAMP.
C             Default: RELAMP=0.1
C Names of the output files:
C     GPAMPR='string'... Name of the output file containing
C             real part of Gaussian packet amplitudes.
C             Default: GPAMPR='gpampr.out'
C     GPAMPI='string'... Name of the output file containing
C             imaginary part of Gaussian packet amplitudes.
C             Default: GPAMPI='gpampi.out'
C     SSC='string'... Name of the output file with composed synthetic
C             seismograms.
C             Default: SSC='ssc.out'
C Optional parameters specifying the form of the quantities
C written in the output formatted files:
C     MINDIG,MAXDIG=positive integers ... See the description in file
C             
C             forms.for.
C     NUMLIN=positive integer ... Number of reals or integers in one
C             line of the output file. See the description in file
C             
C             forms.for.
C Form of the input/output files:
C     FORM='string'... Form of the input/output files: either
C             'FORMATTED' or 'UNFORMATTED' (case insensitive).
C             Default: FORM='FORMATTED'
C     FORMW='string'... Form of the output file. If both FORM and FORMW
C             are specified, FORMW is used for output files.
C             Default: FORMW=FORM
C     FORMR='string'... Form of the input file. If both FORM and FORMR
C             are specified, FORMR is used for output files.
C             Default: FORMR=FORM
C Value of undefined quantities:
C     UNDEF=real... The value to be used for undefined real quantities.
C             Default: UNDEF=undefined value used in forms.for
C
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C.......................................................................
C
C     External functions and subroutines:
      EXTERNAL OARRAY,RARRAY,RARRAI,ERROR,RSEP1,RSEP3T,RSEP3I,RSEP3R
C
C     Filenames and parameters:
      CHARACTER*80 FILSEP,SS,SSC,FGBY22,FGBR22,GPAMPR,GPAMPI,GPSTEP
      INTEGER LU1
      REAL PI
      PARAMETER (LU1=1,PI=3.141592654)
C
C     Other variables:
      INTEGER IRAM(MRAM)
      INTEGER NX,NT,NXT,NXTP,NXTPO
      INTEGER IXR,ITR,IP,IOMEGA,IX,IT,IT0
      INTEGER IGRD,IGRD0,IFGBY0,IFGBR0,IFGBY,IFGBR,IFGB
      INTEGER IGRDFR,IGRDFI,IGRDF0,IGRDF1,IGRDF2,IGRDF3,IXTPM0
      INTEGER NXR,NTR,NP,NOMEGA,NGF
      INTEGER IXMAX,IXMIN,ITMAX,ITMIN,IXR0,IP0,MXR,MP
      REAL OX,OT,DX,DT,XX,T
      REAL DXR,DTR,DP,DOMEGA,OXR,OTR,OP,OOMEGA,XRX,TRT,AX,BT
      REAL XR,TR,P,OMEGA,OMEREF,FR,FI,W,WR,WI,GWM,DXTPO
      REAL DXT,PMAX,ANIMIN,Y0MIN,PKSI,PNY,KAPPA
      REAL AK0,AK0T,AMPLR,AMPLI,AN0R,AN0I
      REAL DXR1,DTR1,DP1,DXRMIN,DTRMIN,DPMIN
      COMPLEX CN,CNT,CN0,CAMPL,CN44,CN44T,CPT
      COMPLEX CK11,CK12,CK22,CKDET,CKI11,CKI12,CKI22
      REAL THETAR,THETAI
      COMPLEX THETA1,THETA4,CNT1,CPT4
      REAL THET2R,THET2I,THET3R,THET3I,THET4R,THET4I
      REAL CNT1R,CNT1I,CPT4R,CPT4I,AUX,COSWW,SINWW
      REAL RELAMP,TMIN,TMAX,WR0,WI0,WR3,WI3,WR4,WI4
      CHARACTER*80 TEXT
      EQUIVALENCE (IRAM,RAM)
C
C.......................................................................
C
C     Reading main input data:
      FILSEP=' '
      WRITE(*,'(A)')'+GPANAL: Enter filename : '
      READ(*,*) FILSEP
      IF(FILSEP.EQ.' ') THEN
C       GPANAL-01
        CALL ERROR('GPANAL-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)')'+GPANAL: 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('SS',SS,'ss.out')
      CALL RSEP3T('FGBY22',FGBY22,'y0m-res.out')
      CALL RSEP3T('FGBR22',FGBR22,'r0m-res.out')
      CALL RSEP3T('GPSTEP',GPSTEP,'gpstep.out')
C
C     Reading input parameters
      CALL RSEP3I('NX',NX,116)
      CALL RSEP3R('DX',DX,25.)
      CALL RSEP3R('OX',OX,0.)
      CALL RSEP3I('NT',NT,750)
      CALL RSEP3R('DT',DT,0.1)
      CALL RSEP3R('OT',OT,0.)
      CALL RSEP3R('PKSI',PKSI,0.5)
      CALL RSEP3R('PNY',PNY,0.5)
      CALL RSEP3R('KAPPA',KAPPA,1.253314137)
      CALL RSEP3R('GWM',GWM,3.)
      CALL RSEP3R('RELAMP',RELAMP,0.1)
      RELAMP=ALOG(RELAMP)
C
C     Reading filenames of the output files
      CALL RSEP3T('GPAMPR',GPAMPR,'gpampr.out')
      CALL RSEP3T('GPAMPI',GPAMPI,'gpampi.out')
      CALL RSEP3T('SSC',SSC,'ssc.out')
      CLOSE(LU1)
C
C     Reading discretization steps
      OPEN(LU1,FILE=GPSTEP,FORM='FORMATTED',STATUS='OLD')
      READ(LU1,'(A)') TEXT
      READ(LU1,*) NXR,DXR,OXR
      READ(LU1,*) NTR,DTR,OTR
      READ(LU1,*) NP, DP, OP
      READ(LU1,*) NOMEGA,DOMEGA,OOMEGA
C     Reading complex-valued constants
      READ(LU1,*) AN0R,AN0I
      READ(LU1,*) AK0,AK0T,Y0MIN
      CLOSE(LU1)
C
      CN0=CMPLX(AN0R,AN0I)
      ANIMIN=Y0MIN
C
C     Writing discretization steps to the screen
      WRITE(*,*) ' '
      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)')          ' NXR   =',NXR
      WRITE(*,'(A,1(G16.6,1X))') ' DXR   =',DXR
      WRITE(*,'(A,1(G16.6,1X))') ' OXR   =',OXR
      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))') ' N0R   =',REAL(CN0)
      WRITE(*,'(A,1(G16.6,1X))') ' N0I   =',AIMAG(CN0)
      WRITE(*,'(A,1(G16.6,1X))') ' K0    =',AK0
      WRITE(*,'(A,1(G16.6,1X))') ' ~K0   =',AK0T
      WRITE(*,'(A,1(G16.6,1X))') ' ANIMIN=',ANIMIN
C
C     Writing constants and grid values to the screen
      WRITE(*,'(A,1(G16.6,1X))') ' PKSI=',PKSI
      WRITE(*,'(A,1(G16.6,1X))') ' PNY =',PNY
      WRITE(*,'(A,1(G16.6,1X))') ' KAPPA =',KAPPA
      WRITE(*,'(A,I6)')          ' NX    =',NX
      WRITE(*,'(A,1(G16.6,1X))') ' DX    =',DX
      WRITE(*,'(A,1(G16.6,1X))') ' OX    =',OX
      WRITE(*,'(A,I6)')          ' NT    =',NT
      WRITE(*,'(A,1(G16.6,1X))') ' DT    =',DT
      WRITE(*,'(A,1(G16.6,1X))') ' OT    =',OT
      WRITE(*,*) ' '
C
C     RAM array dimensions and indices
      NXT=NX*NT
      NXTP=NXR*NTR*NP
      NXTPO=NXR*NTR*NP*NOMEGA
      DXTPO=DXR*DTR*DP*DOMEGA
      DXT=DX*DT
      IFGBY0=NXT
      IFGBR0=IFGBY0+NXTP
      IGRDF0=NXT+2*NXTP
      IXTPM0=NXT+2*NXTP+2*NXTPO
      WRITE(*,*) IXTPM0,NXTPO,NXTP,NXT
C
C     Read synthetic seismograms (time section) to be decomposed
      CALL OARRAY(LU1,SS,'R')
      CALL RARRAY(LU1,' ',.FALSE.,0.,NXT,RAM)
      CLOSE(LU1)
C
C     Read initial parameters of Gaussian beams - imaginary part (Y0)
      CALL OARRAY(LU1,FGBY22,'R')
      CALL RARRAY(LU1,' ',.FALSE.,0.,NXTP,
     *            RAM(NXT+1))
      CLOSE(LU1)
C
C     Read initial parameters of Gaussian beams - real part (R0)
      CALL OARRAY(LU1,FGBR22,'R')
      CALL RARRAY(LU1,' ',.FALSE.,0.,NXTP,
     *            RAM(NXT+NXTP+1))
      CLOSE(LU1)
C
      IF((IGRDF0+NXTPO).GT.MRAM) THEN
C       GPANAL-02
        CALL ERROR('GPANAL-02: Too small array RAM(MRAM)')
C       Too small array RAM(MRAM).  If possible, increase dimension MRAM
C       in include file
C       ram.inc.
      END IF
C
C-----------------------------------------------------------------------
      OPEN(2,FILE='steps1.out')
      DXRMIN=999999.
      DTRMIN=999999.
      DPMIN=999999.
C
C     Computation of minimum values of DXR1,DTR1,DP1
C
C     Loop over positive circular frequency of Gaussian packet
      DO 35 IOMEGA=1,NOMEGA
        WRITE(*,'(A,I6)') '      IOMEGA   =',IOMEGA
        OMEGA=OOMEGA+DOMEGA*(IOMEGA-1)
C
C       Loop over coordinate of intersection of central ray
C       of Gaussian packet with the profile
        DO 50 IXR=1,NXR
          XR=OXR+DXR*(IXR-1)
          IGRDF3=(IXR-1)*NTR
C
C         Loop over corresponding arrival time of Gaussian packet
          DO 45 ITR=1,NTR
            TR=OTR+DTR*(ITR-1)
            IGRDF2=(IGRDF3+ITR-1)*NP
C
C           Loop over component of slowness vector of Gaussian packet
C           along profile
            DO 40 IP=1,NP
              P=OP+DP*(IP-1)
              IGRDF1=IGRDF0+(IGRDF2+IP-1)*NOMEGA
              IFGB=((IXR-1)*NP+IP-1)*NTR+ITR
              IFGBY=IFGBY0+IFGB
              IFGBR=IFGBR0+IFGB
C
C             Complex valued parameter N determines initial shape
C             of Gaussian packets
              CN=CMPLX(RAM(IFGBR),RAM(IFGBY))
C
C             Computation of ~N, equation (16),(Zacek,2003)
C                                         (28),(Zacek,2006)
              CNT=CN*CN0/(2.*CN-CN0)
C
C             Computation of ~p, equation (17),(Zacek,2003)
C                                         (29),(Zacek,2006)
              CPT=-P*CNT/CN
C
                PMAX=PI/(OMEGA*DX)
                IF((P.GT.(-PMAX)).AND.(P.LT.PMAX)) THEN
C
C                 Computation of ~N44, equation (18),(Zacek,2003)
C                                               (30),(Zacek,2006)
                  CN44T=CNT*CMPLX(0.,AK0T)
     *                /(OMEGA*CNT-CMPLX(0.,AK0T)*CPT**2)
C
C                 Computation of N44, equation (19),(Zacek,2003)
C                                              (31),(Zacek,2006)
                  CN44=CN*CMPLX(0.,AK0)
     *               /(OMEGA*CN-CMPLX(0.,AK0)*P**2)
C
C                 Computation of equations (6,10),(Zacek,2003)
C                                          (14,17),(Zacek,2006)
                  CK11=CN+CN44*P*P+CNT+CN44T*CPT*CPT
                  CK12=-CN44*P-CN44T*CPT
                  CK22=CN44+CN44T
                  CKDET=CK11*CK22-CK12*CK12
                  IF(CKDET.EQ.0.) THEN
                    CALL ERROR('GPSTEP-XX: CKDET equal zero')
                  END IF
                  CKI11=CK22/CKDET
                  CKI12=-CK12/CKDET
                  CKI22=CK11/CKDET
C
C                 Computation of equations (25,26),(Zacek,2003)
C                                          (45,46),(Zacek,2006)
                  DXR1=KAPPA*SQRT(-AIMAG(CKI11)*2./OMEGA)
                  DTR1=KAPPA*SQRT(-AIMAG(CKI22)*2./OMEGA)
C
C                 Computation of discretization step DP1 according to
C                 equation (29), (Zacek, 2003)
                  DP1=KAPPA*SQRT(AIMAG(CN0)/OMEGA)
C
                  DXRMIN=AMIN1(DXRMIN,DXR1)
                  DTRMIN=AMIN1(DTRMIN,DTR1)
                  DPMIN=AMIN1(DPMIN,DP1)
                END IF
   40       CONTINUE
   45     CONTINUE
   50   CONTINUE
        WRITE(*,*) DXRMIN,DTRMIN,DPMIN
        WRITE(2,*) DXRMIN,DTRMIN,DPMIN
        WRITE(*,*) IXTPM0,IOMEGA
        RAM(IXTPM0+1+IOMEGA)=DXRMIN
        RAM(IXTPM0+1+IOMEGA+NOMEGA+1)=DPMIN
   35 CONTINUE
      CLOSE(2)
C
C-----------------------------------------------------------------------
C     Decomposition of the time section into 2-D Gaussian packets
C
      WRITE(*,'(A)')'+GPANAL: Decomposition '
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
C     Omission of Gabor functions for lower frequencies
*     OMEREF=2.0*PI*15.0
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
      BT=GWM/SQRT(AK0T)
      OXR=OX
      WRITE(*,'(A,F7.2)') ' OXR=',OXR
C
C     Loop over coordinate of intersection of central ray
C     of Gaussian packet with the profile
      DO 500 IXR=1,NXR
        WRITE(*,'(2(A,I3))') '+GPANAL: Decomposition NXR=',IXR,'/',NXR
        XR=OXR+DXR*(IXR-1)
        IGRDF3=(IXR-1)*NTR
C
C       Loop over corresponding arrival time of Gaussian packet
        DO 450 ITR=1,NTR
          TR=OTR+DTR*(ITR-1)
          IGRDF2=(IGRDF3+ITR-1)*NP
C
C         Loop over component of slowness vector of Gaussian packet
C         along profile
          DO 400 IP=1,NP
            P=OP+DP*(IP-1)
            IGRDF1=IGRDF0+(IGRDF2+IP-1)*NOMEGA
            IFGB=((IXR-1)*NP+IP-1)*NTR+ITR
            IFGBY=IFGBY0+IFGB
            IFGBR=IFGBR0+IFGB
C
C           Complex valued parameter N determines initial shape
C           of Gaussian packets
            CN=CMPLX(RAM(IFGBR),RAM(IFGBY))
C
C           Computation of ~N, equation (16),(Zacek,2003)
C                                       (28),(Zacek,2006)
            CNT=CN*CN0/(2.*CN-CN0)
C
C           Computation of ~p, equation (17),(Zacek,2003)
C                                       (29),(Zacek,2006)
            CPT=-P*CNT/CN
C
C           Loop over positive circular frequency of Gaussian packet
            DO 350 IOMEGA=1,NOMEGA
C             WRITE(*,'(A,I6)') '      IOMEGA   =',IOMEGA
              OMEGA=OOMEGA+DOMEGA*(IOMEGA-1)
              IGRDFR=IGRDF1+IOMEGA
              IGRDFI=IGRDFR+NXTPO
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
C             Omission of Gabor functions for lower frequencies
*             NGF=MAX(INT(SQRT(OMEREF/OMEGA)),1)
*             WRITE(*,'(A,I6)') '      NGF        =',NGF
*             IF(MOD(IXR,NGF).NE.0.OR.MOD(ITR,NGF).NE.0.) THEN
*               RAM(IGRDFR)=0.
*               RAM(IGRDFI)=0.
*             ELSE
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
C             Computation of ~N44, equation (18),(Zacek,2003)
C                                           (30),(Zacek,2006)
              CN44T=CNT*CMPLX(0.,AK0T)
     *            /(OMEGA*CNT-CMPLX(0.,AK0T)*CPT**2)
C
C             Computation of N44, equation (19),(Zacek,2003)
C                                          (31),(Zacek,2006)
              CN44=CN*CMPLX(0.,AK0)
     *           /(OMEGA*CN-CMPLX(0.,AK0)*P**2)
C
C             Computation of equations (6,10),(Zacek,2003)
C                                      (14,17),(Zacek,2006)
C             CK11=CN+CN44*P*P+CNT+CN44T*CPT*CPT
C             CK12=-CN44*P-CN44T*CPT
C             CK22=CN44+CN44T
C             CKDET=CK11*CK22-CK12*CK12
C             IF(CKDET.EQ.0.) THEN
C               CALL ERROR('GPANAL-XX: CKDET equal zero')
C             END IF
C             CKI11=CK22/CKDET
C
C             Computation of equations (25),(Zacek,2003)
C                                      (45),(Zacek,2006)
C             DXR1=KAPPA*SQRT(-AIMAG(CKI11)*2./OMEGA)
              DXR1=RAM(IXTPM0+1+IOMEGA)
C             Computation of discretization step DXR1 according to
C             equation (38), (Zacek, 2003)
C             DXR1=KAPPA*SQRT(-2.*(1.-PNY)/OMEGA*AIMAG(
C    *            (CMPLX(1.,0.)-OOMEGA/OMEGA
C    *            *CMPLX(0.,ANIMIN)/CN*PKSI)/CN))
              MXR=INT(DXR1/DXR)
              IF(MXR.LT.1) THEN
                MXR=1
              END IF
              IXR0=INT(FLOAT(NXR-MXR*INT(FLOAT(NXR-1)
     *            /FLOAT(MXR))-1)/2.)+1
C
C             Computation of discretization step DP1 according to
C             equation (29), (Zacek, 2003)
C             DP1=KAPPA*SQRT(AIMAG(CN0)/OMEGA)
              DP1=RAM(IXTPM0+1+IOMEGA+NOMEGA+1)
              MP=INT(DP1/DP)
              IF(MP.LT.1) THEN
                MP=1
              END IF
              IP0=INT(FLOAT(NP-MP*INT(FLOAT(NP-1)
     *            /FLOAT(MP))-1)/2.)+1
C             IF((MOD((IXR-IXR0),MXR).NE.0).OR.
C    *           (MOD((IP-IP0),MP).NE.0)) THEN
C             WRITE(*,'(A,I6)') '      IOMEGA   =',IOMEGA
C               WRITE(*,*)DXR1,DXR,DP1,DP,MXR,MP,
C    *                  MOD((IXR-IXR0),MXR),MOD((IP-IP0),MP)
C             END IF
C
              IF((MOD((IXR-IXR0),MXR).EQ.0).AND.
     *           (MOD((IP-IP0),MP).EQ.0)) THEN
C
                FR=0.
                FI=0.
                PMAX=PI/(OMEGA*DX)
C
                IF((P.GT.(-PMAX)).AND.(P.LT.PMAX)) THEN
C
C                 Computation of N44, equation (19),(Zacek,2003)
C                                              (31),(Zacek,2006)
C                 CN44=CN*CMPLX(0.,AK0)
C    *               /(OMEGA*CN-CMPLX(0.,AK0)*P**2)
C
C                 Computation of ~N44, equation (18),(Zacek,2003)
C                                               (30),(Zacek,2006)
C                 CN44T=CNT*CMPLX(0.,AK0T)
C    *                /(OMEGA*CNT-CMPLX(0.,AK0T)*CPT**2)
C
C                 Computation of ~a according to equation (38),
C                 (Zacek, 2006)
                  CAMPL=0.5*CSQRT(CMPLX(0.,-1.)*(CN+CNT))
     *                *CSQRT(CMPLX(0.,-1.)*CN44*CN44T
     *                *(OMEGA/CMPLX(0.,AK0)+OMEGA/CMPLX(0.,AK0T)))
C
                  AMPLR=REAL(CAMPL)
                  AMPLI=AIMAG(CAMPL)
C
                  THETA1=0.5*CNT/(CNT-CMPLX(0.,AK0T)/OMEGA*CPT**2)
                  THETA4=THETA1*CMPLX(0.,AK0T)
                  THET4R=REAL (THETA4)
                  THET4I=AIMAG(THETA4)
                  CNT1=CNT*THETA1
                  CNT1R=OMEGA*REAL (CNT1)
                  CNT1I=OMEGA*AIMAG(CNT1)
                  CPT4=CPT*THETA4
                  CPT4R=REAL (CPT4)
                  CPT4I=AIMAG(CPT4)
C
C                 Gaussian beam horizontal halfwidth
                  AX=GWM/SQRT(OMEGA*AIMAG(CNT))
                  IXMIN=INT((XR-AX)/DX)+1
                  IF(IXMIN.LT.1) THEN
                    IXMIN=1
                  END IF
                  IXMAX=INT((XR+AX)/DX)+2
                  IF(IXMAX.GT.NX) THEN
                    IXMAX=NX
                  END IF
                  ITMIN=INT((TR-BT)/DT)+1
                  IF(ITMIN.LT.1) THEN
                    ITMIN=1
                  END IF
                  ITMAX=INT((TR+BT)/DT)+2
                  IF(ITMAX.GT.NT) THEN
                    ITMAX=NT
                  END IF
C
C                 Loop over traces of wavefield (common shot gather)
C                 to be decomposed
                  DO 300 IX=IXMIN,IXMAX
                    XX=OX+DX*(IX-1)
                    XRX=XX-XR
                    IGRD0=(IX-1)*NT
                    AUX=XRX**2
                    THET2R=AUX*CNT1R-OMEGA*P*XRX
                    THET2I=AUX*CNT1I
                    AUX=2.*XRX
                    THET3R=AUX*CPT4R-OMEGA
                    THET3I=AUX*CPT4I
                    AUX=0.5*THET3I/THET4I
                    TMAX=AUX*AUX-(RELAMP+THET2I)/THET4I
                    IF(TMAX.GT.0.) THEN
                      TMAX=SQRT(TMAX)
                      TMIN=(AUX-TMAX)
                      TMAX=(AUX+TMAX)
C                     ITMIN=MAX0(NINT((TMIN+TR-OT)/DT+1.5),ITMIN)
C                     ITMAX=MIN0(NINT((TMAX+TR-OT)/DT+0.5),ITMAX)
                      ITMIN=MAX0(NINT((TMIN+TR-OT)/DT+1.5),1)
                      ITMAX=MIN0(NINT((TMAX+TR-OT)/DT+0.5),NT)
                    ELSE
                      ITMIN=1
                      ITMAX=-1
                    END IF
C
C                   Loop over samples of each trace of wavefield
C                   (common shot gather) to be decomposed
                    IF(ITMIN.LT.ITMAX) THEN
                      IT=(ITMIN+ITMAX)/2
                      T=OT+DT*(IT-1)
                      TRT=T-TR
                      THETAR=THET2R-(THET3R-THET4R*TRT)*TRT
                      THETAI=THET2I-(THET3I-THET4I*TRT)*TRT
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
C                     Omission of Gabor functions for lower frequencies
*                     W=NGF*NGF*RAM(IGRD)*EXP(-OMEGA*AIMAG(THETA))
*                     W=NGF*NGF*EXP(-THETAI)
C
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                      W=EXP(-THETAI)
                      COSWW=W*COS(THETAR)
                      SINWW=W*SIN(THETAR)
                      WR=AMPLR*COSWW-AMPLI*SINWW
                      WI=AMPLR*SINWW+AMPLI*COSWW
                      IGRD=IGRD0+IT
                      AUX=RAM(IGRD)
                      FR=FR+WR*AUX
                      FI=FI+WI*AUX
                      IT0=IT
                      WR0=WR
                      WI0=WI
                      W=EXP( (THET3I-THET4I*(2.*TRT-DT))*DT)
                      AUX=  -(THET3R-THET4R*(2.*TRT-DT))*DT
                      WR3=W*COS(AUX)
                      WI3=W*SIN(AUX)
                      W=EXP(-2.*THET4I*DT*DT)
                      AUX=2.*THET4R*DT*DT
                      WR4=W*COS(AUX)
                      WI4=W*SIN(AUX)
C
                      DO 251 IGRD=IGRD0+IT0+1,IGRD0+ITMAX
                        AUX=WR3*WR4-WI3*WI4
                        WI3=WR3*WI4+WI3*WR4
                        WR3=AUX
                        AUX=WR*WR3-WI*WI3
                        WI =WR*WI3+WI*WR3
                        WR =AUX
C
                        AUX=RAM(IGRD)
                        FR=FR+WR*AUX
                        FI=FI+WI*AUX
  251                 CONTINUE
                      WR=WR0
                      WI=WI0
                      W=EXP(-(THET3I-THET4I*(2.*TRT+DT))*DT)
                      AUX=   (THET3R-THET4R*(2.*TRT+DT))*DT
                      WR3=W*COS(AUX)
                      WI3=W*SIN(AUX)
C
                      DO 252 IGRD=IGRD0+IT0-1,IGRD0+ITMIN,-1
                        AUX=WR3*WR4-WI3*WI4
                        WI3=WR3*WI4+WI3*WR4
                        WR3=AUX
                        AUX=WR*WR3-WI*WI3
                        WI =WR*WI3+WI*WR3
                        WR =AUX
C
                        AUX=RAM(IGRD)
                        FR=FR+WR*AUX
                        FI=FI+WI*AUX
  252                 CONTINUE
                    END IF
  300             CONTINUE
                END IF
C
C               Complex valued amplitudes of GPs
                RAM(IGRDFR)=MXR*MP*DXT*FR*OMEGA**2/(2.*PI**3)
                RAM(IGRDFI)=MXR*MP*DXT*FI*OMEGA**2/(2.*PI**3)
C               RAM(IGRDFR)=MP*DXT*FR*OMEGA**2/(2.*PI**3)
C               RAM(IGRDFI)=MP*DXT*FI*OMEGA**2/(2.*PI**3)
              ELSE
                RAM(IGRDFR)=0.
                RAM(IGRDFI)=0.
              END IF
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
C
C             Omission of Gabor functions for lower frequencies
*             END IF
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  350       CONTINUE
  400     CONTINUE
  450   CONTINUE
  500 CONTINUE
C     Writing output grid values
      IF(GPAMPR.NE.' ') THEN
        CALL OARRAY(LU1,GPAMPR,'W')
        CALL WARRAY(LU1,' ',.FALSE.,0.,.FALSE.,0.,
     *                                    NXTPO,RAM(IGRDF0+1))
        CLOSE(LU1)
      END IF
      IF(GPAMPI.NE.' ') THEN
        CALL OARRAY(LU1,GPAMPI,'W')
        CALL WARRAY(LU1,' ',.FALSE.,0.,.FALSE.,0.,
     *                              NXTPO,RAM(IGRDF0+NXTPO+1))
        CLOSE(LU1)
      END IF
C
      WRITE(*,'(A)') '+GPANAL: 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
C
C=======================================================================
C