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