C
C Program GRDFFT to compute the 1-D, 2-D or 3-D Fourier transformation C of a complex function defined on 1-D, 2-D or 3-D grid of points. C (The dimension of the FFT is less or equal to the grid dimension.) C If the number of gridpoints in any direction of the input grid is not C a power of 2, the input grid is enlarged to the nearest power of 2 C and the functional values in new gridpoints are completed according C to input parameter FFTFIL. C Subroutines from the Numerical Recipes are then used C for the entire FFT. C C Version: 6.00 C Date: 2005, November 12 C C Coded by Petr Bulant C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C http://sw3d.cz/staff/bulant.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 Parameter describing the form of the FFT: C FFT=real ... The Fourier transform F(y) has the form of integral C of f(x)exp(i*FFT*x*y)dx, where f(x) is the input function C to be transformed. The integral is then multiplied by the C factor of (ABS(FFT)/(2.*PI))**(NDIM/2), where NDIM is the C dimension of the part of the input grid to be transformed C (and of the transformed grid). C The mostly used values of FFT are 6.28, -6.28, 1, -1. C If FFT is input as a multiple of PI plus minus 1%, the C value of FFT is rounded to the multiple of PI. C Default: FFT=6.28 (which means 2.*PI) C Data specifying the parameters of the input grid: C O1=real, O2=real, O3=real ... Coordinates of the origin of the C grid. C Default: O1=0. O2=0. O3=0. C D1=real... Grid interval along the X1 axis. C Default: D1=0. C D2=real... Grid interval along the X2 axis. C Default: D2=0. C D3=real... Grid interval along the X3 axis. C Default: D3=0. C N1=positive integer... Number of gridpoints along the X1 axis. C Default: N1=1 C N2=positive integer... Number of gridpoints along the X2 axis. C Default: N2=1 C N3=positive integer... Number of gridpoints along the X3 axis. C Default: N3=1 C N4=positive integer... Number of space grids. The Fourier C transform with respect to X1, X2 and X3 will be repeated C N4 times, for each space grid. C If N4=0, parameter M4 is used to specify the number of C space grids. C Note that each space grid must begin at a new line of C the input file, because each space grid is read by a C a single read statement in free format. C Default: N4=1 C M4='string'... Name of the file containing a single integer number C specifying N4. C Default: M4=' ' means that N4=1. C Data specifying the parameters of the grid for the FFT: C N1FFT=positive integer... Number of gridpoints along the X1 axis. C N2FFT=positive integer... Number of gridpoints along the X2 axis. C N3FFT=positive integer... Number of gridpoints along the X3 axis. C NiFFT must be an integer power of 2, and must be greater C than or equal to the corresponding Ni of the input data. C NiFFT=0 means, that the Fourier transform is not to be C done in the direction of the corresponding axis. In this C case the corresponding Ni number of gridpoints is C considered, this influences the default value of NiOUT. C Default: NiFFT equal to the lowest power of 2, which is C greater or equal to the corresponding Ni. C FFTFIL=real ... Value, which is used at the gridpoints of the grid C for FFT, at which the value is not given by input data. C In present version FFTFIL is used for real part, imaginary C part is zero, with exception of FFTFIL=-999999999. In such C case, input data are linearly interpolated from the values C in the first and in the last gridpoints of the input grid, C to the gridpoints of the grid for the FFT. C Default: FFTFIL=-999999999. (interpolation of the values) C Data specifying the parameters of the output grid: C O1OUT=real, O2OUT=real, O3OUT=real ... Coordinates of the origin C of the output grid. C Default: OiOUT=-1./(2.*Di)*2.*PI/ABS(FFT) C D1OUT=real, D2OUT=real, D3OUT=real ... Grid intervals along the C first, second and third axes of the output grid. DiOUT C must not differ from the default value. C Default: DiOUT=1./(NiFFT*Di)*2.*PI/ABS(FFT) C N1OUT=positive integer, N2OUT=positive integer, C N3OUT=positive integer ... Numbers of gridpoints along the first, C second and third axes of the output grid. C Default: NiOUT=NiFFT C Names of input and output formatted files with the functional values: C FFTINR='string', FFTINI='string' ... real and imaginary parts of C the input function. C Default: FFTINR=' ' FFTINI=' ' C FFTOUTR='string', FFTOUTI='string' ... real and imaginary parts of C the output function. C Default: FFTOUTR=' ' FFTOUTI=' ' C For general description of the files with gridded data refer to C file forms.htm. C Optional parameters specifying the form of the quantities C written in the output formatted files: C MINDIG,MAXDIG=positive integers ... See the description in file C forms.for. C NUMLIN=positive integer ... Number of reals or integers in one C line of the output file. See the description in file C forms.for. C C----------------------------------------------------------------------- C C Subroutines and external functions required: EXTERNAL NFFT,INDRAM,MODF,NCHECK,ERROR,RSEP1,RSEP3T,RSEP3R,RSEP3I, *RARRAY,WARRAY,FOURN,UARRAY REAL UARRAY INTEGER NFFT,INDRAM,MODF C NFFT,INDRAM,MODF,NCHECK ... This file. C ERROR ... File C error.for. C RSEP1,RSEP3T,RSEP3R,RSEP3I ... File C sep.for. C RARRAY,WARRAY,UARRAY ... File C forms.for. C FOURN ... File 'fourn.for' of C Numerical Recipes. C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C Common block /GRDFF/: INTEGER N1FFT,N1N2FF COMMON /GRDFF/ N1N2FF,N1FFT SAVE /GRDFF/ C CHARACTER*80 FILSEP,FM4,FILINR,FILINI,FILOUR,FILOUI INTEGER LU1,LU2,LU3,LU4 PARAMETER (LU1=1,LU2=2,LU3=3,LU4=4) DOUBLE PRECISION PID PARAMETER (PID=3.141592653589793D0) REAL PI PARAMETER (PI=3.141592653589793) INTEGER MODFFT INTEGER N1,N2,N3,N4,N1N2,NN INTEGER N2FFT,N3FFT,NDIFFT(3),NNFFT,NN2FFT,NT1FFT,NT2FFT,NT3FFT INTEGER N1TMP,N2TMP,N3TMP INTEGER N1OUT,N2OUT,N3OUT,N1N2OU,NNOUT REAL UNDEF REAL O1,O2,O3,D1,D2,D3,FFT,FFTFIL, * O1OUT,O2OUT,O3OUT,D1OUT,D2OUT,D3OUT,D1TMP,D2TMP,D3TMP INTEGER IR,II,I1MI,I1MA,I2MI,I2MA,I3MI,I3MA,IO1,IO2,IO3 INTEGER K1MA,K2MA,K3MA,K1,K2,K3,M1,M2,M3 INTEGER IRAM,I1,I2,I3,I4,I,J,K,L,NDIMFF,NFORFF,OFORFF REAL RRA,RRB,RR0,RRD,RIA,RIB,RI0,RID,RRK DOUBLE PRECISION FFTD,O1D,O2D,O3D,D1D,D2D,D3D,M1D,M2D,M3D, * D1OUTD,D2OUTD,D3OUTD,O1OUTD,O2OUTD,O3OUTD, * AUXD,AUXRD,AUXID,RM0RD,RM0ID,RMTRD,RMTID,RMULTD C UNDEF=UARRAY() C C----------------------------------------------------------------------- C C Reading a name of the file with the input data: FILSEP=' ' WRITE(*,'(A)') '+GRDFFT: Enter input filename: ' READ(*,*) FILSEP IF (FILSEP.EQ.' ') THEN C GRDFFT-19 CALL ERROR('GRDFFT-19: No input file specified') C Input file in the form of the SEP (Stanford Exploration Project) C parameter or history file must be specified. C There is no default filename. ENDIF WRITE(*,'(A)') '+GRDFFT: Working... ' C C Reading all the data from the SEP file into the memory: CALL RSEP1(LU1,FILSEP) C C Reading the filenames of the files with the real and imaginary C parts of the input function: CALL RSEP3T('FFTINR',FILINR,' ') CALL RSEP3T('FFTINI',FILINI,' ') IF (FILINR.EQ.' ' .AND. FILINI.EQ.' ') THEN C GRDFFT-01 CALL ERROR('GRDFFT-01: No input files specified.') ENDIF C Reading the filenames of the output files: CALL RSEP3T('FFTOUTR',FILOUR,' ') CALL RSEP3T('FFTOUTI',FILOUI,' ') IF (FILOUR.EQ.' ' .AND. FILOUI.EQ.' ') THEN C GRDFFT-02 CALL ERROR('GRDFFT-02: No output files specified.') ENDIF C Reading the multiplicative constant FFT: CALL RSEP3R('FFT',FFT,2.*PI) FFTD=DBLE(FFT) I=NINT(FFT/PI) IF (I.NE.0) THEN IF (ABS((FFT/PI-FLOAT(I))/FLOAT(I)).LE.0.01) THEN FFT=FLOAT(I)*PI FFTD=DBLE(I)*PID ENDIF ENDIF C Mode of the FFT: IF (FFT.GT.0.) THEN MODFFT=1 ELSEIF (FFT.LT.0.) THEN MODFFT=-1 ELSE C GRDFFT-03 CALL ERROR('GRDFFT-03: Wrong value of FFT.') ENDIF C Reading the values describing the input grid: CALL RSEP3R('O1',O1,0.) CALL RSEP3R('O2',O2,0.) CALL RSEP3R('O3',O3,0.) CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) CALL RSEP3I('N4',N4,1) N1N2=N1*N2 NN=N1N2*N3 IF (N4.LE.0) THEN CALL RSEP3T('M4',FM4,' ') IF (FM4.EQ.' ') THEN N4=1 ELSE OPEN(LU1,FILE=FM4,STATUS='OLD') READ(LU1,*) N4 CLOSE(LU1) ENDIF ENDIF CALL RSEP3R('D1',D1,0.) CALL RSEP3R('D2',D2,0.) CALL RSEP3R('D3',D3,0.) IF (((D1.EQ.0.).AND.(N1.NE.1)).OR. * ((D2.EQ.0.).AND.(N2.NE.1)).OR. * ((D3.EQ.0.).AND.(N3.NE.1))) THEN C GRDFFT-04 CALL ERROR('GRDFFT-04: Wrong input grid.') ENDIF IF ((D1.EQ.0.).AND.(N1.EQ.1)) D1=1. IF ((D2.EQ.0.).AND.(N2.EQ.1)) D2=1. IF ((D3.EQ.0.).AND.(N3.EQ.1)) D3=1. N1TMP=NFFT(N1) N2TMP=NFFT(N2) N3TMP=NFFT(N3) CALL RSEP3I('N1FFT',NT1FFT,N1TMP) CALL RSEP3I('N2FFT',NT2FFT,N2TMP) CALL RSEP3I('N3FFT',NT3FFT,N3TMP) IF (NT1FFT.EQ.0) THEN N1FFT=N1 ELSE N1FFT=NT1FFT CALL NCHECK(N1FFT,N1TMP) ENDIF IF (NT2FFT.EQ.0) THEN N2FFT=N2 ELSE N2FFT=NT2FFT CALL NCHECK(N2FFT,N2TMP) ENDIF IF (NT3FFT.EQ.0) THEN N3FFT=N3 ELSE N3FFT=NT3FFT CALL NCHECK(N3FFT,N3TMP) ENDIF N1N2FF=N1FFT*N2FFT NNFFT=N1N2FF*N3FFT NN2FFT=2*NNFFT C Reading the constant FFTFIL for values missing in input grid: CALL RSEP3R('FFTFIL',FFTFIL,UNDEF) C Dimension of the temporary grid: NFORFF=NNFFT IF (NT1FFT.EQ.0) NFORFF=NFORFF/N1FFT IF (NT2FFT.EQ.0) NFORFF=NFORFF/N2FFT IF (NT3FFT.EQ.0) NFORFF=NFORFF/N3FFT OFORFF=MRAM-2*NFORFF C Reading the values describing the output grid: CALL RSEP3I('N1OUT',N1OUT,N1FFT) CALL RSEP3I('N2OUT',N2OUT,N2FFT) CALL RSEP3I('N3OUT',N3OUT,N3FFT) N1N2OU=N1OUT*N2OUT NNOUT=N1N2OU*N3OUT CALL RSEP3R('O1OUT',O1OUT,-1./(2.*D1)*2.*PI/ABS(FFT)) CALL RSEP3R('O2OUT',O2OUT,-1./(2.*D2)*2.*PI/ABS(FFT)) CALL RSEP3R('O3OUT',O3OUT,-1./(2.*D3)*2.*PI/ABS(FFT)) D1TMP=1./(N1FFT*D1)*2.*PI/ABS(FFT) D2TMP=1./(N2FFT*D2)*2.*PI/ABS(FFT) D3TMP=1./(N3FFT*D3)*2.*PI/ABS(FFT) CALL RSEP3R('D1OUT',D1OUT,D1TMP) CALL RSEP3R('D2OUT',D2OUT,D2TMP) CALL RSEP3R('D3OUT',D3OUT,D3TMP) D1TMP=D1OUT/D1TMP D2TMP=D2OUT/D2TMP D3TMP=D3OUT/D3TMP IF ((ABS(D1TMP-1.).GT.0.001).AND.(N1OUT.NE.1)) THEN IF ((NT1FFT.NE.0).OR.(D1OUT.NE.D1)) THEN C GRDFFT-05 CALL ERROR('GRDFFT-05: Wrong D1OUT.') ENDIF ENDIF IF ((ABS(D2TMP-1.).GT.0.001).AND.(N2OUT.NE.1)) THEN IF ((NT2FFT.NE.0).OR.(D2OUT.NE.D2)) THEN C GRDFFT-16 CALL ERROR('GRDFFT-16: Wrong D2OUT.') ENDIF ENDIF IF ((ABS(D3TMP-1.).GT.0.001).AND.(N3OUT.NE.1)) THEN IF ((NT3FFT.NE.0).OR.(D3OUT.NE.D3)) THEN C GRDFFT-17 CALL ERROR('GRDFFT-17: Wrong D3OUT.') ENDIF ENDIF C IF ((NN2FFT+2*MAX0(NN,NNOUT,NFORFF)).GT.MRAM) THEN C GRDFFT-06 CALL ERROR('GRDFFT-06: Small array RAM.') C The dimension MRAM in the common block RAMC defined in the file C ram.inc should be C at least 2*(N1FFT*N2FFT*N3FFT + C + MAX(N1*N2*N3,N1OUT*N2OUT*N3OUT,N1FFT*N2FFT*N3FFT). ENDIF C C Preparing the quantities needed in double precision: O1D=DBLE(O1) O2D=DBLE(O2) O3D=DBLE(O3) D1D=DBLE(D1) D2D=DBLE(D2) D3D=DBLE(D3) D1OUTD=DBLE(D1OUT) D2OUTD=DBLE(D2OUT) D3OUTD=DBLE(D3OUT) O1OUTD=DBLE(O1OUT) O2OUTD=DBLE(O2OUT) O3OUTD=DBLE(O3OUT) C Preparing number NDIMFF describing the dimension C of the part of the input grid to be transformed. NDIMFF=3 IF ((N1FFT.EQ.1).OR.(NT1FFT.EQ.0)) NDIMFF=NDIMFF-1 IF ((N2FFT.EQ.1).OR.(NT2FFT.EQ.0)) NDIMFF=NDIMFF-1 IF ((N3FFT.EQ.1).OR.(NT3FFT.EQ.0)) NDIMFF=NDIMFF-1 IF (NDIMFF.EQ.0) THEN C GRDFFT-07 CALL ERROR('GRDFFT-07: Input grid is 1*1*1.') ENDIF C Computing the multiplicative factor: RMULTD=1.D0 IF ((N1FFT.NE.1).AND.(NT1FFT.NE.0)) RMULTD=RMULTD*D1D IF ((N2FFT.NE.1).AND.(NT2FFT.NE.0)) RMULTD=RMULTD*D2D IF ((N3FFT.NE.1).AND.(NT3FFT.NE.0)) RMULTD=RMULTD*D3D RMULTD=RMULTD*DSQRT(DABS(FFTD)/(2.D0*PID))**NDIMFF C C Opening files with input and output grids: IF (N4.GT.1) THEN IF (FILINR.NE.' ') THEN OPEN(LU1,FILE=FILINR,FORM='FORMATTED',STATUS='OLD') ENDIF IF (FILINI.NE.' ') THEN OPEN(LU2,FILE=FILINI,FORM='FORMATTED',STATUS='OLD') ENDIF IF (FILOUR.NE.' ') THEN OPEN(LU3,FILE=FILOUR,FORM='FORMATTED') ENDIF IF (FILOUI.NE.' ') THEN OPEN(LU4,FILE=FILOUI,FORM='FORMATTED') ENDIF ENDIF C C Loop over N4 space grids: DO 90, I4=1,N4 C Reading the input function: IR=MRAM-2*NN+1 II=MRAM-NN+1 IF (FILINR.NE.' ') THEN IF (N4.GT.1) THEN CALL RARRAY(LU1,' ' ,'FORMATTED',.TRUE.,0.,NN,RAM(IR)) ELSE CALL RARRAY(LU1,FILINR,'FORMATTED',.TRUE.,0.,NN,RAM(IR)) ENDIF ELSE DO 10, I1=IR,II-1 RAM(I1)=0. 10 CONTINUE ENDIF IF (FILINI.NE.' ') THEN IF (N4.GT.1) THEN CALL RARRAY(LU2,' ' ,'FORMATTED',.TRUE.,0.,NN,RAM(II)) ELSE CALL RARRAY(LU2,FILINI,'FORMATTED',.TRUE.,0.,NN,RAM(II)) ENDIF ELSE DO 11, I1=II,MRAM RAM(I1)=0. 11 CONTINUE ENDIF WRITE(*,'(A)') '+GRDFFT: Working... ' C I=N1FFT-N1 M1=I/2 M1D=DBLE(M1) I1MI=1+M1 I1MA=M1+N1 I=N2FFT-N2 M2=I/2 M2D=DBLE(M2) I2MI=1+M2 I2MA=M2+N2 I=N3FFT-N3 M3=I/2 M3D=DBLE(M3) I3MI=1+M3 I3MA=M3+N3 IF ((I1MI.LT.1).OR.(I2MI.LT.1).OR.(I3MI.LT.1).OR. * (I1MI.GT.N1FFT).OR.(I2MI.GT.N2FFT).OR.(I3MI.GT.N3FFT).OR. * (I1MA.LT.1).OR.(I2MA.LT.1).OR.(I3MA.LT.1).OR. * (I1MA.GT.N1FFT).OR.(I2MA.GT.N2FFT).OR.(I3MA.GT.N3FFT)) THEN C GRDFFT-08 CALL ERROR('GRDFFT-08: Wrong value of IiMA or IiMI.') C This error should not appear. ENDIF C Recording the known values: IRAM=1 DO 24, I3=1,N3FFT DO 23, I2=1,N2FFT DO 22, I1=1,N1FFT IF (IRAM+1.GE.IR) THEN C GRDFFT-10 CALL ERROR('GRDFFT-10: Wrong index of array RAM.') C This error should not appear. ENDIF IF ((I1.GE.I1MI).AND.(I1.LE.I1MA).AND. * (I2.GE.I2MI).AND.(I2.LE.I2MA).AND. * (I3.GE.I3MI).AND.(I3.LE.I3MA)) THEN I=IR + I1-I1MI + (I2-I2MI)*N1 + (I3-I3MI)*N1N2 IF ((I.LT.IR).OR.(I.GE.II)) THEN C C GRDFFT-11 CALL ERROR('GRDFFT-11: Wrong index in the data array') C This error should not appear. ENDIF RAM(IRAM) =RAM(I) RAM(IRAM+1)=RAM(I+NN) ENDIF IRAM=IRAM+2 22 CONTINUE 23 CONTINUE 24 CONTINUE C Interpolating the unknown values in the first axis direction: DO 34, I3=I3MI,I3MA DO 33, I2=I2MI,I2MA RRA=RAM(INDRAM(I1MI,I2,I3)) RRB=RAM(INDRAM(I1MA,I2,I3)) RR0=(RRA+RRB)/2. RRD=(RRA-RRB)/2. RIA=RAM(INDRAM(I1MI,I2,I3)+1) RIB=RAM(INDRAM(I1MA,I2,I3)+1) RI0=(RIA+RIB)/2. RID=(RIA-RIB)/2. DO 31, I1=1,I1MI-1 IF (FFTFIL.EQ.UNDEF) THEN RRK=FLOAT(I1-1)/FLOAT(I1MI-1) RAM(INDRAM(I1,I2,I3))=RR0+RRK*RRD RAM(INDRAM(I1,I2,I3)+1)=RI0+RRK*RID ELSE RAM(INDRAM(I1,I2,I3))=FFTFIL RAM(INDRAM(I1,I2,I3)+1)=0. ENDIF 31 CONTINUE DO 32, I1=I1MA+1,N1FFT IF (FFTFIL.EQ.UNDEF) THEN RRK=FLOAT(I1-N1FFT)/FLOAT(N1FFT-I1MA) RAM(INDRAM(I1,I2,I3))=RR0+RRK*RRD RAM(INDRAM(I1,I2,I3)+1)=RI0+RRK*RID ELSE RAM(INDRAM(I1,I2,I3))=FFTFIL RAM(INDRAM(I1,I2,I3)+1)=0. ENDIF 32 CONTINUE 33 CONTINUE 34 CONTINUE C Interpolating the unknown values in the second axis direction: DO 44, I3=I3MI,I3MA DO 43, I1=1,N1FFT RRA=RAM(INDRAM(I1,I2MI,I3)) RRB=RAM(INDRAM(I1,I2MA,I3)) RR0=(RRA+RRB)/2. RRD=(RRA-RRB)/2. RIA=RAM(INDRAM(I1,I2MI,I3)+1) RIB=RAM(INDRAM(I1,I2MA,I3)+1) RI0=(RIA+RIB)/2. RID=(RIA-RIB)/2. DO 41, I2=1,I2MI-1 IF (FFTFIL.EQ.UNDEF) THEN RRK=FLOAT(I2-1)/FLOAT(I2MI-1) RAM(INDRAM(I1,I2,I3))=RR0+RRK*RRD RAM(INDRAM(I1,I2,I3)+1)=RI0+RRK*RID ELSE RAM(INDRAM(I1,I2,I3))=FFTFIL RAM(INDRAM(I1,I2,I3)+1)=0. ENDIF 41 CONTINUE DO 42, I2=I2MA+1,N2FFT IF (FFTFIL.EQ.UNDEF) THEN RRK=FLOAT(I2-N2FFT)/FLOAT(N2FFT-I2MA) RAM(INDRAM(I1,I2,I3))=RR0+RRK*RRD RAM(INDRAM(I1,I2,I3)+1)=RI0+RRK*RID ELSE RAM(INDRAM(I1,I2,I3))=FFTFIL RAM(INDRAM(I1,I2,I3)+1)=0. ENDIF 42 CONTINUE 43 CONTINUE 44 CONTINUE C Interpolating the unknown values in the third axis direction: DO 54, I1=1,N1FFT DO 53, I2=1,N2FFT RRA=RAM(INDRAM(I1,I2,I3MI)) RRB=RAM(INDRAM(I1,I2,I3MA)) RR0=(RRA+RRB)/2. RRD=(RRA-RRB)/2. RIA=RAM(INDRAM(I1,I2,I3MI)+1) RIB=RAM(INDRAM(I1,I2,I3MA)+1) RI0=(RIA+RIB)/2. RID=(RIA-RIB)/2. DO 51, I3=1,I3MI-1 IF (FFTFIL.EQ.UNDEF) THEN RRK=FLOAT(I3-1)/FLOAT(I3MI-1) RAM(INDRAM(I1,I2,I3))=RR0+RRK*RRD RAM(INDRAM(I1,I2,I3)+1)=RI0+RRK*RID ELSE RAM(INDRAM(I1,I2,I3))=FFTFIL RAM(INDRAM(I1,I2,I3)+1)=0. ENDIF 51 CONTINUE DO 52, I3=I3MA+1,N3FFT IF (FFTFIL.EQ.UNDEF) THEN RRK=FLOAT(I3-N3FFT)/FLOAT(N3FFT-I3MA) RAM(INDRAM(I1,I2,I3))=RR0+RRK*RRD RAM(INDRAM(I1,I2,I3)+1)=RI0+RRK*RID ELSE RAM(INDRAM(I1,I2,I3))=FFTFIL RAM(INDRAM(I1,I2,I3)+1)=0. ENDIF 52 CONTINUE 53 CONTINUE 54 CONTINUE C C Adding the before-FFT multiplicative factor: DO 58, I3=1,N3FFT DO 57, I2=1,N2FFT DO 56, I1=1,N1FFT RMTRD=1.D0 RMTID=0.D0 IF ((N1FFT.NE.1).AND.(NT1FFT.NE.0)) THEN AUXD=FFTD* * (DBLE(I1-1)*D1D*(O1OUTD-DBLE(NINT(O1OUT/D1OUT))*D1OUTD)) AUXRD=DCOS(AUXD) AUXID=DSIN(AUXD) RM0RD=RMTRD RM0ID=RMTID RMTRD=RM0RD*AUXRD - RM0ID*AUXID RMTID=RM0RD*AUXID + RM0ID*AUXRD ENDIF IF ((N2FFT.NE.1).AND.(NT2FFT.NE.0)) THEN AUXD=FFTD* * (DBLE(I2-1)*D2D*(O2OUTD-DBLE(NINT(O2OUT/D2OUT))*D2OUTD)) AUXRD=DCOS(AUXD) AUXID=DSIN(AUXD) RM0RD=RMTRD RM0ID=RMTID RMTRD=RM0RD*AUXRD - RM0ID*AUXID RMTID=RM0RD*AUXID + RM0ID*AUXRD ENDIF IF ((N3FFT.NE.1).AND.(NT3FFT.NE.0)) THEN AUXD=FFTD* * (DBLE(I3-1)*D3D*(O3OUTD-DBLE(NINT(O3OUT/D3OUT))*D3OUTD)) AUXRD=DCOS(AUXD) AUXID=DSIN(AUXD) RM0RD=RMTRD RM0ID=RMTID RMTRD=RM0RD*AUXRD - RM0ID*AUXID RMTID=RM0RD*AUXID + RM0ID*AUXRD ENDIF RM0RD=DBLE(RAM(INDRAM(I1,I2,I3) )) RM0ID=DBLE(RAM(INDRAM(I1,I2,I3)+1)) RAM(INDRAM(I1,I2,I3) )=SNGL(RM0RD*RMTRD - RM0ID*RMTID) RAM(INDRAM(I1,I2,I3)+1)=SNGL(RM0ID*RMTRD + RM0RD*RMTID) 56 CONTINUE 57 CONTINUE 58 CONTINUE C C C Computing the FFT: C Quantities describing the parts of the grid where FFT will C be done: NDIFFT(1)=N1FFT NDIFFT(2)=N2FFT NDIFFT(3)=N3FFT IF (NT1FFT.EQ.0) NDIFFT(1)=1 IF (NT2FFT.EQ.0) NDIFFT(2)=1 IF (NT3FFT.EQ.0) NDIFFT(3)=1 IF (NDIFFT(2).EQ.1) THEN NDIFFT(2)=NDIFFT(3) ENDIF IF (NDIFFT(1).EQ.1) THEN NDIFFT(1)=NDIFFT(2) NDIFFT(2)=NDIFFT(3) ENDIF NFORFF=NNFFT IF (NT1FFT.EQ.0) NFORFF=NFORFF/N1FFT IF (NT2FFT.EQ.0) NFORFF=NFORFF/N2FFT IF (NT3FFT.EQ.0) NFORFF=NFORFF/N3FFT OFORFF=MRAM-2*NFORFF K1MA=1 K2MA=1 K3MA=1 IF (NT1FFT.EQ.0) K1MA=N1FFT IF (NT2FFT.EQ.0) K2MA=N2FFT IF (NT3FFT.EQ.0) K3MA=N3FFT I1MI=1 I1MA=N1FFT I2MI=1 I2MA=N2FFT I3MI=1 I3MA=N3FFT C C FFT loop over subgrids: DO 69, K3=1,K3MA DO 68, K2=1,K2MA DO 67, K1=1,K1MA IF (NT1FFT.EQ.0) THEN I1MI=K1 I1MA=K1 ENDIF IF (NT2FFT.EQ.0) THEN I2MI=K2 I2MA=K2 ENDIF IF (NT3FFT.EQ.0) THEN I3MI=K3 I3MA=K3 ENDIF C Moving the subgrid to the temporary location: K=OFORFF-2 DO 63, I3=I3MI,I3MA DO 62, I2=I2MI,I2MA DO 61, I1=I1MI,I1MA I=INDRAM(I1,I2,I3) J=I+1 K=K+2 L=K+1 IF ((I.LE.0).OR.(I.GT.NN2FFT).OR. * (J.LE.0).OR.(J.GT.NN2FFT).OR. * (K.GT.MRAM).OR.(L.GT.MRAM)) THEN C GRDFFT-15 CALL ERROR('GRDFFT-15: Wrong index of array RAM.') C This error should not appear. ENDIF RAM(K)=RAM(I) RAM(L)=RAM(J) 61 CONTINUE 62 CONTINUE 63 CONTINUE C FFT: CALL FOURN(RAM(OFORFF),NDIFFT,NDIMFF,MODFFT) C Moving the subgrid back from the temporary location: K=OFORFF-2 DO 66, I3=I3MI,I3MA DO 65, I2=I2MI,I2MA DO 64, I1=I1MI,I1MA I=INDRAM(I1,I2,I3) J=I+1 K=K+2 L=K+1 IF ((I.LE.0).OR.(I.GT.NN2FFT).OR. * (J.LE.0).OR.(J.GT.NN2FFT).OR. * (K.GT.MRAM).OR.(L.GT.MRAM)) THEN C GRDFFT-18 CALL ERROR('GRDFFT-18: Wrong index of array RAM.') C This error should not appear. ENDIF RAM(I)=RAM(K) RAM(J)=RAM(L) 64 CONTINUE 65 CONTINUE 66 CONTINUE 67 CONTINUE 68 CONTINUE 69 CONTINUE C End of the FFT loop over subgrids. C C Reordering the computed values to the output field, C adding the multiplicative factor: IO1=MODF(NINT(O1OUT/D1OUT)+1,N1FFT) IF (IO1.LT.0) IO1=IO1+N1FFT IO2=MODF(NINT(O2OUT/D2OUT)+1,N2FFT) IF (IO2.LT.0) IO2=IO2+N2FFT IO3=MODF(NINT(O3OUT/D3OUT)+1,N3FFT) IF (IO3.LT.0) IO3=IO3+N3FFT IF ((IO1.LE.0).OR.(IO2.LE.0).OR.(IO3.LE.0).OR. * (IO1.GT.N1FFT).OR.(IO2.GT.N2FFT).OR.(IO3.GT.N3FFT)) THEN C GRDFFT-12 CALL ERROR('GRDFFT-12: Wrong value of IOi.') C This error should not appear. ENDIF IR=MRAM-2*NNOUT+1 II=MRAM-NNOUT+1 DO 74, I3=1,N3OUT DO 73, I2=1,N2OUT DO 72, I1=1,N1OUT I=INDRAM(MODF(IO1+I1-1,N1FFT),MODF(IO2+I2-1,N2FFT), * MODF(IO3+I3-1,N3FFT)) RMTRD=RMULTD RMTID=DBLE(0.) IF ((N1FFT.NE.1).AND.(NT1FFT.NE.0)) THEN AUXD=FFTD*(O1OUTD+DBLE(I1-1)*D1OUTD)*(O1D-M1D*D1D) AUXRD=DCOS(AUXD) AUXID=DSIN(AUXD) RM0RD=RMTRD RM0ID=RMTID RMTRD=RM0RD*AUXRD - RM0ID*AUXID RMTID=RM0RD*AUXID + RM0ID*AUXRD ENDIF IF ((N2FFT.NE.1).AND.(NT2FFT.NE.0)) THEN AUXD=FFTD*(O2OUTD+DBLE(I2-1)*D2OUTD)*(O2D-M2D*D2D) AUXRD=DCOS(AUXD) AUXID=DSIN(AUXD) RM0RD=RMTRD RM0ID=RMTID RMTRD=RM0RD*AUXRD - RM0ID*AUXID RMTID=RM0RD*AUXID + RM0ID*AUXRD ENDIF IF ((N3FFT.NE.1).AND.(NT3FFT.NE.0)) THEN AUXD=FFTD*(O3OUTD+DBLE(I3-1)*D3OUTD)*(O3D-M3D*D3D) AUXRD=DCOS(AUXD) AUXID=DSIN(AUXD) RM0RD=RMTRD RM0ID=RMTID RMTRD=RM0RD*AUXRD - RM0ID*AUXID RMTID=RM0RD*AUXID + RM0ID*AUXRD ENDIF RM0RD=DBLE(RAM(I)) RM0ID=DBLE(RAM(I+1)) RAM(IR)=SNGL(RM0RD*RMTRD - RM0ID*RMTID) RAM(II)=SNGL(RM0ID*RMTRD + RM0RD*RMTID) IR=IR+1 II=II+1 72 CONTINUE 73 CONTINUE 74 CONTINUE C C C Writing the results of the FFT: IR=MRAM-2*NNOUT+1 II=MRAM-NNOUT+1 IF (N4.GT.1) THEN IF (FILOUR.NE.' ') THEN CALL WARRAY(LU3,' ' ,'FORMATTED',.FALSE.,0.,.FALSE.,0., * NNOUT,RAM(IR)) ENDIF IF (FILOUI.NE.' ') THEN CALL WARRAY(LU4,' ' ,'FORMATTED',.FALSE.,0.,.FALSE.,0., * NNOUT,RAM(II)) ENDIF ELSE IF (FILOUR.NE.' ') THEN CALL WARRAY(LU3,FILOUR,'FORMATTED',.FALSE.,0.,.FALSE.,0., * NNOUT,RAM(IR)) ENDIF IF (FILOUI.NE.' ') THEN CALL WARRAY(LU4,FILOUI,'FORMATTED',.FALSE.,0.,.FALSE.,0., * NNOUT,RAM(II)) ENDIF ENDIF 90 CONTINUE C C Closing files with input and output grids: IF (N4.GT.1) THEN IF (FILINR.NE.' ') THEN CLOSE(LU1) ENDIF IF (FILINI.NE.' ') THEN CLOSE(LU2) ENDIF IF (FILOUR.NE.' ') THEN CLOSE(LU3) ENDIF IF (FILOUI.NE.' ') THEN CLOSE(LU4) ENDIF ENDIF WRITE(*,'(A)') '+GRDFFT: Finished. ' STOP END C C======================================================================= C INTEGER FUNCTION INDRAM(I,J,K) C Common block /GRDFF/: INTEGER N1FFT,N1N2FF COMMON /GRDFF/ N1N2FF,N1FFT SAVE /GRDFF/ INTEGER I,J,K INDRAM=2*((K-1)*N1N2FF + (J-1)*N1FFT + (I-1)) + 1 RETURN END C C======================================================================= C INTEGER FUNCTION MODF(I,J) INTEGER I,J MODF=MOD(I,J) IF (MODF.EQ.0) MODF=J RETURN END C C======================================================================= C INTEGER FUNCTION NFFT(N) INTEGER N IF (N.LE. 1) THEN NFFT= 1 ELSEIF (N.LE. 2) THEN NFFT= 2 ELSEIF (N.LE. 4) THEN NFFT= 4 ELSEIF (N.LE. 8) THEN NFFT= 8 ELSEIF (N.LE. 16) THEN NFFT= 16 ELSEIF (N.LE. 32) THEN NFFT= 32 ELSEIF (N.LE. 64) THEN NFFT= 64 ELSEIF (N.LE. 128) THEN NFFT= 128 ELSEIF (N.LE. 256) THEN NFFT= 256 ELSEIF (N.LE. 512) THEN NFFT= 512 ELSEIF (N.LE. 1024) THEN NFFT= 1024 ELSEIF (N.LE. 2048) THEN NFFT= 2048 ELSEIF (N.LE. 4096) THEN NFFT= 4096 ELSEIF (N.LE. 8192) THEN NFFT= 8192 ELSEIF (N.LE. 16384) THEN NFFT= 16384 ELSEIF (N.LE. 32768) THEN NFFT= 32768 ELSEIF (N.LE. 65536) THEN NFFT= 65536 ELSEIF (N.LE. 131072) THEN NFFT= 131072 ELSEIF (N.LE. 262144) THEN NFFT= 262144 ELSEIF (N.LE. 524288) THEN NFFT= 524288 ELSEIF (N.LE. 1048576) THEN NFFT= 1048576 ELSEIF (N.LE. 2097152) THEN NFFT= 2097152 ELSEIF (N.LE. 4194304) THEN NFFT= 4194304 ELSEIF (N.LE. 8388608) THEN NFFT= 8388608 ELSEIF (N.LE. 16777216) THEN NFFT= 16777216 ELSEIF (N.LE. 33554432) THEN NFFT= 33554432 ELSEIF (N.LE. 67108864) THEN NFFT= 67108864 ELSEIF (N.LE. 134217728) THEN NFFT= 134217728 ELSEIF (N.LE. 268435456) THEN NFFT= 268435456 ELSEIF (N.LE. 536870912) THEN NFFT= 536870912 ELSEIF (N.LE.1073741824) THEN NFFT=1073741824 ELSE C GRDFFT-09 CALL ERROR ('GRDFFT-09: Too large N') C One of the N1, N2, N3 or N4 is greater than 2**31. ENDIF cc Old: cc REAL AUX cc AUX=LOG10(FLOAT(N))/LOG10(2.) cc NFFT=2**INT(AUX+0.999999) RETURN END C C======================================================================= C SUBROUTINE NCHECK(N1,N2) INTEGER N1,N2 EXTERNAL NFFT INTEGER NFFT IF (N1.LT.N2) THEN C GRDFFT-13 CALL ERROR * ('GRDFFT-13: Wrong specification of the grid for FFT') C Number of gridpoints for FFT must be greater than or equal to C the number of gridpoints in data. ENDIF IF (N1.NE.NFFT(N1)) THEN C GRDFFT-14 CALL ERROR * ('GRDFFT-14: Wrong specification of the grid for FFT') C Number of gridpoints for FFT must be an integer power of 2. ENDIF RETURN END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'forms.for' C forms.for INCLUDE 'length.for' C length.for INCLUDE 'fourn.for' C fourn.for of Numerical Recipes C C======================================================================= C