C
C Program SGFHOM to generate the structural Gabor functions which shape
C is optimized for a zero-offset surface seismic reflection survey
C in a homogeneous 2-D velocity model
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 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 output files:
C     SGF='string'... Name of the output file with structural Gabor
C             functions.
C             Input for programs
C             sgfmat.for
C             and sgpini.for.
C             The file is not generated if SGF=' '.
C             Description of file SGF.
C             Default: SGF='sgf.out'
C     PTS='string'... Name of the output file with the central points of
C             structural Gabor functions.  If several Gabor functions
C             have equal central point, this central point appears only
C             once in file PTS.
C             Input for program
C             mtt.for
C             which should calculate the the Green tensors corresponding
C             to the waves incident from the given sources to these
C             points.  The Green tensors are then read by program
C             sgpini.for.
C             The file is not generated if PTS=' '.
C             Description of file PTS.
C             Default: PTS='pts.out'
C     SGFPTS='string'... Name of the output file linking the structural
C             Gabor functions of file SGF to central points of file PTS.
C             Input for program
C             sgpini.for.
C             The file is not generated if SGFPTS=' '.
C             Description of file SGFPTS.
C             Default: SGFPTS='sgfpts.out'
C Data specifying the phase-space region covered by Gabor functions:
C     X1MIN=real, X1MAX=real, X2MIN=real, X2MAX=real, X3MIN=real,
C     X3MAX=real... Cartesian coordinates specifying the boundaries of
C             the rectangular spatial volume inside which the central
C             points of structural Gabor functions are located.
C             In the current 2-D version, this volume should be reduced
C             to a 2-D rectangle by specifying X1MIN=X1MAX or
C             X2MIN=X2MAX or X3MIN=X3MAX.  For a 2-D rectangle, only 2-D
C             wavenumber vector is considered.
C             It is also possible to reduce this volume to a 1-D line
C             segment for testing purposes.  For a 1-D line segment,
C             only 1-D wavenumber vector is considered.
C             Defaults: X1MIN=0.0, X1MAX=0.0, X2MIN=0.0, X2MAX=0.0,
C               X3MIN=0.0, X3MAX=0.0
C     WKMIN=real... Minimum absolute value of the projection of the
C             structural wavenumber vector onto the normal to the
C             surface at which the seismic sources and receivers are
C             situated.  WKMIN must be non-negative.
C             Default: WKMIN=0.
C     WKMAX=real... Maximum norm of the structural wavenumber vector.
C             No default: WKMAX must be specified and must be greater
C               than WKMIN.
C Data specifying the phase-space density of Gabor functions:
C     REDUNDANCY=real... Redundancy ratio with respect to one spatial
C             direction and the corresponding wavenumber.
C             The redundancy ratio corresponds to 2*PI/(DX*DK),
C             where DX*DK is the product of the spatial and wavenumber
C             discretization steps.
C             The redundancy ratio should be greater than 1.
C             Default: REDUNDANCY=1.1
C     LATTICEV=integer... Structure of the lattice of the central points
C             of Gabor functions in the coordinate plane WX3, WK3, where
C             coordinate WX3 and wavenumber component WK3 are measured
C             perpendicularly to the surface at which the seismic
C             sources and receivers are situated.
C             Parameter LATTICEV specifies the position of hyperbolae
C             WX3*WK3=constant containing the chains of central points
C             of Gabor functions.
C             LATTICEV=0:  The hyperbolae are positioned with respect to
C               WX3MIN and WKMIN.  Point (WX3,WK3)=(WX3MIN,WKMIN) is
C               situated in the middle between two hyperbolae.
C               There is a gap (-WKMIN,+WKMIN) between positive and
C               negative vertical wavenumbers WK3.
C             LATTICEV=1:  The hyperbolae are independent of the given
C               phase-space domain.  There is a rough continuation
C               between positive and negative vertical wavenumbers WK3,
C               except for the central points removed from interval
C               (-WKMIN,+WKMIN).  If WKMIN=0 (or WX3MIN=0), options
C               LATTICEV=0 and LATTICEV=1 are equal.
C             LATTICEV=2:  The hyperbolae are independent of the given
C               phase-space domain, and are located strictly according
C               to Klimes(2008), with at least one chain of points
C               missing between positive and negative vertical
C               wavenumbers WK3 even for WKMIN=0.
C             Default: LATTICEV=0
C     LATTICEH=integer... Structure of the lattice of the central points
C             WX3=constant, WK3=constant of the phase space, where
C             coordinate WX3 and wavenumber component WK3 are measured
C             perpendicularly to the surface at which the seismic
C             sources and receivers are situated.
C             LATTICEH=0:  Oblique lattice with very regular coverage
C               according to equations (77)-(80) in 2-D.
C             LATTICEH=1:  Rectangular lattice according to equations
C               (72)-(74).
C             LATTICEH=2:  Rectangular lattice according to equations
C               (72), (75) and (76).
C             Default: LATTICEH=0
C Data specifying the surface with seismic sources and receivers:
C     SRF1=real, SRF2=real, SRF3=real... Unit vector perpendicular to
C             the surface at which the seismic sources and receivers are
C             situated.  If the specified vector is not unit, it is
C             normalized by this program.
C             The surface should not intersect the rectangular spatial
C             volume inside which the central points of structural Gabor
C             functions are located.
C             The vector should point from the surface towards the
C             rectangular spatial volume inside which the central points
C             of structural Gabor functions are located.
C             In the current 2-D version, the specified vector should
C             be parallel with the specified 2-D rectangle.
C             Defaults: SRF1=0.0, SRF2=0.0, SRF3=-1.0
C     SRF0=real... Value of the scalar product of the coordinates
C             of the seismic sources and receivers with unit vector
C             (SRF1,SRF2,SRF3).
C             Default: SRF0=0.0
C Output format:
C     MAXDIG=positive integer... Minimum number of digits of a positive
C             number in the output format.
C             MAXDIG must be less than 10.
C             Default: MAXDIG=6
C     MINDIG=positive integer... Number of digits to change editing F
C             to editing G.
C             MINDIG should be less than MAXDIG.
C             Default: MINDIG=4
C
C                                                     
C File SGF with structural Gabor functions:
C     (version of form PTS):
C (1) None to several strings terminated by / (a slash)
C (2) For each couple of Gabor functions data (2.1):
C (2.1) 'IGF',X1,X2,X3,K1,K2,K3,RK11,RK12,RK22,RK13,RK23,RK33,
C       YK11,YK12,YK22,YK13,YK23,YK33,/
C     'IGF'... Name of the data line.  The names are odd
C       integers.
C     X1,X2,X3... Coordinates of the central point of the
C       structural Gabor function.
C     K1,K2,K3... Real-valued structural wavenumber.
C     RK11,RK12,RK22,RK13,RK23,RK33... Real part of the symmetric 3*3
C       matrix K describing the envelope of the structural Gabor
C       function.
C     YK11,YK12,YK22,YK13,YK23,YK33... Imaginary part of the symmetric
C       3*3 matrix K describing the envelope of the structural Gabor
C       function.
C     One data line describes two Gabor functions:
C     Gabor function IGF with wavenumber (K1,K2,K3), matrix K,
C       and unknown complex-valued model coefficients;
C     Gabor function IGF+1 with wavenumber (-K1,-K2,-K3),
C       complex-conjugate matrix K, and complex-conjugate model
C       coefficients.
C (3) / or end of file.
C
C                                                     
C File PTS with central points of the structural Gabor functions
C     for calculating quantities corresponding to the incident wave
C     by program mtt.for
C     (version of form PTS):
C (1) None to several strings terminated by / (a slash)
C (2) For each point data (2.1):
C (2.1) 'IPT',X1,X2,X3,/
C     'IPT'... Name of the point corresponding to the smallest index IGF
C       of Gabor functions centred at the point.
C     X1,X2,X3... Coordinates of the central point common to several
C       structural Gabor functions.
C (3) / or end of file.
C
C                                                  
C File SGFPTS linking the structural Gabor functions of file SGF to the
C     central points of file PTS:
C (1) None to several strings terminated by / (a slash)
C (2) For each Gabor function of file SGF data (2.1):
C (2.1) 'IGF','IPT',/
C     'IGF'.. Name of the structural Gabor function of file SGF.
C             The names are odd integers.
C     'IPT'.. Name of the central point of file PTS.
C             The name corresponds to the smallest index of Gabor
C             functions centred at the point.
C (3) / or end of file.
C
C=======================================================================
C
C     External functions and subroutines:
      EXTERNAL RSEP1,RSEP3T,RSEP3R,ERROR,FORM2
C
C     Filenames and parameters:
      CHARACTER*80 FSEP,FSGF,FPTS,FSGFPT
      INTEGER LU1,LU2,LU3
      PARAMETER (LU1=1,LU2=2,LU3=3)
C
C     Input data:
      INTEGER LATICV,LATICH
      REAL X1MIN,X1MAX,X2MIN,X2MAX,X3MIN,X3MAX,WKMIN,WKMAX
      REAL REDUND,SRF0,SRF1,SRF2,SRF3
C
C     Output data:
      INTEGER NOUT,IGF,IGFMAX,IPT
      PARAMETER (NOUT=18)
      REAL X1,X2,X3,OUT(NOUT)
      CHARACTER*(11+8*NOUT) FORMAT
C
C     Transformation matrix:
C     REAL E10,E20,E30,E11,E21,E31,E12,E22,E32,E13,E23,E33
      REAL E10,E20,E30,E11,E21,E31,E13,E23,E33
      REAL E1111,E1211,E2211,E1311,E2311,E3311
C     REAL E1113,E1213,E2213,E1313,E2313,E3313
      REAL E1133,E1233,E2233,E1333,E2333,E3333
C
C     Constants describing the Gabor functions:
      REAL PI,STEP0,Y0,R0,ALPHA,DW4
      PARAMETER (PI=3.141592654)
C
C     Gabor functions in working coordinates:
      REAL WX1,WX3,WK1,WK3,Y11,Y33,R11,R33
      REAL WX1MIN,WX1MAX,WX1INI,WX3MIN,WX3MAX,DWX1,DWK1,DWK1I2,WK1INI
C
C     Other variables:
      INTEGER I1MIN,I1MAX,I2MIN,I2MAX,I3MIN,I3MAX,I4MIN,I4MAX
      INTEGER I1,I2,I3,I4,I
      REAL W4INI,W4,AUX
C
C.......................................................................
C
C     Reading main input data:
      WRITE(*,'(A)') '+SGFHOM: Enter input filename: '
      FSEP=' '
      READ (*,*) FSEP
      IF(FSEP.EQ.' ') THEN
C       SGFHOM-01
        CALL ERROR('SGFHOM-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)') '+SGFHOM: Working...            '
C
C     Reading output filenames:
      CALL RSEP3T('SGF'   ,FSGF  ,'sgf.out')
      CALL RSEP3T('PTS'   ,FPTS  ,'pts.out' )
      CALL RSEP3T('SGFPTS',FSGFPT,'sgfpts.out')
C
C     Reading other input SEP parameters:
      CALL RSEP3R('X1MIN' ,X1MIN ,0.0)
      CALL RSEP3R('X1MAX' ,X1MAX ,0.0)
      CALL RSEP3R('X2MIN' ,X2MIN ,0.0)
      CALL RSEP3R('X2MAX' ,X2MAX ,0.0)
      CALL RSEP3R('X3MIN' ,X3MIN ,0.0)
      CALL RSEP3R('X3MAX' ,X3MAX ,0.0)
      CALL RSEP3R('WKMIN' ,WKMIN ,0.0)
      CALL RSEP3R('WKMAX' ,WKMAX ,0.0)
      CALL RSEP3R('REDUNDANCY',REDUND,1.1)
      CALL RSEP3I('LATTICEV',LATICV,0)
      CALL RSEP3I('LATTICEH',LATICH,0)
      CALL RSEP3R('SRF0'  ,SRF0  ,0.0)
      CALL RSEP3R('SRF1'  ,SRF1  ,0.0)
      CALL RSEP3R('SRF2'  ,SRF2  ,0.0)
      CALL RSEP3R('SRF3'  ,SRF3  ,-1.0)
      IF(X1MIN.GT.X1MAX.OR.X2MIN.GT.X2MAX.OR.X3MIN.GT.X3MAX) THEN
C       SGFHOM-02
        CALL ERROR('SGFHOM-02: Wrong boundaries of the spatial volume')
C       Minimum coordinate bounding the spatial volume is greater than
C       the corresponding maximum coordinate.
      END IF
      IF(WKMIN.LT.0.0.OR.WKMIN.GE.WKMAX) THEN
C       SGFHOM-03
        CALL ERROR('SGFHOM-03: Wrong boundaries of the wavenumbers')
C       The minimum wavenumber must be non-negative and smaller than the
C       maximum wavenumber.
      END IF
      IF(X1MIN.LT.X1MAX.AND.X2MIN.LT.X2MAX.AND.X3MIN.LT.X3MAX) THEN
C       SGFHOM-04
        CALL ERROR('SGFHOM-04: 3-D spatial volume')
C       In the current 2-D version, the spatial volume specified by
C       X1MIN, X1MAX, X2MIN, X2MAX, X3MIN, X3MAX must be 2-D.
      END IF
C
C     Opening the output files and writing their beginnings:
      IF(FSGF.NE.' ') THEN
        OPEN (LU1,FILE=FSGF,FORM='FORMATTED')
        WRITE(LU1,'(A)') '/'
      END IF
      IF(FPTS.NE.' ') THEN
        OPEN (LU2,FILE=FPTS,FORM='FORMATTED')
        WRITE(LU2,'(A)') '/'
      END IF
      IF(FSGFPT.NE.' ') THEN
        OPEN (LU3,FILE=FSGFPT,FORM='FORMATTED')
        WRITE(LU3,'(A)') '/'
      END IF
C
C     Transformation matrix from working to output coordinates:
      AUX=SQRT(SRF1*SRF1+SRF2*SRF2+SRF3*SRF3)
      IF(AUX.LE.0.0) THEN
C       SGFHOM-05
        CALL ERROR('SGFHOM-05: Zero vector (SRF1,SRF2,SRF3)')
      END IF
C     E12=0.0
C     E22=0.0
C     E32=0.0
      E13=SRF1/AUX
      E23=SRF2/AUX
      E33=SRF3/AUX
      E10=E13*SRF0
      E20=E23*SRF0
      E30=E33*SRF0
      IF(X1MIN.EQ.X1MAX) THEN
        IF(SRF1.NE.0.0) THEN
C         SGFHOM-06
          CALL ERROR('SGFHOM-06: Non-zero SRF1')
C         In the current 2-D version, specified vector (SRF1,SRF2,SRF3)
C         must be parallel with the 2-D rectangle specified by
C         X1MIN, X1MAX, X2MIN, X2MAX, X3MIN, X3MAX.
        END IF
        E11= 0.0
        E21=-E33
        E31= E23
C       Disabling the check of position affected by rounding errors
        X1MIN=X1MIN-ABS(X1MIN)
        X1MAX=X1MAX-ABS(X1MAX)
      END IF
      IF(X2MIN.EQ.X2MAX) THEN
        IF(SRF2.NE.0.0) THEN
C         SGFHOM-07
          CALL ERROR('SGFHOM-07: Non-zero SRF1')
C         In the current 2-D version, specified vector (SRF1,SRF2,SRF3)
C         must be parallel with the 2-D rectangle specified by
C         X1MIN, X1MAX, X2MIN, X2MAX, X3MIN, X3MAX.
        END IF
        E11=-E33
        E21= 0.0
        E31= E13
C       Disabling the check of position affected by rounding errors
        X2MIN=X2MIN-ABS(X2MIN)
        X2MAX=X2MAX-ABS(X2MAX)
      END IF
      IF(X3MIN.EQ.X3MAX) THEN
        IF(SRF3.NE.0.0) THEN
C         SGFHOM-08
          CALL ERROR('SGFHOM-08: Non-zero SRF1')
C         In the current 2-D version, specified vector (SRF1,SRF2,SRF3)
C         must be parallel with the 2-D rectangle specified by
C         X1MIN, X1MAX, X2MIN, X2MAX, X3MIN, X3MAX.
        END IF
        E11=-E23
        E21= E13
        E31= 0.0
        IF(X1MIN.EQ.X1MAX.OR.X2MIN.EQ.X2MAX) THEN
          E11=0.0
          E21=0.0
        END IF
C       Disabling the check of position affected by rounding errors
        X3MIN=X3MIN-ABS(X3MIN)
        X3MAX=X3MAX-ABS(X3MAX)
      END IF
      WX1MIN=AMIN1(E11*X1MIN,E11*X1MAX)
     *      +AMIN1(E21*X2MIN,E21*X2MAX)
     *      +AMIN1(E31*X3MIN,E31*X3MAX)
      WX1MAX=AMAX1(E11*X1MIN,E11*X1MAX)
     *      +AMAX1(E21*X2MIN,E21*X2MAX)
     *      +AMAX1(E31*X3MIN,E31*X3MAX)
      WX3MIN=AMIN1(E13*X1MIN,E13*X1MAX)
     *      +AMIN1(E23*X2MIN,E23*X2MAX)
     *      +AMIN1(E33*X3MIN,E33*X3MAX)
      WX3MAX=AMAX1(E13*X1MIN,E13*X1MAX)
     *      +AMAX1(E23*X2MIN,E23*X2MAX)
     *     +AMAX1(E33*X3MIN,E33*X3MAX)
      WX1INI=0.5*(WX1MIN+WX1MAX)
      IF(WX3MIN.LT.0.0) THEN
        IF(WX3MAX.LE.0.0) THEN
C         SGFHOM-09
          CALL ERROR('SGFHOM-09: Wrong vector (SRF1,SRF2,SRF3)')
C         Vector (SRF1,SRF2,SRF3) is not oriented towards the spatial
C         volume specified by X1MIN, X1MAX, X2MIN, X2MAX, X3MIN, X3MAX.
        ELSE
C         SGFHOM-10
          CALL ERROR('SGFHOM-10: Wrong position of the surface')
C         The surface with seismic sources and receivers intersects
C         the spatial volume specified by X1MIN, X1MAX, X2MIN, X2MAX,
C         X3MIN, X3MAX.
        END IF
      END IF
      E1111=E11*E11
      E1211=E11*E21
      E2211=E21*E21
      E1311=E11*E31
      E2311=E21*E31
      E3311=E31*E31
C     E1113=E11*E13*2.0
C     E1213=E11*E23+E13*E21
C     E2213=E21*E23*2.0
C     E1313=E11*E33+E13*E21
C     E2313=E21*E33+E13*E21
C     E3313=E31*E33*2.0
      E1133=E13*E13
      E1233=E13*E23
      E2233=E23*E23
      E1333=E13*E33
      E2333=E23*E33
      E3333=E33*E33
C
C     Initial values for output:
      IGF=1
      FORMAT(1:11)='(A,I6.6,2A,'
C     Maximum integer fitting into the above format
      IGFMAX=999999
C
C     Constants describing the Gabor functions:
      STEP0=SQRT(2.0*PI/REDUND)
      Y0=SQRT(3.0)/4.0
      R0=5.0/4.0
      AUX=(1.0+R0)**2/Y0+Y0
      ALPHA=2.0/AUX
      DW4=0.5*STEP0*SQRT(AUX)
C
C     Positioning hyperbolae containing the chains of lattice points
      IF(LATICV.EQ.0) THEN
C       Hyperbolae positioned with respect to WX3MIN and WKMIN
        W4INI=SQRT(WX3MIN*WKMIN)+0.5*DW4
      ELSE IF(LATICV.EQ.1) THEN
C       Hyperbolae independent of the phase-space domain
        W4INI=0.5*DW4
      ELSE IF(LATICV.EQ.2) THEN
C       Hyperbolae according to Klimes(2008)
        W4INI=0.0
      ELSE
C       SGFHOM-11
        CALL ERROR('SGFHOM-11: Wrong SEP parameter LATTICEV')
C       Refer to the input data.
      END IF
C
C     Loops over Gabor functions:
      I4MIN=NINT((SQRT(WX3MIN*WKMIN)-W4INI)/DW4+0.499)
      I4MAX=NINT((SQRT(WX3MAX*WKMAX)-W4INI)/DW4-0.499)
      DO 14 I4=I4MIN,I4MAX
        W4=W4INI+DW4*FLOAT(I4)
        I3MIN=NINT(ALOG(W4/WKMAX)*W4/DW4/ALPHA+0.499)
        IF(WX3MIN.GT.0.0) THEN
          I3MIN=MAX0(NINT(ALOG(WX3MIN/W4)*W4/DW4/ALPHA+0.499),I3MIN)
        END IF
        I3MAX=NINT(ALOG(WX3MAX/W4)*W4/DW4/ALPHA-0.499)
        IF(WKMIN.GT.0.0) THEN
          I3MAX=MIN0(NINT(ALOG(W4/WKMIN)*W4/DW4/ALPHA-0.499),I3MAX)
        END IF
        DO 13 I3=I3MIN,I3MAX
          WX3=W4*EXP(ALPHA*FLOAT(I3)*DW4/W4)
          WK3=W4*W4/WX3
          Y33=Y0*WK3/WX3
          R33=R0*WK3/WX3
          Y11=Y33
          R11=R33
          IF(LATICH.EQ.0) THEN
C           Basic option according to equations (77)-(80) in 2-D
            AUX=SQRT(Y11*2.0/SQRT(3.0))
            DWX1=STEP0/AUX
            DWK1=STEP0*AUX
            DWK1I2=0.5*DWK1+R11*DWX1
          ELSE IF(LATICH.EQ.1) THEN
C           First option according to equations (72)-(74)
            AUX=SQRT(Y11)
            DWX1=STEP0/AUX
            DWK1=STEP0*AUX
            DWK1I2=0.0
          ELSE IF(LATICH.EQ.2) THEN
C           Second option according to equations (72), (75) and (76)
            AUX=SQRT(R11**2/Y11+Y11)
            DWX1=STEP0/AUX
            DWK1=STEP0*AUX
            DWK1I2=0.0
          ELSE
C           SGFHOM-12
            CALL ERROR('SGFHOM-12: Wrong SEP parameter LATTICEH')
C           Refer to the input data.
          END IF
C
C         Loop over the horizontal coordinate:
          I2MAX=INT((WX1MAX-WX1INI)/DWX1)
          I2MIN=-I2MAX
          DO 12 I2=I2MIN,I2MAX
            WX1=WX1INI+DWX1*FLOAT(I2)
            WK1INI=DWK1I2*FLOAT(I2)
C
C           Transformation to output coordinates:
C           Coordinates of the Gabor function
            X1=E10+E11*WX1+E13*WX3
            X2=E20+E21*WX1+E23*WX3
            X3=E30+E31*WX1+E33*WX3
C           Checking the phase-space position of the Gabor function
            IF(X1MIN.LE.X1.AND.X1.LE.X1MAX.AND.
     *         X2MIN.LE.X2.AND.X2.LE.X2MAX.AND.
     *         X3MIN.LE.X3.AND.X3.LE.X3MAX.AND.
     *         WKMIN.LE.WK3.AND.WK3.LE.WKMAX) THEN
C             Writing the central coordinates of the Gabor functions
              OUT(1)=X1
              OUT(2)=X2
              OUT(3)=X3
              IF(FPTS.NE.' ') THEN
                IF(IGF.GE.IGFMAX) THEN
C                 
C                 SGFHOM-13
                  CALL ERROR('SGFHOM-13: Too many Gabor functions')
C                 The number of Gabor functions exceeds the maximum
C                 number which fits into the string describing each
C                 Gabor function, see IGFMAX.
                END IF
                IPT=IGF
                CALL FORM2(3,OUT,OUT,FORMAT(12:11+8*3))
                WRITE(LU2,FORMAT) '''',IPT,'''',(' ',OUT(I),I=1,3),' /'
              END IF
C
C             Loop over the horizontal component of wavenumber vector:
              I1MAX=NINT(( SQRT(WKMAX*WKMAX-WK3*WK3)-WK1INI)/DWK1-0.5)
              I1MIN=NINT((-SQRT(WKMAX*WKMAX-WK3*WK3)-WK1INI)/DWK1+0.5)
              IF(WX1MIN.EQ.WX1MAX) THEN
                I1MIN=0
                I1MAX=0
                WK1INI=0.0
              END IF
              DO 11 I1=I1MIN,I1MAX
                WK1=WK1INI+DWK1*FLOAT(I1)
C
C               Transformation to output coordinates:
                IF(FSGF.NE.' ') THEN
C                 Wavenumber vector with opposite sign
                  OUT( 4)=-E11*WK1-E13*WK3
                  OUT( 5)=-E21*WK1-E23*WK3
                  OUT( 6)=-E31*WK1-E33*WK3
C                 Real part of matrix K
                  OUT( 7)=E1111*Y11+E1133*Y33
                  OUT( 8)=E1211*Y11+E1233*Y33
                  OUT( 9)=E2211*Y11+E2233*Y33
                  OUT(10)=E1311*Y11+E1333*Y33
                  OUT(11)=E2311*Y11+E2333*Y33
                  OUT(12)=E3311*Y11+E3333*Y33
C                 Imaginary part of matrix K with opposite sign
                  OUT(13)=E1111*R11+E1133*R33
                  OUT(14)=E1211*R11+E1233*R33
                  OUT(15)=E2211*R11+E2233*R33
                  OUT(16)=E1311*R11+E1333*R33
                  OUT(17)=E2311*R11+E2333*R33
                  OUT(18)=E3311*R11+E3333*R33
C                 Writing
                  IF(IGF.GE.IGFMAX) THEN
C                   
C                   SGFHOM-14
                    CALL ERROR('SGFHOM-14: Too many Gabor functions')
C                   The number of Gabor functions exceeds the maximum
C                   number which fits into the string describing each
C                   Gabor function, see IGFMAX.
                  END IF
                  CALL FORM2(NOUT,OUT,OUT,FORMAT(12:11+8*NOUT))
                  WRITE(LU1,FORMAT)
     *                      '''',IGF,'''',(' ',OUT(I),I=1,NOUT),' /'
                END IF
                IF(FSGFPT.NE.' ') THEN
                  WRITE(LU3,'(3(A,I6.6))') '''',IGF,''' ''',IPT,''' /'
                END IF
                IGF=IGF+2
   11         CONTINUE
            END IF
   12     CONTINUE
   13   CONTINUE
   14 CONTINUE
C
C     Closing the output files:
      IF(FSGF.NE.' ') THEN
        WRITE(LU1,'(A)') '/'
        CLOSE(LU1)
      END IF
      IF(FPTS.NE.' ') THEN
        WRITE(LU2,'(A)') '/'
        CLOSE(LU2)
      END IF
      IF(FSGFPT.NE.' ') THEN
        WRITE(LU3,'(A)') '/'
        CLOSE(LU3)
      END IF
      WRITE(*,'(A,I9,A)') '+SGFHOM: Done.',IGF-1,' Gabor functions'
      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
C
C=======================================================================
C