C
C Program SGFGRD to calculate the grid values of a real-valued quantity
C decomposed into the structural Gabor functions
C
C Version: 6.20
C Date: 2008, June 10
C
C Coded by: Ludek Klimes
C     Department of Geophysics, Charles University Prague,
C     Ke Karlovu 3, 121 16 Praha 2, Czech Republic,
C     http://sw3d.cz/staff/klimes.htm
C
C.......................................................................
C
C Two real-valued coefficients stand in for the complex-valued
C coefficient of each odd Gabor function.  The odd real-valued
C coefficient is the real part, and the even real-valued
C coefficient is the imaginary part.  These two real-valued coefficients
C simultaneously represent the complex-valued coefficient of the
C successive even Gabor function because the two Gabor functions and
C the corresponding complex-valued coefficients are complex-conjugate.
C Each odd real-valued coefficient thus corresponds to twice the
C real part of the Gabor function of the same odd index, and each even
C real-valued coefficient corresponds to twice the imaginary
C part of the Gabor function of the same even index.
C
C.......................................................................
C                                                    
C Description of data files:
C
C Input data read from the standard input device (*):
C     The data are read by the list directed input (free format) and
C     consist of a single string 'SEP':
C     'SEP'...String in apostrophes containing the name of the input
C             SEP parameter or history file with the input data.
C     No default, 'SEP' must be specified and cannot be blank.
C
C                                                     
C Input data file 'SEP':
C     File 'SEP' has the form of the SEP
C     parameter file.  The parameters, which do not differ from their
C     defaults, need not be specified in file 'SEP'.
C Data specifying input and output files:
C     SGF='string'... Name of the input file with structural Gabor
C             functions.  Input parameters N1, N2, N3 determine whether
C             the structural Gabor functions are 1-D, 2-D or 3-D.
C             The wavemumber components and matrix elements
C             corresponding to the direction in which the number of
C             gridpoints is Ni=1 should equal zero.
C             Description of file SGF.
C             Default: SGF='sgf.out'
C     SGFAMP='string'... Name of the header file of the input
C             real-valued vector of the amplitudes of Gabor functions,
C             corresponding to a real-valued quantity to be gridded.
C             For general description of the files with matrices refer
C             to file forms.htm.
C             Default: SGFAMP='sgfamp.out'
C     GRD='string'... Name of the formatted output file with the gridded
C             real-valued quantity.  The file contains N1*N2*N3 values.
C             For general description of files with gridded data refer
C             to file forms.htm.
C             Default: GRD='sgfgrd.out'
C Data specifying dimensions of the input grid:
C     N1,N2,N3=integers... Numbers of gridpoints along the X1,X2,X3
C             axes, respectively.  These numbers also determine whether
C             the structural Gabor functions are 1-D, 2-D or 3-D.
C             Defaults: N1=1, N2=1, N3=1
C     O1,O2,O3=reals... Coordinates of the origin of the grid, i.e.,
C             of the first gridpoint.
C             Defaults: O1=0, O2=0, O3=0
C     D1,D2,D3=reals... Grid intervals along the X1,X2,X3 axes,
C             respectively.
C             Defaults: D1=1, D2=1, D3=1
C Numerical parameters:
C     RELAMP=positive real... Relative decay of the Gaussian envelope
C             at which the loop over the points of the input grid is
C             terminated.
C             The relative error due to this economizing roughly
C             corresponds to the value of RELAMP.
C             Default: RELAMP=0.001
C Form of the output files with matrices:
C     FORMM='string' ... Form of the output files with matrices. Allowed
C             values are FORMM='formatted' and FORMM='unformatted'.
C             Default: FORMM='formatted'
C Optional parameters specifying the form of the quantities
C written to the output formatted files:
C     MINDIG,MAXDIG=positive integers ... See the description in file
C             forms.for.
C     NUMLIN=positive integer ... Number of reals in one line of the
C             output file. See the description in file
C             forms.for.
C
C=======================================================================
C
C Common block /RAMC/:
      INCLUDE 'ram.inc'
C     ram.inc
C
C.......................................................................
C
C     External functions and subroutines:
      EXTERNAL ERROR,RSEP1,RSEP3T,RSEP3I,RSEP3R,RMATH,RMATD
C
C     Constants:
      REAL COEF1D,COEF2D,COEF3D
      PARAMETER (COEF1D=1.772453851)
      PARAMETER (COEF2D=3.141592654)
      PARAMETER (COEF3D=5.568327997)
C
C     Filenames and parameters:
      CHARACTER*80 FSEP,FSGF,FAMP,FDAT,FGRD
      INTEGER LU1
      PARAMETER (LU1=1)
C
C     Input data:
      INTEGER N1,N2,N3,N1N2N3,NQ,NZ,NSGF,NQNSGF,NDIM
      INTEGER NELEM,M1,M2
      REAL O1,O2,O3,D1,D2,D3,RELAMP,RELLOG,ZERO(12)
      CHARACTER*3 SPARSE
      CHARACTER*4 SYMM
      CHARACTER*11 FORM
      CHARACTER*1 TEXT
C
C     Gabor function b (beta):
      REAL BX1,BX2,BX3,BK1,BK2,BK3
      REAL BY11,BY12,BY22,BY13,BY23,BY33,BYDET
      REAL BR11,BR12,BR22,BR13,BR23,BR33,BAMPR,BAMPI
C     2-D projection of 3-D Gabor function:
      REAL AX1,AX2,AK1,AK2
      REAL AY13,AY23,AY33
      REAL AR33
C     Coordinate differences:
      REAL DX1,DX2,DX3
C     Calculation of a Gabor function:
      REAL EXPR,EXPI,C,S
C     Calculation of scalar product with the given grid:
      REAL EXP0R,EXP0I,EXP1R,EXP1I,EXP1MR,EXP1MI,EXP2R,EXP2I
      INTEGER J1,J2,J3,K1,K2,K3
C
C     Auxiliary variables:
      INTEGER ISGF,IQ
      INTEGER I1,I2,I3,I
      REAL DET,AUX
C
C.......................................................................
C
C     Reading main input data:
      WRITE(*,'(A)') '+SGFGRD: Enter input filename: '
      FSEP=' '
      READ (*,*) FSEP
      IF(FSEP.EQ.' ') THEN
C       SGFGRD-01
        CALL ERROR('SGFGRD-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.
      END IF
      CALL RSEP1(LU1,FSEP)
      WRITE(*,'(A)') '+SGFGRD: Working...            '
C
C     Reading input and output filenames:
      CALL RSEP3T('SGF'   ,FSGF,'sgf.out')
      CALL RSEP3T('SGFAMP',FAMP,'sgfamp.out')
      CALL RSEP3T('GRD'   ,FGRD,'sgfgrd.out')
      IF(FSGF.EQ.' ') THEN
C       SGFGRD-02
        CALL ERROR('SGFGRD-02: Blank name of input file SGF')
      END IF
      IF(FAMP.EQ.' ') THEN
C       SGFGRD-03
        CALL ERROR('SGFGRD-03: Blank name of input file AMP')
      END IF
      IF(FGRD.EQ.' ') THEN
C       SGFGRD-04
        CALL ERROR('SGFGRD-04: Blank name of input file GRD')
      END IF
C
C     Reading other input SEP parameters:
      CALL RSEP3I('N1',N1,1)
      CALL RSEP3I('N2',N2,1)
      CALL RSEP3I('N3',N3,1)
      CALL RSEP3R('O1',O1,0.0)
      CALL RSEP3R('O2',O2,0.0)
      CALL RSEP3R('O3',O3,0.0)
      CALL RSEP3R('D1',D1,1.0)
      CALL RSEP3R('D2',D2,1.0)
      CALL RSEP3R('D3',D3,1.0)
      N1N2N3=N1*N2*N3
      CALL RSEP3R('RELAMP',RELAMP,0.001)
      RELLOG=-ALOG(RELAMP)
C
C     Determination of NDIM:
C     1-D: NDIM=1,2,3
C     2-D: NDIM=4,5,6
C     3-D: NDIM=7
      NDIM=-1
      IF(N1.GT.1) NDIM=NDIM+2
      IF(N2.GT.1) NDIM=NDIM+3
      IF(N3.GT.1) NDIM=NDIM+4
      IF(NDIM.EQ.-1) THEN
C       SGFGRD-05
        CALL ERROR('SGFGRD-05: N1=N2=N3=1 is not allowed')
      END IF
      IF(NDIM.EQ.2) THEN
        N1=N2
        N2=1
        O1=O2
        D1=D2
      ELSE IF(NDIM.EQ.3) THEN
        N1=N3
        N3=1
        O1=O3
        D1=D3
      ELSE IF(NDIM.EQ.5) THEN
        N2=N3
        N3=1
        O2=O3
        D2=D3
      ELSE IF(NDIM.EQ.6) THEN
        N1=N2
        N2=N3
        N3=1
        O1=O2
        O2=O3
        D1=D2
        D2=D3
      ELSE IF(NDIM.EQ.8) THEN
        NDIM=7
      END IF
C
C     Determination of NQ and NZ:
C     NQ is the number of reals stored for each Gabor function
C     NZ is the number of input zeros for each Gabor function
      IF(NDIM.LE.3) THEN
C       1-D
        NQ=6
        NZ=12
      ELSE IF(NDIM.LE.6) THEN
C       2-D
        NQ=12
        NZ=7
      ELSE
C       3-D:
        NQ=20
        NZ=0
      END IF
C
C     Reading the input vector of the amplitudes of Gabor functions:
      CALL RMATH(LU1,FAMP,FDAT,M1,M2,SPARSE,NELEM,SYMM,FORM)
      IF(SPARSE.NE.' ') THEN
C       SGFGRD-06
        CALL ERROR('SGFGRD-06: Input vector is sparse')
      END IF
      IF(M2.NE.1.OR.SYMM.NE.' '.OR.NELEM.NE.M1) THEN
C       SGFGRD-07
        CALL ERROR('SGFGRD-07: Input vector is not a vector')
      END IF
      IF(NQ*(M1/2+1).GT.MRAM) THEN
C       SGFGRD-08
        CALL ERROR
     *        ('SGFGRD-08: Too small array RAM for Gabor functions')
C       The input parameters of Gabor functions do not fit into array
C       RAM(MRAM).
C       The number of parameters of each Gabor function is
C       1-D: NQ=6
C       2-D: NQ=12
C       3-D: NQ=20
      END IF
      CALL RMATD(LU1,FDAT,1,SPARSE,NELEM,FORM,1)
      DO 1 ISGF=M1/2,1,-1
        I1=2*ISGF
        I2=NQ*ISGF
        RAM(I2-1)=RAM(I1-1)
        RAM(I2)  =RAM(I1)
    1 CONTINUE
C
C     Reading the Gabor functions:
      OPEN(LU1,FILE=FSGF,FORM='FORMATTED',STATUS='OLD')
      READ(LU1,*) (TEXT,I=1,20)
      NQNSGF=0
   10 CONTINUE
        IF(NQNSGF+NQ.GT.MRAM) THEN
C         SGFGRD-09
          CALL ERROR
     *          ('SGFGRD-09: Too small array RAM for Gabor functions')
C         The input parameters of Gabor functions do not fit into array
C         RAM(MRAM).
C         The number of parameters of each Gabor function is
C         1-D: NQ=6
C         2-D: NQ=14
C         3-D: NQ=24
        END IF
        TEXT='$'
        GO TO (11,12,13,14,15,16,17) NDIM
C         SGFGRD-10
          CALL ERROR('SGFGRD-10: Wrong value of NDIM')
   11   CONTINUE
          READ(LU1,*,END=20) TEXT,RAM(NQNSGF+1),AUX,AUX,
     *            RAM(NQNSGF+2),ZERO(1),ZERO(2),
     *            RAM(NQNSGF+3),(ZERO(I),I=3,7),
     *            RAM(NQNSGF+4),(ZERO(I),I=8,12)
          GO TO 18
   12   CONTINUE
          READ(LU1,*,END=20) TEXT,AUX,RAM(NQNSGF+1),AUX,
     *            ZERO(1),RAM(NQNSGF+2),(ZERO(I),I=2,4),
     *            RAM(NQNSGF+3),(ZERO(I),I=5,9),
     *            RAM(NQNSGF+4),(ZERO(I),I=10,12)
          GO TO 18
   13   CONTINUE
          READ(LU1,*,END=20) TEXT,AUX,AUX,RAM(NQNSGF+1),
     *            ZERO(1),ZERO(2),RAM(NQNSGF+2),(ZERO(I),I=3,7),
     *            RAM(NQNSGF+3),(ZERO(I),I=8,12),
     *            RAM(NQNSGF+4)
          GO TO 18
   14   CONTINUE
          READ(LU1,*,END=20) TEXT,RAM(NQNSGF+1),RAM(NQNSGF+2),AUX,
     *            RAM(NQNSGF+3),RAM(NQNSGF+4),ZERO(1),
     *            RAM(NQNSGF+5),RAM(NQNSGF+6),RAM(NQNSGF+7),
     *            (ZERO(I),I=2,4),
     *            RAM(NQNSGF+8),RAM(NQNSGF+9),RAM(NQNSGF+10),
     *            (ZERO(I),I=5,7)
          GO TO 18
   15   CONTINUE
          READ(LU1,*,END=20) TEXT,RAM(NQNSGF+1),AUX,RAM(NQNSGF+2),
     *            RAM(NQNSGF+3),ZERO(1),RAM(NQNSGF+4),
     *            RAM(NQNSGF+5),ZERO(2),ZERO(3),
     *            RAM(NQNSGF+6),ZERO(4),RAM(NQNSGF+7),
     *            RAM(NQNSGF+8),ZERO(5),ZERO(6),
     *            RAM(NQNSGF+9),ZERO(7),RAM(NQNSGF+10)
          GO TO 18
   16   CONTINUE
          READ(LU1,*,END=20) TEXT,AUX,RAM(NQNSGF+1),RAM(NQNSGF+2),
     *            ZERO(1),RAM(NQNSGF+3),RAM(NQNSGF+4),
     *            ZERO(2),ZERO(3),RAM(NQNSGF+5),
     *            ZERO(4),RAM(NQNSGF+6),RAM(NQNSGF+7),
     *            ZERO(5),ZERO(6),RAM(NQNSGF+8),
     *            ZERO(7),RAM(NQNSGF+9),RAM(NQNSGF+10)
          GO TO 18
   17   CONTINUE
          READ(LU1,*,END=20) TEXT,(RAM(I),I=NQNSGF+1,NQNSGF+18)
          GO TO 18
   18   CONTINUE
        IF(TEXT.EQ.'$') GO TO 20
        NQNSGF=NQNSGF+NQ
        DO 19 I=1,NZ
          IF(ZERO(I).NE.0.0) THEN
C           SGFGRD-11
            CALL ERROR ('SGFGRD-11: Non-zero input quantity')
C           The input wavemumber components and matrix elements
C           corresponding to the direction in which the number of
C           gridpoints is Ni=1 should equal zero.
          END IF
   19   CONTINUE
      GO TO 10
   20 CONTINUE
      CLOSE(LU1)
      NSGF=NQNSGF/NQ
      IF(M1.NE.2*NSGF) THEN
C       SGFGRD-12
        CALL ERROR('SGFGRD-12: Wrong length of the input vector')
      END IF
C     NSGF is the number of odd Gabor functions.
C     RAM(1:NQNSGF) contain the parameters of odd Gabor functions.
C
C     Checking the positive definiteness of the real part of matrix K,
C     calculating the squares of maximum differences between coordinates
C     and wavenumber components, and halving the components of matrix K:
      IF(NDIM.LE.3) THEN
C       1-D:
        DO 32 I=0,NQNSGF-NQ,NQ
          IF(RAM(I+3).LE.0.0) THEN
C           SGFGRD-13
            CALL ERROR('SGFGRD-13: Indefinite real part of matrix K')
          END IF
C         Halving the components of matrix K
          DO 31 I1=I+3,I+4
            RAM(I1)=0.5*RAM(I1)
   31     CONTINUE
   32   CONTINUE
      ELSE IF(NDIM.LE.6) THEN
C       2-D:
        DO 42 I=0,NQNSGF-NQ,NQ
          DET=RAM(I+5)*RAM(I+7)-RAM(I+6)**2
          IF(RAM(I+5).LE.0.0.OR.DET.LE.0.0) THEN
C           SGFGRD-14
            CALL ERROR('SGFGRD-14: Indefinite real part of matrix K')
          END IF
C         Halving the components of matrix K
          DO 41 I1=I+5,I+10
            RAM(I1)=0.5*RAM(I1)
   41     CONTINUE
   42   CONTINUE
      ELSE
C       3-D:
        DO 52 I=0,NQNSGF-NQ,NQ
C         Matrix Y
          BY11=RAM(I+7)
          BY12=RAM(I+8)
          BY22=RAM(I+9)
          BY13=RAM(I+10)
          BY23=RAM(I+11)
          BY33=RAM(I+12)
C         Matrix R
          BR11=-RAM(I+13)
          BR12=-RAM(I+14)
          BR22=-RAM(I+15)
          BR13=-RAM(I+16)
          BR23=-RAM(I+17)
          BR33=-RAM(I+18)
C         Last column of inverse matrix to Y multiplied by det(Y)
          AY13=BY12*BY23-BY13*BY22
          AY23=BY12*BY13-BY23*BY11
          AY33=BY11*BY22-BY12*BY12
C         det(Y)
          BYDET=BY13*AY13+BY23*AY23+BY33*AY33
          IF(BY11.LE.0.0.OR.AY33.LE.0.0.OR.BYDET.LE.0.0) THEN
C           SGFGRD-15
            CALL ERROR('SGFGRD-15: Indefinite real part of matrix K')
          END IF
C         Halving the components of matrix K
          DO 51 I1=I+7,I+18
            RAM(I1)=0.5*RAM(I1)
   51     CONTINUE
   52   CONTINUE
      END IF
C
C     Check for the memory required for the calculation of the grid:
      IF(NQNSGF+N1N2N3.GT.MRAM) THEN
C       SGFGRD-16
        CALL ERROR('SGFGRD-16: Too small array RAM for the grid')
C       The input parameters of Gabor functions and the input grid do
C       not fit together into array
C       RAM(MRAM).
C       The number of grid values is N1*N2*N3.
C       The number of parameters of each Gabor function is
C       1-D: NQ=6
C       2-D: NQ=12
C       3-D: NQ=20
      END IF
C
C     Initializing the grid values:
      I=NQNSGF
      DO 103 I3=1,N3
        DO 102 I2=1,N2
          DO 101 I1=1,N1
            I=I+1
            RAM(I)=0.0
  101     CONTINUE
  102   CONTINUE
  103 CONTINUE
C
C.......................................................................
C
C     Calculating the grid values:
      IF(NDIM.LE.3) THEN
C       1-D:
        DO 214 ISGF=0,NSGF-1
          IF(MOD(ISGF,50).EQ.0) THEN
            WRITE(*,'(2(A,I7))')
     *        '+SGFGRD: Gabor function',2*ISGF,' /',2*NSGF
          END IF
          IQ=NQ*ISGF
          IF(RAM(IQ+5).NE.0.0.OR.RAM(IQ+5).NE.0.0) THEN
C           Quantities describing the Gabor function
            BX1=RAM(IQ+1)
            BK1=RAM(IQ+2)
C           Half matrix Y
            BY11=RAM(IQ+3)
C           Half matrix R
            BR11=-RAM(IQ+4)
C           Determinant of half matrix Y
            BYDET=BY11
C           Amplitude coefficients
            AUX=2.0*SQRT(SQRT(4.0*BYDET)/COEF2D)
            BAMPR=AUX*RAM(IQ+5)
            BAMPI=AUX*RAM(IQ+6)
C           Extent of the Gabor function along axis X1
            AUX=SQRT(RELLOG/BY11)
            J1=MAX0(INT((BX1-AUX)/D1+0.999),0)
            K1=MIN0(INT((BX1+AUX)/D1+0.001),N1-1)
            IF(J1.LE.K1) THEN
C             Index of the central point along the gridline
              M1=MAX0(J1,MIN0(NINT(BX1/D1),K1))
C             Relative coordinate of the central point
              DX1=O1+D1*FLOAT(M1)-BX1
C             Exponent at the central point
              EXP0R=DX1*BY11*DX1
              EXP0I=DX1*(BK1+BR11*DX1)
C             The first derivative of the exponent
              EXP1R=BY11*DX1
              EXP1I=BR11*DX1
              EXP1R=D1*(EXP1R+EXP1R)
              EXP1I=D1*(EXP1I+EXP1I+BK1)
C             Half the second derivative of the exponent
              EXP2R=D1*BY11*D1
              EXP2I=D1*BR11*D1
C             Half the second derivative minus the first derivative
              EXP1MR=EXP2R-EXP1R
              EXP1MI=EXP2I-EXP1I
C             Half the second derivative plus the first derivative
              EXP1R=EXP2R+EXP1R
              EXP1I=EXP2I+EXP1I
C             Second derivative of the exponent
              EXP2R=EXP2R+EXP2R
              EXP2I=EXP2I+EXP2I
C             Exponential function at the central point
              AUX=EXP(-EXP0R)
              C=COS(EXP0I)
              S=SIN(EXP0I)
              EXP0R=C*AUX
              EXP0I=S*AUX
C             Exponential correction at the central point (+)
              AUX=EXP(-EXP1R)
              C=COS(EXP1I)
              S=SIN(EXP1I)
              EXP1R=C*AUX
              EXP1I=S*AUX
C             Exponential correction at the central point (-)
              AUX=EXP(-EXP1MR)
              C=COS(EXP1MI)
              S=SIN(EXP1MI)
              EXP1MR=C*AUX
              EXP1MI=S*AUX
C             Constant correction to the correction
              AUX=EXP(-EXP2R)
              C=COS(EXP2I)
              S=SIN(EXP2I)
              EXP2R=C*AUX
              EXP2I=S*AUX
C
C             Contribution to the integral at the central point
              I=NQNSGF+1+M1
              RAM(I)=RAM(I)+BAMPR*EXP0R-BAMPI*EXP0I
C
C             Exponential function at the central point
              EXPR=EXP0R
              EXPI=EXP0I
C             Loop over the gridpoints
              DO 210 I1=I+1,I+K1-M1
                AUX =EXPR*EXP1R-EXPI*EXP1I
                EXPI=EXPR*EXP1I+EXPI*EXP1R
                EXPR=AUX
                AUX  =EXP1R*EXP2R-EXP1I*EXP2I
                EXP1I=EXP1R*EXP2I+EXP1I*EXP2R
                EXP1R=AUX
                RAM(I1)=RAM(I1)+BAMPR*EXPR-BAMPI*EXPI
  210         CONTINUE
C
C             Exponential function at the central point
              EXPR=EXP0R
              EXPI=EXP0I
C             Exponential correction at the central point
              EXP1R=EXP1MR
              EXP1I=EXP1MI
C             Loop over the gridpoints
              DO 211 I1=I-1,I+J1-M1,-1
                AUX =EXPR*EXP1R-EXPI*EXP1I
                EXPI=EXPR*EXP1I+EXPI*EXP1R
                EXPR=AUX
                AUX  =EXP1R*EXP2R-EXP1I*EXP2I
                EXP1I=EXP1R*EXP2I+EXP1I*EXP2R
                EXP1R=AUX
                RAM(I1)=RAM(I1)+BAMPR*EXPR-BAMPI*EXPI
  211         CONTINUE
            END IF
          END IF
  214   CONTINUE
      ELSE IF(NDIM.LE.6) THEN
C       2-D:
        DO 224 ISGF=0,NSGF-1
          IF(MOD(ISGF,50).EQ.0) THEN
            WRITE(*,'(2(A,I7))')
     *        '+SGFGRD: Gabor function',2*ISGF,' /',2*NSGF
          END IF
          IQ=NQ*ISGF
          IF(RAM(IQ+11).NE.0.0.OR.RAM(IQ+12).NE.0.0) THEN
C           Quantities describing the Gabor function
            BX1=RAM(IQ+1)
            BX2=RAM(IQ+2)
            BK1=RAM(IQ+3)
            BK2=RAM(IQ+4)
C           Half matrix Y
            BY11=RAM(IQ+5)
            BY12=RAM(IQ+6)
            BY22=RAM(IQ+7)
C           Half matrix R
            BR11=-RAM(IQ+8)
            BR12=-RAM(IQ+9)
            BR22=-RAM(IQ+10)
C           Determinant of half matrix Y
            BYDET=BY11*BY22-BY12*BY12
C           Amplitude coefficients
            AUX=2.0*SQRT(SQRT(4.0*BYDET)/COEF2D)
            BAMPR=AUX*RAM(IQ+11)
            BAMPI=AUX*RAM(IQ+12)
C           Extent of the Gabor function along axis X2
            AUX=SQRT(RELLOG/(BY22-BY12*BY12/BY11))
            J2=MAX0(INT((BX2-AUX-O2)/D2+0.999),0)
            K2=MIN0(INT((BX2+AUX-O2)/D2+0.001),N2-1)
            DO 222 I2=J2,K2
              DX2=O2+D2*FLOAT(I2)-BX2
C             Extent of the Gabor function along axis X1
              AUX=RELLOG-BY22*DX2*DX2
              IF(AUX.GE.0.0) THEN
C               Halfwidth of the Gabor function along axis X1
                AUX=SQRT(AUX/BY11)
C               Central point along the gridline
                DX1=BX1-BY12*DX2/BY11-O1
                J1=MAX0(INT((DX1-AUX)/D1+0.999),0)
                K1=MIN0(INT((DX1+AUX)/D1+0.001),N1-1)
                IF(J1.LE.K1) THEN
C                 Index of the central point along the gridline
                  M1=MAX0(J1,MIN0(NINT(DX1/D1),K1))
C                 Relative coordinate of the central point
                  DX1=O1+D1*FLOAT(M1)-BX1
C                 Exponent at the central point
                  EXP0R=DX1*(BY11*DX1+2.0*BY12*DX2)+DX2*BY22*DX2
                  EXP0I=      DX1*(BK1+BR11*DX1+2.0*BR12*DX2)
                  EXP0I=EXP0I+DX2*(BK2+BR22*DX2)
C                 The first derivative of the exponent
                  EXP1R=BY11*DX1+BY12*DX2
                  EXP1I=BR11*DX1+BR12*DX2
                  EXP1R=D1*(EXP1R+EXP1R)
                  EXP1I=D1*(EXP1I+EXP1I+BK1)
C                 Half the second derivative of the exponent
                  EXP2R=D1*BY11*D1
                  EXP2I=D1*BR11*D1
C                 Half the second derivative minus the first derivative
                  EXP1MR=EXP2R-EXP1R
                  EXP1MI=EXP2I-EXP1I
C                 Half the second derivative plus the first derivative
                  EXP1R=EXP2R+EXP1R
                  EXP1I=EXP2I+EXP1I
C                 Second derivative of the exponent
                  EXP2R=EXP2R+EXP2R
                  EXP2I=EXP2I+EXP2I
C                 Exponential function at the central point
                  AUX=EXP(-EXP0R)
                  C=COS(EXP0I)
                  S=SIN(EXP0I)
                  EXP0R=C*AUX
                  EXP0I=S*AUX
C                 Exponential correction at the central point (+)
                  AUX=EXP(-EXP1R)
                  C=COS(EXP1I)
                  S=SIN(EXP1I)
                  EXP1R=C*AUX
                  EXP1I=S*AUX
C                 Exponential correction at the central point (-)
                  AUX=EXP(-EXP1MR)
                  C=COS(EXP1MI)
                  S=SIN(EXP1MI)
                  EXP1MR=C*AUX
                  EXP1MI=S*AUX
C                 Constant correction to the correction
                  AUX=EXP(-EXP2R)
                  C=COS(EXP2I)
                  S=SIN(EXP2I)
                  EXP2R=C*AUX
                  EXP2I=S*AUX
C
C                 Contribution to the integral at the central point
                  I=NQNSGF+1+M1+N1*I2
                  RAM(I)=RAM(I)+BAMPR*EXP0R-BAMPI*EXP0I
C
C                 Exponential function at the central point
                  EXPR=EXP0R
                  EXPI=EXP0I
C                 Loop over the gridpoints
                  DO 220 I1=I+1,I+K1-M1
                    AUX =EXPR*EXP1R-EXPI*EXP1I
                    EXPI=EXPR*EXP1I+EXPI*EXP1R
                    EXPR=AUX
                    AUX  =EXP1R*EXP2R-EXP1I*EXP2I
                    EXP1I=EXP1R*EXP2I+EXP1I*EXP2R
                    EXP1R=AUX
                    RAM(I1)=RAM(I1)+BAMPR*EXPR-BAMPI*EXPI
  220             CONTINUE
C
C                 Exponential function at the central point
                  EXPR=EXP0R
                  EXPI=EXP0I
C                 Exponential correction at the central point
                  EXP1R=EXP1MR
                  EXP1I=EXP1MI
C                 Loop over the gridpoints
                  DO 221 I1=I-1,I+J1-M1,-1
                    AUX =EXPR*EXP1R-EXPI*EXP1I
                    EXPI=EXPR*EXP1I+EXPI*EXP1R
                    EXPR=AUX
                    AUX  =EXP1R*EXP2R-EXP1I*EXP2I
                    EXP1I=EXP1R*EXP2I+EXP1I*EXP2R
                    EXP1R=AUX
                    RAM(I1)=RAM(I1)+BAMPR*EXPR-BAMPI*EXPI
  221             CONTINUE
                END IF
              END IF
  222       CONTINUE
          END IF
  224   CONTINUE
      ELSE
C       3-D:
        DO 234 ISGF=0,NSGF-1
          IF(MOD(ISGF,50).EQ.0) THEN
            WRITE(*,'(2(A,I7))')
     *        '+SGFGRD: Gabor function',2*ISGF,' /',2*NSGF
          END IF
          IQ=NQ*ISGF
          IF(RAM(IQ+19).NE.0.0.OR.RAM(IQ+20).NE.0.0) THEN
C           Quantities describing the Gabor function
            BX1=RAM(IQ+1)
            BX2=RAM(IQ+2)
            BX3=RAM(IQ+3)
            BK1=RAM(IQ+4)
            BK2=RAM(IQ+5)
            BK3=RAM(IQ+6)
C           Half matrix Y
            BY11=RAM(IQ+7)
            BY12=RAM(IQ+8)
            BY22=RAM(IQ+9)
            BY13=RAM(IQ+10)
            BY23=RAM(IQ+11)
            BY33=RAM(IQ+12)
C           Half matrix R
            BR11=-RAM(IQ+13)
            BR12=-RAM(IQ+14)
            BR22=-RAM(IQ+15)
            BR13=-RAM(IQ+16)
            BR23=-RAM(IQ+17)
            BR33=-RAM(IQ+18)
C           Determinant of the 2*2 submatrix of half Y
            DET=BY11*BY22-BY12*BY12
C           Determinant of half matrix Y
            BYDET=BY33*DET-BY11*BY23*BY23
     *           +BY13*(2.0*BY12*BY23-BY22*BY13)
C           Amplitude coefficients
            AUX=2.0*SQRT(SQRT(8.0*BYDET)/COEF3D)
            BAMPR=AUX*RAM(IQ+19)
            BAMPI=AUX*RAM(IQ+20)
C           Matrix inverse to the 2*2 submatrix of Y times (BY13,BY23)
            AY13=( BY22*BY13-BY12*BY23)/DET
            AY23=(-BY12*BY13+BY11*BY23)/DET
C           Extent of the Gabor function along axis X3
            AUX=SQRT(RELLOG*DET/BYDET)
            J3=MAX0(INT((BX3-AUX-O3)/D3+0.999),0)
            K3=MIN0(INT((BX3+AUX-O3)/D3+0.001),N3-1)
            DO 233 I3=J3,K3
              DX3=O3+D3*FLOAT(I3)-BX3
C
C             Transforming 3-D Gabor packet to 2-D Gabor packet
C             Shift of the central point
              DX1=-AY13*DX3
              DX2=-AY23*DX3
C             Central point of the 2-D Gabor packet
              AX1=BX1+DX1
              AX2=BX2+DX2
C             Wavenumber vetor of the 2-D Gabor packet
              AK1=BK1+2.0*(BR11*DX1+BR12*DX2+BR13*DX3)
              AK2=BK2+2.0*(BR12*DX1+BR22*DX2+BR23*DX3)
C             Exponent of the 2-D Gabor packet at its central point
              AR33=     DX1*(BK1+BR11*DX1+2.0*(BR12*DX2+BR13*DX3))
              AR33=AR33+DX2*(BK2+BR22*DX2+2.0* BR23*DX3)
              AR33=AR33+DX3*(BK3+BR33*DX3)
              AY33=(BY33-BY13*AY13-BY23*AY23)*DX3*DX3
C             Matrices BY11,BY12,BY22 and BR11,BR12,BR22 keep unchanged.
C
C             Extent of the Gabor function along axis X2
              AUX=RELLOG-AY33
              IF(AUX.GE.0.0) THEN
                AUX=SQRT(AUX*BY11/DET)
                J2=MAX0(INT((AX2-AUX-O2)/D2+0.999),0)
                K2=MIN0(INT((AX2+AUX-O2)/D2+0.001),N2-1)
                DO 232 I2=J2,K2
                  DX2=O2+D2*FLOAT(I2)-AX2
C                 Extent of the Gabor function along axis X1
                  AUX=RELLOG-AY33-(DET/BY11)*DX2*DX2
                  IF(AUX.GE.0.0) THEN
                    AUX=SQRT(AUX/BY11)
                    DX1=AX1-BY12*DX2/BY11-O1
                    J1=MAX0(INT((DX1-AUX)/D1+0.999),0)
                    K1=MIN0(INT((DX1+AUX)/D1+0.001),N1-1)
                    IF(J1.LE.K1) THEN
C                     Index of the central point along the gridline
                      M1=MAX0(J1,MIN0(NINT(DX1/D1),K1))
C                     Relative coordinate of the central point
                      DX1=O1+D1*FLOAT(M1)-AX1
C                     Exponent at the central point
                      EXP0R=AY33
                      EXP0R=EXP0R+DX1*(BY11*DX1+2.0*BY12*DX2)
                      EXP0R=EXP0R+DX2* BY22*DX2
                      EXP0I=AR33
                      EXP0I=EXP0I+DX1*(AK1+BR11*DX1+2.0*BR12*DX2)
                      EXP0I=EXP0I+DX2*(AK2+BR22*DX2)
C                     The first derivative of the exponent
                      EXP1R=BY11*DX1+BY12*DX2
                      EXP1I=BR11*DX1+BR12*DX2
                      EXP1R=D1*(EXP1R+EXP1R)
                      EXP1I=D1*(EXP1I+EXP1I+AK1)
C                     Half the second derivative of the exponent
                      EXP2R=D1*BY11*D1
                      EXP2I=D1*BR11*D1
C                     Half the second derivative minus first derivative
                      EXP1MR=EXP2R-EXP1R
                      EXP1MI=EXP2I-EXP1I
C                     Half the second derivative plus first derivative
                      EXP1R=EXP2R+EXP1R
                      EXP1I=EXP2I+EXP1I
C                     Second derivative of the exponent
                      EXP2R=EXP2R+EXP2R
                      EXP2I=EXP2I+EXP2I
C                     Exponential function at the central point
                      AUX=EXP(-EXP0R)
                      C=COS(EXP0I)
                      S=SIN(EXP0I)
                      EXP0R=C*AUX
                      EXP0I=S*AUX
C                     Exponential correction at the central point (+)
                      AUX=EXP(-EXP1R)
                      C=COS(EXP1I)
                      S=SIN(EXP1I)
                      EXP1R=C*AUX
                      EXP1I=S*AUX
C                     Exponential correction at the central point (-)
                      AUX=EXP(-EXP1MR)
                      C=COS(EXP1MI)
                      S=SIN(EXP1MI)
                      EXP1MR=C*AUX
                      EXP1MI=S*AUX
C                     Constant correction to the correction
                      AUX=EXP(-EXP2R)
                      C=COS(EXP2I)
                      S=SIN(EXP2I)
                      EXP2R=C*AUX
                      EXP2I=S*AUX
C
C                     Contribution to the integral at the central point
                      I=NQNSGF+1+M1+N1*(I2+N2*I3)
                      RAM(I)=RAM(I)+BAMPR*EXP0R-BAMPI*EXP0I
C
C                     Exponential function at the central point
                      EXPR=EXP0R
                      EXPI=EXP0I
C                     Loop over the gridpoints
                      DO 230 I1=I+1,I+K1-M1
                        AUX =EXPR*EXP1R-EXPI*EXP1I
                        EXPI=EXPR*EXP1I+EXPI*EXP1R
                        EXPR=AUX
                        AUX  =EXP1R*EXP2R-EXP1I*EXP2I
                        EXP1I=EXP1R*EXP2I+EXP1I*EXP2R
                        EXP1R=AUX
                        RAM(I1)=RAM(I1)+BAMPR*EXPR-BAMPI*EXPI
  230                 CONTINUE
C
C                     Exponential function at the central point
                      EXPR=EXP0R
                      EXPI=EXP0I
C                     Exponential correction at the central point
                      EXP1R=EXP1MR
                      EXP1I=EXP1MI
C                     Loop over the gridpoints
                      DO 231 I1=I-1,I+J1-M1,-1
                        AUX =EXPR*EXP1R-EXPI*EXP1I
                        EXPI=EXPR*EXP1I+EXPI*EXP1R
                        EXPR=AUX
                        AUX  =EXP1R*EXP2R-EXP1I*EXP2I
                        EXP1I=EXP1R*EXP2I+EXP1I*EXP2R
                        EXP1R=AUX
                        RAM(I1)=RAM(I1)+BAMPR*EXPR-BAMPI*EXPI
  231                 CONTINUE
                    END IF
                  END IF
  232           CONTINUE
              END IF
  233       CONTINUE
          END IF
  234   CONTINUE
      END IF
C
C     Writing the output grid:
      FORM='FORMATTED'
      CALL WARRAY(LU1,FGRD,FORM,.FALSE.,0.0,.FALSE.,0.0,
     *                                           N1N2N3,RAM(NQNSGF+1))
C
      WRITE(*,'(A)') '+SGFGRD: Done.                 '
      STOP
      END
C
C=======================================================================
C
      INCLUDE 'error.for'
C     error.for
      INCLUDE 'sep.for'
C     sep.for
      INCLUDE 'length.for'
C     length.for
      INCLUDE 'forms.for'
C     forms.for
      INCLUDE 'mat.for'
C     mat.for
C
C=======================================================================
C