diff --git a/Make_include_template b/Make_include_template
index 8353d6e971869c7337e17e7fb5bcfe4d60f88049..09dbb7fa49ee24a9e81c9d057696fcda6595c193 100644
--- a/Make_include_template
+++ b/Make_include_template
@@ -82,8 +82,7 @@ OPTC += -fopenmp
#LDFLAGS = -static -Ofast
#############################################################################
-# BLAS and LAPACK libraries
-# and FFT LIBRARIES
+# BLAS and LAPACK libraries and FFT LIBRARIES
#MKLROOT=/opt/intel/mkl/
MKLLIB=${MKLROOT}/lib
OPTC += -DMKL -I$(MKLROOT)/include -I$(MKLROOT)/include/fftw
@@ -98,7 +97,8 @@ BLAS = -L${MKLROOT}/lib/ -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm
#for intel compilers
#BLAS = -mkl
#
-#If you do not want to use MKL
+#If you do not want or can not use MKL you can try the basic lapack and blas libaries that are default installed on most Linux
+#systems
#BLAS = -llapack -lblas
#############################################################################
@@ -108,7 +108,7 @@ BLAS = -L${MKLROOT}/lib/ -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm
#BLAS = -L$(AMDROOT)/lib -lacml -lfortran -lffio -lrt -lm
#LIBARIES
-LIBS= -L$L -lgenfft $(BLAS)
+LIBS= -L$L -lgenfft $(BLAS)
########################################################################
diff --git a/README b/README
index 4587443e9c43046de6df1731342f2d43c36d0874..3c02538e179d943feb3c34633eeccc89210fd476 100644
--- a/README
+++ b/README
@@ -13,19 +13,20 @@ For more information, go to http://software.seg.org/2017/00XX .
You must read and accept usage terms at:
http://software.seg.org/disclaimer.txt before use.
-DOI reference of release
-1.4.0 https://zenodo.org/badge/latestdoi/23060862
REFERENCES
---------
--1- If the Finite Difference code has helped you in your research please refer to our paper in your publications:
+-0- DOI reference of this software release
+https://zenodo.org/badge/latestdoi/23060862
+
+-1- If the Finite Difference code has helped you in your research please refer to the following paper in your publications:
Finite-difference modeling experiments for seismic interferometry
Jan Thorbecke and Deyan Draganov
2011, Geophysics, Vol. 76, no. 6 (November-December); p H1--H18, doi: 10.1190/GEO2010-0039.1
Download: https://janth.home.xs4all.nl/Publications/Articles/ThorbeckeDraganov2012.pdf
--2- If the Machenko code has helped you in your research please refer to our paper in your publications:
+-2- If the Machenko code has helped you in your research please refer to this paper in your publications:
Implementation of the Marchenko method
Jan Thorbecke, Evert Slob, Joeri Brackenhoff, Joost van der Neut, and Kees Wapenaar
@@ -39,6 +40,15 @@ Lele Zhang1 and Evert Slob2
2019, Geophysics, Vol. 84, no. 1 (January-February); p. A7-A11, doi: 10.1190/GEO2018-0548.1
Download: http://homepage.tudelft.nl/t4n4v/BeyondInterferometry/geo_19h.pdf
+-4- If you use the fdacrtmc code of Max Holicki please refer to the following paper:
+
+Acoustic directional snapshot wavefield decomposition
+Holicki, M., Drijkoningen, G., and Wapenaar, K., 2019, Acoustic directional snapshot wavefield decomposition: Geophysical
+Prospecting, Vol. 67, 32-51
+Download: http://homepage.tudelft.nl/t4n4v/4_Journals/Geophys.Prosp/GP_19a.pdf
+
+-5-
+
INSTALLATION
-------------
@@ -99,19 +109,22 @@ The MDD kernels depend on BLAS and LAPACK calls. Free downloads of these librari
https://www.netlib.org/blas/index.html
https://www.netlib.org/lapack/index.html
+
MKL libraries
-------------
If you are running on x64_86 processors you can download (for free) Intel's highly optimised MKL package:
https://software.intel.com/en-us/mkl/choose-download
-Usually the MKL libraries are installed in $MKL_ROOT. If that variable is not set, you can try to find the correct path by searching for one of the libraries:
+These Libraries include highly optimised libraries of BLAS, LAPACK, and FFT(W).
+Usually the MKL libraries are installed in $MKLROOT. If that variable is not set, you can try to find the correct path by searching for one of the libraries:
find /opt/intel -name libmkl_gf_lp64.so
and adjust MKLROOT in Make_include accordingly.
You can also completely disable the use of MKL by commenting out the MKL parts in Make_include.
+
SEISMIC UNIX
-----------
If you want to use the .su files with SU from CWP:
@@ -126,12 +139,28 @@ so the XDRFLAG is empty. If you have made that change you have to remake SU with
make remake
make xtremake
+
ZFP
---
-The fdacrtmc code makes use of ZFP compression to store the snaphots in CPU memory. This package is included in this repository for your convenience. The latest package and detailed explanation can be found on:
+The fdacrtmc and the 3D Marchenko code makes use of ZFP compression to store snaphots in CPU memory or do effient IO. This package is included in this repository for your convenience. The latest package and detailed explanation can be found on:
https://github.com/LLNL/zfp
+FDACRTMC
+--------
+fdacrtmc uses FFTW and the wisdom computations are stored on disk for re-usage. This directory is defined in fdacrtmc.h
+
+#ifndef WISDOMDIR
+#define WISDOMDIR "/tmp/fftw/"
+#endif
+
+For the code to run properly the directory /tmp/fftw/ must exist. It is also possible to compile the fdacrtmc code that defines
+WISDOMDIR. In fdacrtmc/Makefile add:
+
+CFLAGS += -DWISDOMDIR="/directory/that/exists"
+
+or you can change the name of WISDOMDIR in fdacrtmc.h
+
MISC
----
@@ -147,7 +176,7 @@ The latest version of the source code and manual can be found at:
http://www.xs4all.nl/~janth/Software/Software.html
- or at github:
+or at github:
git clone https://github.com/JanThorbecke/OpenSource.git
git clone git://github.com/JanThorbecke/OpenSource.git
diff --git a/marchenko/Makefile b/marchenko/Makefile
index 848b6a8624741950e794ba9e25da1e43e5ab8213..cd1c5416851cc624e7455be23ff1bf08c6a6a00a 100644
--- a/marchenko/Makefile
+++ b/marchenko/Makefile
@@ -42,6 +42,7 @@ SRCP = marchenko_primaries.c \
readData.c \
readShotData.c \
readTinvData.c \
+ findFirstBreak.c \
writeData.c \
writeDataIter.c \
wallclock_time.c \
diff --git a/marchenko/demo/oneD/primaries.scr b/marchenko/demo/oneD/primaries.scr
index e11a01563aa2872d7b539f9b7694c2b64fe31a7f..b5398efee445d97e3c0a7ca5e7a326c9b3fb01c5 100755
--- a/marchenko/demo/oneD/primaries.scr
+++ b/marchenko/demo/oneD/primaries.scr
@@ -1,6 +1,6 @@
#!/bin/bash -x
#SBATCH -J marchenko_primaries
-#SBATCH --cpus-per-task=40
+#SBATCH --cpus-per-task=10
#SBATCH --ntasks=1
#SBATCH --time=2:40:00
@@ -8,7 +8,7 @@
#cd $SLURM_SUBMIT_DIR
export PATH=$HOME/src/OpenSource/bin:$PATH:
-export OMP_NUM_THREADS=40
+export OMP_NUM_THREADS=10
#shot record to remove internal multiples
select=451
@@ -17,7 +17,7 @@ makewave fp=20 dt=0.004 file_out=wave.su nt=2048 t0=0.0 scale=0 scfft=1
../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
nshots=901 verbose=2 istart=40 iend=500 fmax=90 pad=1024 \
- niter=31 smooth=10 niterec=2 niterskip=1 shift=20 file_rr=pred_rr.su T=0 file_update=update.su
+ niter=31 smooth=10 niterec=2 niterskip=50 shift=20 file_rr=pred_rr.su T=0 file_update=update.su
exit;
diff --git a/marchenko/demo/oneD/primariesPlane.scr b/marchenko/demo/oneD/primariesPlane.scr
index 0ab5c64d70e0541f87459979f760c0386452840e..510bb9ab88871334f05dc8c4d39106be42ad73ed 100755
--- a/marchenko/demo/oneD/primariesPlane.scr
+++ b/marchenko/demo/oneD/primariesPlane.scr
@@ -1,6 +1,6 @@
#!/bin/bash -x
#SBATCH -J marchenko_primaries
-#SBATCH --cpus-per-task=40
+#SBATCH --cpus-per-task=10
#SBATCH --ntasks=1
#SBATCH --time=2:40:00
@@ -8,20 +8,20 @@
#cd $SLURM_SUBMIT_DIR
export PATH=$HOME/src/OpenSource/bin:$PATH:
-export OMP_NUM_THREADS=40
+export OMP_NUM_THREADS=10
#shot record to remove internal multiples
select=451
makewave fp=20 dt=0.004 file_out=wave.su nt=1024 t0=0.0 scale=0 scfft=1
-#../../marchenko_primaries file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=10 \
+#../../marchenko_primaries file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=30 \
# nshots=901 verbose=2 istart=276 iend=277 fmax=90 pad=1024 \
-# xorig=900 niter=2 smooth=10 niterskip=1 file_iter=iterplane.su shift=20 file_rr=plane2_10_rr.su T=0
+# niter=2 smooth=10 niterskip=1 file_iter=iterplane.su shift=20 file_rr=plane2_10_rr.su T=0
#exit
-../../marchenko_primaries file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=-10 \
+../../marchenko_primaries file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=30.0 \
nshots=901 verbose=2 istart=40 iend=1000 fmax=90 pad=1024 \
- niter=20 smooth=10 niterskip=40 niterec=0 shift=20 file_rr=plane2_10_rr.su T=0
+ niter=20 smooth=10 niterskip=60 niterec=2 shift=20 file_rr=plane2_10_rr.su T=0
exit
diff --git a/marchenko/findFirstBreak.c b/marchenko/findFirstBreak.c
new file mode 100644
index 0000000000000000000000000000000000000000..100974d81a47eb914d4c9d3e197bcd0e1b40e697
--- /dev/null
+++ b/marchenko/findFirstBreak.c
@@ -0,0 +1,109 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <assert.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+int findFirstBreak(float *shot, int nx, int nt, int ishot, float *maxval, int tr, int hw, int verbose)
+{
+ int i, j;
+ int jmax, tstart, tend;
+ float xmax, tmax, lmax;
+
+ /* find consistent (one event) maximum related to maximum value */
+
+ /* find global maximum
+ xmax=0.0;
+ for (i = 0; i < nx; i++) {
+ tmax=0.0;
+ jmax = 0;
+ for (j = 0; j < nt; j++) {
+ lmax = fabs(shot[i*nt+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ if (lmax > xmax) {
+ ishot = i;
+ xmax=lmax;
+ }
+ }
+ }
+ maxval[i] = jmax;
+ }
+ */
+
+ /* alternative find maximum at source position ishot */
+ tmax=0.0;
+ jmax = 0;
+ if (tr == 0) {
+ for (j = 0; j < nt; j++) {
+ lmax = fabs(shot[ishot*nt+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ if (lmax > xmax) {
+ xmax=lmax;
+ }
+ }
+ }
+ }
+ else {
+ for (j = nt-1; j >= 0; j--) {
+ lmax = fabs(shot[ishot*nt+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ if (lmax > xmax) {
+ xmax=lmax;
+ }
+ }
+ }
+ }
+ maxval[ishot] = jmax;
+ if (verbose >= 3) vmess("Mute max at src-trace %d is sample %d", ishot, maxval[ishot]);
+
+ /* search forward in trace direction from maximum in file */
+ for (i = ishot+1; i < nx; i++) {
+ tstart = MAX(0, (maxval[i-1]-hw));
+ tend = MIN(nt-1, (maxval[i-1]+hw));
+ jmax=tstart;
+ tmax=0.0;
+ for(j = tstart; j <= tend; j++) {
+ lmax = fabs(shot[i*nt+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ }
+ }
+ maxval[i] = jmax;
+ }
+ /* search backward in trace direction from maximum in file */
+ for (i = ishot-1; i >=0; i--) {
+ tstart = MAX(0, (maxval[i+1]-hw));
+ tend = MIN(nt-1, (maxval[i+1]+hw));
+ jmax=tstart;
+ tmax=0.0;
+ for(j = tstart; j <= tend; j++) {
+ lmax = fabs(shot[i*nt+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ }
+ }
+ maxval[i] = jmax;
+ }
+ if (tr != 0) {
+ for (i = 0; i < nx; i++) maxval[i] = nt-1 - maxval[i];
+ }
+
+ return 0;
+}
+
diff --git a/marchenko/marchenko_primaries.c b/marchenko/marchenko_primaries.c
index cc98f668ddf199dcc92923c998d2a4bdf480b6a6..d893a330715027cf0e5069297341b693a8a82bc0 100644
--- a/marchenko/marchenko_primaries.c
+++ b/marchenko/marchenko_primaries.c
@@ -29,6 +29,7 @@ int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx
int readTinvData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nfoc, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
+int findFirstBreak(float *shot, int nx, int nt, int ishot, float *maxval, int tr, int hw, int verbose);
int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int t0shift, int iter);
int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
int readData(FILE *fp, float *data, segy *hdrs, int n1);
@@ -110,7 +111,7 @@ int main (int argc, char **argv)
FILE *fp_out, *fp_rr, *fp_w, *fp_up;
size_t nread;
int i, j, k, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath;
- int size, n1, n2, ntap, tap, di, ntraces;
+ int size, n1, n2, ntap, tap, di, ntraces, tr;
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
int reci, countmin, mode, ixa, ixb, n2out, verbose, ntfft;
int iter, niter, niterec, recur, niterskip, niterrun, tracf, *muteW;
@@ -121,9 +122,9 @@ int main (int argc, char **argv)
double t0, t1, t2, t3, t4, tsyn, tread, tfft, tcopy, tii;
double energyMi, *energyM0;
float tt0, d1, d2, f1, f2, fxsb, fxse, ft, fx, *xsyn, dxsrc;
- float *M0, *DD, *RR, *SRC, dt, dx, dxs, scl, mem;
+ float *M0, *DD, *RR, *SRC, dt, dx, dxs, scl, mem, scltap;
float *rtrace, *tmpdata, *f1min, *f1plus, *RMi, *Mi, *trace;
- float *Mup, *Msp;
+ float *Mup, *Msp, *maxval;
float xmin, xmax, scale, tsq;
float Q, f0, *ixmask, *costaper;
float src_velo, src_angle, grad2rad, p, *twplane;
@@ -228,20 +229,20 @@ int main (int argc, char **argv)
//if(!getparint("xorig", &xorig)) xorig=-(nxs-1)/2;
/* compute time delay for plane-wave responses */
- twplane = (float *) calloc(nxs,sizeof(float));
+ twplane = (float *) calloc(nxs,sizeof(float)); /* initialize with zeros */
if (plane_wave==1) {
grad2rad = 17.453292e-3;
p = sin(src_angle*grad2rad)/src_velo;
if (p < 0.0) {
for (i=0; i<nxs; i++) {
twplane[i] = fabsf((nxs-i-1)*dxs*p)+tt0;
+ if (verbose >=3) vmess("plane-wave delay-time i=%d x=%f t=%f", i, dxs*(nxs-i-1), twplane[i]);
}
}
else {
for (i=0; i<nxs; i++) {
twplane[i] = (i)*dxs*p+tt0;
- //twplane[i] = dxs*(xorig+i)*tan(src_angle*grad2rad)/src_velo;
-// fprintf(stderr,"plane-wave i=%d x=%f t=%f %f\n", i, dxs*(xorig+i), twplane[i], (xorig+i)*dxs*p);
+ if (verbose >=3) vmess("plane-wave delay-time i=%d x=%f t=%f", i, dxs*(i), twplane[i]);
}
}
}
@@ -367,9 +368,9 @@ int main (int argc, char **argv)
if (tap == 2 || tap == 3) {
tapersh = (float *)malloc(nx*sizeof(float));
- for (j = 0; j < ntap; j++)
- tapersh[j] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
- for (j = ntap; j < nx-ntap; j++)
+ for (j = 1; j < ntap; j++)
+ tapersh[j-1] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
+ for (j = ntap-1; j < nx-ntap; j++)
tapersh[j] = 1.0;
for (j = nx-ntap; j < nx; j++)
tapersh[j] =(cos(PI*(j-(nx-ntap))/ntap)+1)/2.0;
@@ -382,7 +383,6 @@ int main (int argc, char **argv)
}
}
}
- free(tapersh);
}
/* check consistency of header values */
@@ -413,6 +413,7 @@ int main (int argc, char **argv)
scl = 1.0/((float)2.0*ntfft);
rtrace = (float *)calloc(ntfft,sizeof(float));
ctrace = (complex *)calloc(nfreq+1,sizeof(complex));
+
for (i = 0; i < xnx[ishot]; i++) {
for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
ctrace[j].r = Refl[ishot*nw*nx+m*nx+i].r*cwave[j].r + Refl[ishot*nw*nx+m*nx+i].i*cwave[j].i;
@@ -423,9 +424,22 @@ int main (int argc, char **argv)
for (j = 0; j < nts; j++) {
DD[0*nxs*nts+i*nts+j] = -1.0*scl*rtrace[j];
}
+ /* remove taper at edges for selected shot record */
+ //if (tap == 2 || tap == 3) scltap= scl/tapersh[i];
+ //for (j = 0; j < nts; j++) {
+ // DD[0*nxs*nts+i*nts+j] = -1.0*scltap*rtrace[j];
+ //}
}
-
- /* construct plane wave from all shot records */
+ /* set timereversal for searching First Break */
+ /* for experimenting with non flat truncation windows */
+ //maxval = (float *)calloc(nxs,sizeof(float));
+ //tr = 1;
+ //findFirstBreak(DD, xnx[ishot], nts, ishot, maxval, tr, hw, verbose);
+ //twplane = (float *) calloc(nxs,sizeof(float));
+ //for (i=0; i<nxs; i++) twplane[i] = dt*(maxval[i]-maxval[ishot]);
+ //for (i = 0; i < xnx[ishot]; i++) fprintf(stderr,"maxval[%d] = %f tw=%f\n", i, dt*maxval[i], twplane[i]);
+
+ /* construct plane wave (time reversed and multiplid with -1) from all shot records */
if (plane_wave==1) {
for (l=0; l<nshots; l++) {
memset(ctrace, 0, sizeof(complex)*(nfreq+1));
@@ -451,6 +465,7 @@ int main (int argc, char **argv)
tom = j*deltom*twplane[l];
csum.r = cos(-tom);
csum.i = sin(-tom);
+ /* Optional add wavelet to SRC-field */
//cwav.r = csum.r*cwave[j].r - csum.i*cwave[j].i;
//cwav.i = csum.i*cwave[j].r + csum.r*cwave[j].i;
ctrace[j] = csum;
@@ -523,8 +538,21 @@ int main (int argc, char **argv)
vmess("number of time samples fft nt nts = %d %d %d", ntfft, nt, nts);
vmess("time sampling = %e ", dt);
vmess("smoothing taper for time-window = %d ", smooth);
+ if (plane_wave) {
+ vmess("*** Plane-wave processing selected *** ");
+ vmess("plane-wave angle = %f ", src_angle);
+ vmess("plane-wave velocity = %f ", src_velo);
+ vmess("plane-wave t0-shift = %f ", tt0);
+ vmess("plane-wave DD input vector = DDplane%03d.su ", src_angle);
+ vmess("plane-wave SRC for migration = SRCplane%03d.su ", src_angle);
+ }
if (file_rr != NULL) vmess("RR output file = %s ", file_rr);
- if (file_iter != NULL) vmess("Iterations output file = %s ", file_iter);
+ if (file_iter != NULL) {
+ vmess("Iterations output file = %s ", file_iter);
+ vmess("Initialisation input data = M0_%06d.su ", 1000*istart);
+ vmess("f1min intermediate array = f1min_%03d'iter'.su ", istart);
+ vmess("Mi intermediate array = Mi_%03d'iter'.su ", istart);
+ }
}
/*================ initializations ================*/
@@ -589,7 +617,15 @@ int main (int argc, char **argv)
if (plane_wave) {
writeDataIter("SRCplane.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, NINT(src_angle));
- writeDataIter("DDplane.su", DD, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, NINT(src_angle));
+ /* make DD causal again and undo the -1 multplication */
+ for (l=0; l<nshots; l++) {
+ j=0;
+ SRC[l*nts+j] = -1.0*DD[l*nts+j];
+ for (j = 1; j < nts; j++) {
+ SRC[l*nts+j] = -1.0*DD[l*nts+nts-j];
+ }
+ }
+ writeDataIter("DDplane.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, NINT(src_angle));
}
free(SRC);
@@ -813,8 +849,10 @@ int main (int argc, char **argv)
for (i = 0; i < npos; i++) {
ix = ixpos[i];
iw = NINT((ii*dt+twplane[ix])/dt);
- RR[l*nxs*nts+i*nts+iw] = f1min[l*nxs*nts+i*nts+iw];
- if (file_update != NULL) Msp[l*nxs*nts+i*nts+iw] = Mup[l*nxs*nts+i*nts+iw];
+ if ( iw<nts && iw>=0 ) {
+ RR[l*nxs*nts+i*nts+iw] = f1min[l*nxs*nts+i*nts+iw];
+ if (file_update != NULL) Msp[l*nxs*nts+i*nts+iw] = Mup[l*nxs*nts+i*nts+iw];
+ }
}
}
diff --git a/utils/fconv.c b/utils/fconv.c
index 42686ef4e70c2cb3ef144109161ec7f8af2cf7d6..d44c92aaac7372422ee41f86260f09a1a82d57f3 100644
--- a/utils/fconv.c
+++ b/utils/fconv.c
@@ -507,7 +507,7 @@ void corr(float *data1, float *data2, float *cov, int nrec, int nsam, float dt,
/* inverse frequency-time FFT and scale result */
sign = 1;
- scl = 1.0/(float)optn;
+ scl = dt/(float)optn;
crmfft(&ccov[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
scl_data(rdata1,optn,nrec,scl,cov,nsam);
@@ -578,7 +578,7 @@ void corr3(float *data1, float *data2, float *cov, int nrec, int nsam, float dt)
/* inverse frequency-time FFT and scale result */
sign = 1;
- scl = 1.0/(float)optn;
+ scl = dt/(float)optn;
crmfft(&ccov[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
scl_data(rdata1,optn,nrec,scl,cov,ntout);
@@ -689,7 +689,7 @@ void deconv(float *data1, float *data2, float *decon, int nrec, int nsam,
/* inverse frequency-time FFT and scale result */
sign = 1;
- scl = 1.0/(float)optn;
+ scl = dt/(float)optn;
crmfft(&cdec[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
scl_data(rdata1,optn,nrec,scl,decon,nsam);
@@ -777,7 +777,7 @@ void power(float *data1, float *data2, float *cov, int nrec, int nsam, float dt,
/* inverse frequency-time FFT and scale result */
sign = 1;
- scl = 1.0/(float)optn;
+ scl = dt/(float)optn;
crmfft(&acov[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
scl_data(rdata1,optn,nrec,scl,cov,nsam);
@@ -866,7 +866,7 @@ void convol(float *data1, float *data2, float *con, int nrec, int nsam, float dt
/* inverse frequency-time FFT and scale result */
sign = 1;
- scl = 1.0/((float)(optn));
+ scl = dt/((float)(optn));
crmfft(&ccon[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
scl_data(rdata1,optn,nrec,scl,con,nsam);
@@ -979,7 +979,7 @@ void cohr(float *data1, float *data2, float *cov, int nrec, int nsam, float dt,
/* inverse frequency-time FFT and scale result */
sign = 1;
- scl = 1.0/(float)optn;
+ scl = dt/(float)optn;
crmfft(&ccov[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
scl_data(rdata1,optn,nrec,scl,cov,nsam);