diff --git a/Make_include_OSX b/Make_include_OSX
index ffcf0c7f947965e4fb6c9be8d77a80cd9fd22bd5..f84793d29350a7f7e0d20d47ca0f8bd42571af07 100644
--- a/Make_include_OSX
+++ b/Make_include_OSX
@@ -19,12 +19,12 @@ L = $(ROOT)/lib
#GNU
#CC = gcc-mp-5
-CC = gcc
-FC = gfortran
+#CC = gcc
+#FC = gfortran
# Linux gcc version 4.x
-OPTC = -O3 -ffast-math
+#OPTC = -O3 -ffast-math
#to include parallelisation with OpenMP
-OPTC += -fopenmp
+#OPTC += -fopenmp
# for better performing code you can try:
#OPTC = -O3 -fno-trapping-math -ffast-math -funroll-all-loops -mavx -fomit-frame-pointer -mfpmath=sse -ftree-vectorizer-verbose=1
@@ -83,7 +83,9 @@ LDFLAGS = -Ofast
#############################################################################
# BLAS and LAPACK libraries
+# and FFT LIBRARIES
MKLLIB=${MKLROOT}/lib
+OPTC += -DMKL -I$(MKLROOT)/include -I$(MKLROOT)/include/fftw
#for GNU compilers
#you might need to add intel64 to : ${MKLROOT}/lib/intel64
#BLAS = -L${MKLROOT}/lib/ -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
@@ -96,27 +98,14 @@ BLAS = -Wl,-rpath ${MKLLIB} -L${MKLROOT}/lib/ -lmkl_intel_lp64 -lmkl_sequential
#BLAS = -mkl
#############################################################################
-# FOR FFT LIBRARIES
-#AMD ACML 4.4.0
+# AMD ACML 4.4.0
#AMDROOT = /home/thorbcke/amdsdk/v1.0/acml/open64_64
#OPTC += -DACML440 -I$(AMDROOT)/include
#BLAS = -L$(AMDROOT)/lib -lacml -lfortran -lffio -lrt -lm
#Intel MKL
-MKLLIB=${MKLROOT}/lib
-OPTC += -DMKL -I$(MKLROOT)/include -I$(MKLROOT)/include/fftw
-#for GNU compilers
-#FFT = -Wl,-rpath ${MKLLIB} -L${MKLLIB} -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
-#for clang on OSX
-FFT = -L${MKLROOT}/lib/ -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
-#on linux you want to use groups and MKL is in lib/intel64
-#MKLLIB=${MKLROOT}/lib/intel64
-#FFT = -Wl,-rpath ${MKLLIB} -Wl,--start-group -L${MKLLIB} -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -Wl,--end-group -lpthread -lm -ldl
-#for Intel compilers
-#FFT = -Wl,-rpath ${MKLLIB} -L${MKLLIB} -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
#LIBARIES
-LIBS= -L$L -lgenfft $(FFT) $(BLAS)
-
+LIBS= -L$L -lgenfft $(BLAS)
########################################################################
# standard CFLAGS
diff --git a/Make_include_template b/Make_include_template
index ad3a45193b691d93f2508603b6bac3cd4f6b9ff6..002c45dc914840e33369eeefdd8cf28af97a87e4 100644
--- a/Make_include_template
+++ b/Make_include_template
@@ -83,13 +83,15 @@ OPTC += -fopenmp
#############################################################################
# BLAS and LAPACK libraries
-MKLROOT=/opt/intel/mkl/
+# and FFT LIBRARIES
+#MKLROOT=/opt/intel/mkl/
MKLLIB=${MKLROOT}/lib
+OPTC += -DMKL -I$(MKLROOT)/include -I$(MKLROOT)/include/fftw
#for GNU compilers
#you might need to add intel64 to : ${MKLROOT}/lib/intel64
BLAS = -L${MKLROOT}/lib/ -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
-#for GNU on OSX
-#BLAS = -L${MKLROOT}/lib/ -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
+#for GNU/clang on OSX
+#BLAS = -Wl,-rpath ${MKLLIB} -L${MKLROOT}/lib/ -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
#on linux you want to use groups and MKL is in lib/intel64
#MKLLIB=${MKLROOT}/lib/intel64
#BLAS = -Wl,-rpath ${MKLLIB} -Wl,--start-group -L${MKLLIB} -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -Wl,--end-group -lpthread -lm -ldl
@@ -97,27 +99,13 @@ BLAS = -L${MKLROOT}/lib/ -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm
#BLAS = -mkl
#############################################################################
-# FOR FFT LIBRARIES
-#AMD ACML 4.4.0
+# AMD ACML 4.4.0
#AMDROOT = /home/thorbcke/amdsdk/v1.0/acml/open64_64
#OPTC += -DACML440 -I$(AMDROOT)/include
#BLAS = -L$(AMDROOT)/lib -lacml -lfortran -lffio -lrt -lm
-#Intel MKL
-MKLROOT=/opt/intel/mkl/
-MKLLIB=${MKLROOT}/lib
-OPTC += -DMKL -I$(MKLROOT)/include
-#for GNU compilers
-FFT = -Wl,-rpath ${MKLLIB} -L${MKLLIB} -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
-#for GNU on OSX
-#FFT = -L${MKLROOT}/lib/ -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
-#on linux you want to use groups and MKL is in lib/intel64
-#MKLLIB=${MKLROOT}/lib/intel64
-#FFT = -Wl,-rpath ${MKLLIB} -Wl,--start-group -L${MKLLIB} -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -Wl,--end-group -lpthread -lm -ldl
-#for Intel compilers
-#FFT = -Wl,-rpath ${MKLLIB} -L${MKLLIB} -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
#LIBARIES
-LIBS= -L$L -lgenfft $(FFT) $(BLAS)
+LIBS= -L$L -lgenfft $(BLAS)
########################################################################
diff --git a/README b/README
index 8dd89ff767db5b607815a3ab35059e82eb80ebef..4ec806b5c552611cf6831df4f69b9b6eb459f0a2 100644
--- a/README
+++ b/README
@@ -117,12 +117,18 @@ SEISMIC UNIX
If you want to use the .su files with SU from CWP:
git clone https://github.com/JohnWStockwellJr/SeisUnix
-==> Please make sure that SU is compiled without XDR (in $CWPROOT/Makefile.config make sure that XDRFLAG is NOT set). The SU output files of fdelmodc are all based on local IEEE data.
+==> Please make sure that SU is compiled without XDR (in $CWPROOT/Makefile.config make sure that XDRFLAG is NOT enabled). The SU output files of fdelmodc are all based on local IEEE data.
+To exclude the XDRFLAG in SU you have to use the following line in $CWPROOT/src/Makefile.config around line 35:
+
+XDRFLAG =
+
+so the XDRFLAG is empty. If you have made that change you have to remake SU with the commands:
+make remake
+make xtremane
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 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:
https://github.com/LLNL/zfp
diff --git a/fdelmodc/applySource.c b/fdelmodc/applySource.c
index ecd1bd13fef06ba5627f6118d17f7d3c72b9ef26..3aeb972c1404aaa91a8e9b3edf17685edea4e40a 100644
--- a/fdelmodc/applySource.c
+++ b/fdelmodc/applySource.c
@@ -143,6 +143,8 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
/* Force source */
+#pragma omp critical
+{
if (src.type == 6) {
vx[ix*n1+iz] += src_ampl*rox[ix*n1+iz]/(l2m[ix*n1+iz]);
/* stable implementation from "Numerical Techniques for Conservation Laws with Source Terms" by Justin Hudson */
@@ -338,6 +340,7 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
txz[ix*n1+iz] -= Mxz*src_ampl;
} /* src.type */
} /* ischeme */
+}
} /* loop over isrc */
return 0;
diff --git a/fdelmodc/getParameters.c b/fdelmodc/getParameters.c
index edf157883b26f007055e567fec5225bc541c9747..d882ba97c05a56c85c4314bc10ba87480ddcdb25 100644
--- a/fdelmodc/getParameters.c
+++ b/fdelmodc/getParameters.c
@@ -283,7 +283,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
/* Check stability and dispersion setting */
- if (cp_max > dx*stabfactor/dt) {
+ if (cp_max > (dx*stabfactor)/dt) {
vwarn("************ ! Stability ! ****************");
vwarn("From the input file maximum P-wave velocity");
vwarn("in the current model is %f !!", cp_max);
diff --git a/fdelmodc/getWaveletInfo.c b/fdelmodc/getWaveletInfo.c
index 2f3734aae6c38e54653fab909ec5e936a157d8ce..efeee9b4e9bfb8904f0b403f2881c6e0f4ed96ba 100644
--- a/fdelmodc/getWaveletInfo.c
+++ b/fdelmodc/getWaveletInfo.c
@@ -41,7 +41,7 @@ int getWaveletInfo(char *file_src, int *n1, int *n2, float *d1, float *d2, float
off_t bytes;
int ret, one_shot, ntraces;
int optn, nfreq, i, iwmax;
- float *trace;
+ float *trace, df;
float ampl, amplmax, tampl, tamplmax;
complex *ctrace;
segy hdr;
@@ -73,8 +73,9 @@ int getWaveletInfo(char *file_src, int *n1, int *n2, float *d1, float *d2, float
/* check to find out number of traces in shot gather */
- optn = optncr(*n1);
+ optn = optncr(*n1);
nfreq = optn/2 + 1;
+ df = 1.0/(optn*(*d1));
ctrace = (complex *)malloc(nfreq*sizeof(complex));
one_shot = 1;
trace = (float *)malloc(optn*sizeof(float));
@@ -107,12 +108,13 @@ int getWaveletInfo(char *file_src, int *n1, int *n2, float *d1, float *d2, float
iwmax = i;
}
}
+
/* from the maximum amplitude position look for the largest frequency
* which has an amplitude 400 times weaker than the maximum amplitude */
for (i=iwmax;i<nfreq;i++) {
ampl = sqrt(ctrace[i].r*ctrace[i].r+ctrace[i].i*ctrace[i].i);
if (400*ampl < amplmax) {
- *fmax = (i-1)*(1.0/(optn*(*d1)));
+ *fmax = (i-1)*df;
break;
}
}
diff --git a/marchenko/Makefile b/marchenko/Makefile
index b125699273c1db2562f3da2eac5d0090bb8e6229..7fadd2cca106c94b3ae012eb6aaa78128c749097 100644
--- a/marchenko/Makefile
+++ b/marchenko/Makefile
@@ -6,7 +6,7 @@ include ../Make_include
#OPTC += -g -O0 -Wall
-ALL: fmute marchenko marchenko_primaries
+ALL: fmute marchenko #marchenko_primaries
SRCJ = fmute.c \
getFileInfo.c \
@@ -35,20 +35,20 @@ SRCH = marchenko.c \
docpkge.c \
getpars.c
-SRCP = marchenko_primaries.c \
- synthesis.c \
- getFileInfo.c \
- readData.c \
- readShotData.c \
- readTinvData.c \
- writeData.c \
- writeDataIter.c \
- wallclock_time.c \
- name_ext.c \
- verbosepkg.c \
- atopkge.c \
- docpkge.c \
- getpars.c
+#SRCP = marchenko_primaries.c \
+# synthesis.c \
+# getFileInfo.c \
+# readData.c \
+# readShotData.c \
+# readTinvData.c \
+# writeData.c \
+# writeDataIter.c \
+# wallclock_time.c \
+# name_ext.c \
+# verbosepkg.c \
+# atopkge.c \
+# docpkge.c \
+# getpars.c
OBJJ = $(SRCJ:%.c=%.o)
diff --git a/marchenko/demo/oneD/epsModel.scr b/marchenko/demo/oneD/epsModel.scr
index 5ae0b460f468bf00cb8804d8882d6fa35a4f7885..6fddac7469d8998493c1c41577775adf32709f5b 100755
--- a/marchenko/demo/oneD/epsModel.scr
+++ b/marchenko/demo/oneD/epsModel.scr
@@ -43,26 +43,14 @@ sufft < wavefw.su | suamp | sugain scale=$dt | supsgraph \
label1="frequency (1/s)" label2="amplitude" \
d1num=10 wbox=6 hbox=3 x1end=125 x2end=1.1 > wavefw_freq.eps
-
#shot record
file=shot5_rp.su
sumax < $file mode=abs outpar=nep
-clipref=`cat nep | awk '{print $1/3}'`
+clipref=`cat nep | awk '{print $1/15}'`
suwind key=gx min=-2250000 max=2250000 < $file | \
- supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
- label1="time (s)" label2="lateral distance (m)" \
- n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+ label1="time sample number" label2="lateral distance (m)" \
+ n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=500 d1=1 f1=0 d1num=100 \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > shotx0_rp.eps
-#Initial focusing operator
-file=iniFocus_rp.su
-file_base=${file%.su}
-sumax < $file mode=abs outpar=nep
-clipref=`cat nep | awk '{print $1/3}'`
-suwind key=gx min=-2250000 max=2250000 < $file | \
- supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
- label1="time (s)" label2="lateral distance (m)" \
- n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
- f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
-
rm nep
diff --git a/marchenko/demo/oneD/iterations.scr b/marchenko/demo/oneD/iterations.scr
index 07839ed9456ba7c0b8d3284183b59f4a4194fff2..f8245068c6d7c759ec5cee25c1a862202dd924ea 100755
--- a/marchenko/demo/oneD/iterations.scr
+++ b/marchenko/demo/oneD/iterations.scr
@@ -3,22 +3,23 @@
#second reflector at time:
# 800/1800+600/2300
# .70531400966183574878 => sample 176
-#thrids reflector at time:
+#third reflector at time:
# 800/1800+600/2300+800/2000
# 1.10531400966183574878 sample 276
select=451
+for istart in 276
#for istart in 176 200 276 400
-for istart in 176
+#for istart in 246 256 266 276 286 296 306 316
do
(( iend = istart + 1 ))
../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
- nshots=601 verbose=1 istart=$istart iend=$iend fmax=90 \
- pad=1024 niter=45 smooth=10 niterskip=600 shift=20 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su
+ nshots=901 verbose=2 istart=$istart iend=$iend fmax=90 \
+ pad=0 niter=4 smooth=10 niterskip=600 shift=20 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su
-suximage < f1min_${istart}043.su clip=1 title="${istart}" &
+#suximage < f1min_${istart}043.su clip=1 title="${istart}" &
done
diff --git a/marchenko/demo/oneD/primariesFrame.scr b/marchenko/demo/oneD/primariesFrame.scr
index 5017c4b908a441e7cec67ef68cdb183adeb513a5..30e656b6fbaeec0891a805583324a619aeb2b82d 100755
--- a/marchenko/demo/oneD/primariesFrame.scr
+++ b/marchenko/demo/oneD/primariesFrame.scr
@@ -23,19 +23,19 @@ fconv file_in1=shotsx0.su file_in2=wave.su file_out=shotw.su
basop file_in=shotw.su choice=5 | sugain scale=-1 > DD.su
nts=1024
-itime=300 # select time sample to work on 1.2 seconds
+itime=276 # select time sample to work on 1.2 seconds
shift=20
#mute time nts-ii+shift to compute G_d
(( itmax = nts-itime+shift ))
-suzero itmax=$itmax < DD.su > G_d.su
+#suzero itmax=$itmax < DD.su > G_d.su
#f1min = -DD(-t) = shotw.su
-cp shotw.su f1min0.su
+#cp shotw.su f1min0.su
#first iteration
-cp G_d.su N0.su
+cp G_d_276000.su N0.su
#compute R*N0
-fconv file_in1=shotsdx5_rp.su file_in2=N0.su file_out=fconvN1.su verbose=1 fmax=90
+fconv file_in1=shotsdx5_rp.su file_in2=N0.su file_out=fconvN1.su verbose=1 fmax=120
suwind key=fldr min=451 max=451 < fconvN1.su | suximage
@@ -49,6 +49,7 @@ suzero itmax=$itmax < N1.su > N1m.su
exit
+
#display
file=fconvN1.su
file_base=${file%.su}
diff --git a/marchenko/marchenko_primaries.c b/marchenko/marchenko_primaries.c
deleted file mode 100644
index 80c23983eca7e283ff611ec2e8fa2e149f2f72dc..0000000000000000000000000000000000000000
--- a/marchenko/marchenko_primaries.c
+++ /dev/null
@@ -1,756 +0,0 @@
-#include "par.h"
-#include "segy.h"
-#include <time.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <math.h>
-#include <assert.h>
-#include <genfft.h>
-
-int omp_get_max_threads(void);
-int omp_get_num_threads(void);
-void omp_set_num_threads(int num_threads);
-
-#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))
-
-#ifndef COMPLEX
-typedef struct _complexStruct { /* complex number */
- float r,i;
-} complex;
-#endif/* complex */
-
-int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int nshots, int nx, int nxs, float fxsb, float dxs, int ntfft, int mode, float scale, float tsq, float Q, float f0, int reci, int *nshots_r, int *isxcount, int *reci_xsrc, int *reci_xrcv, float *ixmask, int verbose);
-
-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 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);
-int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
-int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
-double wallclock_time(void);
-
-void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int
-Nfoc, float *xrcv, float *xsrc, int *xnx, float fxse, float fxsb, float dxs, float dxsrc, float dx, int ntfft, int
-nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int npos, double *tfft, int *isxcount, int
-*reci_xsrc, int *reci_xrcv, float *ixmask, int verbose);
-
-void synthesisPositions(int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nfoc, float *xrcv, float *xsrc, int *xnx, float fxse, float fxsb, float dxs, float dxsrc, float dx, int nshots, int *ixpos, int *npos, int *isxcount, int countmin, int reci, int verbose);
-
-
-
-/*********************** self documentation **********************/
-char *sdoc[] = {
-" ",
-" MARCHENKO_primaries - Iterative primary reflections retrieval",
-" ",
-" marchenko_primaries file_tinv= file_shot= [optional parameters]",
-" ",
-" Required parameters: ",
-" ",
-" file_tinv= ............... shot-record from R to remove internal multiples",
-" file_shot= ............... Reflection response: R",
-" ",
-" Optional parameters: ",
-" ",
-" INTEGRATION ",
-" ishot=nshots/2 ........... shot number(s) to remove internal multiples ",
-" file_src=spike ........... convolve ishot(s) with source wavelet",
-" file_tinv= ............... use file_tinv to remove internal multiples",
-" COMPUTATION",
-" tap=0 .................... lateral taper focusing(1), shot(2) or both(3)",
-" ntap=0 ................... number of taper points at boundaries",
-" fmin=0 ................... minimum frequency in the Fourier transform",
-" fmax=70 .................. maximum frequency in the Fourier transform",
-" file_src= ................ optional source wavelet to convolve selected ishot's",
-" MARCHENKO ITERATIONS ",
-" niter=22 ................. number of iterations to inialize and restart",
-" niterec=2 ................ number of iterations in recursive part of the time-samples",
-" niterskip=50 ............. restart scheme each niterskip time-samples with niter iterations",
-" istart=20 ................ start sample of iterations for primaries",
-" iend=nt .................. end sample of iterations for primaries",
-" MUTE-WINDOW ",
-" shift=20 ................. number of points to account for wavelet (epsilon in papers)",
-" smooth=5 ................. number of points to smooth mute with cosine window",
-" REFLECTION RESPONSE CORRECTION ",
-" tsq=0.0 .................. scale factor n for t^n for true amplitude recovery",
-" Q=0.0 .......,............ Q correction factor",
-" f0=0.0 ................... ... for Q correction factor",
-" scale=2 .................. scale factor of R for summation of Ni with G_d",
-" pad=0 .................... amount of samples to pad the reflection series",
-" reci=0 ................... 1; add receivers as shots 2; only use receivers as shot positions",
-" countmin=0 ............... 0.3*nxrcv; minumum number of reciprocal traces for a contribution",
-" OUTPUT DEFINITION ",
-" file_rr= ................. output file with primary only shot record",
-" file_iter= ............... output file with -Ni(-t) for each iteration",
-" T=0 ...................... :1 compute transmission-losses compensated primaries ",
-" verbose=0 ................ silent option; >0 displays info",
-" ",
-" ",
-" author : Lele Zhang & Jan Thorbecke : 2019 ",
-" ",
-NULL};
-/**************** end self doc ***********************************/
-
-int main (int argc, char **argv)
-{
- FILE *fp_out, *fp_rr, *fp_w;
- 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 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;
- int hw, ii, ishot, istart, iend;
- int smooth, *ixpos, npos, ix, m, pad, T, perc;
- int nshots_r, *isxcount, *reci_xsrc, *reci_xrcv, shift;
- float fmin, fmax, *tapersh, *tapersy, fxf, dxf, *xsrc, *xrcv, *zsyn, *zsrc, *xrcvsyn;
- double t0, t1, t2, t3, t4, tsyn, tread, tfft, tcopy, tii;
- float d1, d2, f1, f2, fxsb, fxse, ft, fx, *xsyn, dxsrc;
- float *G_d, *DD, *RR, dt, dx, dxs, scl, mem;
- float *rtrace, *tmpdata, *f1min, *f1plus, *iRN, *Ni, *trace;
- float xmin, xmax, scale, tsq;
- float Q, f0, *ixmask, *costaper;
- complex *Refl, *Fop, *ctrace, *cwave;
- char *file_tinv, *file_shot, *file_rr, *file_src, *file_iter;
- segy *hdrs_out, hdr;
-
- initargs(argc, argv);
- requestdoc(1);
-
- tsyn = tread = tfft = tcopy = tii = 0.0;
- t0 = wallclock_time();
-
- if (!getparstring("file_shot", &file_shot)) file_shot = NULL;
- if (!getparstring("file_tinv", &file_tinv)) file_tinv = NULL;
- if(!getparstring("file_src", &file_src)) file_src = NULL;
- if (!getparstring("file_rr", &file_rr)) verr("parameter file_rr not found");
- if (!getparstring("file_iter", &file_iter)) file_iter = NULL;
-
- if (!getparint("verbose", &verbose)) verbose = 0;
- if (!getparfloat("fmin", &fmin)) fmin = 0.0;
- if (!getparfloat("fmax", &fmax)) fmax = 70.0;
- if (!getparint("reci", &reci)) reci = 0;
- reci=0; // source-receiver reciprocity is not yet fully build into the code
- if (!getparfloat("scale", &scale)) scale = 2.0;
- if (!getparfloat("Q", &Q)) Q = 0.0;
- if (!getparfloat("tsq", &tsq)) tsq = 0.0;
- if (!getparfloat("f0", &f0)) f0 = 0.0;
- if (!getparint("tap", &tap)) tap = 0;
- if (!getparint("ntap", &ntap)) ntap = 0;
- if (!getparint("pad", &pad)) pad = 0;
- if (!getparint("T", &T)) T = 0;
- if (T>0) T=-1;
- else T=1;
-
-
- if(!getparint("niter", &niter)) niter = 22;
- if(!getparint("niterec", &niterec)) niterec = 2;
- if(!getparint("niterskip", &niterskip)) niterskip = 50;
- if(!getparint("hw", &hw)) hw = 15;
- if(!getparint("shift", &shift)) shift=20;
- if(!getparint("smooth", &smooth)) smooth = 5;
- if(!getparint("ishot", &ishot)) ishot=300;
- if (reci && ntap) vwarn("tapering influences the reciprocal result");
-
- smooth = MIN(smooth, shift);
- if (smooth) {
- costaper = (float *)malloc(smooth*sizeof(float));
- scl = M_PI/((float)smooth);
- for (i=0; i<smooth; i++) {
- costaper[i] = 0.5*(1.0+cos((i+1)*scl));
- }
- }
-
-/*================ Reading info about shot and initial operator sizes ================*/
-
- if (file_tinv != NULL) { /* G_d is read from file_tinv */
- ngath = 0; /* setting ngath=0 scans all traces; n2 contains maximum traces/gather */
- ret = getFileInfo(file_tinv, &n1, &n2, &ngath, &d1, &d2, &f1, &f2, &xmin, &xmax, &scl, &ntraces);
- Nfoc = ngath;
- nxs = n2;
- nts = n1;
- dxs = d2;
- fxsb = f2;
- }
-
- ngath = 0; /* setting ngath=0 scans all traces; nx contains maximum traces/gather */
- ret = getFileInfo(file_shot, &nt, &nx, &ngath, &d1, &dx, &ft, &fx, &xmin, &xmax, &scl, &ntraces);
- nshots = ngath;
-
- if (!getparfloat("dt", &dt)) dt = d1;
- if(!getparint("istart", &istart)) istart=20;
- if(!getparint("iend", &iend)) iend=nt;
-
- if (file_tinv == NULL) {/* 'G_d' is one of the shot records */
- if(!getparint("ishot", &ishot)) ishot=1+(nshots-1)/2;
- ishot -= 1; /* shot numbering starts at 0 */
- Nfoc = 1;
- nxs = nx;
- nts = nt;
- dxs = dx;
- fxsb = fx;
- }
- assert (nxs >= nshots);
-
- ntfft = optncr(MAX(nt+pad, nts+pad));
- nfreq = ntfft/2+1;
- nw_low = (int)MIN((fmin*ntfft*dt), nfreq-1);
- nw_low = MAX(nw_low, 1);
- nw_high = MIN((int)(fmax*ntfft*dt), nfreq-1);
- nw = nw_high - nw_low + 1;
- scl = 1.0/((float)ntfft);
- if (!getparint("countmin", &countmin)) countmin = 0.3*nx;
-
-/*================ Allocating all data arrays ================*/
-
- Fop = (complex *)calloc(nxs*nw*Nfoc,sizeof(complex));
- f1min = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
- f1plus = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
- iRN = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
- Ni = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
- G_d = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
- DD = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
- RR = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
- trace = (float *)malloc(ntfft*sizeof(float));
- ixpos = (int *)malloc(nxs*sizeof(int));
- xrcvsyn = (float *)calloc(Nfoc*nxs,sizeof(float));
- xsyn = (float *)malloc(Nfoc*sizeof(float));
- zsyn = (float *)malloc(Nfoc*sizeof(float));
- xnxsyn = (int *)calloc(Nfoc,sizeof(int));
-
- Refl = (complex *)malloc(nw*nx*nshots*sizeof(complex));
- xsrc = (float *)calloc(nshots,sizeof(float));
- zsrc = (float *)calloc(nshots,sizeof(float));
- xrcv = (float *)calloc(nshots*nx,sizeof(float));
- xnx = (int *)calloc(nshots,sizeof(int));
-
- if (reci!=0) {
- reci_xsrc = (int *)malloc((nxs*nxs)*sizeof(int));
- reci_xrcv = (int *)malloc((nxs*nxs)*sizeof(int));
- isxcount = (int *)calloc(nxs,sizeof(int));
- ixmask = (float *)calloc(nxs,sizeof(float));
- }
-
-/*================ Read focusing operator(s) ================*/
-/* G_d = p_0^+ = G_d (-t) ~ Tinv */
-
- if (file_tinv != NULL) { /* G_d is named DD */
- muteW = (int *)calloc(Nfoc*nxs,sizeof(int));
- mode=-1; /* apply complex conjugate to read in data */
- readTinvData(file_tinv, xrcvsyn, xsyn, zsyn, xnxsyn, Nfoc, nxs, ntfft,
- mode, muteW, DD, hw, verbose);
- /* reading data added zero's to the number of time samples to be the same as ntfft */
- nts = ntfft;
- free(muteW);
-
- /* check consistency of header values */
- if (xrcvsyn[0] != 0 || xrcvsyn[1] != 0 ) fxsb = xrcvsyn[0];
- fxse = fxsb + (float)(nxs-1)*dxs;
- dxf = (xrcvsyn[nxs-1] - xrcvsyn[0])/(float)(nxs-1);
- if (NINT(dxs*1e3) != NINT(fabs(dxf)*1e3)) {
- vmess("dx in hdr.d1 (%.3f) and hdr.gx (%.3f) not equal",d2, dxf);
- if (dxf != 0) dxs = fabs(dxf);
- vmess("dx in operator => %f", dxs);
- }
- }
-
-/* ========================= Opening optional wavelet file ====================== */
-
- cwave = (complex *)calloc(ntfft,sizeof(complex));
- if (file_src != NULL){
- if (verbose) vmess("Reading wavelet from file %s.", file_src);
-
- fp_w = fopen(file_src, "r");
- if (fp_w == NULL) verr("error on opening input file_src=%s", file_src);
- nread = fread(&hdr, 1, TRCBYTES, fp_w);
- assert (nread == TRCBYTES);
- tmpdata = (float *)malloc(hdr.ns*sizeof(float));
- nread = fread(tmpdata, sizeof(float), hdr.ns, fp_w);
- assert (nread == hdr.ns);
- fclose(fp_w);
-
- /* Check dt wavelet is the same as reflection data */
- if ( rint(dt*1000) != rint(hdr.dt/1000.0) ) {
- vwarn("file_src dt != dt of file_shot %e != %e", hdr.dt/1e6, dt);
- }
- rtrace = (float *)calloc(ntfft,sizeof(float));
-
- /* add zero's to the number of time samples to be the same as ntfft */
- /* Add zero-valued samples in middle */
- if (hdr.ns <= ntfft) {
- for (i = 0; i < hdr.ns/2; i++) {
- rtrace[i] = tmpdata[i];
- rtrace[ntfft-1-i] = tmpdata[hdr.ns-1-i];
- }
- }
- else {
- vwarn("file_src has more samples than reflection data: truncated to ntfft = %d samples in middle are removed ", ntfft);
- for (i = 0; i < ntfft/2; i++) {
- rtrace[i] = tmpdata[i];
- rtrace[ntfft-1-i] = tmpdata[hdr.ns-1-i];
- }
- }
-
- rc1fft(rtrace, cwave, ntfft, -1);
- free(tmpdata);
- free(rtrace);
- }
- else {
- for (i = 0; i < nfreq; i++) cwave[i].r = 1.0;
- }
-
-
-/*================ Reading shot records ================*/
-
- mode=1;
- readShotData(file_shot, xrcv, xsrc, zsrc, xnx, Refl, nw, nw_low, ngath, nx, nx, fxsb, dxs, ntfft,
- mode, scale, tsq, Q, f0, reci, &nshots_r, isxcount, reci_xsrc, reci_xrcv, ixmask,verbose);
-
- 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++)
- tapersh[j] = 1.0;
- for (j = nx-ntap; j < nx; j++)
- tapersh[j] =(cos(PI*(j-(nx-ntap))/ntap)+1)/2.0;
- if (verbose) vmess("Taper for shots applied ntap=%d", ntap);
- for (l = 0; l < nshots; l++) {
- for (j = 1; j < nw; j++) {
- for (i = 0; i < nx; i++) {
- Refl[l*nx*nw+j*nx+i].r *= tapersh[i];
- Refl[l*nx*nw+j*nx+i].i *= tapersh[i];
- }
- }
- }
- free(tapersh);
- }
-
- /* check consistency of header values */
- fxf = xsrc[0];
- if (nx > 1) dxf = (xrcv[nx-1] - xrcv[0])/(float)(nx-1);
- else dxf = d2;
- if (NINT(dx*1e3) != NINT(fabs(dxf)*1e3)) {
- vmess("dx in hdr.d1 (%.3f) and hdr.gx (%.3f) not equal",dx, dxf);
- if (dxf != 0) dx = fabs(dxf);
- else verr("gx hdrs not set");
- vmess("dx used => %f", dx);
- }
-
- dxsrc = (float)xsrc[1] - xsrc[0];
- if (dxsrc == 0) {
- vwarn("sx hdrs are not filled in!!");
- dxsrc = dx;
- }
-
-/*================ Defining focusing operator(s) from R ================*/
-/* G_d = -R(ishot,-t)*/
-
- /* use ishot from Refl, complex-conjugate(time reverse), scale with -1 and convolve with wavelet */
- if (file_tinv == NULL) {
- if (verbose) vmess("Selecting G_d from Refl of %s", file_shot);
- nts = ntfft;
-
- 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;
- ctrace[j].i = -Refl[ishot*nw*nx+m*nx+i].i*cwave[j].r + Refl[ishot*nw*nx+m*nx+i].r*cwave[j].i;;
- }
- /* transfrom result back to time domain */
- cr1fft(ctrace, rtrace, ntfft, 1);
- for (j = 0; j < nts; j++) {
- DD[0*nxs*nts+i*nts+j] = -1.0*scl*rtrace[j];
- }
- }
- free(ctrace);
- free(rtrace);
-
- xsyn[0] = xsrc[ishot];
- zsyn[0] = zsrc[ishot];
- xnxsyn[0] = xnx[ishot];
- fxse = fxsb = xrcv[0];
- /* check consistency of header values */
- for (l=0; l<nshots; l++) {
- for (i = 0; i < nx; i++) {
- fxsb = MIN(xrcv[l*nx+i],fxsb);
- fxse = MAX(xrcv[l*nx+i],fxse);
- }
- }
- dxf = dx;
- dxs = dx;
- }
-
- /* define tapers to taper edges of acquisition */
- if (tap == 1 || tap == 3) {
- tapersy = (float *)malloc(nxs*sizeof(float));
- for (j = 0; j < ntap; j++)
- tapersy[j] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
- for (j = ntap; j < nxs-ntap; j++)
- tapersy[j] = 1.0;
- for (j = nxs-ntap; j < nxs; j++)
- tapersy[j] =(cos(PI*(j-(nxs-ntap))/ntap)+1)/2.0;
- if (verbose) vmess("Taper for operator applied ntap=%d", ntap);
- for (l = 0; l < Nfoc; l++) {
- for (i = 0; i < nxs; i++) {
- for (j = 0; j < nts; j++) {
- DD[l*nxs*nts+i*nts+j] *= tapersy[i];
- }
- }
- }
- free(tapersy);
- }
-
-/*================ Check the size of the files ================*/
-
- if (NINT(dxsrc/dx)*dx != NINT(dxsrc)) {
- vwarn("source (%.2f) and receiver step (%.2f) don't match",dxsrc,dx);
- if (reci == 2) vwarn("step used from operator (%.2f) ",dxs);
- }
- di = NINT(dxf/dxs);
- if ((NINT(di*dxs) != NINT(dxf)) && verbose)
- vwarn("dx in receiver (%.2f) and operator (%.2f) don't match",dx,dxs);
- if (verbose) {
- vmess("Number of receivers in focusop = %d", nxs);
- vmess("number of shots = %d", nshots);
- vmess("number of receiver/shot = %d", nx);
- vmess("first model position = %.2f", fxsb);
- vmess("last model position = %.2f", fxse);
- vmess("first source position fxf = %.2f", fxf);
- vmess("source distance dxsrc = %.2f", dxsrc);
- vmess("last source position = %.2f", fxf+(nshots-1)*dxsrc);
- vmess("receiver distance dxf = %.2f", dxf);
- vmess("direction of increasing traces = %d", di);
- 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 (file_rr != NULL) vmess("RR output file = %s ", file_rr);
- if (file_iter != NULL) vmess("Iterations output file = %s ", file_iter);
- }
-
-/*================ initializations ================*/
-
- if (reci) n2out = nxs;
- else n2out = nshots;
- mem = Nfoc*n2out*ntfft*sizeof(float)/1048576.0;
- if (verbose) {
- vmess("Time-sample range processed = %d:%d", istart, iend);
- vmess("number of output traces = %d", n2out);
- vmess("number of output samples = %d", ntfft);
- vmess("Size of output data/file = %.1f MB", mem);
- }
-
- ixa=0; ixb=0;
-
- synthesisPositions(nx, nt, nxs, nts, dt, xsyn, Nfoc, xrcv, xsrc, xnx, fxse, fxsb,
- dxs, dxsrc, dx, nshots, ixpos, &npos, isxcount, countmin, reci, verbose);
-
- if (verbose) {
- vmess("synthesisPositions: nshots=%d npos=%d", nshots, npos);
- }
-
-/*================ set variables for output data ================*/
-
- n1 = nts;
- n2 = n2out;
- f1 = ft;
- f2 = fxsb+dxs*ixpos[0];
- d1 = dt;
- if (reci == 0) d2 = dxsrc;
- else if (reci == 1) d2 = dxs;
- else if (reci == 2) d2 = dx;
-
- hdrs_out = (segy *) calloc(n2,sizeof(segy));
- if (hdrs_out == NULL) verr("allocation for hdrs_out");
- size = nxs*nts;
-
- for (i = 0; i < n2; i++) {
- hdrs_out[i].ns = n1;
- hdrs_out[i].trid = 1;
- hdrs_out[i].dt = dt*1000000;
- hdrs_out[i].f1 = f1;
- hdrs_out[i].f2 = f2;
- hdrs_out[i].d1 = d1;
- hdrs_out[i].d2 = d2;
- hdrs_out[i].trwf = n2out;
- hdrs_out[i].scalco = -1000;
- hdrs_out[i].gx = NINT(1000*(f2+i*d2));
- hdrs_out[i].scalel = -1000;
- hdrs_out[i].tracl = i+1;
- }
- t1 = wallclock_time();
- tread = t1-t0;
- if(verbose) {
- vmess("*******************************************");
- vmess("***** Computing Marchenko for all steps****");
- vmess("*******************************************");
- fprintf(stderr," %s: Progress: %3d%%",xargv[0],0);
- }
- perc=(iend-istart)/100;if(!perc)perc=1;
-
-/*================ start loop over number of time-samples ================*/
-
- for (ii=istart; ii<iend; ii++) {
-
-/*================ initialization ================*/
-
- /* once every 'niterskip' time-steps start from fresh G_d and do niter (~20) iterations */
- if ( ((ii-istart)%niterskip==0) || (ii==istart) ) {
- niterrun=niter;
- recur=0;
- if (verbose>2) vmess("Doing %d iterations to reset recursion at time-sample %d\n",niterrun,ii);
- for (l = 0; l < Nfoc; l++) {
- for (i = 0; i < nxs; i++) {
- for (j = 0; j < nts; j++) {
- G_d[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+i*nts+j];
- }
- /* apply mute window for samples above nts-ii */
- for (j = 0; j < nts-ii+T*shift-smooth; j++) {
- G_d[l*nxs*nts+i*nts+j] = 0.0;
- }
- for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) {
- G_d[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
- }
- //for (j = 0; j < nts-ii+T*shift; j++) {
- // G_d[l*nxs*nts+i*nts+j] = 0.0;
- //}
- }
- for (i = 0; i < npos; i++) {
- ix = ixpos[i];
- j = 0;
- f1min[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+j];
- for (j = 1; j < nts; j++) {
- f1min[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j];
- }
- /* apply mute window for samples above nts-ii */
- //for (j = 0; j < nts-ii+T*shift; j++) {
- // f1min[l*nxs*nts+i*nts+j] = 0.0;
- //}
- }
- }
- }
- else { /* use f1min from previous iteration as starting point and do niterec iterations */
- niterrun=niterec;
- recur=1;
- if (verbose>2) vmess("Doing %d iterations using previous result at time-sample %d",niterrun,ii);
- for (l = 0; l < Nfoc; l++) {
- for (i = 0; i < npos; i++) {
- j=0;
- ix = ixpos[i];
- G_d[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts] - f1min[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- G_d[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j] - f1min[l*nxs*nts+i*nts+nts-j];
- }
- /* apply mute window for samples above nts-ii */
- for (j = 0; j < nts-ii+T*shift-smooth; j++) {
- G_d[l*nxs*nts+i*nts+j] = 0.0;
- }
- for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) {
- G_d[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
- }
- //for (j = 0; j < nts-ii+T*shift; j++) {
- // G_d[l*nxs*nts+i*nts+j] = 0.0;
- //}
- }
- }
- }
-/*================ initialization ================*/
-
- memcpy(Ni, G_d, Nfoc*nxs*ntfft*sizeof(float));
-
-/*================ number of Marchenko iterations ================*/
-
- for (iter=0; iter<niterrun; iter++) {
-
- t2 = wallclock_time();
-
-/*================ construction of Ni(-t) = - \int R(x,t) Ni(t) ================*/
-
- synthesis(Refl, Fop, Ni, iRN, nx, nt, nxs, nts, dt, xsyn, Nfoc,
- xrcv, xsrc, xnx, fxse, fxsb, dxs, dxsrc, dx, ntfft, nw, nw_low, nw_high, mode,
- reci, nshots, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv, ixmask, verbose);
-
- if (file_iter != NULL) {
- writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 1, 1000*ii+iter);
- }
-
- t3 = wallclock_time();
- tsyn += t3 - t2;
-
- /* N_k(x,t) = -N_(k-1)(x,-t) */
- for (l = 0; l < Nfoc; l++) {
- for (i = 0; i < npos; i++) {
- j = 0;
- ix = ixpos[i];
- Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+ix*nts+j];
- for (j = 1; j < nts; j++) {
- Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+ix*nts+nts-j];
- }
- }
- }
-
- if (iter % 2 == 0) { /* even iterations: => f_1^+(t) */
- /* apply muting for the acausal part */
- for (l = 0; l < Nfoc; l++) {
- for (i = 0; i < npos; i++) {
- /* apply mute window for samples after ii */
- for (j = ii-T*shift; j < nts; j++) {
- Ni[l*nxs*nts+i*nts+j] = 0.0;
- }
- for (j = ii-T*shift+smooth, k=0; j < ii-T*shift; j++, k++) {
- Ni[l*nxs*nts+i*nts+j] *= costaper[k];
- }
- /* apply mute window for delta function at t=0*/
- for (j = 0; j < shift-smooth; j++) {
- Ni[l*nxs*nts+i*nts+j] = 0.0;
- }
- for (j = shift-smooth, k=1; j < shift; j++, k++) {
- Ni[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
- }
- f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
- }
- }
- }
- if (file_iter != NULL) {
- writeDataIter("f1plus.su", f1plus, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter);
- }
- }
- else {/* odd iterations: => f_1^-(t) */
- for (l = 0; l < Nfoc; l++) {
- for (i = 0; i < npos; i++) {
- ix = ixpos[i];
- if (recur==1) { /* use f1min from previous iteration */
- for (j = 0; j < nts; j++) {
- Ni[l*nxs*nts+i*nts+j] += DD[l*nxs*nts+ix*nts+j];
- }
- j = 0;
- f1min[l*nxs*nts+i*nts+j] = -Ni[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- f1min[l*nxs*nts+i*nts+j] = -Ni[l*nxs*nts+i*nts+nts-j];
- }
- }
- else {
- j = 0;
- f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j];
- }
- }
- /* apply mute window for delta function at t=0*/
- for (j = nts-shift+smooth; j < nts; j++) {
- Ni[l*nxs*nts+i*nts+j] = 0.0;
- }
- for (j = nts-shift, k=0; j < nts-shift+smooth; j++, k++) {
- Ni[l*nxs*nts+i*nts+j] *= costaper[k];
- }
- /* apply mute window for samples above nts-ii */
- /* important to apply this mute after updating f1min */
- //for (j = 0; j < nts-ii+T*shift; j++) {
- // Ni[l*nxs*nts+i*nts+j] = 0.0;
- //}
- /* apply mute window for samples above nts-ii */
- for (j = 0; j < nts-ii+T*shift-smooth; j++) {
- Ni[l*nxs*nts+i*nts+j] = 0.0;
- }
- for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) {
- Ni[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
- }
- }
- }
- if (file_iter != NULL) {
- writeDataIter("f1min.su", f1min, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter);
- }
- } /* end else (iter) branch */
-
- t2 = wallclock_time();
- tcopy += t2 - t3;
- if (verbose>2) vmess("*** Iteration %d finished ***", iter);
-
- } /* end of iterations */
-
- for (l = 0; l < Nfoc; l++) {
- for (i = 0; i < npos; i++) {
- RR[l*nxs*nts+i*nts+ii] = f1min[l*nxs*nts+i*nts+ii];
- }
- }
-
- /* To Do? optional write intermediate RR results to file */
-
- if (verbose) {
- if(!((iend-ii-istart)%perc)) fprintf(stderr,"\b\b\b\b%3d%%",(ii-istart)*100/(iend-istart));
- if((ii-istart)==10)t4=wallclock_time();
- if((ii-istart)==20){
- t4=(wallclock_time()-t4)*((iend-istart)/10.0);
- fprintf(stderr,"\r %s: Estimated total compute time = %.2fs.\n %s: Progress: %3d%%",xargv[0],t4,xargv[0],(ii-istart)/((iend-istart)/100));
- }
- //t4=wallclock_time();
- tii=(t4-t1)*((float)(iend-istart)/(ii-istart+1.0))-(t4-t1);
- //vmess("Remaining compute time at time-sample %d = %.2f s.",ii, tii);
- }
-
- } /* end of time iterations ii */
-
- free(Ni);
- free(G_d);
- free(f1min);
- free(f1plus);
-
- t2 = wallclock_time();
- if (verbose) {
- fprintf(stderr,"\b\b\b\b%3d%%\n",100);
- vmess("Total CPU-time marchenko = %.3f", t2-t0);
- vmess("with CPU-time synthesis = %.3f", tsyn);
- vmess("with CPU-time copy array = %.3f", tcopy);
- vmess(" CPU-time fft data = %.3f", tfft);
- vmess("and CPU-time read data = %.3f", tread);
- }
-
-/*================ write output files ================*/
-
- fp_rr = fopen(file_rr, "w+");
- if (fp_rr==NULL) verr("error on creating output file %s", file_rr);
-
- tracf = 1;
- for (l = 0; l < Nfoc; l++) {
- if (reci) f2 = fxsb;
- else f2 = fxf;
-
- for (i = 0; i < n2; i++) {
- hdrs_out[i].fldr = l+1;
- hdrs_out[i].sx = NINT(xsyn[l]*1000);
- hdrs_out[i].offset = (long)NINT((f2+i*d2) - xsyn[l]);
- hdrs_out[i].tracf = tracf++;
- hdrs_out[i].selev = NINT(zsyn[l]*1000);
- hdrs_out[i].sdepth = NINT(-zsyn[l]*1000);
- hdrs_out[i].f1 = f1;
- }
- ret = writeData(fp_rr, (float *)&RR[l*size], hdrs_out, n1, n2);
- if (ret < 0 ) verr("error on writing output file.");
- }
- ret = fclose(fp_rr);
- if (ret < 0) verr("err %d on closing output file",ret);
-
- if (verbose) {
- t1 = wallclock_time();
- vmess("and CPU-time write data = %.3f", t1-t2);
- }
-
-/*================ free memory ================*/
-
- free(hdrs_out);
-
- exit(0);
-}
-
-
diff --git a/marchenko/synthesis.c b/marchenko/synthesis.c
index f149dbadfec328da28c94ac602c836bfc2c9f75f..828da0749d5d9158ad5f14314e999a6437e4c5e6 100644
--- a/marchenko/synthesis.c
+++ b/marchenko/synthesis.c
@@ -35,11 +35,6 @@ typedef struct _complexStruct { /* complex number */
double wallclock_time(void);
-void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int
-Nfoc, float *xrcv, float *xsrc, int *xnx, float fxse, float fxsb, float dxs, float dxsrc, float dx, int ntfft, int
-nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int npos, double *tfft, int *isxcount, int
-*reci_xsrc, int *reci_xrcv, float *ixmask, int verbose);
-
void synthesisPositions(int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nfoc, float *xrcv, float *xsrc, int *xnx,
float fxse, float fxsb, float dxs, float dxsrc, float dx, int nshots, int *ixpos, int *npos, int *isxcount, int countmin, int reci, int verbose);
@@ -49,9 +44,9 @@ int linearsearch(int *array, size_t N, int value);
/* Refl has the full acquisition grid R(x_r, x_s)
* Fop has the acquisition grid of the operator, ideally this should be equal to the acquisition grid of Refl,
* so all traces can be used to compute R*Fop.
- * The output iRN has the traces in the grid of Fop, these are the x_s positions of R(x_r,x_s) */
+ * The output RNi has the traces in the grid of Fop, these are the x_s positions of R(x_r,x_s) */
-void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int
+void synthesis(complex *Refl, complex *Fop, float *Top, float *RNi, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int
Nfoc, float *xrcv, float *xsrc, int *xnx, float fxse, float fxsb, float dxs, float dxsrc, float dx, int ntfft, int
nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int npos, double *tfft, int *isxcount, int
*reci_xsrc, int *reci_xrcv, float *ixmask, int verbose)
@@ -90,7 +85,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
t0 = wallclock_time();
/* reset output data to zero */
- memset(&iRN[0], 0, Nfoc*nxs*nts*sizeof(float));
+ memset(&RNi[0], 0, Nfoc*nxs*nts*sizeof(float));
ctrace = (complex *)calloc(ntfft,sizeof(complex));
/* this first check is done to support an acquisition geometry that has more receiver than source
@@ -142,7 +137,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
/*================ SYNTHESIS ================*/
#pragma omp parallel default(none) \
- shared(iRN, dx, npe, nw, verbose, nshots, xnx) \
+ shared(RNi, dx, npe, nw, verbose, nshots, xnx) \
shared(Refl, Nfoc, reci, xsrc, xsyn, fxsb, fxse, nxs, dxs) \
shared(nx, dxsrc, nfreq, nw_low, nw_high, fxb, fxe) \
shared(Fop, size, nts, ntfft, scl, ixrcv) \
@@ -177,7 +172,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
/* place result at source position ixsrc; dx = receiver distance */
for (j = 0; j < nts; j++)
- iRN[l*size+ixsrc*nts+j] += rtrace[j]*scl*dx;
+ RNi[l*size+ixsrc*nts+j] += rtrace[j]*scl*dx;
} /* end of Nfoc loop */
@@ -196,7 +191,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
if (reci != 0) {
#pragma omp parallel default(none) \
- shared(iRN, dx, nw, verbose) \
+ shared(RNi, dx, nw, verbose) \
shared(Refl, Nfoc, reci, xsrc, xsyn, fxsb, fxse, nxs, dxs) \
shared(nx, dxsrc, nfreq, nw_low, nw_high, fxb, fxe) \
shared(reci_xrcv, reci_xsrc, ixmask, isxcount) \
@@ -233,7 +228,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
/* place result at source position ixsrc; dxsrc = shot distance */
for (j = 0; j < nts; j++)
- iRN[l*size+ixsrc*nts+j] = ixmask[ixsrc]*(iRN[l*size+ixsrc*nts+j]+rtrace[j]*scl*dxsrc);
+ RNi[l*size+ixsrc*nts+j] = ixmask[ixsrc]*(RNi[l*size+ixsrc*nts+j]+rtrace[j]*scl*dxsrc);
} /* end of Nfoc loop */
diff --git a/utils/freqwave.c b/utils/freqwave.c
index 591537747f614f86337fb3c25144ad32cbf1b318..5879c329abfc1470e0f870fb393674c273a912ee 100644
--- a/utils/freqwave.c
+++ b/utils/freqwave.c
@@ -43,9 +43,9 @@ void freqwave(float *wave, int nt, float dt, float fp, float fmin, float flef, f
optn = optncr(nt);
nfreq = 1+(optn/2);
df = 1.0/(dt*optn);
- iof = MAX(NINT(fmax/df), NINT(fp/df));
att = pow(10.0, db/20.0);
-
+ if (fp < 0) iof = NINT(fmax/df);
+ else iof = NINT(fp/df);
if (iof > nfreq) verr("characterizing frequency aliased");
cwave = (complex *)malloc(nfreq*sizeof(complex));