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