diff --git a/Make_include_template b/Make_include_template
index 09dbb7fa49ee24a9e81c9d057696fcda6595c193..50d0665da5ee49792d976dfd52f09e7c17f73204 100644
--- a/Make_include_template
+++ b/Make_include_template
@@ -27,11 +27,11 @@ OPTC = -O3 -ffast-math
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
-# Linux gcc version 3.x
-#OPTC = -O3 -ffast-math -funroll-all-loops -mfpmath=sse
+#OPTC = -O3 -fno-trapping-math -ffast-math -funroll-all-loops -fomit-frame-pointer -ftree-vectorizer-verbose=1
+# Linux gcc very old version 3.x
+#OPTC = -O3 -ffast-math -funroll-all-loops
-#for double precision use FFTlib and emmod
+#enable double precision only for use in FFTlib and emmod
#OPTC += -DDOUBLE
#OPTF += -fdefault-double-8 -fdefault-real-8
diff --git a/Makefile b/Makefile
index 1351be294e29137f4d3ecbd503ac264307f3279f..86415c08131202fcd87cf41f2ef7f23211bd89d6 100644
--- a/Makefile
+++ b/Makefile
@@ -2,6 +2,7 @@
all: mkdirs
cd FFTlib ; $(MAKE)
+ cd zfp ; $(MAKE) install
cd fdelmodc ; $(MAKE) install
cd fdelmodc3D ; $(MAKE) install
cd utils ; $(MAKE) install
@@ -10,7 +11,6 @@ all: mkdirs
cd corrvir ; $(MAKE) install
cd raytime ; $(MAKE) install
cd MDD ; $(MAKE) install
- cd zfp ; $(MAKE)
cd fdacrtmc ; $(MAKE) install
@@ -21,6 +21,7 @@ mkdirs:
clean:
cd FFTlib ; $(MAKE) $@
+ cd zfp ; $(MAKE) $@
cd fdelmodc ; $(MAKE) $@
cd fdelmodc3D ; $(MAKE) $@
cd utils ; $(MAKE) $@
@@ -29,11 +30,11 @@ clean:
cd corrvir ; $(MAKE) $@
cd raytime ; $(MAKE) $@
cd MDD ; $(MAKE) $@
- cd zfp ; $(MAKE) $@
cd fdacrtmc ; $(MAKE) $@
realclean:
cd FFTlib ; $(MAKE) $@
+ cd zfp ; $(MAKE) $@
cd fdelmodc ; $(MAKE) $@
cd fdelmodc3D ; $(MAKE) $@
cd fdelrtmc ; $(MAKE) $@
@@ -43,7 +44,6 @@ realclean:
cd corrvir ; $(MAKE) $@
cd raytime ; $(MAKE) $@
cd MDD ; $(MAKE) $@
- cd zfp ; $(MAKE) $@
cd fdacrtmc ; $(MAKE) $@
rm -f lib/*
rm -f include/*
diff --git a/README b/README.md
similarity index 90%
rename from README
rename to README.md
index 3c02538e179d943feb3c34633eeccc89210fd476..e31ebf72ba99e49263e9ff91331cf8454f532f6d 100644
--- a/README
+++ b/README.md
@@ -1,19 +1,25 @@
-%%%%%%%
+
+ACKNOWLEDGEMENT
+---------
+This work received funding from the European Research Council (grant 742703) and the NWO Domain Applied and Engineering Sciences (grant 13939).
+
+LICENSE
+---------
THE ACCOMPANYING PROGRAM IS PROVIDED UNDER THE TERMS OF THIS COMMON PUBLIC LICENSE ("AGREEMENT"). ANY USE, REPRODUCTION OR DISTRIBUTION OF THE PROGRAM CONSTITUTES RECIPIENT'S ACCEPTANCE OF THIS AGREEMENT.
A copy of this license can be found in the file 'Common_Public_License.txt' in the directory where you have found this README.
http://www.opensource.org/licenses/cpl1.0.php
+SU
+--
Some routines are from Seismic Unix and include the SU LEGAL_STATEMENT in the source code.
-%%%%%%%
Copyright (c) 2017 by the Society of Exploration Geophysicists.
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.
-
REFERENCES
---------
-0- DOI reference of this software release
@@ -36,7 +42,7 @@ Download: https://janth.home.xs4all.nl/Publications/Articles/ThorbeckeGPY2017.pd
-3- When you are using the marchenko_primaries algorithm developed by Lele Zhang please refer to the following paper
Free-surface and internal multiple elimination in one step without adaptive subtraction
-Lele Zhang1 and Evert Slob2
+Lele Zhang and Evert Slob
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
@@ -49,37 +55,40 @@ Download: http://homepage.tudelft.nl/t4n4v/4_Journals/Geophys.Prosp/GP_19a.pdf
-5-
-
INSTALLATION
-------------
1) To compile and link the code you first have to set the ROOT variable in the Make_include file which can be found in the directory where you have found this README.
-You can use Make_include_template as a first start: cp Make_include_template Make_include
+You can use Make_include_template as a first start:
+
+```
+cp Make_include_template Make_include
+```
2) Check the compiler and CFLAGS options in the file Make_include and adapt to the system you are using. The default options are set for a the GNU C-compiler on a Linux system. A Fortran or C++ compiler are not needed to compile the code. The Makefile has been tested with GNU make.
3) If the compiler options are set in the Make_include file you can type
-
-> make
-
+```
+make clean
+make
+```
and the Makefile will make:
- FFT library
-- fdelmodc
-- marchenko
+- fdelmodc / 3D
+- marchenko / 3D
- utils
The libraries will be placed in the lib/ directory and the executables in the bin/ directory.
To use the executables don't forget to include the pathname in your PATH:
+```
bash/sh:
export PATH='path_to_this_directory'/bin:$PATH:
csh:
setenv PATH 'path_to_this_directory'/bin:$PATH:
-
-
-> make clean
+```
Finite Difference Modeling: FDELMODC
------------------------------------
@@ -89,19 +98,21 @@ If the compilation has finished without errors and produced an executable called
in the directory fdelmodc/demo/
-The demo directory contains scripts which demonstrate the different possibilities of the modeling program.
+The demo directory contains many scripts which demonstrate the different possibilities of the modeling program.
To reproduce the Figures shown in the GEOPHYICS manuscript "Finite-difference modeling experiments for seismic interferometry" the scripts in FiguresPaper directory can be used. Please read the README in the FiguresPaper directory for more instructions and guidelines.
Marchenko method : MARCHENKO
----------------------------
-If the compilation has finished without errors and produced an executable called bin/marchenko you can run one of the demo programs by running a set of scripts that are explained in a README in one of the directories marchenko/demo/*
+If the compilation has finished without errors and produced an executable called bin/marchenko you can run one of the demo programs by running a set of scripts that are explained in a README in one of the directories marchenko/demo/oneD or demo/twoD
To reproduce the Figures shown in the GEOPHYICS paper "Implementation of the Marchenko method" the scripts in marchenko/demo/oneD directory can be used. The README in this directory gives more instructions and guidelines.
To reproduce the Figures shown in the Scientific Reports paper "Virtual acoustics in inhomogeneous media with single-sided access" the scripts in marchenko/demo/ScientificReports directory can be used. The README in this directory gives more instructions and guidelines.
+To reproduce the Figures shown in the GEOPHYICS paper "Implementation of the Marchenko Multiple Elimination algorithm" the scripts in marchenko/demo/oneD directory can be used. The README_PRIMARIES in this directory gives more instructions and guidelines.
+
MDD
---
The MDD kernels depend on BLAS and LAPACK calls. Free downloads of these libraries can be found on
diff --git a/fdelmodc/Makefile b/fdelmodc/Makefile
index f214df9f50540585f5d0f3aa2949d3731d1441b8..37f88801d8c6e93aee11b2f6825a91da971293a4 100644
--- a/fdelmodc/Makefile
+++ b/fdelmodc/Makefile
@@ -8,6 +8,7 @@ ALLINC = -I.
#LIBS += -L$L -lgenfft -lm $(LIBSM)
#OPTC = -g -Wall -fsignaling-nans -O0 -fopenmp
#OPTC += -fopenmp -Waddress
+#OPTC += -qopt-report=5 -g
#OPTC := $(subst -O3 -ffast-math, -O1 -g ,$(OPTC))
#PGI options for compiler feedback
#OPTC += -Mprof=lines
diff --git a/fdelmodc/applySource.c b/fdelmodc/applySource.c
index 3aeb972c1404aaa91a8e9b3edf17685edea4e40a..3b9e0377fc3b9052a65bfa0a58479bc647812e91 100644
--- a/fdelmodc/applySource.c
+++ b/fdelmodc/applySource.c
@@ -89,7 +89,10 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
ix = src.x[isrc] + ibndx;
iz = src.z[isrc] + ibndz;
}
- else { /* plane wave and point sources */
+ else if (src.plane) {/* plane wave sources */
+ ix = ixsrc + ibndx + src.x[isrc];
+ iz = izsrc + ibndz + src.z[isrc];
+ else { /* point sources */
ix = ixsrc + ibndx + is0 + isrc;
iz = izsrc + ibndz;
}
@@ -156,6 +159,10 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
/* stable implementation changes amplitude and more work is needed */
//vz[ix*n1+iz] = 0.5*(vz[ix*n1+iz-1]+vz[ix*n1+iz+1])+src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
//vz[ix*n1+iz] = 0.25*(vz[ix*n1+iz-2]+vz[ix*n1+iz-1]+vz[ix*n1+iz]+vz[ix*n1+iz+1])+src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
+ }
+ else if (src.type == 10) {
+ tzz[ix*n1+iz-1] -= src_ampl;
+ tzz[ix*n1+iz+1] += src_ampl;
} /* src.type */
diff --git a/fdelmodc/demo/modelOilGas.scr b/fdelmodc/demo/modelOilGas.scr
index 5fa26904388a3d67e07cc72a6bc7861d1b5bec16..6227c3739cffebeb4815b8ba8315ac5dd01980d2 100755
--- a/fdelmodc/demo/modelOilGas.scr
+++ b/fdelmodc/demo/modelOilGas.scr
@@ -16,19 +16,19 @@ makemod sizex=5000 sizez=2500 dx=$dx dz=$dx cp0=$cp ro0=$rho \
makewave w=g1 fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
-export OMP_NUM_THREADS=1
+export OMP_NUM_THREADS=4
zsrc=1100
zsrc=0
which fdelmodc
-fdelmodc \
+../fdelmodc \
file_cp=syncl_cp.su ischeme=1 iorder=4 \
file_den=syncl_ro.su \
file_src=wave.su \
file_rcv=shot_fd.su \
- src_type=7 \
+ src_type=1 \
src_orient=1 \
src_injectionrate=0 \
rec_type_vz=1 \
@@ -37,13 +37,13 @@ fdelmodc \
dtrcv=0.004 \
rec_delay=0.1 \
verbose=2 \
- tmod=2.01 \
+ tmod=0.01 \
dxrcv=10.0 \
xrcv1=-2250 xrcv2=2250 \
zrcv1=0 zrcv2=0 \
xsrc=0 zsrc=$zsrc \
file_snap=snapF_$zsrc \
- tsnap1=0.1 tsnap2=2.0 dtsnap=0.05 dxsnap=$dx dzsnap=$dx \
+ tsnap1=0.1 tsnap2=0.0 dtsnap=0.05 dxsnap=$dx dzsnap=$dx \
ntaper=101 \
snapwithbnd=1 \
left=2 right=2 top=2 bottom=2
diff --git a/fdelmodc/demo/modelOilGas.scr.ok b/fdelmodc/demo/modelOilGas.scr.ok
deleted file mode 100755
index 89907c0d5c8891a821c3b06d27cdbd1f25b94fdf..0000000000000000000000000000000000000000
--- a/fdelmodc/demo/modelOilGas.scr.ok
+++ /dev/null
@@ -1,95 +0,0 @@
-#!/bin/bash
-
-#export PATH=$HOME/bin:$HOME/src/OpenSource/utils:$PATH:
-
-cp=2000
-rho=2500
-dx=2.5
-dt=0.0005
-
-
-makemod sizex=5000 sizez=2500 dx=$dx dz=$dx cp0=$cp ro0=$rho \
- orig=-2500,0 file_base=syncl.su \
- intt=def x=-2500,0,2500 z=250,250,250 poly=0 cp=2300 ro=2000 \
- intt=def x=-2500,-2000,-1000,-800,0,800,2500 z=650,650,700,750,900,750,600 poly=2 cp=2600 ro=2500 \
- intt=def x=-2500,0,2500 z=1390,1390,1390 poly=0 cp=2000 ro=2000
-
-makewave w=g1 fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
-
-export OMP_NUM_THREADS=8
-
-zsrc=1100
-zsrc=0
-
-../fdelmodc \
- file_cp=syncl_cp.su ischeme=1 iorder=4 \
- file_den=syncl_ro.su \
- file_src=wave.su \
- file_rcv=shot_fd.su \
- src_type=7 \
- src_orient=1 \
- src_injectionrate=0 \
- rec_type_vz=1 \
- rec_type_p=1 \
- rec_int_vz=2 \
- dtrcv=0.004 \
- rec_delay=0.1 \
- verbose=2 \
- tmod=2.01 \
- dxrcv=10.0 \
- xrcv1=-2250 xrcv2=2250 \
- zrcv1=0 zrcv2=0 \
- xsrc=0 zsrc=$zsrc \
- file_snap=snapF_$zsrc \
- tsnap1=0.1 tsnap2=2.0 dtsnap=0.05 dxsnap=$dx dzsnap=$dx \
- ntaper=101 \
- snapwithbnd=1 \
- left=2 right=2 top=2 bottom=2
-
-suxmovie < snapF_${zsrc}_svz.su loop=1 clip=1e-13
-
-exit
-
-makemod sizex=5000 sizez=2500 dx=$dx dz=$dx cp0=$cp ro0=$rho \
- orig=-2500,0 file_base=hom.su
-
-fdelmodc \
- file_cp=hom_cp.su ischeme=1 iorder=4 \
- file_den=hom_ro.su \
- file_src=wave.su \
- file_rcv=shot_hom_fd.su \
- src_type=7 \
- src_orient=1 \
- src_injectionrate=1 \
- rec_type_vz=1 \
- rec_type_p=1 \
- rec_int_vz=2 \
- dtrcv=0.004 \
- rec_delay=0.1 \
- verbose=2 \
- tmod=4.195 \
- dxrcv=10.0 \
- xrcv1=-2250 xrcv2=2250 \
- zrcv1=0 zrcv2=0 \
- xsrc=0 zsrc=0 \
- ntaper=250 \
- left=4 right=4 top=1 bottom=4
-
-cp=2000
-rho=1000
-dx=10
-
-makemod sizex=5000 sizez=2500 dx=$dx dz=2.5 cp0=$cp ro0=$rho \
- orig=-2500,0 file_base=syncl_migr.su \
- intt=def x=-2500,0,2500 z=250,250,250 poly=0 cp=2300 ro=5000 \
- intt=def x=-2500,-2000,-1000,-800,0,800,2500 z=650,650,700,750,900,750,600 poly=2 cp=2600 ro=1000 \
- intt=def x=-2500,0,2500 z=1390,1390,1390 poly=0 cp=2000 ro=5000
-
-sudiff shot_fd_rvz.su shot_hom_fd_rvz.su > diff_rvz.su
-
-makewave fp=20 dt=0.004 file_out=wavedt.su nt=1024 t0=0.0
-
-migr file_shot=diff_rvz.su file_src=wavedt.su file_vel=syncl_migr_cp.su nshots=1 \
- file_image=migr0.su verbose=3 imc=0
-
-
diff --git a/fdelmodc/demo/testFreeSurface.scr b/fdelmodc/demo/testFreeSurface.scr
index d4591319a6c5cf0453ab6937488bc7787e3e9cad..b824c9d8939ee0fe951da29a75b86e0e9acc50d1 100755
--- a/fdelmodc/demo/testFreeSurface.scr
+++ b/fdelmodc/demo/testFreeSurface.scr
@@ -12,17 +12,19 @@ cd /Users/jan/src/OpenSource/fdelmodc/demo
cp=2000
rho=1000
-dx=4
-dt=0.0010
+dx=1
+dt=0.00010
makemod sizex=3000 sizez=1000 dx=$dx dz=$dx cp0=$cp ro0=$rho cs0=1500 \
orig=-1500,0 file_base=freesurf.su \
intt=def x=-1500,0,1500 z=20,20,20 poly=0 cp=2300 ro=2000 cs=1600 \
intt=def x=-1500,0,1500 z=50,50,50 poly=0 cp=2100 ro=1400 cs=1300
+export OMP_NUM_THREADS=1
+
makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
-../fdelmodc \
+../../bin/fdelmodc \
file_cp=freesurf_cp.su ischeme=3 \
file_cs=freesurf_cs.su \
file_den=freesurf_ro.su \
diff --git a/fdelmodc/fdelmodc.c b/fdelmodc/fdelmodc.c
index 878f6b61edb35ca53052e38db12355786b7a8e06..6f6242a26899474661f8ce1ddb8249b0bea4b364 100644
--- a/fdelmodc/fdelmodc.c
+++ b/fdelmodc/fdelmodc.c
@@ -127,13 +127,13 @@ char *sdoc[] = {
" verbose=0 ......... silent mode; =1: display info",
" ",
" SHOT AND GENERAL SOURCE DEFINITION:",
-" src_type=1 ........ 1=P 2=Txz 3=Tzz 4=Txx 5=S-pot 6=Fx 7=Fz 8=P-pot 9=double-couple",
+" src_type=1 ........ 1=P 2=Txz 3=Tzz 4=Txx 5=S-pot 6=Fx 7=Fz 8=P-pot 9=double-couple 10=Fz by P",
" src_orient=1 ...... orientation of the source",
" - 1=monopole",
" - 2=dipole +/- vertical oriented",
" - 3=dipole - + horizontal oriented",
-//" - 4=dipole +/0/-",
-//" - 5=dipole + -",
+" - 4=dipole +/0/-",
+" - 5=dipole + -",
" dip=0 ............. dip for double-couple source",
" strike=0 .......... strike for double-couple source",
" xsrc=middle ....... x-position of (first) shot ",
@@ -749,6 +749,23 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, izsrc, it, src_nwav, verbose)
free(p);
free(q);
}
+ if (bnd.ntap) {
+ free(bnd.tapx);
+ free(bnd.tapz);
+ free(bnd.tapxz);
+ }
+ free(bnd.surface);
+ free(shot.x);
+ free(shot.z);
+ free(src.x);
+ free(src.z);
+ free(src.tbeg);
+ free(src.tend);
+ free(rec.x);
+ free(rec.z);
+ free(rec.xr);
+ free(rec.zr);
+ if(wav.nsamp!=NULL) free(wav.nsamp);
#ifdef MPI
MPI_Finalize();
diff --git a/fdelmodc/getParameters.c b/fdelmodc/getParameters.c
index d882ba97c05a56c85c4314bc10ba87480ddcdb25..a316a1f53635357c73247199b83486bd47b12810 100644
--- a/fdelmodc/getParameters.c
+++ b/fdelmodc/getParameters.c
@@ -173,7 +173,6 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
/* check if receiver delays is defined; option inactive: add delay time to total modeling time */
if (!getparfloat("rec_delay",&rdelay)) rdelay=0.0;
rec->delay=NINT(rdelay/mod->dt);
-// mod->tmod += rdelay;
mod->nt = NINT(mod->tmod/mod->dt);
dt = mod->dt;
@@ -925,6 +924,9 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
wav->nst = nsamp; /* put total number of samples in nst part */
}
}
+ if (src->type==7) { /* set also src_injectionrate=1 */
+ src->injectionrate=1;
+ }
if (src->multiwav) {
if (wav->nx != nsrc) {
@@ -956,6 +958,8 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
case 6 : fprintf(stderr,"Fx "); break;
case 7 : fprintf(stderr,"Fz "); break;
case 8 : fprintf(stderr,"P-potential"); break;
+ case 9 : fprintf(stderr,"double-couple"); break;
+ case 10 : fprintf(stderr,"Fz on P grid with +/-"); break;
}
fprintf(stderr,"\n");
if (wav->random) vmess("Wavelet has a random signature with fmax=%.2f", wav->fmax);
@@ -1110,7 +1114,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
rec->skipdt=NINT(dtrcv/dt);
dtrcv = mod->dt*rec->skipdt;
if (!getparfloat("rec_delay",&rdelay)) rdelay=0.0;
- if (!getparint("rec_ntsam",&rec->nt)) rec->nt=NINT((mod->tmod-rdelay)/dtrcv)+1;
+ if (!getparint("rec_ntsam",&rec->nt)) rec->nt=(int)round((mod->tmod-rdelay+0.01*mod->dt)/dtrcv)+1;
if (!getparint("rec_int_p",&rec->int_p)) rec->int_p=0;
if (!getparint("rec_int_vx",&rec->int_vx)) rec->int_vx=0;
if (!getparint("rec_int_vz",&rec->int_vz)) rec->int_vz=0;
@@ -1118,7 +1122,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
if (!getparint("scale",&rec->scale)) rec->scale=0;
if (!getparfloat("dxspread",&dxspread)) dxspread=0;
if (!getparfloat("dzspread",&dzspread)) dzspread=0;
- rec->nt=MIN(rec->nt, NINT((mod->tmod-rdelay)/dtrcv)+1);
+ rec->nt=MIN(rec->nt, NINT((mod->tmod-rdelay+0.01*mod->dt)/dtrcv)+1);
/* allocation of receiver arrays is done in recvPar */
/*
diff --git a/marchenko/demo/oneD/README_PRIMARIES b/marchenko/demo/oneD/README_PRIMARIES
index 5eda7197ef806033b2f1286d28a433d9d53b2c64..55ecc68005042fc0d3ff7d4bcd66dc62dacfc690 100644
--- a/marchenko/demo/oneD/README_PRIMARIES
+++ b/marchenko/demo/oneD/README_PRIMARIES
@@ -2,28 +2,30 @@ Description of files:
1) model.scr computes the model and the 'basis' shot of R => shot5_rp.su
-2) p5all.scr create from basis shot full Reflection response matrix => shotsdx5_rp.su (3.3 GB)
-3) primaries.scr computes the internal multiple attenuated (middle) shot record.
-4) itertions.scr computes the intermediate results of the multiple attenutation scheme and produces all output files that are used in the manuscript.
-5) epsPrimaries.scr selected output form 4) are converted to .eps pictures that are used in the Figures to explain the method.
- reproduce the postscript files of the manuscript using SU postscript plotting programs.
-
-extra scripts
+ - runtime on 4 cores is 4-5 minutes and produces a 3.3 GB data file with 901 shots
+2) itertions.scr computes the intermediate results of the multiple attenutation scheme and produces all output files that are used in the manuscript.
+ - runtime on 4 cores is
+3) epsPrimaries.scr selected output from step 2) are converted to .eps pictures that are used in the Figures to explain the method.
+ To reproduce the postscript files of the manuscript SU postscript plotting programs are required.
+-3) epsModel.scr to generate the postscript files for the numerical model
+
+optional scripts not needed to reproduce the figures:
++) primaries.scr computes the internal multiple attenuated (middle) shot record.
+ - runtime on 4 cores is ~500 s.
+) primariesPlane.scr: computes the internal moval scheme for plane-waves (see Meles 2020)
+) clean: remove all produced files and start with a clean directory
+--------------------------
To reproduce the Figures in the Manuscript:
-
--------------------------
* Figure 2: Model + Initial wavefield
-==> run model.scr to generate the data .su files: this will take 3-4 minutes. The files generate are:
+==> run model.scr to generate the data .su files: this will take 4-5 minutes. The files generate are:
- hom_cp.su, hom_ro.su
- model10_cp.su, model10_ro.su
- shot5_fd_rp.su
- shot5_hom_fd_rp.su
- shot5_rp.su
- - wave.su
- wavefw.su
==> run './epsPrimaries.scr Figure2' to generate the postscript files of Figure 2
@@ -32,24 +34,23 @@ model_cp_line.eps => Figure 2a
model_ro_line.eps => Figure 2b
shotx0_rp.eps => Figure 2c
+It also produces two extra pictures of the wavelet used in the FD modelling:
+wavefw_freq.eps
+wavefw.eps
+
--------------------------
* Figure 3: First Iteration
-The full R matrix is build up from the the shot record computed with model.scr
-
-==> run p5all.scr to generate the full R matrix for a fixed spread geometry. This will take less than one minute. The file generated is
- - shotsdx5_rp.su this file has a size of 3.3 GB
-
-This R is the only input of the Marchenko Primaries algorithm.
-
-==> run './iterations.scr ./iterations.scr Figure3467' to compute the intermediate results of the first iterations of the Marchenko Primaries algorithm.
-This will take 10 seconds. The generated files are:
+==> run './iterations.scr Figure34910' to compute the intermediate results of the first iterations of the Marchenko Primaries algorithm.
+This will take 15 seconds. The generated files are:
- M0_276000.su
- - iter_2760##.su (not used)
- Mi_2760##.su
- - f1min_2760##.su
+ - k1min_2760##.su
+ - v1plus_2760##.su
+ - iter_2760##.su (not used)
- pred_rr_276.su (not used)
+ - DDshot_450.su (not used): selected shot record convolved with file_src=wave.su
where ## ranges from 01 to 33
@@ -66,117 +67,123 @@ fconvN0flip.eps => Figure 3d
Mi_276001.eps => Figure 3e
+--------------------------
+* Figure 4 second iteration
To generate the postscript files for Figure 4:
==> run './epsPrimaries.scr Figure4'
+This will produce the following files:
+
+fconvN1fulltime.eps => Figure 4c
+fconvN1flip.eps => Figure 4d
+Mi_276002.eps => Figure 4e
+
+The window time function in Figure 5 is not reproduced.
+
+--------------------------
+* Figure 6 v1plus and convergence
+==> run './iterations.scr Figure6' to compute the marchenko results with ii=200
+
To generate the postscript files for Figure 6:
==> run './epsPrimaries.scr Figure6'
-To generate the postscript files for Figure 7:
-==> run './epsPrimaries.scr Figure7'
-
-
-
-shotx0_rp.eps => Figure 5 R == Figure 3c
-f1min_001.eps => Figure 5 f^-_1,0
-iter_002.eps => Figure 5 -N_1
-f1plus_002.eps => Figure 5 f^+_1,0
-
--- Figure 6 column 1
-iter_001.eps
-iter_002.eps
-iter_003.eps
-iter_004.eps
--- Figure 6 column 2
-f1min_001.eps
-f1min_002.eps
-f1min_003.eps
-f1min_004.eps
--- Figure 6 column 3
-p0plus_flip.eps
-f1plus_002.eps
-f1plus_003.eps
-f1plus_004.eps
--- Figure 6 column 4
-pgreen_001.eps
-pgreen_002.eps
-pgreen_003.eps
-pgreen_004.eps
--- Figure 6 column 5
-compare_001.eps
-compare_002.eps
-compare_003.eps
-compare_004.eps
-
-
-Note that the script epsIterwithLabels.scr produces the same figures, but with axis-labels.
+This will produce the following files:
+v1plus_200001.eps => Figure 6a
+v1plus_max.eps => Figure 6b
+k1min_200030.eps => Figure 6b
--------------------------
-* Figure 7: Comparison of Marchenko result with reference
+* Figure 8 To compute the convergence for a strong contrast medium:
+cd strongContrast
+==> run ./model.scr
+==> run ./iterations.scr Figure8
-To compute the marchenko results for 8 iterations.
+To generate the postscript files for Figure 8:
+==> run './epsPrimaries.scr Figure8'
-==> run marchenko.scr This will take less than 1 minute. The generated files are:
- - pgreen.su, pgreen512.su
- - diffref.su
- - Gplus0.su
- - Gmin0.su
- - f1plus0.su
- - f1min0.su
- - f2.su
+This will produce the following files:
+v1plusStrong_max.eps => Figure 8
+Don't forget to go back to the main directory with the results
+cd ../
-At the end of the run the script will display in X11 a comparison of the middle trace.
-To make the postscript figure
+--------------------------
+* Figure 9 iterations M_i
-==> run epsCompare.scr
+To generate the postscript files for Figure 9:
+==> run './epsPrimaries.scr Figure9'
+
+This will produce the following files:
+Mi_276002.eps => Figure 9b
+Mi_276004.eps => Figure 9d
+Mi_276012.eps => Figure 9f
+Mi_276020.eps => Figure 9h
+Mi_276001.eps => Figure 9a
+Mi_276003.eps => Figure 9c
+Mi_276011.eps => Figure 9e
+Mi_276019.eps => Figure 9g
-mergeGreenRef.eps => Figure 7
--------------------------
-* Figure 8: snapshots of back propagating f2 in actual medium
-
-To compute the snapshots
-
-==> run backpropf2.scr This will take about 1 minute. The generated output file is
- - backpropf2_sp.su
-
-The postscript files of Figure 8 are generated with
-
-==> run epsBackprop.scr
-
--- Figure 8 column 1
-backpropf2_-0.30.eps
-backpropf2_-0.15.eps
-backpropf2_-0.03.eps
-backpropf2_-0.02.eps
-backpropf2_0.00.eps
--- Figure 8 column 2
-backpropf2_0.30.eps
-backpropf2_0.15.eps
-backpropf2_0.03.eps
-backpropf2_0.02.eps
-backpropf2_0.00.eps
--- Figure 8 column 3
-backpropf2sum_0.30.eps
-backpropf2sum_0.15.eps
-backpropf2sum_0.03.eps
-backpropf2sum_0.02.eps
-backpropf2_0.00.eps
-
-
-The figures in the appendix, to explain the different options in the programs, are reproduced by
-
-==> run figAppendi.scr
-
--- Figure A-1
-noise_above0.eps
-noise_above1.eps
-noise_above-1.eps
-noise_above2.eps
-noise_above4.eps
-
--- Figure A-2
-iniFocus_shifts.eps
+* Figure 10 iterations M_i and k_1^-
+
+To generate the postscript files for Figure 10:
+==> run './epsPrimaries.scr Figure10'
+
+This will produce the following files:
+Mi_276002flip.eps => Figure 10a
+k1min_276002.eps => Figure 10b
+Mi_276004flip.eps => Figure 10c
+k1min_276004.eps => Figure 10d
+Mi_276010flip.eps => Figure 10e
+k1min_276010.eps => Figure 10f
+Mi_276020flip.eps => Figure 10g
+k1min_276020.eps => Figure 10h
+
+
+--------------------------
+* Figure 11 iterations k_1^- for different ii 246:316:10
+
+To generate the data
+==> run ./iterations.scr Figure11
+this will take ~2 minutes and generate a lot of files
+
+To generate the postscript files for Figure 11:
+==> run './epsPrimaries.scr Figure11'
+
+This will produce the following files:
+k1min_246032.eps => Figure 11a
+k1min_256032.eps => Figure 11b
+k1min_266032.eps => Figure 11c
+k1min_276032.eps => Figure 11d
+k1min_286032.eps => Figure 11e
+k1min_296032.eps => Figure 11f
+k1min_306032.eps => Figure 11g
+k1min_316032.eps => Figure 11h
+
+
+--------------------------
+* Figure 13 iterations M_i and k_1^- for ii-276 T-MME scheme
+
+To generate the data
+==> run ./iterations.scr Figure13
+this will take ~15 seconds
+****NOTE this will overwrite the results of the MME-scheme !
+
+To generate the postscript files for Figure 13:
+==> run './epsPrimaries.scr Figure13'
+
+This will produce the following files:
+Mi_276002T.eps => Figure 13a
+Mi_276004T.eps => Figure 13b
+Mi_276001T.eps => Figure 13c
+Mi_276003T.eps => Figure 13d
+k1min_276002T.eps => Figure 13e
+k1min_276004T.eps => Figure 13f
+k1min_276010T.eps => Figure 13g
+k1min_276020T.eps => Figure 13h
+
+
+
diff --git a/marchenko/demo/oneD/epsPrimaries.scr b/marchenko/demo/oneD/epsPrimaries.scr
index 00949d8fe8d62007ba0bae6ec5dda71b5bf1ef44..9b444ea14762d86382557b214b73fa24fc3382d7 100755
--- a/marchenko/demo/oneD/epsPrimaries.scr
+++ b/marchenko/demo/oneD/epsPrimaries.scr
@@ -19,7 +19,7 @@ ns=`surange < Mi_276014.su | grep ns | awk '{print $2}'`
nsd=400 #number of samples to display
(( nstart = ns - nsd ))
-file=f1min_${istart}002.su
+file=k1min_${istart}002.su
sumax < $file mode=abs outpar=nep
clipf1=`cat nep | awk '{print $1/11}'`
clipm1=`cat nep | awk '{print $1/22}'`
@@ -30,9 +30,9 @@ file_base=${file%.su}
sumax < $file mode=abs outpar=nep
clipref=`cat nep | awk '{print $1/15}'`
-echo "clipiter="$clipiter "clipref="$clipref "clipf1="$clipf1
-clipiter=$clipref
-clipf1=$clipref
+echo "clipiter="$clipiter "clipref="$clipref "clipf1="$clipf1 "clipm1="$clipm1
+#clipiter=$clipref
+#clipf1=$clipref
# Initialisation and First iteration
if [[ "$1" == "Figure3" ]];
@@ -55,23 +55,27 @@ select=451
#original shot record from Reflection matrix
suwind key=fldr min=$select max=$select < shotsdx5_rp.su > shotsx0.su
#first iteration
-cp M0_276000.su N0.su
+suwind itmin=516 itmax=1539 < M0_276000.su > N0.su
#compute R*N0
fconv file_in1=shotsx0.su file_in2=N0.su file_out=fconvN0.su verbose=1 fmax=90
+file=fconvN0.su
+sumax < $file mode=abs outpar=nep
+clipiter=`cat nep | awk '{print $1/15}'`
+
file=fconvN0.su
file_base=${file%.su}
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
n1tic=2 d2=5 x1beg=0 d1=1 f1=0 f1num=0 d1num=200 \
label1="time sample number" label2="lateral distance" \
- f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}fulltime.eps
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipiter > ${file_base}fulltime.eps
file=fconvN0.su
file_base=${file%.su}
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
- n1tic=2 d2=5 x1beg=$nstart x1end=$ns d1=1 f1=0 f1num=700 d1num=100 \
+ n1tic=2 d2=5 x1beg=624 x1end=1024 d1=1 f1=0 f1num=700 d1num=100 \
label1="time sample number" label2="lateral distance" \
- f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}zoom.eps
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipiter > ${file_base}zoom.eps
suflip < fconvN0.su flip=3 | sugain scale=-1 > fconvN0flip.su
file=fconvN0flip.su
@@ -79,7 +83,7 @@ file_base=${file%.su}
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 f1=0 f1num=0 d1num=100 \
label1="time sample number" label2="lateral distance" \
- f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipiter > $file_base.eps
iter=1
piter=$(printf %03d $iter)
@@ -88,7 +92,7 @@ file_base=${file%.su}
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \
- f2=-2250 f2num=-2000 d2num=1000 clip=$clipiter > $file_base.eps
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
fi
@@ -100,26 +104,76 @@ suwind itmax=1023 < Mi_276001.su > N0.su
#compute R*N0
fconv file_in1=shotsx0.su file_in2=N0.su file_out=fconvN1.su verbose=1 fmax=90
+file=fconvN1.su
+sumax < $file mode=abs outpar=nep
+clipiter=`cat nep | awk '{print $1/15}'`
suflip < fconvN1.su flip=3 | sugain scale=-1 > fconvN1flip.su
+
file=fconvN1flip.su
file_base=${file%.su}
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
n1tic=2 d2=5 x1beg=624 x1end=1024 d1=1 f1=0 f1num=600 d1num=100 \
label1="time sample number" label2="lateral distance" \
- f2=-2250 f2num=-2000 d2num=1000 clip=$clipm1 > $file_base.eps
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipiter > $file_base.eps
file=fconvN1.su
file_base=${file%.su}
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
n1tic=2 d2=5 x1beg=0 d1=1 f1=0 f1num=0 d1num=200 \
label1="time sample number" label2="lateral distance" \
- f2=-2250 f2num=-2000 d2num=1000 clip=$clipm1 > ${file_base}fulltime.eps
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipiter > ${file_base}fulltime.eps
+
+iter=2
+piter=$(printf %03d $iter)
+file=Mi_${istart}${piter}.su
+file_base=${file%.su}
+suwind itmin=516 itmax=1539 < $file > N2.su
+supsimage < N2.su hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+ n1tic=2 d2=5 x1beg=624 x1end=1024 d1=1 f1=0 d1num=100 \
+ label1="time sample number" label2="lateral distance" \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+fi
+
+if [[ "$1" == "Figure6" ]];
+then
+#convergence plot
+rm vplus.su
+istart=200
+for (( iter=1; iter<=29; iter+=2 ))
+do
+piter=$(printf %03d $iter)
+echo $piter
+file=v1plus_${istart}${piter}.su
+file_base=${file%.su}
+suwind < $file key=offset min=0 max=0 >> vplus.su
+done
+
+basop file_in=vplus.su file_out=vplus_env.su choice=4
+sumax < vplus_env.su verbose=1 mode=max
+supsmax < vplus_env.su n=15 verbose=1 mode=max \
+ style=normal linewidth=2.0 f1=1 labelsize=14 label1="iteration number" label2="amplitude" \
+ d1num=2 d2num=0.2 wbox=6 hbox=3 grid1=dot grid2=dot n1tic=2 n2tic=2 x2end=2.9 x2beg=1.8 > v1plus_max.eps
+
+file=v1plus_${istart}001.su
+file_base=${file%.su}
+supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+ n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
+ label1="time sample number" label2="lateral distance" \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > $file_base.eps
+
+file=k1min_${istart}030.su
+file_base=${file%.su}
+supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 < $file\
+ n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
+ label1="time sample number" label2="lateral distance" \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > $file_base.eps
fi
#Iterations
-if [[ "$1" == "Figure5" ]];
+if [[ "$1" == "Figure9" ]];
then
-for (( iter=0; iter<=21; iter+=2 ))
+for iter in 2 4 12 20
do
piter=$(printf %03d $iter)
echo $piter
@@ -129,13 +183,13 @@ file_base=${file%.su}
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
n1tic=2 d2=5 x1beg=$nstart x1end=$ns d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \
- f2=-2250 f2num=-2000 d2num=1000 clip=$clipiter > $file_base.eps
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
done
fi
-if [[ "$1" == "Figure5" ]];
+if [[ "$1" == "Figure9" ]];
then
-for (( iter=1; iter<=21; iter+=2 ))
+for iter in 1 3 11 19
do
piter=$(printf %03d $iter)
echo $piter
@@ -144,13 +198,14 @@ file_base=${file%.su}
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \
- f2=-2250 f2num=-2000 d2num=1000 clip=$clipiter > $file_base.eps
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
done
fi
-exit;
+if [[ "$1" == "Figure10" ]];
+then
#iterations
-for (( iter=0; iter<=31; iter+=2 ))
+for iter in 2 4 10 20
do
piter=$(printf %03d $iter)
echo $piter
@@ -166,7 +221,7 @@ suflip flip=3 < $file | supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > ${file_base}flip.eps
-file=f1min_${istart}${piter}.su
+file=k1min_${istart}${piter}.su
file_base=${file%.su}
supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 < $file\
n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
@@ -174,19 +229,63 @@ supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 < $file\
f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > $file_base.eps
done
+fi
+if [[ "$1" == "Figure11" ]];
+then
iter=32
piter=$(printf %03d $iter)
echo $piter
for (( istart=246; istart<=316; istart+=10 ))
do
-file=f1min_${istart}${piter}.su
+file=k1min_${istart}${piter}.su
file_base=${file%.su}
supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 < $file\
n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > $file_base.eps
done
+fi
+
+if [[ "$1" == "Figure13" ]];
+then
+istart=276
+for iter in 2 4
+do
+piter=$(printf %03d $iter)
+echo $piter
+file=Mi_${istart}${piter}.su
+file_base=${file%.su}
+#possibly add suflip flip=3 to flip the time-axis
+supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+ n1tic=2 d2=5 x1beg=$nstart x1end=$ns d1=1 d1num=100 \
+ label1="time sample number" label2="lateral distance" \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}T.eps
+done
+for iter in 1 3
+do
+piter=$(printf %03d $iter)
+echo $piter
+file=Mi_${istart}${piter}.su
+file_base=${file%.su}
+supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+ n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
+ label1="time sample number" label2="lateral distance" \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}T.eps
+done
+for iter in 2 4 10 20
+do
+piter=$(printf %03d $iter)
+echo $piter
+file=k1min_${istart}${piter}.su
+file_base=${file%.su}
+supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 < $file\
+ n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
+ label1="time sample number" label2="lateral distance" \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > ${file_base}T.eps
+done
+
+fi
exit;
diff --git a/marchenko/demo/oneD/iterations.scr b/marchenko/demo/oneD/iterations.scr
index 3cc1ec7981379ceab761f5aed388890c6c2246fa..671376345f4a18dd8218288e2170095d6695ae54 100755
--- a/marchenko/demo/oneD/iterations.scr
+++ b/marchenko/demo/oneD/iterations.scr
@@ -3,41 +3,57 @@
#second reflector at time:
# 800/1800+600/2300
# .70531400966183574878 => sample 176
-#third reflector at time:
+#third reflector at time model.scr:
# 800/1800+600/2300+800/2000
# 1.10531400966183574878 sample 276
+#third reflector at time modepm.scr
+#800/1800+600/3200+520/2000
+#.96531400966183574878 sample 241
+
+R=shotsdx5_rp.su
+makewave fp=20 dt=0.004 file_out=wave.su nt=2048 t0=0.0 scale=0 scfft=1
select=451
-if [[ "$1" == "Figure3467" ]];
+if [[ "$1" == "Figure34910" ]];
then
for istart in 276
do
(( iend = istart + 1 ))
-../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
+../../marchenko_primaries file_shot=$R ishot=$select file_src=wave.su \
nshots=901 verbose=2 istart=$istart iend=$iend fmax=90 \
pad=512 niter=33 smooth=8 niterskip=600 shift=16 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su
done
fi
-if [[ "$1" == "Figure8" ]];
+if [[ "$1" == "Figure6" ]];
+then
+istart=200
+(( iend = istart + 1 ))
+../../marchenko_primaries file_shot=$R ishot=$select file_src=wave.su \
+ nshots=901 verbose=2 istart=$istart iend=$iend fmax=90 \
+ pad=512 niter=30 smooth=8 niterskip=600 shift=16 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su
+fi
+
+if [[ "$1" == "Figure11" ]];
then
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 \
+../../marchenko_primaries file_shot=$R ishot=$select file_src=wave.su \
nshots=901 verbose=2 istart=$istart iend=$iend fmax=90 \
pad=512 niter=33 smooth=8 niterskip=600 shift=16 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su
done
fi
-if [[ "$1" == "Figure10" ]];
+if [[ "$1" == "Figure13" ]];
then
istart=276
(( iend = istart + 1 ))
-../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
+../../marchenko_primaries file_shot=$R ishot=$select file_src=wave.su \
nshots=901 verbose=2 istart=$istart iend=$iend fmax=90 \
pad=512 niter=33 smooth=8 niterskip=600 shift=16 file_rr=pred_rr_${istart}.su T=1 file_iter=iterT.su
fi
+
diff --git a/marchenko/demo/oneD/model.scr b/marchenko/demo/oneD/model.scr
index 9165e8758862d4cef0b497ca585bea46336e9d13..b97d41800605d7af489dac030e90265174f5b17b 100755
--- a/marchenko/demo/oneD/model.scr
+++ b/marchenko/demo/oneD/model.scr
@@ -1,4 +1,9 @@
#!/bin/bash
+#SBATCH -J model_1.5D
+#SBATCH --cpus-per-task=8
+#SBATCH --ntasks=1
+#SBATCH --time=0:15:00
+
#adjust this PATH to where the code is installed
export PATH=$HOME/src/OpenSource/bin:$PATH:
@@ -19,7 +24,7 @@ makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096
#define wavelet for reference and intial focusing field.
makewave fp=25 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
-export OMP_NUM_THREADS=4
+export OMP_NUM_THREADS=8
#Model shot record in middle of model
fdelmodc \
@@ -74,4 +79,28 @@ fdelmodc \
#subtract direct wave from full model shot record: this defines R
sudiff shot5_fd_rp.su shot5_hom_fd_rp.su > shot5_rp.su
+#########################################################
+# Generate the full R matrix for a fixed spread geometry.
+
+dxshot=5000 # with scalco factor of 1000
+ishot=0
+nshots=901
+
+rm shotsdx5_rp.su
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -2250000 + ${ishot}*${dxshot} ))
+ (( tr1 = $nshots - ${ishot} ))
+ (( tr2 = ${tr1} + $nshots - 1 ))
+ echo xsrc=$xsrc tr1=$tr1 tr2=$tr2
+
+ (( ishot = $ishot + 1))
+
+ suwind < shot5_rp.su key=tracl min=$tr1 max=$tr2 | \
+ sushw key=sx,gx,fldr,trwf \
+ a=$xsrc,-2250000,$ishot,$nshots b=0,$dxshot,0,0 j=0,$nshots,0,0 | \
+ suchw key1=offset key2=gx key3=sx c=-1 d=1000 >> shotsdx5_rp.su
+done
diff --git a/marchenko/demo/oneD/primaries.scr b/marchenko/demo/oneD/primaries.scr
index b5398efee445d97e3c0a7ca5e7a326c9b3fb01c5..94b37e1ca0081efa4b5b1581d24b6e444478d768 100755
--- a/marchenko/demo/oneD/primaries.scr
+++ b/marchenko/demo/oneD/primaries.scr
@@ -1,28 +1,32 @@
#!/bin/bash -x
#SBATCH -J marchenko_primaries
-#SBATCH --cpus-per-task=10
+#SBATCH --cpus-per-task=20
#SBATCH --ntasks=1
-#SBATCH --time=2:40:00
+##SBATCH --time=1:29:00
+#SBATCH -p max2h
+# Generate the model and Reflection data
+#modelpm.scr
#cd $SLURM_SUBMIT_DIR
export PATH=$HOME/src/OpenSource/bin:$PATH:
-export OMP_NUM_THREADS=10
+export OMP_NUM_THREADS=20
#shot record to remove internal multiples
+R=shotsdx5_rp.su
select=451
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 \
+../../marchenko_primaries file_shot=$R 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=50 shift=20 file_rr=pred_rr.su T=0 file_update=update.su
+ niter=31 shift=20 smooth=10 niterec=2 niterskip=50 file_rr=pred_rr.su T=0 file_update=update.su
exit;
#for reference original shot record from Reflection matrix
-suwind key=fldr min=$select max=$select < shotsdx5_rp.su itmax=2047 > shotsx0.su
+suwind key=fldr min=$select max=$select < $R itmax=2047 > shotsx0.su
fconv file_in1=shotsx0.su file_in2=wave.su file_out=shotw.su
# for displaying results
diff --git a/marchenko/demo/oneD/primariesPlane.scr b/marchenko/demo/oneD/primariesPlane.scr
index 510bb9ab88871334f05dc8c4d39106be42ad73ed..94fe87958e2e20a687b2b4c6c6b2af0e50a74f85 100755
--- a/marchenko/demo/oneD/primariesPlane.scr
+++ b/marchenko/demo/oneD/primariesPlane.scr
@@ -2,7 +2,7 @@
#SBATCH -J marchenko_primaries
#SBATCH --cpus-per-task=10
#SBATCH --ntasks=1
-#SBATCH --time=2:40:00
+#SBATCH --time=0:30:00
#cd $SLURM_SUBMIT_DIR
@@ -10,31 +10,9 @@
export PATH=$HOME/src/OpenSource/bin:$PATH:
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=30 \
-# nshots=901 verbose=2 istart=276 iend=277 fmax=90 pad=1024 \
-# 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=30.0 \
- nshots=901 verbose=2 istart=40 iend=1000 fmax=90 pad=1024 \
- niter=20 smooth=10 niterskip=60 niterec=2 shift=20 file_rr=plane2_10_rr.su T=0
-
-exit
-
-#for reference original shot record from Reflection matrix
-suwind key=fldr min=$select max=$select < shotsdx5_rp.su itmax=2047 > shotsx0.su
-fconv file_in1=shotsx0.su file_in2=wave.su file_out=shotw.su
-
-# for displaying results
-
-(suwind key=offset min=0 max=0 < pred_rr.su ; suwind key=offset min=0 max=0 < shotw.su) | suxgraph &
-
-sudiff shotw.su pred_rr.su > diff.su
-suximage < shotw.su x1end=2.5 clip=1 title="original shot"&
-suximage < pred_rr.su x1end=2.5 clip=1 title="shot with multiples removed"&
-suximage < diff.su x1end=2.5 clip=1 title="removed multiples"&
+../../marchenko_primaries1 file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=10.0 \
+ nshots=901 verbose=4 istart=40 iend=500 fmax=90 pad=1024 \
+ niter=30 smooth=10 niterskip=60 niterec=2 shift=20 file_rr=plane1_10_rr.su T=0
diff --git a/marchenko/demo/oneD/strongContrast/epsPrimaries.scr b/marchenko/demo/oneD/strongContrast/epsPrimaries.scr
new file mode 100755
index 0000000000000000000000000000000000000000..ee99f2234d6d916097aef29b45628e1eb8f233c3
--- /dev/null
+++ b/marchenko/demo/oneD/strongContrast/epsPrimaries.scr
@@ -0,0 +1,25 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin/:$PATH:
+
+if [[ "$1" == "Figure8" ]];
+then
+rm vplus.su
+istart=200
+for (( iter=1; iter<=29; iter+=2 ))
+do
+piter=$(printf %03d $iter)
+echo $piter
+file=v1plus_${istart}${piter}.su
+file_base=${file%.su}
+suwind < $file key=offset min=0 max=0 >> vplus.su
+done
+
+basop file_in=vplus.su file_out=vplus_env.su choice=4
+sumax < vplus_env.su verbose=1 mode=max
+supsmax < vplus_env.su n=15 verbose=1 mode=max \
+ style=normal linewidth=2.0 f1=1 labelsize=14 label1="iteration number" label2="amplitude" \
+ d1num=2 d2num=0.5 wbox=6 hbox=3 grid1=dot grid2=dot n1tic=2 n2tic=2 x2end=6.9 x2beg=1.8 > v1plusStrong_max.eps
+
+fi
+
diff --git a/marchenko/demo/oneD/strongContrast/iterations.scr b/marchenko/demo/oneD/strongContrast/iterations.scr
new file mode 100755
index 0000000000000000000000000000000000000000..6a71918049bde156839a1f247db1a5ce10e1a63a
--- /dev/null
+++ b/marchenko/demo/oneD/strongContrast/iterations.scr
@@ -0,0 +1,25 @@
+#!/bin/bash
+
+#second reflector at time:
+# 800/1800+600/2300
+# .70531400966183574878 => sample 176
+#third reflector at time model.scr:
+# 800/1800+600/2300+800/2000
+# 1.10531400966183574878 sample 276
+#third reflector at time modepm.scr
+#800/1800+600/3200+520/2000
+#.96531400966183574878 sample 241
+
+makewave fp=20 dt=0.004 file_out=wave.su nt=2048 t0=0.0 scale=0 scfft=1
+R=shotsdx5_rp.su
+
+select=451
+
+if [[ "$1" == "Figure8" ]];
+then
+istart=200
+(( iend = istart + 1 ))
+../../../marchenko_primaries file_shot=$R ishot=$select file_src=wave.su \
+ nshots=901 verbose=2 istart=$istart iend=$iend fmax=90 \
+ pad=512 niter=30 smooth=8 niterskip=600 shift=16 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su
+fi
diff --git a/marchenko/demo/oneD/strongContrast/model.scr b/marchenko/demo/oneD/strongContrast/model.scr
new file mode 100755
index 0000000000000000000000000000000000000000..bd651affb43e888d6fcaf5eea912f624ed124238
--- /dev/null
+++ b/marchenko/demo/oneD/strongContrast/model.scr
@@ -0,0 +1,106 @@
+#!/bin/bash
+#SBATCH -J marchenko_primaries
+#SBATCH --cpus-per-task=10
+#SBATCH --ntasks=1
+#SBATCH --time=1:59:00
+#SBATCH -p max2h
+
+#adjust this PATH to where the code is installed
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+dx=2.5
+dt=0.0005
+
+#define gridded model for FD computations
+makemod sizex=10000 sizez=1400 dx=$dx dz=$dx cp0=1800 ro0=500 \
+ orig=-5000,0 file_base=model10.su verbose=2 \
+ intt=def x=-5000,5000 z=400,400 poly=0 cp=2300 ro=5000 \
+ intt=def x=-5000,5000 z=700,700 poly=0 cp=2000 ro=500 \
+ intt=def x=-5000,5000 z=1100,1100 poly=0 cp=2500 ro=5000
+
+#define wavelet for modeling R
+makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3 scale=0 scfft=1
+
+#define wavelet for reference and intial focusing field.
+makewave fp=25 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+export OMP_NUM_THREADS=10
+
+#Model shot record in middle of model
+fdelmodc \
+ file_cp=model10_cp.su ischeme=1 iorder=4 \
+ file_den=model10_ro.su \
+ file_src=wavefw.su \
+ file_rcv=shot5_fd.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.3 \
+ verbose=2 \
+ tmod=4.392 \
+ dxrcv=5.0 \
+ xrcv1=-4500 xrcv2=4500 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=0 \
+ npml=101 \
+ left=2 right=2 top=2 bottom=2
+
+#define homogenoeus model to compute direct wave only
+makemod sizex=10000 sizez=1200 dx=$dx dz=$dx cp0=1800 ro0=1000 \
+ orig=-5000,0 file_base=hom.su verbose=2
+
+#Model direct wave only in middle of model
+fdelmodc \
+ file_cp=hom_cp.su ischeme=1 iorder=4 \
+ file_den=hom_ro.su \
+ file_src=wavefw.su \
+ file_rcv=shot5_hom_fd.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.3 \
+ verbose=2 \
+ tmod=4.392 \
+ dxrcv=5.0 \
+ xrcv1=-4500 xrcv2=4500 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=0 \
+ npml=101 \
+ left=2 right=2 top=2 bottom=2
+
+#subtract direct wave from full model shot record: this defines R
+sudiff shot5_fd_rp.su shot5_hom_fd_rp.su > shot5_rp.su
+
+
+# Generate the full R matrix for a fixed spread geometry.
+
+dxshot=5000 # with scalco factor of 1000
+ishot=0
+nshots=901
+
+rm shotsdx5_rp.su
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -2250000 + ${ishot}*${dxshot} ))
+ (( tr1 = $nshots - ${ishot} ))
+ (( tr2 = ${tr1} + $nshots - 1 ))
+ echo xsrc=$xsrc tr1=$tr1 tr2=$tr2
+
+ (( ishot = $ishot + 1))
+
+ suwind < shot5_rp.su key=tracl min=$tr1 max=$tr2 | \
+ sushw key=sx,gx,fldr,trwf \
+ a=$xsrc,-2250000,$ishot,$nshots b=0,$dxshot,0,0 j=0,$nshots,0,0 | \
+ suchw key1=offset key2=gx key3=sx c=-1 d=1000 >> shotsdx5_rp.su
+
+done
diff --git a/marchenko/marchenko_primaries.c b/marchenko/marchenko_primaries.c
index d893a330715027cf0e5069297341b693a8a82bc0..6c572a4111c8bfe0bbe6a0c88afcfd4febd2fcb3 100644
--- a/marchenko/marchenko_primaries.c
+++ b/marchenko/marchenko_primaries.c
@@ -18,6 +18,7 @@ void omp_set_num_threads(int num_threads);
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#define ISODD(n) ((n) & 01)
#ifndef COMPLEX
typedef struct _complexStruct { /* complex number */
@@ -37,9 +38,7 @@ 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 *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
+void synthesisp(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);
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);
@@ -72,7 +71,6 @@ char *sdoc[] = {
" src_angle=0 .............. angle with horizontal of plane source array",
" src_velo=1500 ............ velocity to use in src_angle definition",
" t0=0.1 ................... time shift in plane-wave source wavelet for migration",
-//" xorig=nrecv/2 ............ center of plane-wave",
" MARCHENKO ITERATIONS ",
" niter=22 ................. number of iterations to initialize and restart",
" niterec=2 ................ number of iterations in recursive part of the time-samples",
@@ -92,16 +90,17 @@ char *sdoc[] = {
//" countmin=0 ............... 0.3*nxrcv; minimum number of reciprocal traces for a contribution",
" OUTPUT DEFINITION ",
" file_rr= ................. output file with primary only shot record",
+" file_dd= ................. output file with input of the algorithm ",
" file_iter= ............... output file with -Mi(-t) for each iteration: writes",
" ............... M0.su=M0 : initialisation of algorithm",
" ............... RMi: iterative terms ",
-" ............... f1min.su: f1min terms ",
+" ............... k1min.su: k1min terms ",
" file_update= ............. output file with updates only => removed internal multiples",
" T=0 ...................... :1 compute transmission-losses compensated primaries ",
" verbose=0 ................ silent option; >0 displays info",
" ",
" ",
-" author : Lele Zhang & Jan Thorbecke : 2019 ",
+" author : Lele Zhang & Jan Thorbecke : 2020 ",
" ",
NULL};
/**************** end self doc ***********************************/
@@ -109,27 +108,28 @@ NULL};
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, tr;
+ size_t nread, size;
+ int i, j, k, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath, nacq;
+ int 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 reci, countmin, mode, n2out, verbose, ntfft;
int iter, niter, niterec, recur, niterskip, niterrun, tracf, *muteW;
int hw, ii, iw, ishot, istart, iend;
- int smooth, *ixpos, npos, ix, m, pad, T, isms, isme, perc;
- int nshots_r, *isxcount, *reci_xsrc, *reci_xrcv, shift, plane_wave, xorig;
+ int smooth, *ixpos, *ixp, npos, ix, ixrcv, m, pad, T, isms, isme, perc;
+ int nshots_r, *isxcount, *reci_xsrc, *reci_xrcv, shift, plane_wave;
float fmin, fmax, tom, deltom, *tapersh, *tapersy, fxf, dxf, *xsrc, *xrcv, *zsyn, *zsrc, *xrcvsyn;
- double t0, t1, t2, t3, t4, tsyn, tread, tfft, tcopy, tii;
+ double t0, t1, t2, t3, t4, t5, ttime, 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, scltap;
- float *rtrace, *tmpdata, *f1min, *f1plus, *RMi, *Mi, *trace;
+ float *rtrace, *tmpdata, *k1min, *v1plus, *RMi, *Mi, *trace;
float *Mup, *Msp, *maxval;
float xmin, xmax, scale, tsq;
float Q, f0, *ixmask, *costaper;
float src_velo, src_angle, grad2rad, p, *twplane;
complex *Refl, *Fop, *ctrace, *cwave, csum, cwav;
char *file_tinv, *file_shot, *file_rr, *file_src, *file_iter, *file_update;
+ char *file_dd;
segy *hdrs_out, hdr;
initargs(argc, argv);
@@ -142,6 +142,7 @@ int main (int argc, char **argv)
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_dd", &file_rr)) file_dd = NULL;
if (!getparstring("file_iter", &file_iter)) file_iter = NULL;
if (!getparstring("file_update", &file_update)) file_update = NULL;
@@ -169,8 +170,9 @@ int main (int argc, char **argv)
if(!getparint("ishot", &ishot)) ishot=300;
if(!getparint("plane_wave", &plane_wave)) plane_wave=0;
if(!getparfloat("src_angle", &src_angle)) src_angle = 0.0;
- if (!getparfloat("src_velo",&src_velo)) src_velo=1500.;
- if (!getparfloat("t0",&tt0)) tt0=0.1;
+ if(!getparfloat("src_velo",&src_velo)) src_velo=1500.;
+ if(!getparfloat("t0",&tt0)) tt0=0.1;
+ if( (niterskip>1) && ISODD(niter) ) niter++;
if (T>0) {
T=-1;
@@ -188,7 +190,7 @@ int main (int argc, char **argv)
smooth = MIN(smooth, shift);
if (smooth) {
costaper = (float *)malloc(smooth*sizeof(float));
- scl = M_PI/((float)smooth);
+ scl = M_PI/((float)smooth+1);
for (i=0; i<smooth; i++) {
costaper[i] = 0.5*(1.0+cos((i+1)*scl));
}
@@ -200,6 +202,7 @@ int main (int argc, char **argv)
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;
+ nacq = n2;
nxs = n2;
nts = n1;
dxs = d2;
@@ -220,27 +223,28 @@ int main (int argc, char **argv)
if(!getparint("ishot", &ishot)) ishot=1+(nshots-1)/2;
ishot -= 1; /* shot numbering starts at 0 */
Nfoc = 1;
- nxs = nx;
+ nacq = NINT((xmax-xmin)/dx)+1;
+ nxs = nshots;
nts = nt;
dxs = dx;
- fxsb = fx;
+ fxsb = xmin;
+ fxse = xmax;
}
- assert (nxs >= nshots); /* ToDo allow other geometries */
- //if(!getparint("xorig", &xorig)) xorig=-(nxs-1)/2;
+ assert (nacq >= nshots); /* ToDo allow other geometries */
/* compute time delay for plane-wave responses */
- twplane = (float *) calloc(nxs,sizeof(float)); /* initialize with zeros */
+ twplane = (float *) calloc(nacq,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]);
+ for (i=0; i<nacq; i++) {
+ twplane[i] = fabsf((nacq-i-1)*dxs*p)+tt0;
+ if (verbose >=3) vmess("plane-wave delay-time i=%d x=%f t=%f", i, dxs*(nacq-i-1), twplane[i]);
}
}
else {
- for (i=0; i<nxs; i++) {
+ for (i=0; i<nacq; i++) {
twplane[i] = (i)*dxs*p+tt0;
if (verbose >=3) vmess("plane-wave delay-time i=%d x=%f t=%f", i, dxs*(i), twplane[i]);
}
@@ -260,19 +264,19 @@ int main (int argc, char **argv)
/*================ 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));
- RMi = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ Fop = (complex *)calloc(nacq*nw*Nfoc,sizeof(complex));
+ RMi = (float *)calloc(Nfoc*nacq*ntfft,sizeof(float));
+ DD = (float *)calloc(Nfoc*nacq*ntfft,sizeof(float));
Mi = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
- M0 = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
- DD = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ M0 = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ k1min = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ v1plus = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
SRC = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
RR = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
Mup = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
Msp = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
trace = (float *)malloc(ntfft*sizeof(float));
- ixpos = (int *)malloc(nxs*sizeof(int));
+ ixpos = (int *)malloc(nacq*sizeof(int));
energyM0= (double *)malloc(Nfoc*sizeof(double));
xrcvsyn = (float *)calloc(Nfoc*nxs,sizeof(float));
xsyn = (float *)malloc(Nfoc*sizeof(float));
@@ -363,8 +367,8 @@ int main (int argc, char **argv)
/*================ 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);
+ readShotData(file_shot, xrcv, xsrc, zsrc, xnx, Refl, nw, nw_low, nshots, 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));
@@ -402,19 +406,28 @@ int main (int argc, char **argv)
dxsrc = dx;
}
+/*================ set focusing postions in R ================*/
+
+ synthesisPositions(nx, nt, nacq, 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);
+ }
+
/*================ Defining focusing operator(s) from R ================*/
/* M0 = -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 M0 from Refl of %s", file_shot);
- nts = ntfft;
-
- scl = 1.0/((float)2.0*ntfft);
+ 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++) {
+ ixrcv = NINT((xrcv[ishot*nx+i]-fxsb)/dxs);
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;;
@@ -422,13 +435,8 @@ int main (int argc, char **argv)
/* 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];
+ DD[0*nxs*nts+ixrcv*nts+j] = 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];
- //}
}
/* set timereversal for searching First Break */
/* for experimenting with non flat truncation windows */
@@ -439,14 +447,15 @@ int main (int argc, char **argv)
//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 */
+ /* construct plane wave (time reversed and multiplied with -1) from all shot records */
if (plane_wave==1) {
for (l=0; l<nshots; l++) {
+ ix = ixpos[l];
memset(ctrace, 0, sizeof(complex)*(nfreq+1));
for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
- tom = j*deltom*twplane[l];
csum.r=0.0; csum.i=0.0;
for (i = 0; i < xnx[l]; i++) {
+ // ToDo for general acquisitons use ix at receiver position */
tom = j*deltom*twplane[i];
csum.r += Refl[l*nw*nx+m*nx+i].r*cos(-tom) - Refl[l*nw*nx+m*nx+i].i*sin(-tom);
csum.i += Refl[l*nw*nx+m*nx+i].i*cos(-tom) + Refl[l*nw*nx+m*nx+i].r*sin(-tom);
@@ -458,7 +467,7 @@ int main (int argc, char **argv)
/* transfrom result back to time domain */
cr1fft(ctrace, rtrace, ntfft, 1);
for (j = 0; j < nts; j++) {
- DD[0*nxs*nts+l*nts+j] = -1.0*scl*rtrace[j];
+ DD[0*nxs*nts+ix*nts+j] = 1.0*scl*rtrace[j];
}
/* compute Source wavelet for plane-wave imaging */
for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
@@ -507,13 +516,16 @@ int main (int argc, char **argv)
if (verbose) vmess("Taper for operator applied ntap=%d", ntap);
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < nxs; i++) {
+ ix = ixpos[i];
for (j = 0; j < nts; j++) {
- DD[l*nxs*nts+i*nts+j] *= tapersy[i];
+ DD[l*nxs*nts+ix*nts+j] *= tapersy[i];
}
}
}
free(tapersy);
}
+ ixp = (int *)malloc(nxs*sizeof(int));
+ for (i=0; i<nxs; i++) ixp[i] = i;
/*================ Check the size of the files ================*/
@@ -525,9 +537,10 @@ int main (int argc, char **argv)
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 receivers in tinv = %d", nxs);
vmess("number of shots = %d", nshots);
vmess("number of receiver/shot = %d", nx);
+ vmess("number of receiver acquisition = %d", nacq);
vmess("first model position = %.2f", fxsb);
vmess("last model position = %.2f", fxse);
vmess("first source position fxf = %.2f", fxf);
@@ -550,14 +563,14 @@ int main (int argc, char **argv)
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("k1min intermediate array = k1min_%03d'iter'.su ", istart);
vmess("Mi intermediate array = Mi_%03d'iter'.su ", istart);
}
}
/*================ initializations ================*/
- if (reci) n2out = nxs;
+ if (reci) n2out = nacq;
else n2out = nshots;
mem = Nfoc*n2out*ntfft*sizeof(float)/1048576.0;
if (verbose) {
@@ -567,15 +580,6 @@ int main (int argc, char **argv)
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;
@@ -616,16 +620,22 @@ int main (int argc, char **argv)
perc=(iend-istart)/100;if(!perc)perc=1;
if (plane_wave) {
- writeDataIter("SRCplane.su", SRC, 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("SRCplane.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, NINT(src_angle));
+ }
+ /* make DD causal again */
+ for (l=0; l<npos; l++) {
+ ix = ixpos[l];
+ j=0;
+ SRC[l*nts+j] = DD[ix*nts+j];
+ for (j = 1; j < nts; j++) {
+ SRC[l*nts+j] = DD[ix*nts+nts-j];
}
- writeDataIter("DDplane.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, NINT(src_angle));
+ }
+ if (plane_wave) {
+ writeDataIter("DDplane.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, nxs, 0, NINT(src_angle));
+ }
+ else {
+ writeDataIter("DDshot.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, nxs, 0, ishot+1);
}
free(SRC);
@@ -635,19 +645,22 @@ int main (int argc, char **argv)
/*================ initialization ================*/
+ t5 = wallclock_time();
+ memset(M0, 0, Nfoc*nxs*ntfft*sizeof(float));
+ memset(v1plus, 0, Nfoc*nxs*ntfft*sizeof(float));
/* once every 'niterskip' time-steps start from fresh M0 and do niter (~20) iterations */
if ( ((ii-istart)%niterskip==0) || (ii==istart) ) {
niterrun=niter;
recur=0;
- if (verbose>1) vmess("Doing %d iterations to reset recursion at time-sample %d\n",niterrun,ii);
+ 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++) {
- iw = NINT((ii*dt+twplane[i])/dt);
- //iw = ii;
+ for (i = 0; i < npos; i++) {
+ ix = ixpos[i];
+ iw = NINT((ii*dt+twplane[i])/dt);
for (j = 0; j < nts; j++) {
- M0[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+i*nts+j];
+ M0[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+j];
}
- /* apply mute window for samples above nts-ii */
+ /* apply mute window for samples above nts-ii */
for (j = 0; j < MIN(nts, nts-iw+isms); j++) {
M0[l*nxs*nts+i*nts+j] = 0.0;
}
@@ -658,30 +671,26 @@ int main (int argc, char **argv)
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];
+ k1min[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];
+ k1min[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+ix*nts+nts-j];
}
}
}
- if (file_iter != NULL) {
- writeDataIter("M0.su", M0, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii);
- }
}
- else { /* use f1min from previous iteration as starting point and do niterec iterations */
+ else { /* use k1min from previous iteration as starting point and do niterec iterations */
niterrun=niterec;
recur=1;
if (verbose>1) 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];
iw = NINT((ii*dt+twplane[ix])/dt);
- //iw = ii;
- M0[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts] - f1min[l*nxs*nts+i*nts+j];
+ M0[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+ix*nts] - k1min[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- M0[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j] - f1min[l*nxs*nts+i*nts+nts-j];
+ M0[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+ix*nts+nts-j] - k1min[l*nxs*nts+i*nts+nts-j];
}
+
/* apply mute window for samples above nts-ii */
for (j = 0; j < MIN(nts,nts-iw+isms); j++) {
M0[l*nxs*nts+i*nts+j] = 0.0;
@@ -692,6 +701,10 @@ int main (int argc, char **argv)
}
}
}
+ if (file_iter != NULL) {
+ writeDataIter("M0.su", M0, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii);
+ }
+
/*================ initialization ================*/
memcpy(Mi, M0, Nfoc*nxs*ntfft*sizeof(float));
@@ -705,25 +718,36 @@ int main (int argc, char **argv)
/*================ construction of Mi(-t) = - \int R(x,t) Mi(t) ================*/
- synthesis(Refl, Fop, Mi, RMi, nx, nt, nxs, nts, dt, xsyn, Nfoc,
+ synthesisp(Refl, Fop, Mi, RMi, nx, nt, nacq, 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, RMi, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter+1);
+ writeDataIter(file_iter, RMi, hdrs_out, ntfft, nacq, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter+1);
}
if (verbose >=2) {
for (l = 0; l < Nfoc; l++) {
- energyMi = 0.0;
- for (i = 0; i < npos; i++) {
- for (j = 0; j < nts; j++) {
- energyMi += RMi[l*nxs*nts+ix*nts+j]*RMi[l*nxs*nts+ix*nts+j];
- }
- }
- if (iter % 2 == 0) {
- if (iter==0) energyM0[Nfoc] = energyMi;
- vmess(" - iSyn %d: Mi at iteration %d has energy %e; relative to M0 %e", l, iter, sqrt(energyMi), sqrt(energyMi/energyM0[Nfoc]));
+ if (iter % 2 == 0) {
+ energyMi = 0.0;
+ for (i = 0; i < npos; i++) {
+ ix = ixpos[i];
+ if (recur==0) {
+ for (j = 0; j < nts; j++) {
+ energyMi += RMi[l*nacq*nts+ix*nts+j]*RMi[l*nacq*nts+ix*nts+j];
+ }
+ }
+ else {
+ mem = RMi[l*nacq*nts+ix*nts+nts-1]+DD[l*nacq*nts+ix*nts];
+ energyMi += mem*mem;
+ for (j = 1; j < nts; j++) {
+ mem = RMi[l*nacq*nts+ix*nts+nts-j]+DD[l*nacq*nts+ix*nts+j];
+ energyMi += mem*mem;
+ }
+ }
+ }
+ if ( (iter==0) && (recur==0) ) energyM0[l] = energyMi;
+ vmess(" - ii %d: Mi at iteration %d has energy %e; relative to M0 %e", ii, iter, sqrt(energyMi), sqrt(energyMi/energyM0[l]));
}
}
}
@@ -736,20 +760,19 @@ int main (int argc, char **argv)
for (i = 0; i < npos; i++) {
j = 0;
ix = ixpos[i];
- Mi[l*nxs*nts+i*nts+j] = -RMi[l*nxs*nts+ix*nts+j];
+ Mi[l*nxs*nts+i*nts+j] = RMi[l*nacq*nts+ix*nts+j];
for (j = 1; j < nts; j++) {
- Mi[l*nxs*nts+i*nts+j] = -RMi[l*nxs*nts+ix*nts+nts-j];
+ Mi[l*nxs*nts+i*nts+j] = RMi[l*nacq*nts+ix*nts+nts-j];
}
}
}
- if (iter % 2 == 0) { /* even iterations: => f_1^+(t) */
+ if (iter % 2 == 0) { /* even iterations, correlation => v_1^+(t) */
/* apply muting for the acausal part */
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
ix = ixpos[i];
iw = NINT((ii*dt+twplane[ix])/dt);
- //iw = ii;
/* apply mute window for samples after ii */
for (j = MAX(0,iw-isme); j < nts; j++) {
Mi[l*nxs*nts+i*nts+j] = 0.0;
@@ -765,64 +788,57 @@ int main (int argc, char **argv)
for (j = MAX(0,iw+shift-smooth), k=1; j < MAX(0,iw+shift); j++, k++) {
Mi[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
}
- f1plus[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- f1plus[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+j];
+ for (j = 0; j < nts; j++) {
+ v1plus[l*nxs*nts+i*nts+j] += Mi[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);
- }
-*/
+ if (file_iter != NULL) {
+ writeDataIter("v1plus.su", v1plus, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii+iter+1);
+ }
}
- else {/* odd iterations: => f_1^-(t) */
+ else {/* odd iterations, convolution => k_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 */
+ if (recur==1) { /* use k1min from previous iteration */
for (j = 0; j < nts; j++) {
- Mi[l*nxs*nts+i*nts+j] += DD[l*nxs*nts+ix*nts+j];
+ Mi[l*nxs*nts+i*nts+j] += -DD[l*nxs*nts+ix*nts+j];
}
j = 0;
- f1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+j];
+ k1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- f1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+nts-j];
+ k1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+nts-j];
}
+ //j = 0;
+ //k1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+j];
+ //for (j = 1; j < nts; j++) {
+ // k1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+nts-j];
+ //}
if (file_update != NULL) {
j=0;
- Mup[l*nxs*nts+i*nts+j] += f1min[l*nxs*nts+i*nts+j]+DD[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- Mup[l*nxs*nts+i*nts+j] += f1min[l*nxs*nts+i*nts+j]+DD[l*nxs*nts+i*nts+nts-j];
- }
+ Mup[l*nxs*nts+i*nts+j] += DD[l*nxs*nts+ix*nts+j] - k1min[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ Mup[l*nxs*nts+i*nts+j] += DD[l*nxs*nts+ix*nts+nts-j] - k1min[l*nxs*nts+i*nts+j];
+ }
}
}
else {
j = 0;
- f1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+j];
+ k1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- f1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+nts-j];
+ k1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+nts-j];
}
if (file_update != NULL) {
j=0;
- Mup[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+j];
+ Mup[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- Mup[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+nts-j];
+ Mup[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+nts-j];
}
}
}
- /* apply mute window for delta function at t=0*/
- iw = NINT((twplane[ix])/dt);
- for (j = nts-shift+smooth+iw; j < nts; j++) {
- Mi[l*nxs*nts+i*nts+j] = 0.0;
- }
- for (j = nts-shift+iw, k=0; j < MIN(nts, nts-shift+smooth+iw); j++, k++) {
- Mi[l*nxs*nts+i*nts+j] *= costaper[k];
- }
/* apply mute window for samples above nts-ii */
iw = NINT((ii*dt+twplane[ix])/dt);
- //iw = ii;
for (j = 0; j < MIN(nts,nts-iw+isms); j++) {
Mi[l*nxs*nts+i*nts+j] = 0.0;
}
@@ -832,11 +848,12 @@ int main (int argc, char **argv)
}
}
if (file_iter != NULL) {
- writeDataIter("f1min.su", f1min, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter+1);
+ writeDataIter("k1min.su", k1min, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii+iter+1);
}
} /* end else (iter) branch */
+
if (file_iter != NULL) {
- writeDataIter("Mi.su", Mi, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter+1);
+ writeDataIter("Mi.su", Mi, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii+iter+1);
}
t2 = wallclock_time();
@@ -850,7 +867,7 @@ int main (int argc, char **argv)
ix = ixpos[i];
iw = NINT((ii*dt+twplane[ix])/dt);
if ( iw<nts && iw>=0 ) {
- RR[l*nxs*nts+i*nts+iw] = f1min[l*nxs*nts+i*nts+iw];
+ RR[l*nxs*nts+i*nts+iw] = k1min[l*nxs*nts+i*nts+iw];
if (file_update != NULL) Msp[l*nxs*nts+i*nts+iw] = Mup[l*nxs*nts+i*nts+iw];
}
}
@@ -870,13 +887,16 @@ int main (int argc, char **argv)
//vmess("Remaining compute time at time-sample %d = %.2f s.",ii, tii);
}
+ ttime = wallclock_time()-t5;
+ if (verbose>2) vmess("Compute time at time-sample %d = %.3f s.",ii, ttime);
+
} /* end of time iterations ii */
free(Mi);
free(energyM0);
free(M0);
- free(f1min);
- free(f1plus);
+ free(k1min);
+ free(v1plus);
t2 = wallclock_time();
if (verbose) {
@@ -936,5 +956,3 @@ int main (int argc, char **argv)
exit(0);
}
-
-
diff --git a/marchenko/synthesis.c b/marchenko/synthesis.c
index e975d8746ffa1b44fcf34a15f9b3fa0db03772fb..7531c25a1ad7977e7b44190a6bcd164ce15cea00 100644
--- a/marchenko/synthesis.c
+++ b/marchenko/synthesis.c
@@ -89,16 +89,16 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
ctrace = (complex *)calloc(ntfft,sizeof(complex));
/* this first check is done to support an acquisition geometry that has more receiver than source
- * postions. In the first iteration the int R(x_r,x_s) Fop(x_r) d x_r results in a grid on x_s.
- * so for the next interations onlt x_s traces have to be computed on Fop */
+ * positions. In the first iteration the int R(x_r,x_s) Fop(x_r) d x_r results in a grid on x_s.
+ * so for the next interations only x_s traces have to be computed on Fop */
if (!first) {
- /* transform muted Ni (Top) to frequency domain, input for next iteration */
+ /* transform muted Ni (Top) to frequency domain, input for next iteration */
for (l = 0; l < Nfoc; l++) {
/* set Fop to zero, so new operator can be defined within ixpos points */
memset(&Fop[l*nxs*nw], 0, nxs*nw*sizeof(complex));
for (i = 0; i < npos; i++) {
- rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
ix = ixpos[i];
+ rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
for (iw=0; iw<nw; iw++) {
Fop[l*nxs*nw+iw*nxs+ix].r = ctrace[nw_low+iw].r;
Fop[l*nxs*nw+iw*nxs+ix].i = mode*ctrace[nw_low+iw].i;
@@ -107,7 +107,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
}
}
else { /* only for first call to synthesis using all nxs traces in G_d */
- /* transform G_d to frequency domain, over all nxs traces */
+ /* transform G_d to frequency domain, over all nxs traces */
first=0;
for (l = 0; l < Nfoc; l++) {
/* set Fop to zero, so new operator can be defined within all ix points */
@@ -248,6 +248,202 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
return;
}
+/* same routine as synthesis, but without the special case for the first time called to include receiver positions
+ * that are not covered by a shot position */
+
+/*================ Convolution and Integration ================*/
+/* 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 RNi has the traces in the grid of Fop, these are the x_s positions of R(x_r,x_s) */
+
+void synthesisp(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)
+{
+ int nfreq, size, inx;
+ float scl;
+ int i, j, l, m, iw, ix, k, ixsrc, il, ik;
+ float *rtrace, idxs;
+ float fxb, fxe;
+ complex *sum, *ctrace;
+ int npe;
+ static int first=1, *ixrcv;
+ static double t0, t1, t;
+
+ if (fxsb < 0) fxb = 1.001*fxsb;
+ else fxb = 0.999*fxsb;
+ if (fxse > 0) fxe = 1.001*fxse;
+ else fxe = 0.999*fxse;
+
+ size = nxs*nts;
+ nfreq = ntfft/2+1;
+ /* scale factor 1/N for backward FFT,
+ * scale dt for correlation/convolution along time,
+ * scale dx (or dxsrc) for integration over receiver (or shot) coordinates */
+ scl = 1.0*dt/((float)ntfft);
+
+#ifdef _OPENMP
+ npe = omp_get_max_threads();
+ /* parallelisation is over number of shot positions (nshots) */
+ if (npe > nshots) {
+ vmess("Number of OpenMP threads set to %d (was %d)", nshots, npe);
+ omp_set_num_threads(nshots);
+ }
+#endif
+
+ t0 = wallclock_time();
+
+ /* reset output data to zero */
+ memset(&RNi[0], 0, Nfoc*nxs*nts*sizeof(float));
+ ctrace = (complex *)calloc(ntfft,sizeof(complex));
+
+ /* transform muted Ni (Top) to frequency domain, input for next iteration */
+ for (l = 0; l < Nfoc; l++) {
+ /* set Fop to zero, so new operator can be defined within ixpos points */
+ memset(&Fop[l*nxs*nw], 0, nxs*nw*sizeof(complex));
+ for (i = 0; i < npos; i++) {
+ ix = ixpos[i];
+ rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
+ for (iw=0; iw<nw; iw++) {
+ Fop[l*nxs*nw+iw*nxs+ix].r = ctrace[nw_low+iw].r;
+ Fop[l*nxs*nw+iw*nxs+ix].i = mode*ctrace[nw_low+iw].i;
+ }
+ }
+ }
+ free(ctrace);
+
+ if (first) {
+ first=0;
+ idxs = 1.0/dxs;
+ ixrcv = (int *)malloc(nshots*nx*sizeof(int));
+ for (k=0; k<nshots; k++) {
+ for (i = 0; i < nx; i++) {
+ ixrcv[k*nx+i] = NINT((xrcv[k*nx+i]-fxsb)*idxs);
+ }
+ }
+ }
+ t1 = wallclock_time();
+ *tfft += t1 - t0;
+
+ if (reci == 0 || reci == 1) {
+
+/*================ SYNTHESIS ================*/
+
+#pragma omp parallel default(none) \
+ 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) \
+ private(l, ix, j, m, i, sum, rtrace, k, ixsrc, inx)
+{ /* start of parallel region */
+ sum = (complex *)malloc(nfreq*sizeof(complex));
+ rtrace = (float *)calloc(ntfft,sizeof(float));
+
+/* Loop over total number of shots */
+#pragma omp for schedule(guided,1)
+ for (k=0; k<nshots; k++) {
+ if ((xsrc[k] < fxb) || (xsrc[k] > fxe)) continue;
+ ixsrc = NINT((xsrc[k] - fxsb)/dxs);
+ inx = xnx[k]; /* number of traces per shot */
+
+ for (l = 0; l < Nfoc; l++) {
+ /* compute integral over receiver positions */
+ /* multiply R with Fop and sum over nx */
+ memset(&sum[0].r,0,nfreq*2*sizeof(float));
+ for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
+ for (i = 0; i < inx; i++) {
+ ix = ixrcv[k*nx+i];
+ sum[j].r += Refl[k*nw*nx+m*nx+i].r*Fop[l*nw*nxs+m*nxs+ix].r -
+ Refl[k*nw*nx+m*nx+i].i*Fop[l*nw*nxs+m*nxs+ix].i;
+ sum[j].i += Refl[k*nw*nx+m*nx+i].i*Fop[l*nw*nxs+m*nxs+ix].r +
+ Refl[k*nw*nx+m*nx+i].r*Fop[l*nw*nxs+m*nxs+ix].i;
+ }
+ }
+
+ /* transfrom result back to time domain */
+ cr1fft(sum, rtrace, ntfft, 1);
+
+ /* place result at source position ixsrc; dx = receiver distance */
+ for (j = 0; j < nts; j++)
+ RNi[l*size+ixsrc*nts+j] += rtrace[j]*scl*dx;
+
+ } /* end of Nfoc loop */
+
+ if (verbose>4) vmess("*** Shot gather %d processed ***", k);
+
+ } /* end of nparallel shots (k) loop */
+ free(sum);
+ free(rtrace);
+
+} /* end of parallel region */
+
+
+ } /* end of if reci */
+
+/* if reciprocal traces are enabled start a new loop over reciprocal shot positions */
+ if (reci != 0) {
+
+#pragma omp parallel default(none) \
+ 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) \
+ shared(Fop, size, nts, ntfft, scl) \
+ private(l, ix, j, m, i, k, sum, rtrace, ik, il, ixsrc, inx)
+{ /* start of parallel region */
+ sum = (complex *)malloc(nfreq*sizeof(complex));
+ rtrace = (float *)calloc(ntfft,sizeof(float));
+
+#pragma omp for schedule(guided,1)
+ for (k=0; k<nxs; k++) {
+ if (isxcount[k] == 0) continue;
+ ixsrc = k;
+ inx = isxcount[ixsrc]; /* number of traces per reciprocal shot */
+
+ for (l = 0; l < Nfoc; l++) {
+ /* compute integral over (reciprocal) source positions */
+ /* multiply R with Fop and sum over nx */
+ memset(&sum[0].r,0,nfreq*2*sizeof(float));
+ for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
+ for (i = 0; i < inx; i++) {
+ il = reci_xrcv[ixsrc*nxs+i];
+ ik = reci_xsrc[ixsrc*nxs+i];
+ ix = NINT((xsrc[il] - fxsb)/dxs);
+ sum[j].r += Refl[il*nw*nx+m*nx+ik].r*Fop[l*nw*nxs+m*nxs+ix].r -
+ Refl[il*nw*nx+m*nx+ik].i*Fop[l*nw*nxs+m*nxs+ix].i;
+ sum[j].i += Refl[il*nw*nx+m*nx+ik].i*Fop[l*nw*nxs+m*nxs+ix].r +
+ Refl[il*nw*nx+m*nx+ik].r*Fop[l*nw*nxs+m*nxs+ix].i;
+ }
+ }
+
+ /* transfrom result back to time domain */
+ cr1fft(sum, rtrace, ntfft, 1);
+
+ /* place result at source position ixsrc; dxsrc = shot distance */
+ for (j = 0; j < nts; j++)
+ RNi[l*size+ixsrc*nts+j] = ixmask[ixsrc]*(RNi[l*size+ixsrc*nts+j]+rtrace[j]*scl*dxsrc);
+
+ } /* end of Nfoc loop */
+
+ } /* end of parallel reciprocal shots (k) loop */
+ free(sum);
+ free(rtrace);
+
+ } /* end of parallel region */
+
+ } /* end of if reci */
+
+ t = wallclock_time() - t0;
+ if (verbose>2) {
+ vmess("OMP: parallel region = %f seconds (%d threads)", t, npe);
+ }
+
+ return;
+}
+
+
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)
{
diff --git a/marchenko3D/Makefile b/marchenko3D/Makefile
index 3afe49a69b1b3c2430dde18f665ddec80cbd616d..0bd671993adda8c7504fbf9be2b5794d33e94644 100644
--- a/marchenko3D/Makefile
+++ b/marchenko3D/Makefile
@@ -2,7 +2,8 @@
include ../Make_include
-LIBS += -L$L -lgenfft -lm $(LIBSM)
+LIBS += -L$L -lgenfft -lzfp -lm $(LIBSM)
+CFLAGS +=
#OPTC += -g -O0 -Wall
#ALL: fmute marchenko marchenko2
diff --git a/utils/getFileInfo.c b/utils/getFileInfo.c
index 5aec1a4def8a94cfdfe14f1d60a03ea434e19cdd..d93010a4f54f1023a0aafc64a57ea7c8d9731860 100644
--- a/utils/getFileInfo.c
+++ b/utils/getFileInfo.c
@@ -43,6 +43,8 @@ int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *
ret = fseeko( fp, 0, SEEK_END );
if (ret<0) perror("fseeko");
bytes = ftello( fp );
+ *xmax = hdr.gx;
+ *xmin = hdr.gx;
*n1 = hdr.ns;
if ( (hdr.trid == 1) && (hdr.dt != 0) ) {
@@ -112,9 +114,7 @@ int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *
one_shot = 1;
igath = 0;
fseeko( fp, 0, SEEK_SET );
- offset = (NINT((hdr.gx-hdr.sx)*scl*100)/100.0);
- *xmax = offset;
- *xmin = offset;
+ //offset = (NINT((hdr.gx-hdr.sx)*scl*100)/100.0);
dxrcv = *d2;
while (!end_of_file) {
@@ -125,6 +125,8 @@ int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *
sx_shot = hdr.sx;
gx_start = hdr.gx;
gx_end = hdr.gx;
+ *xmax = MAX(MAX(hdr.gx,*xmax),hdr.sx);
+ *xmin = MIN(MIN(hdr.gx,*xmin),hdr.sx);
itrace = 0;
while (one_shot) {
@@ -132,9 +134,11 @@ int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *
itrace++;
if (hdr.gx != gx_end) dxrcv = MIN(dxrcv,abs(hdr.gx-gx_end));
gx_end = hdr.gx;
- offset = (NINT((hdr.gx-hdr.sx)*scl*100)/100.0);
- *xmax = MAX(*xmax,offset);
- *xmin = MIN(*xmin,offset);
+ //offset = (NINT((hdr.gx-hdr.sx)*scl*100)/100.0);
+ //*xmax = MAX(*xmax,offset);
+ //*xmin = MIN(*xmin,offset);
+ *xmax = MAX(MAX(hdr.gx,*xmax),hdr.sx);
+ *xmin = MIN(MIN(hdr.gx,*xmin),hdr.sx);
nread = fread( &hdr, 1, TRCBYTES, fp );
if (nread != TRCBYTES) {
one_shot = 0;
@@ -167,14 +171,18 @@ int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *
fseeko( fp, -trace_sz, SEEK_END );
nread = fread( &hdr, 1, TRCBYTES, fp );
- offset = (NINT((hdr.gx-hdr.sx)*scl*100)/100.0);
- *xmax = MAX(*xmax,offset);
- *xmin = MIN(*xmin,offset);
- *xmin = MIN(sx_shot,hdr.sx*scl);
+ //offset = (NINT((hdr.gx-hdr.sx)*scl*100)/100.0);
+ //*xmax = MAX(*xmax,offset);
+ //*xmin = MIN(*xmin,offset);
+ //*xmin = MIN(sx_shot,hdr.sx*scl);
+ *xmax = MAX(MAX(hdr.gx,*xmax),hdr.sx);
+ *xmin = MIN(MIN(hdr.gx,*xmin),hdr.sx);
*ngath = ntraces/(*n2);
}
// *nxm = NINT((*xmax-*xmin)/dxrcv)+1;
*nxm = (int)ntraces;
+ *xmax *= scl;
+ *xmin *= scl;
fclose( fp );
free(trace);
diff --git a/zfp/Makefile b/zfp/Makefile
index cc696640144c360aa1126111ed47430a17795dbf..1f73d1047f2d720d0fa03fb8136b22248467a9f4 100644
--- a/zfp/Makefile
+++ b/zfp/Makefile
@@ -44,6 +44,9 @@ ifneq ($(BUILD_EXAMPLES),0)
@cd examples; $(MAKE) clean all
endif
+install:
+ -cp -rp include/* $I
+ -cp -rp lib/* $L
# run basic regression tests
test: