diff --git a/fdelmodc3D/applySource3D.c b/fdelmodc3D/applySource3D.c
index 51c21f35b62421ced0b6d98d15da04dec6dee9de..74e88f32e5039c99875b56218f1213c3c9dc3c9e 100644
--- a/fdelmodc3D/applySource3D.c
+++ b/fdelmodc3D/applySource3D.c
@@ -77,7 +77,7 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l
}
}
}
-
+
/*
* for plane wave sources the sources are placed
* around the central shot position
@@ -119,7 +119,6 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l
else { /* multi-wavelet sources */
src_ampl = src_nwav[isrc][id1]*(id2-time/dt) + src_nwav[isrc][id2]*(time/dt-id1);
}
-
if (src_ampl==0.0) continue;
if ( ((ix-ibndx)<0) || ((ix-ibndx)>mod.nx) ) continue; /* source outside grid */
if ( ((iy-ibndy)<0) || ((iy-ibndy)>mod.ny) ) continue; /* source outside grid */
diff --git a/fdelmodc3D/defineSource3D.c b/fdelmodc3D/defineSource3D.c
index 64564a85f0a1ee34732594ebee698bb77b2d90e8..30e5192bd8fd4b6de1fc9b8ce1a3f3672d9ccf82 100644
--- a/fdelmodc3D/defineSource3D.c
+++ b/fdelmodc3D/defineSource3D.c
@@ -115,7 +115,7 @@ long defineSource3D(wavPar wav, srcPar src, modPar mod, recPar rec, float **src_
optnscale = optn;
nfreqscale = optnscale/2 + 1;
}
-// fprintf(stderr,"define S optn=%li ns=%li %e nt=%li %e\n", optn, wav.ns, wav.ds, optnscale, wav.dt);
+ // fprintf(stderr,"define S optn=%li ns=%li %e nt=%li %e\n", optn, wav.ns, wav.ds, optnscale, wav.dt);
ctrace = (complex *)calloc(nfreqscale,sizeof(complex));
trace = (float *)calloc(optnscale,sizeof(float));
diff --git a/fdelmodc3D/fdelmodc3D.c b/fdelmodc3D/fdelmodc3D.c
index c9b31306fde496f3f4a108de882f6016caeaa386..57901b8aba1c1be61ff7713173ae855a57f263ad 100644
--- a/fdelmodc3D/fdelmodc3D.c
+++ b/fdelmodc3D/fdelmodc3D.c
@@ -141,6 +141,7 @@ char *sdoc[] = {
" - 3=dipole - + horizontal oriented",
" dip=0 ............. dip for double-couple source",
" strike=0 .......... strike for double-couple source",
+" rake=0 ............ rake for double-couple source",
" xsrc=middle ....... x-position of (first) shot ",
" ysrc=middle ....... y-position of (first) shot ",
" zsrc=zmin ......... z-position of (first) shot ",
@@ -270,9 +271,14 @@ char *sdoc[] = {
" NOTES: For viscoelastic media dispersion and stability are not always",
" guaranteed by the calculated criteria, especially for Q values smaller than 13",
"",
+" Original fdelmodc by",
" Jan Thorbecke 2011",
" TU Delft",
" E-mail: janth@xs4all.nl ",
+" 3D extension by",
+" Joeri Brackenhoff 2019",
+" TU Delft",
+" E-mail: jabrackenhoff@tudelft.nl ",
" 2015 Contributions from Max Holicki",
"",
NULL};
diff --git a/marchenko3D/TWtransform.c b/marchenko3D/TWtransform.c
index 8895271e4a3884b8c1a8951ae9abe6725530dc21..12bcb7a790be00fb2e16d6c8eb7b4b64f759facc 100644
--- a/marchenko3D/TWtransform.c
+++ b/marchenko3D/TWtransform.c
@@ -28,21 +28,22 @@ char *sdoc[] = {
" ",
" TWtransform - Transform data from uncompressed time domain to compressed frequency domain",
" ",
-" TWtransform file_T= file_W= [optional parameters]",
+" TWtransform file_in= file_out= [optional parameters]",
" ",
" Required parameters: ",
" ",
-" file_T= .................. File containing the uncompressed time domain data",
-" file_W= .................. Output for the compressed frequency domain data",
+" file_in= ................. File containing the uncompressed time domain data",
+" file_out= ................ Output for the (compressed) frequency domain data",
" ",
" Optional parameters: ",
" ",
-" verbose=0 ................ silent option; >0 displays info",
+" verbose=1 ................ silent option; >0 displays info",
" fmin=0 ................... minimum frequency in the output",
" fmax=70 .................. maximum frequency in the output",
" mode=1 ................... sign of the frequency transform",
-" zfp=0 .................... compress the transformed data using zfp",
-" tolerance=1e-3 ........... accuracy of the zfp compression",
+" zfp=0 .................... (=1) compress the transformed data using zfp",
+" tolerance=1e-3 ........... accuracy of the zfp compression,",
+" smaller values give more accuracy to the compressed data but will decrease the compression rate",
" weight=2.0 ............... scaling of the reflection data",
" ",
" ",
@@ -63,7 +64,7 @@ int main (int argc, char **argv)
long inx, iny, gy;
float scl, *trace, dt, dx, dy, ft, fx, fy, fmin, fmax, *cdata, scale;
double tolerance;
- char *file_T, *file_W;
+ char *file_in, *file_out;
complex *ctrace;
zfptop zfpt;
zfpmar zfpm;
@@ -71,10 +72,10 @@ int main (int argc, char **argv)
initargs(argc, argv);
requestdoc(1);
- if (!getparstring("file_T", &file_T)) file_T = "in.su";
- if (file_T==NULL) verr("No file_in is given");
- if (!getparstring("file_W", &file_W)) file_W = NULL;
- if (file_W==NULL) verr("No file_W is given");
+ if (!getparstring("file_in", &file_in)) file_in = NULL;
+ if (file_in==NULL) verr("No file_in is given");
+ if (!getparstring("file_out", &file_out)) file_out = NULL;
+ if (file_out==NULL) verr("No file_out is given");
if (!getparfloat("fmin", &fmin)) fmin = 0;
if (!getparfloat("fmax", &fmax)) fmax = 70;
if (!getpardouble("tolerance", &tolerance)) tolerance = 1e-3;
@@ -84,7 +85,7 @@ int main (int argc, char **argv)
if (!getparlong("zfp", &zfp)) zfp = 0;
nshots = 0;
- ret = getFileInfo3D(file_T, &nt, &nx, &ny, &nshots, &dt, &dx, &dy, &ft, &fx, &fy, &scl, &ntraces);
+ ret = getFileInfo3D(file_in, &nt, &nx, &ny, &nshots, &dt, &dx, &dy, &ft, &fx, &fy, &scl, &ntraces);
ntfft = loptncr(nt);
nfreq = ntfft/2+1;
@@ -97,19 +98,19 @@ int main (int argc, char **argv)
/* Reading first header */
- if (file_T == NULL) fp = stdin;
- else fp = fopen( file_T, "r" );
+ if (file_in == NULL) fp = stdin;
+ else fp = fopen( file_in, "r" );
if ( fp == NULL ) {
- fprintf(stderr,"input file %s has an error\n", file_T);
+ fprintf(stderr,"input file %s has an error\n", file_in);
perror("error in opening file: ");
fflush(stderr);
return -1;
}
- if (file_W == NULL) fp_out = stdin;
- else fp_out = fopen( file_W, "w+" );
+ if (file_out == NULL) fp_out = stdin;
+ else fp_out = fopen( file_out, "w+" );
if ( fp_out == NULL ) {
- fprintf(stderr,"input file %s has an error\n", file_W);
+ fprintf(stderr,"input file %s has an error\n", file_out);
perror("error in opening file: ");
fflush(stderr);
return -1;
@@ -324,4 +325,4 @@ long zfpcompress(float* data, long nx, long ny, long nz, double tolerance, zfpma
// fclose(file);
return 1;
-}
\ No newline at end of file
+}
diff --git a/marchenko3D/demo/README b/marchenko3D/demo/marchenko2D/README
similarity index 100%
rename from marchenko3D/demo/README
rename to marchenko3D/demo/marchenko2D/README
diff --git a/marchenko3D/demo/ScientificReports/NatureSnapshots.tex b/marchenko3D/demo/marchenko2D/ScientificReports/NatureSnapshots.tex
similarity index 100%
rename from marchenko3D/demo/ScientificReports/NatureSnapshots.tex
rename to marchenko3D/demo/marchenko2D/ScientificReports/NatureSnapshots.tex
diff --git a/marchenko3D/demo/ScientificReports/README b/marchenko3D/demo/marchenko2D/ScientificReports/README
similarity index 100%
rename from marchenko3D/demo/ScientificReports/README
rename to marchenko3D/demo/marchenko2D/ScientificReports/README
diff --git a/marchenko3D/demo/ScientificReports/back_injrate_planes.scr b/marchenko3D/demo/marchenko2D/ScientificReports/back_injrate_planes.scr
similarity index 100%
rename from marchenko3D/demo/ScientificReports/back_injrate_planes.scr
rename to marchenko3D/demo/marchenko2D/ScientificReports/back_injrate_planes.scr
diff --git a/marchenko3D/demo/ScientificReports/backpropf2.scr b/marchenko3D/demo/marchenko2D/ScientificReports/backpropf2.scr
similarity index 100%
rename from marchenko3D/demo/ScientificReports/backpropf2.scr
rename to marchenko3D/demo/marchenko2D/ScientificReports/backpropf2.scr
diff --git a/marchenko3D/demo/ScientificReports/check.scr b/marchenko3D/demo/marchenko2D/ScientificReports/check.scr
similarity index 100%
rename from marchenko3D/demo/ScientificReports/check.scr
rename to marchenko3D/demo/marchenko2D/ScientificReports/check.scr
diff --git a/marchenko3D/demo/ScientificReports/clean b/marchenko3D/demo/marchenko2D/ScientificReports/clean
similarity index 100%
rename from marchenko3D/demo/ScientificReports/clean
rename to marchenko3D/demo/marchenko2D/ScientificReports/clean
diff --git a/marchenko3D/demo/ScientificReports/direct.scr b/marchenko3D/demo/marchenko2D/ScientificReports/direct.scr
similarity index 100%
rename from marchenko3D/demo/ScientificReports/direct.scr
rename to marchenko3D/demo/marchenko2D/ScientificReports/direct.scr
diff --git a/marchenko3D/demo/ScientificReports/epsBack.scr b/marchenko3D/demo/marchenko2D/ScientificReports/epsBack.scr
similarity index 100%
rename from marchenko3D/demo/ScientificReports/epsBack.scr
rename to marchenko3D/demo/marchenko2D/ScientificReports/epsBack.scr
diff --git a/marchenko3D/demo/ScientificReports/initialFocus.scr b/marchenko3D/demo/marchenko2D/ScientificReports/initialFocus.scr
similarity index 100%
rename from marchenko3D/demo/ScientificReports/initialFocus.scr
rename to marchenko3D/demo/marchenko2D/ScientificReports/initialFocus.scr
diff --git a/marchenko3D/demo/ScientificReports/marchenko.scr b/marchenko3D/demo/marchenko2D/ScientificReports/marchenko.scr
similarity index 100%
rename from marchenko3D/demo/ScientificReports/marchenko.scr
rename to marchenko3D/demo/marchenko2D/ScientificReports/marchenko.scr
diff --git a/marchenko3D/demo/ScientificReports/model.scr b/marchenko3D/demo/marchenko2D/ScientificReports/model.scr
similarity index 100%
rename from marchenko3D/demo/ScientificReports/model.scr
rename to marchenko3D/demo/marchenko2D/ScientificReports/model.scr
diff --git a/marchenko3D/demo/ScientificReports/remove_direct.scr b/marchenko3D/demo/marchenko2D/ScientificReports/remove_direct.scr
similarity index 100%
rename from marchenko3D/demo/ScientificReports/remove_direct.scr
rename to marchenko3D/demo/marchenko2D/ScientificReports/remove_direct.scr
diff --git a/marchenko3D/demo/ScientificReports/shotslurm.scr b/marchenko3D/demo/marchenko2D/ScientificReports/shotslurm.scr
similarity index 100%
rename from marchenko3D/demo/ScientificReports/shotslurm.scr
rename to marchenko3D/demo/marchenko2D/ScientificReports/shotslurm.scr
diff --git a/marchenko3D/demo/WS15/MarchenkoWorkshop.pdf b/marchenko3D/demo/marchenko2D/WS15/MarchenkoWorkshop.pdf
similarity index 100%
rename from marchenko3D/demo/WS15/MarchenkoWorkshop.pdf
rename to marchenko3D/demo/marchenko2D/WS15/MarchenkoWorkshop.pdf
diff --git a/marchenko3D/demo/WS15/README.1 b/marchenko3D/demo/marchenko2D/WS15/README.1
similarity index 100%
rename from marchenko3D/demo/WS15/README.1
rename to marchenko3D/demo/marchenko2D/WS15/README.1
diff --git a/marchenko3D/demo/WS15/README.2 b/marchenko3D/demo/marchenko2D/WS15/README.2
similarity index 100%
rename from marchenko3D/demo/WS15/README.2
rename to marchenko3D/demo/marchenko2D/WS15/README.2
diff --git a/marchenko3D/demo/WS15/README.3 b/marchenko3D/demo/marchenko2D/WS15/README.3
similarity index 100%
rename from marchenko3D/demo/WS15/README.3
rename to marchenko3D/demo/marchenko2D/WS15/README.3
diff --git a/marchenko3D/demo/WS15/README.4 b/marchenko3D/demo/marchenko2D/WS15/README.4
similarity index 100%
rename from marchenko3D/demo/WS15/README.4
rename to marchenko3D/demo/marchenko2D/WS15/README.4
diff --git a/marchenko3D/demo/WS15/README.5 b/marchenko3D/demo/marchenko2D/WS15/README.5
similarity index 100%
rename from marchenko3D/demo/WS15/README.5
rename to marchenko3D/demo/marchenko2D/WS15/README.5
diff --git a/marchenko3D/demo/WS15/job.pbs b/marchenko3D/demo/marchenko2D/WS15/job.pbs
similarity index 100%
rename from marchenko3D/demo/WS15/job.pbs
rename to marchenko3D/demo/marchenko2D/WS15/job.pbs
diff --git a/marchenko3D/demo/WS15/setup.sh b/marchenko3D/demo/marchenko2D/WS15/setup.sh
similarity index 100%
rename from marchenko3D/demo/WS15/setup.sh
rename to marchenko3D/demo/marchenko2D/WS15/setup.sh
diff --git a/marchenko3D/demo/oneD/README b/marchenko3D/demo/marchenko2D/oneD/README
similarity index 100%
rename from marchenko3D/demo/oneD/README
rename to marchenko3D/demo/marchenko2D/oneD/README
diff --git a/marchenko3D/demo/oneD/backProp_f2sum_movie.scr b/marchenko3D/demo/marchenko2D/oneD/backProp_f2sum_movie.scr
similarity index 100%
rename from marchenko3D/demo/oneD/backProp_f2sum_movie.scr
rename to marchenko3D/demo/marchenko2D/oneD/backProp_f2sum_movie.scr
diff --git a/marchenko3D/demo/oneD/backpropf2.scr b/marchenko3D/demo/marchenko2D/oneD/backpropf2.scr
similarity index 100%
rename from marchenko3D/demo/oneD/backpropf2.scr
rename to marchenko3D/demo/marchenko2D/oneD/backpropf2.scr
diff --git a/marchenko3D/demo/oneD/clean b/marchenko3D/demo/marchenko2D/oneD/clean
similarity index 100%
rename from marchenko3D/demo/oneD/clean
rename to marchenko3D/demo/marchenko2D/oneD/clean
diff --git a/marchenko3D/demo/oneD/conv.gnp b/marchenko3D/demo/marchenko2D/oneD/conv.gnp
similarity index 100%
rename from marchenko3D/demo/oneD/conv.gnp
rename to marchenko3D/demo/marchenko2D/oneD/conv.gnp
diff --git a/marchenko3D/demo/oneD/conv.txt b/marchenko3D/demo/marchenko2D/oneD/conv.txt
similarity index 100%
rename from marchenko3D/demo/oneD/conv.txt
rename to marchenko3D/demo/marchenko2D/oneD/conv.txt
diff --git a/marchenko3D/demo/oneD/epsBackprop.scr b/marchenko3D/demo/marchenko2D/oneD/epsBackprop.scr
similarity index 100%
rename from marchenko3D/demo/oneD/epsBackprop.scr
rename to marchenko3D/demo/marchenko2D/oneD/epsBackprop.scr
diff --git a/marchenko3D/demo/oneD/epsCompare.scr b/marchenko3D/demo/marchenko2D/oneD/epsCompare.scr
similarity index 100%
rename from marchenko3D/demo/oneD/epsCompare.scr
rename to marchenko3D/demo/marchenko2D/oneD/epsCompare.scr
diff --git a/marchenko3D/demo/oneD/epsIterwithLabels.scr b/marchenko3D/demo/marchenko2D/oneD/epsIterwithLabels.scr
similarity index 100%
rename from marchenko3D/demo/oneD/epsIterwithLabels.scr
rename to marchenko3D/demo/marchenko2D/oneD/epsIterwithLabels.scr
diff --git a/marchenko3D/demo/oneD/epsMarchenkoIter.scr b/marchenko3D/demo/marchenko2D/oneD/epsMarchenkoIter.scr
similarity index 100%
rename from marchenko3D/demo/oneD/epsMarchenkoIter.scr
rename to marchenko3D/demo/marchenko2D/oneD/epsMarchenkoIter.scr
diff --git a/marchenko3D/demo/oneD/epsModel.scr b/marchenko3D/demo/marchenko2D/oneD/epsModel.scr
similarity index 100%
rename from marchenko3D/demo/oneD/epsModel.scr
rename to marchenko3D/demo/marchenko2D/oneD/epsModel.scr
diff --git a/marchenko3D/demo/oneD/figAppendix.scr b/marchenko3D/demo/marchenko2D/oneD/figAppendix.scr
similarity index 100%
rename from marchenko3D/demo/oneD/figAppendix.scr
rename to marchenko3D/demo/marchenko2D/oneD/figAppendix.scr
diff --git a/marchenko3D/demo/oneD/initialFocus.scr b/marchenko3D/demo/marchenko2D/oneD/initialFocus.scr
similarity index 100%
rename from marchenko3D/demo/oneD/initialFocus.scr
rename to marchenko3D/demo/marchenko2D/oneD/initialFocus.scr
diff --git a/marchenko3D/demo/oneD/initialFocus1300.scr b/marchenko3D/demo/marchenko2D/oneD/initialFocus1300.scr
similarity index 100%
rename from marchenko3D/demo/oneD/initialFocus1300.scr
rename to marchenko3D/demo/marchenko2D/oneD/initialFocus1300.scr
diff --git a/marchenko3D/demo/oneD/initialFocusPlane.scr b/marchenko3D/demo/marchenko2D/oneD/initialFocusPlane.scr
similarity index 100%
rename from marchenko3D/demo/oneD/initialFocusPlane.scr
rename to marchenko3D/demo/marchenko2D/oneD/initialFocusPlane.scr
diff --git a/marchenko3D/demo/oneD/iterations.scr b/marchenko3D/demo/marchenko2D/oneD/iterations.scr
similarity index 100%
rename from marchenko3D/demo/oneD/iterations.scr
rename to marchenko3D/demo/marchenko2D/oneD/iterations.scr
diff --git a/marchenko3D/demo/oneD/marchenko.scr b/marchenko3D/demo/marchenko2D/oneD/marchenko.scr
similarity index 100%
rename from marchenko3D/demo/oneD/marchenko.scr
rename to marchenko3D/demo/marchenko2D/oneD/marchenko.scr
diff --git a/marchenko3D/demo/oneD/marchenkoIter.scr b/marchenko3D/demo/marchenko2D/oneD/marchenkoIter.scr
similarity index 100%
rename from marchenko3D/demo/oneD/marchenkoIter.scr
rename to marchenko3D/demo/marchenko2D/oneD/marchenkoIter.scr
diff --git a/marchenko3D/demo/oneD/marchenkoPlane.scr b/marchenko3D/demo/marchenko2D/oneD/marchenkoPlane.scr
similarity index 100%
rename from marchenko3D/demo/oneD/marchenkoPlane.scr
rename to marchenko3D/demo/marchenko2D/oneD/marchenkoPlane.scr
diff --git a/marchenko3D/demo/oneD/marchenkodt.scr b/marchenko3D/demo/marchenko2D/oneD/marchenkodt.scr
similarity index 100%
rename from marchenko3D/demo/oneD/marchenkodt.scr
rename to marchenko3D/demo/marchenko2D/oneD/marchenkodt.scr
diff --git a/marchenko3D/demo/oneD/migr.scr b/marchenko3D/demo/marchenko2D/oneD/migr.scr
similarity index 100%
rename from marchenko3D/demo/oneD/migr.scr
rename to marchenko3D/demo/marchenko2D/oneD/migr.scr
diff --git a/marchenko3D/demo/oneD/model.scr b/marchenko3D/demo/marchenko2D/oneD/model.scr
similarity index 100%
rename from marchenko3D/demo/oneD/model.scr
rename to marchenko3D/demo/marchenko2D/oneD/model.scr
diff --git a/marchenko3D/demo/oneD/model2.scr b/marchenko3D/demo/marchenko2D/oneD/model2.scr
similarity index 100%
rename from marchenko3D/demo/oneD/model2.scr
rename to marchenko3D/demo/marchenko2D/oneD/model2.scr
diff --git a/marchenko3D/demo/oneD/p5all.scr b/marchenko3D/demo/marchenko2D/oneD/p5all.scr
similarity index 100%
rename from marchenko3D/demo/oneD/p5all.scr
rename to marchenko3D/demo/marchenko2D/oneD/p5all.scr
diff --git a/marchenko3D/demo/oneD/primaries.scr b/marchenko3D/demo/marchenko2D/oneD/primaries.scr
similarity index 100%
rename from marchenko3D/demo/oneD/primaries.scr
rename to marchenko3D/demo/marchenko2D/oneD/primaries.scr
diff --git a/marchenko3D/demo/oneD/primariesFrame.scr b/marchenko3D/demo/marchenko2D/oneD/primariesFrame.scr
similarity index 100%
rename from marchenko3D/demo/oneD/primariesFrame.scr
rename to marchenko3D/demo/marchenko2D/oneD/primariesFrame.scr
diff --git a/marchenko3D/demo/oneD/referenceShot.scr b/marchenko3D/demo/marchenko2D/oneD/referenceShot.scr
similarity index 100%
rename from marchenko3D/demo/oneD/referenceShot.scr
rename to marchenko3D/demo/marchenko2D/oneD/referenceShot.scr
diff --git a/marchenko3D/demo/test/test.scr b/marchenko3D/demo/marchenko2D/test/test.scr
similarity index 100%
rename from marchenko3D/demo/test/test.scr
rename to marchenko3D/demo/marchenko2D/test/test.scr
diff --git a/marchenko3D/demo/twoD/README b/marchenko3D/demo/marchenko2D/twoD/README
similarity index 100%
rename from marchenko3D/demo/twoD/README
rename to marchenko3D/demo/marchenko2D/twoD/README
diff --git a/marchenko3D/demo/twoD/backProp_makemovie.scr b/marchenko3D/demo/marchenko2D/twoD/backProp_makemovie.scr
similarity index 100%
rename from marchenko3D/demo/twoD/backProp_makemovie.scr
rename to marchenko3D/demo/marchenko2D/twoD/backProp_makemovie.scr
diff --git a/marchenko3D/demo/twoD/backpropf2.scr b/marchenko3D/demo/marchenko2D/twoD/backpropf2.scr
similarity index 100%
rename from marchenko3D/demo/twoD/backpropf2.scr
rename to marchenko3D/demo/marchenko2D/twoD/backpropf2.scr
diff --git a/marchenko3D/demo/twoD/check.scr b/marchenko3D/demo/marchenko2D/twoD/check.scr
similarity index 100%
rename from marchenko3D/demo/twoD/check.scr
rename to marchenko3D/demo/marchenko2D/twoD/check.scr
diff --git a/marchenko3D/demo/twoD/clean b/marchenko3D/demo/marchenko2D/twoD/clean
similarity index 100%
rename from marchenko3D/demo/twoD/clean
rename to marchenko3D/demo/marchenko2D/twoD/clean
diff --git a/marchenko3D/demo/twoD/direct.scr b/marchenko3D/demo/marchenko2D/twoD/direct.scr
similarity index 100%
rename from marchenko3D/demo/twoD/direct.scr
rename to marchenko3D/demo/marchenko2D/twoD/direct.scr
diff --git a/marchenko3D/demo/twoD/eps.scr b/marchenko3D/demo/marchenko2D/twoD/eps.scr
similarity index 100%
rename from marchenko3D/demo/twoD/eps.scr
rename to marchenko3D/demo/marchenko2D/twoD/eps.scr
diff --git a/marchenko3D/demo/twoD/epsPlane.scr b/marchenko3D/demo/marchenko2D/twoD/epsPlane.scr
similarity index 100%
rename from marchenko3D/demo/twoD/epsPlane.scr
rename to marchenko3D/demo/marchenko2D/twoD/epsPlane.scr
diff --git a/marchenko3D/demo/twoD/epsPrimaries.scr b/marchenko3D/demo/marchenko2D/twoD/epsPrimaries.scr
similarity index 100%
rename from marchenko3D/demo/twoD/epsPrimaries.scr
rename to marchenko3D/demo/marchenko2D/twoD/epsPrimaries.scr
diff --git a/marchenko3D/demo/twoD/homg_reference.scr b/marchenko3D/demo/marchenko2D/twoD/homg_reference.scr
similarity index 100%
rename from marchenko3D/demo/twoD/homg_reference.scr
rename to marchenko3D/demo/marchenko2D/twoD/homg_reference.scr
diff --git a/marchenko3D/demo/twoD/homgpng.scr b/marchenko3D/demo/marchenko2D/twoD/homgpng.scr
similarity index 100%
rename from marchenko3D/demo/twoD/homgpng.scr
rename to marchenko3D/demo/marchenko2D/twoD/homgpng.scr
diff --git a/marchenko3D/demo/twoD/homgview.scr b/marchenko3D/demo/marchenko2D/twoD/homgview.scr
similarity index 100%
rename from marchenko3D/demo/twoD/homgview.scr
rename to marchenko3D/demo/marchenko2D/twoD/homgview.scr
diff --git a/marchenko3D/demo/twoD/initialFocus.scr b/marchenko3D/demo/marchenko2D/twoD/initialFocus.scr
similarity index 100%
rename from marchenko3D/demo/twoD/initialFocus.scr
rename to marchenko3D/demo/marchenko2D/twoD/initialFocus.scr
diff --git a/marchenko3D/demo/twoD/initialFocus_pbs.scr b/marchenko3D/demo/marchenko2D/twoD/initialFocus_pbs.scr
similarity index 100%
rename from marchenko3D/demo/twoD/initialFocus_pbs.scr
rename to marchenko3D/demo/marchenko2D/twoD/initialFocus_pbs.scr
diff --git a/marchenko3D/demo/twoD/initialFocus_slurm.scr b/marchenko3D/demo/marchenko2D/twoD/initialFocus_slurm.scr
similarity index 100%
rename from marchenko3D/demo/twoD/initialFocus_slurm.scr
rename to marchenko3D/demo/marchenko2D/twoD/initialFocus_slurm.scr
diff --git a/marchenko3D/demo/twoD/initialPlane.scr b/marchenko3D/demo/marchenko2D/twoD/initialPlane.scr
similarity index 100%
rename from marchenko3D/demo/twoD/initialPlane.scr
rename to marchenko3D/demo/marchenko2D/twoD/initialPlane.scr
diff --git a/marchenko3D/demo/twoD/marchenko.scr b/marchenko3D/demo/marchenko2D/twoD/marchenko.scr
similarity index 100%
rename from marchenko3D/demo/twoD/marchenko.scr
rename to marchenko3D/demo/marchenko2D/twoD/marchenko.scr
diff --git a/marchenko3D/demo/twoD/marchenkoPlane.scr b/marchenko3D/demo/marchenko2D/twoD/marchenkoPlane.scr
similarity index 100%
rename from marchenko3D/demo/twoD/marchenkoPlane.scr
rename to marchenko3D/demo/marchenko2D/twoD/marchenkoPlane.scr
diff --git a/marchenko3D/demo/twoD/marchenko_ray.pbs b/marchenko3D/demo/marchenko2D/twoD/marchenko_ray.pbs
similarity index 100%
rename from marchenko3D/demo/twoD/marchenko_ray.pbs
rename to marchenko3D/demo/marchenko2D/twoD/marchenko_ray.pbs
diff --git a/marchenko3D/demo/twoD/marchenko_ray.scr b/marchenko3D/demo/marchenko2D/twoD/marchenko_ray.scr
similarity index 100%
rename from marchenko3D/demo/twoD/marchenko_ray.scr
rename to marchenko3D/demo/marchenko2D/twoD/marchenko_ray.scr
diff --git a/marchenko3D/demo/twoD/model.scr b/marchenko3D/demo/marchenko2D/twoD/model.scr
similarity index 100%
rename from marchenko3D/demo/twoD/model.scr
rename to marchenko3D/demo/marchenko2D/twoD/model.scr
diff --git a/marchenko3D/demo/twoD/primaries.scr b/marchenko3D/demo/marchenko2D/twoD/primaries.scr
similarity index 100%
rename from marchenko3D/demo/twoD/primaries.scr
rename to marchenko3D/demo/marchenko2D/twoD/primaries.scr
diff --git a/marchenko3D/demo/twoD/rayvsp.scr b/marchenko3D/demo/marchenko2D/twoD/rayvsp.scr
similarity index 100%
rename from marchenko3D/demo/twoD/rayvsp.scr
rename to marchenko3D/demo/marchenko2D/twoD/rayvsp.scr
diff --git a/marchenko3D/demo/twoD/referenceShot.scr b/marchenko3D/demo/marchenko2D/twoD/referenceShot.scr
similarity index 100%
rename from marchenko3D/demo/twoD/referenceShot.scr
rename to marchenko3D/demo/marchenko2D/twoD/referenceShot.scr
diff --git a/marchenko3D/demo/twoD/remove_direct.scr b/marchenko3D/demo/marchenko2D/twoD/remove_direct.scr
similarity index 100%
rename from marchenko3D/demo/twoD/remove_direct.scr
rename to marchenko3D/demo/marchenko2D/twoD/remove_direct.scr
diff --git a/marchenko3D/demo/twoD/shots_pbs.scr b/marchenko3D/demo/marchenko2D/twoD/shots_pbs.scr
similarity index 100%
rename from marchenko3D/demo/twoD/shots_pbs.scr
rename to marchenko3D/demo/marchenko2D/twoD/shots_pbs.scr
diff --git a/marchenko3D/demo/twoD/shots_slurm.scr b/marchenko3D/demo/marchenko2D/twoD/shots_slurm.scr
similarity index 100%
rename from marchenko3D/demo/twoD/shots_slurm.scr
rename to marchenko3D/demo/marchenko2D/twoD/shots_slurm.scr
diff --git a/marchenko3D/demo/marchenko3D/oneD/README b/marchenko3D/demo/marchenko3D/oneD/README
new file mode 100644
index 0000000000000000000000000000000000000000..52d17c72cb350b8b428f3c3944492402f14d5a17
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/README
@@ -0,0 +1,6 @@
+In this directory, you can run the 3D marchenko code.
+!!!!!!!!!!!!!!!!!!!!!!!!!!WARNING!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+Due to the size of the data, one requires around 256 GB of RAM. Keep this limitation in mind
+
+To avoid the need for large storage space on your directory, we reccomend using the ZFP compressed option.
+This is the standard option included in the data. You can just execute execute.sh and everything will be done for you. It will take a considerable amount of time, depending on your machine.
diff --git a/marchenko3D/demo/marchenko3D/oneD/execute.sh b/marchenko3D/demo/marchenko3D/oneD/execute.sh
new file mode 100755
index 0000000000000000000000000000000000000000..f183f9e9dc334595fd842b125d6a07c7f710b8df
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/execute.sh
@@ -0,0 +1,15 @@
+#!/bin/bash
+
+./model.sh
+
+./shotmodel.sh
+
+./makeR.sh
+
+./farrmod.sh
+
+./fmute3D.sh
+
+./zfp.sh
+
+./marchenko3d.sh
diff --git a/marchenko3D/demo/marchenko3D/oneD/farrmod.sh b/marchenko3D/demo/marchenko3D/oneD/farrmod.sh
new file mode 100755
index 0000000000000000000000000000000000000000..4cfbf5211d87618c33e6987dff124e4250fc669b
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/farrmod.sh
@@ -0,0 +1,32 @@
+#!/bin/bash
+
+dt=0.001
+
+makewave fp=15 fmax=22 dt=${dt} file_out=wavefpmod.su nt=8192 t0=0.2 scale=1 scfft=0
+
+export OMP_NUM_THREADS=8
+
+#Model shot record in middle of model
+fdelmodc3D \
+ file_cp=cp3d.su ischeme=1 iorder=4 \
+ file_den=ro3d.su \
+ file_src=wavefpmod.su \
+ file_rcv=farrmod.su \
+ src_type=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.2 \
+ verbose=2 \
+ tmod=1.2 \
+ dxrcv=10.0 \
+ xrcv1=-1000 xrcv2=1000 \
+ dyrcv=10.0 \
+ yrcv1=-300 yrcv2=300 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 ysrc=0 zsrc=900 \
+ ntaper=61 \
+ left=4 right=4 top=4 bottom=4 front=4 back=4 \
diff --git a/marchenko3D/demo/marchenko3D/oneD/fmute3D.sh b/marchenko3D/demo/marchenko3D/oneD/fmute3D.sh
new file mode 100755
index 0000000000000000000000000000000000000000..8b89c7354eaef93cff259c4df8e0240660626654
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/fmute3D.sh
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+#Mute the coda from the modeled first arrival
+fmute3D file_mute=farrmod_rp.su file_shot=farrmod_rp.su file_out=farr.su above=2 shift=18 smooth=10 hw=5
diff --git a/marchenko3D/demo/marchenko3D/oneD/makeR.sh b/marchenko3D/demo/marchenko3D/oneD/makeR.sh
new file mode 100755
index 0000000000000000000000000000000000000000..e31cf348249a47faec8c10db9635619c0c2a7c69
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/makeR.sh
@@ -0,0 +1,5 @@
+#!/bin/bash
+
+#create the reflection data that was modeled in 1 1D medium
+makeR1D file_in=shotx10y10.su file_out=reflx10y10.su verbose=2 nxrcv=201 nyrcv=61 nxsrc=201 nysrc=61
+
diff --git a/marchenko3D/demo/marchenko3D/oneD/marchenko3d.sh b/marchenko3D/demo/marchenko3D/oneD/marchenko3d.sh
new file mode 100755
index 0000000000000000000000000000000000000000..68439928f64b33cb290aa83fdbebce4712b8c16f
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/marchenko3d.sh
@@ -0,0 +1,19 @@
+#!/bin/bash
+
+#Marchenko using time data
+#marchenko3D file_shot=reflx10y10.su verbose=2 \
+# file_tinv=farr.su \
+# niter=10 shift=15 smooth=10 hw=5 \
+# file_green=green_timeshot.su file_file_f2=f2_timeshot.su
+
+#Marchenko using frequency data
+#marchenko3D file_shotw=reflx10y10_W.bin verbose=2 \
+# file_tinv=farr.su \
+# niter=10 shift=15 smooth=10 hw=5 \
+# file_green=green_freqshot.su file_file_f2=f2_freqshot.su
+
+#Marchenko using zfp data
+marchenko3D file_shotzfp=reflx10y10_zfp.bin verbose=2 \
+ file_tinv=farr.su \
+ niter=10 shift=15 smooth=10 hw=5 \
+ file_green=green_zfpshot.su file_file_f2=f2_zfpshot.su
diff --git a/marchenko3D/demo/marchenko3D/oneD/model.sh b/marchenko3D/demo/marchenko3D/oneD/model.sh
new file mode 100755
index 0000000000000000000000000000000000000000..3725e0aff7c937946560dccc3f1a64d6720f4379
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/model.sh
@@ -0,0 +1,40 @@
+#!/bin/bash
+
+#adjust this PATH to where the code is installed
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+dx=10
+ny=181
+dy=10000
+
+#define gridded model for FD computations
+makemod sizex=5000 sizez=1400 dx=$dx dz=$dx cp0=1800 ro0=1000 \
+ orig=-2500,0 file_base=model5.su verbose=2 \
+ intt=def x=-2500,2500 z=400,400 poly=0 cp=2300 ro=3000 \
+ intt=def x=-2500,2500 z=700,700 poly=0 cp=2000 ro=1100 \
+ intt=def x=-2500,2500 z=1100,1100 poly=0 cp=2500 ro=4000
+
+#define homogenoeus model to compute direct wave only
+makemod sizex=5000 sizez=300 dx=$dx dz=$dx cp0=1800 ro0=1000 \
+ orig=-2500,0 file_base=hom.su verbose=2
+
+rm cp3d.su
+rm ro3d.su
+rm cphom3d.su
+rm rohom3d.su
+
+for (( i=0; i<$ny; i++ ));do
+
+ (( ypos = -1*($ny-1)/2 + $i ))
+ (( yloc = $dy*$ypos ))
+ echo $yloc
+
+ sushw < model5_cp.su key=gy,scalel a=$yloc,-1000 >> cp3d.su
+ sushw < model5_ro.su key=gy,scalel a=$yloc,-1000 >> ro3d.su
+ sushw < hom_cp.su key=gy,scalel a=$yloc,-1000 >> cphom3d.su
+ sushw < hom_ro.su key=gy,scalel a=$yloc,-1000 >> rohom3d.su
+
+done
+
+
+
diff --git a/marchenko3D/demo/marchenko3D/oneD/shotmodel.sh b/marchenko3D/demo/marchenko3D/oneD/shotmodel.sh
new file mode 100755
index 0000000000000000000000000000000000000000..cd6c947964220e7a2bbfd5573a8d9adcae7eb9bc
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/shotmodel.sh
@@ -0,0 +1,62 @@
+#!/bin/bash
+
+dt=0.001
+
+#define wavelet for modeling R
+makewave w=fw fmin=0 flef=5 frig=25 fmax=30 dt=$dt file_out=wavefw.su nt=4096 t0=0.3 scale=0 scfft=1
+
+#Set the number of threads for parallel jobs. More is not always better
+export OMP_NUM_THREADS=8
+
+#Model shot record in middle of model
+fdelmodc3D \
+ file_cp=cp3d.su ischeme=1 iorder=4 \
+ file_den=ro3d.su \
+ file_src=wavefw.su \
+ file_rcv=shotx10y10.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=10.0 \
+ xrcv1=-2000 xrcv2=2000 \
+ dyrcv=10.0 \
+ yrcv1=-600 yrcv2=600 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 ysrc=0 zsrc=0 \
+ ntaper=41 \
+ left=4 right=4 top=4 bottom=4 front=4 back=4
+
+#Model direct wave only in middle of model
+fdelmodc3D \
+ file_cp=cphom3d.su ischeme=1 iorder=4 \
+ file_den=rohom3d.su \
+ file_src=wavefw.su \
+ file_rcv=shotx10y10hom.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=10.0 \
+ xrcv1=-2000 xrcv2=2000 \
+ dyrcv=10.0 \
+ yrcv1=-600 yrcv2=600 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 ysrc=0 zsrc=0 \
+ ntaper=41 \
+ left=4 right=4 top=4 bottom=4 front=4 back=4
+
+#subtract direct wave from full model shot record: this defines R
+sudiff shotx10y10_rp.su shotx10y10hom_rp.su > shotx10y10.su
diff --git a/marchenko3D/demo/marchenko3D/oneD/zfp.sh b/marchenko3D/demo/marchenko3D/oneD/zfp.sh
new file mode 100755
index 0000000000000000000000000000000000000000..45f436827b4c75bfecd7b5c535053002b2930443
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/zfp.sh
@@ -0,0 +1,9 @@
+#!/bin/bash
+
+# Pre-transform the data to the frequency domain
+#TWtransform file_in=reflx10y10.su file_out=reflx10y10_W.bin verbose=1 fmin=0 fmax=30 zfp=0 tolerance=1e-7
+
+# Pre-transform the data to the frequency domain and apply zfp compression
+TWtransform file_in=reflx10y10.su file_out=reflx10y10_zfp.bin verbose=1 fmin=0 fmax=30 zfp=1 tolerance=1e-7
+
+rm reflx10y10.su
diff --git a/marchenko3D/fmute3D.c b/marchenko3D/fmute3D.c
index 56152cf9098553d6db33077e67b9011624cb10d3..b8fea43b5865f93957d3506d7c5b8f94e50999b0 100644
--- a/marchenko3D/fmute3D.c
+++ b/marchenko3D/fmute3D.c
@@ -302,9 +302,11 @@ int main (int argc, char **argv)
/* search forward in y-trace direction from maximum in file */
for (i = iymax+1; i < ny1; i++) {
+ tstart = MAX(0, (maxval[(i-1)*nx1+ixmax]-hw));
+ tend = MIN(nt1-1, (maxval[(i-1)*nx1+ixmax]+hw));
tmax=0.0;
- jmax = 0;
- for (j = 0; j < nt1; j++) {
+ jmax = tstart;
+ for (j = tstart; j < tend; j++) {
lmax = fabs(tmpdata[i*nx1*nt1+ixmax*nt1+j]);
if (lmax > tmax) {
jmax = j;
@@ -347,9 +349,11 @@ int main (int argc, char **argv)
/* search backward in y-trace direction from maximum in file */
for (i = iymax-1; i >= 0; i--) {
+ tstart = MAX(0, (maxval[(i+1)*nx1+ixmax]-hw));
+ tend = MIN(nt1-1, (maxval[(i+1)*nx1+ixmax]+hw));
tmax=0.0;
- jmax = 0;
- for (j = 0; j < nt1; j++) {
+ jmax = tstart;
+ for (j = tstart; j < tend; j++) {
lmax = fabs(tmpdata[i*nx1*nt1+ixmax*nt1+j]);
if (lmax > tmax) {
jmax = j;
diff --git a/marchenko3D/homogeneousg3D.c b/marchenko3D/homogeneousg3D.c
index 5da29e05ba74f8b67185dba4bf7c9eedffbf70cc..6f3d0f7c7ad064236c2e194db84a1a187e0fad08 100644
--- a/marchenko3D/homogeneousg3D.c
+++ b/marchenko3D/homogeneousg3D.c
@@ -7,6 +7,25 @@
#include <assert.h>
#include <genfft.h>
+/*
+The schemes in this module use a variety of retrieval representations
+For more information about Green's function retrieval see:
+Brackenhoff, J., Thorbecke, J., & Wapenaar, K. (2019).
+Virtual sources and receivers in the real Earth: Considerations for practical applications.
+Journal of Geophysical Research: Solid Earth, 124, 11802– 11821.
+https://doi.org/10.1029/2019JB018485
+
+Brackenhoff, J., Thorbecke, J., and Wapenaar, K.:
+Monitoring of induced distributed double-couple sources using Marchenko-based virtual receivers.
+Solid Earth, 10, 1301–1319,
+https://doi.org/10.5194/se-10-1301-2019, 2019.
+
+Wapenaar, K., Brackenhoff, J., Thorbecke, J. et al.
+Virtual acoustics in inhomogeneous media with single-sided access.
+Sci Rep 8, 2497 (2018).
+https://doi.org/10.1038/s41598-018-20924-x
+*/
+
#ifndef MAX
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#endif
@@ -29,7 +48,6 @@ void convol2(float *data1, float *data2, float *con, long nrec, long nsam, float
long readSnapData3D(char *filename, float *data, segy *hdrs, long nsnaps, long nx, long ny, long nz, long sx, long ex, long sy, long ey, long sz, long ez);
long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, float *d1, float *d2, float *d3, float *f1, float *f2, float *f3, float *sclsxgxsygy, long *nxm);
-//void kxwfilter(float *data, long nt, long nx, float dt, float dx, float fmin, float fmax, float angle, float cp, float perc);
void timeDiff(float *data, long nsam, long nrec, float dt, float fmin, float fmax, long opt);
void depthDiff(float *data, long nsam, long nrec, float dt, float dx, float fmin, float fmax, float c, long opt);
@@ -38,22 +56,19 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *f1p, float *f1
{
char *file_inp;
long i, j, k, l, ret, scheme, count=0, n_source;
- long is, kxwfilt, n1, n2, n3, ngath, ntraces, zsrc;
+ long is, n1, n2, n3, ngath, ntraces, zsrc;
float d1, d2, d3, f1, f2, f3, scl;
- float *conv, fmin, fmax, alpha, cp, rho, perc;
+ float *conv, fmin, fmax, cp, rho;
float *greenjkz, *input, *inputjkz, *tmp1, *tmp2;
double t0, t2, tfft;
segy *hdr_inp;
if (!getparstring("file_inp", &file_inp)) file_inp = NULL;
if (!getparlong("scheme", &scheme)) scheme = 0;
- if (!getparlong("kxwfilt", &kxwfilt)) kxwfilt = 0;
if (!getparfloat("fmin", &fmin)) fmin = 0.0;
if (!getparfloat("fmax", &fmax)) fmax = 100.0;
- if (!getparfloat("alpha", &alpha)) alpha = 65.0;
if (!getparfloat("cp", &cp)) cp = 1000.0;
if (!getparfloat("rho", &rho)) rho = 1000.0;
- if (!getparfloat("perc", &perc)) perc = 0.15;
tfft = 0.0;
ret = 0;
@@ -70,6 +85,8 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *f1p, float *f1
if (NINT(d2*1000.0)!=NINT(dx*1000.0)) verr("dx sampling between input (%.3e) and retrieved (%.3e) is not equal",d2,dx);
if (NINT(d3*1000.0)!=NINT(dy*1000.0)) verr("dy sampling between input (%.3e) and retrieved (%.3e) is not equal",d3,dy);
+ scl = dx*dy*dt;
+
input = (float *)calloc(n_source*nx*ny*nt,sizeof(float));
hdr_inp = (segy *)calloc(n_source*nx*ny,sizeof(segy));
readSnapData3D(file_inp, input, hdr_inp, n_source, nx, ny, nt, 0, nx, 0, ny, 0, nt);
@@ -86,7 +103,22 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *f1p, float *f1
sz[l] = hdr_inp[l*nx*ny].sdepth;
}
- if (scheme==1) { // Scale the Green's functions if the classical scheme is used
+ if (scheme==0) {
+ if (verbose) vmess("Marchenko Homogeneous Green's function retrieval with G source");
+ }
+ else if (scheme==1) {
+ if (verbose) vmess("Marchenko Homogeneous Green's function retrieval with f2 source");
+ }
+ else if (scheme==2) {
+ if (verbose) vmess("Marchenko Green's function retrieval with source depending on position");
+ }
+ else if (scheme==3) {
+ if (verbose) vmess("Marchenko Green's function retrieval with G source");
+ }
+ else if (scheme==4) {
+ if (verbose) vmess("Marchenko Green's function retrieval with f2 source");
+ }
+ else if (scheme==5) { // Scale the Green's functions if the classical scheme is used
if (verbose) vmess("Classical Homogeneous Green's function retrieval");
greenjkz = (float *)calloc(Nfoc*nx*ny*nt,sizeof(float));
inputjkz = (float *)calloc(n_source*nx*ny*nt,sizeof(float));
@@ -108,26 +140,14 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *f1p, float *f1
}
}
}
- else if (scheme==2) {
- if (verbose) vmess("Marchenko Green's function retrieval with G source");
- }
- else if (scheme==6) {
- if (verbose) vmess("Marchenko Green's function retrieval with f2 source");
- }
- else if (scheme==7) {
- if (verbose) vmess("Marchenko Green's function retrieval with source depending on position");
- }
- else if (scheme==3) {
+ else if (scheme==6) {
if (verbose) vmess("Marchenko Homogeneous Green's function retrieval with multiple sources");
if (verbose) vmess("Looping over %li source positions",n_source);
}
- else if (scheme==4) {
+ else if (scheme==7) {
if (verbose) vmess("Back propagation with multiple sources");
if (verbose) vmess("Looping over %li source positions",n_source);
}
- else if (scheme==5) {
- if (verbose) vmess("Marchenko Homogeneous Green's function retrieval with f2 source");
- }
else if (scheme==8) { // 0=f1p 1=f1m
if (verbose) vmess("f1+ redatuming");
if (n_source<2) verr("Not enough input for the homogeneous Green's function");
@@ -150,12 +170,11 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *f1p, float *f1
if (verbose) vmess("2i IM(f1) redatuming");
inputjkz = (float *)calloc(n_source*nx*ny*nt,sizeof(float));
for (k = 0; k < ny; k++) {
+ depthDiff(&input[k*nx*nt] , nt, nx, dt, dx, fmin, fmax, cp, 1);
for (l = 0; l < nx*nt; l++) {
inputjkz[k*nx*nt+l] = input[k*nx*nt+l];
}
conjugate(&inputjkz[k*nx*nt], nt, nx, dt);
- depthDiff(&inputjkz[k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
- depthDiff(&input[k*nx*nt] , nt, nx, dt, dx, fmin, fmax, cp, 1);
}
}
else {
@@ -166,11 +185,11 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *f1p, float *f1
private(i,j,k,is,conv,tmp1,tmp2,zsrc)
{
conv = (float *)calloc(nx*nt,sizeof(float));
- if (scheme==1) {
+ if (scheme==5) {
tmp1 = (float *)calloc(nx*nt,sizeof(float));
tmp2 = (float *)calloc(nx*nt,sizeof(float));
}
- if (scheme==3 || scheme==8 || scheme==9 || scheme==10) tmp1 = (float *)calloc(nx*nt,sizeof(float));
+ if (scheme==6 || scheme==8 || scheme==9 || scheme==10) tmp1 = (float *)calloc(nx*nt,sizeof(float));
#pragma omp for schedule(guided,1)
for (l = 0; l < Nfoc; l++) {
@@ -185,70 +204,77 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *f1p, float *f1
depthDiff(&f2p[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
convol(&input[k*nx*nt], &f2p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
timeDiff(conv, nt, nx, dt, fmin, fmax, -3);
- if (kxwfilt) {
- //kxwfilter(conv, nt, nx, dt, dx, fmin, fmax, alpha, cp, perc);
- }
for (i=0; i<nx; i++) {
for (j=0; j<nt/2; j++) {
- HomG[(j+nt/2)*Nfoc+l] += conv[i*nt+j]/rho;
- HomG[j*Nfoc+l] += conv[i*nt+(j+nt/2)]/rho;
+ HomG[(j+nt/2)*Nfoc+l] += scl*conv[i*nt+j]/rho;
+ HomG[j*Nfoc+l] += scl*conv[i*nt+(j+nt/2)]/rho;
}
}
}
}
- else if (scheme==5) { //Marchenko representation with f2 source
+ else if (scheme==1) { //Marchenko representation with f2 source
for (k = 0; k < ny; k++) {
depthDiff(&green[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
convol(&input[k*nx*nt], &green[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
timeDiff(conv, nt, nx, dt, fmin, fmax, -3);
- if (kxwfilt) {
- //kxwfilter(conv, nt, nx, dt, dx, fmin, fmax, alpha, cp, perc);
- }
for (i=0; i<nx; i++) {
for (j=0; j<nt/2; j++) {
- HomG[(j+nt/2)*Nfoc+l] += conv[i*nt+j]/rho;
- HomG[j*Nfoc+l] += conv[i*nt+(j+nt/2)]/rho;
+ HomG[(j+nt/2)*Nfoc+l] += scl*conv[i*nt+j]/rho;
+ HomG[j*Nfoc+l] += scl*conv[i*nt+(j+nt/2)]/rho;
}
}
}
}
- else if (scheme==8) { // f1+ redatuming 0=f1p 1=f1m
+ else if (scheme==2) { //Marchenko representation without time-reversal using varying sources
for (k = 0; k < ny; k++) {
- convol2(&input[0*ny*nx*nt+k*nx*nt], &f1p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, fmin, fmax, 1);
- convol2(&input[1*ny*nx*nt+k*nx*nt], &f1m[l*ny*nx*nt+k*nx*nt], tmp1, nx, nt, dt, fmin, fmax, 1);
- for (i=0; i<nx; i++) {
- for (j=0; j<nt/2; j++) {
- HomG[(j+nt/2)*Nfoc+l] -= 2.0*(conv[i*nt+j] + tmp1[i*nt+j])/rho;
- HomG[j*Nfoc+l] -= 2.0*(conv[i*nt+(j+nt/2)] + tmp1[i*nt+(j+nt/2)])/rho;
+ depthDiff(&f2p[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ depthDiff(&green[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ for (is=0; is<(n_source/2); is+=2) {
+ zsrc = labs(hdr_inp[is].selev);
+ if (zsrc > NINT(1000.0*zsyn[l])) {
+ if (verbose > 1) vmess("Homogeneous Green's function at %li uses G source (zsrc=%li)",NINT(1000.0*zsyn[l]));
+ convol(&input[is*ny*nx*nt+k*nx*nt], &f2p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
+ }
+ else {
+ if (verbose > 1) vmess("Homogeneous Green's function at %li uses f2 source (zsrc=%li)",NINT(1000.0*zsyn[l]));
+ convol(&input[(is+1)*ny*nx*nt+k*nx*nt], &green[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
+ }
+ timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
+ for (i=0; i<nx; i++) {
+ for (j=0; j<nt/2; j++) {
+ HomG[(j+nt/2)*Nfoc+l] += 2.0*scl*conv[i*nt+j]/rho;
+ HomG[j*Nfoc+l] += 2.0*scl*conv[i*nt+(j+nt/2)]/rho;
+ }
}
}
}
}
- else if (scheme==9) { // f1- redatuming 0=f1p 1=f1m
+ else if (scheme==3) { //Marchenko representation without time-reversal G source
for (k = 0; k < ny; k++) {
- convol2(&input[0*ny*nx*nt+k*nx*nt], &f1m[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, fmin, fmax, 1);
- convol2(&input[1*ny*nx*nt+k*nx*nt], &f1p[l*ny*nx*nt+k*nx*nt], tmp1, nx, nt, dt, fmin, fmax, 1);
+ depthDiff(&f2p[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ convol2(&input[k*nx*nt], &f2p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, fmin, fmax, 1);
for (i=0; i<nx; i++) {
for (j=0; j<nt/2; j++) {
- HomG[(j+nt/2)*Nfoc+l] += 2.0*(conv[i*nt+j] + tmp1[i*nt+j])/rho;
- HomG[j*Nfoc+l] += 2.0*(conv[i*nt+(j+nt/2)] + tmp1[i*nt+(j+nt/2)])/rho;
+ HomG[(j+nt/2)*Nfoc+l] += 2.0*scl*conv[i*nt+j]/rho;
+ HomG[j*Nfoc+l] += 2.0*scl*conv[i*nt+(j+nt/2)]/rho;
}
}
}
}
- else if (scheme==10) { // 2i IM(f1) redatuming
+ else if (scheme==4) { //Marchenko representation without time-reversal f2 source
for (k = 0; k < ny; k++) {
- convol2(&input[k*nx*nt] , &f1m[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, fmin, fmax, 2);
- convol2(&inputjkz[k*nx*nt], &f1p[l*ny*nx*nt+k*nx*nt], tmp1, nx, nt, dt, fmin, fmax, 2);
+ depthDiff(&green[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
+ convol(&input[k*nx*nt], &green[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
+ timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
for (i=0; i<nx; i++) {
for (j=0; j<nt/2; j++) {
- HomG[(j+nt/2)*Nfoc+l] += 4.0*(conv[i*nt+j] - tmp1[i*nt+j])/rho;
- HomG[j*Nfoc+l] += 4.0*(conv[i*nt+(j+nt/2)] - tmp1[i*nt+(j+nt/2)])/rho;
+ HomG[(j+nt/2)*Nfoc+l] += 2.0*scl*conv[i*nt+j]/rho;
+ HomG[j*Nfoc+l] += 2.0*scl*conv[i*nt+(j+nt/2)]/rho;
}
}
}
}
- else if (scheme==1) { //classical representation
+ else if (scheme==5) { //classical representation
for (k = 0; k < ny; k++) {
convol(&greenjkz[l*ny*nx*nt+k*nx*nt], &input[k*nx*nt], tmp1, nx, nt, dt, 0);
convol(&green[l*ny*nx*nt+k*nx*nt], &inputjkz[k*nx*nt], tmp2, nx, nt, dt, 0);
@@ -260,108 +286,89 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *f1p, float *f1
timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
for (i=0; i<nx; i++) {
for (j=0; j<nt/2; j++) {
- HomG[(j+nt/2)*Nfoc+l] += conv[i*nt+j]/rho;
- HomG[j*Nfoc+l] += conv[i*nt+(j+nt/2)]/rho;
- }
- }
- }
- }
- else if (scheme==2) { //Marchenko representation without time-reversal G source
- for (k = 0; k < ny; k++) {
- depthDiff(&f2p[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
- convol2(&input[k*nx*nt], &f2p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, fmin, fmax, 1);
- for (i=0; i<nx; i++) {
- for (j=0; j<nt/2; j++) {
- HomG[(j+nt/2)*Nfoc+l] += 2.0*conv[i*nt+j]/rho;
- HomG[j*Nfoc+l] += 2.0*conv[i*nt+(j+nt/2)]/rho;
- }
- }
- }
- }
- else if (scheme==6) { //Marchenko representation without time-reversal f2 source
- for (k = 0; k < ny; k++) {
- depthDiff(&green[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
- convol(&input[k*nx*nt], &green[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
- timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
- for (i=0; i<nx; i++) {
- for (j=0; j<nt/2; j++) {
- HomG[(j+nt/2)*Nfoc+l] += 2.0*conv[i*nt+j]/rho;
- HomG[j*Nfoc+l] += 2.0*conv[i*nt+(j+nt/2)]/rho;
+ HomG[(j+nt/2)*Nfoc+l] += scl*conv[i*nt+j]/rho;
+ HomG[j*Nfoc+l] += scl*conv[i*nt+(j+nt/2)]/rho;
}
}
}
}
- else if (scheme==7) { //Marchenko representation without time-reversal using varying sources
- for (k = 0; k < ny; k++) {
- depthDiff(&f2p[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
- depthDiff(&green[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
- for (is=0; is<(n_source/2); is+=2) {
- zsrc = labs(hdr_inp[is].selev);
- if (zsrc > NINT(1000.0*zsyn[l])) {
- if (verbose > 1) vmess("Homogeneous Green's function at %li uses G source (zsrc=%li)",NINT(1000.0*zsyn[l]));
- convol(&input[is*ny*nx*nt+k*nx*nt], &f2p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
- }
- else {
- if (verbose > 1) vmess("Homogeneous Green's function at %li uses f2 source (zsrc=%li)",NINT(1000.0*zsyn[l]));
- convol(&input[(is+1)*ny*nx*nt+k*nx*nt], &green[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
- }
- timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
- for (i=0; i<nx; i++) {
- for (j=0; j<nt/2; j++) {
- HomG[(j+nt/2)*Nfoc+l] += 2.0*conv[i*nt+j]/rho;
- HomG[j*Nfoc+l] += 2.0*conv[i*nt+(j+nt/2)]/rho;
- }
- }
- }
- }
- }
- else if (scheme==3) { //Marchenko representation with multiple shot gathers
+ else if (scheme==6) { //Marchenko representation with multiple shot gathers
for (k = 0; k < ny; k++) {
depthDiff(&f2p[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
for (is=0; is<n_source; is++) {
convol(&input[is*ny*nx*nt+k*nx*nt], &f2p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
timeDiff(conv, nt, nx, dt, fmin, fmax, -3);
- if (kxwfilt) {
- //kxwfilter(conv, nt, nx, dt, dx, fmin, fmax, alpha, cp, perc);
- }
for (i=0; i<nx; i++) {
for (j=0; j<nt/2; j++) {
- HomG[is*nt*Nfoc+(j+nt/2)*Nfoc+l] += conv[i*nt+j]/rho;
- HomG[is*nt*Nfoc+j*Nfoc+l] += conv[i*nt+(j+nt/2)]/rho;
+ HomG[is*nt*Nfoc+(j+nt/2)*Nfoc+l] += scl*conv[i*nt+j]/rho;
+ HomG[is*nt*Nfoc+j*Nfoc+l] += scl*conv[i*nt+(j+nt/2)]/rho;
}
}
}
}
}
- else if (scheme==4) { //Marchenko representation with multiple shot gathers without time-reversal
+ else if (scheme==7) { //Marchenko representation with multiple shot gathers without time-reversal
for (k = 0; k < ny; k++) {
depthDiff(&f2p[l*ny*nx*nt+k*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
for (is=0; is<n_source; is++) {
convol(&input[is*ny*nx*nt+k*nx*nt], &f2p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, 0);
timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
- if (kxwfilt) {
- //kxwfilter(conv, nt, nx, dt, dx, fmin, fmax, alpha, cp, perc);
- }
for (i=0; i<nx; i++) {
for (j=0; j<nt/2; j++) {
- HomG[is*nt*Nfoc+(j+nt/2)*Nfoc+l] += conv[i*nt+j]/rho;
- HomG[is*nt*Nfoc+j*Nfoc+l] += conv[i*nt+(j+nt/2)]/rho;
+ HomG[is*nt*Nfoc+(j+nt/2)*Nfoc+l] += scl*conv[i*nt+j]/rho;
+ HomG[is*nt*Nfoc+j*Nfoc+l] += scl*conv[i*nt+(j+nt/2)]/rho;
}
}
}
}
}
+ else if (scheme==8) { // f1+ redatuming 0=f1p 1=f1m
+ for (k = 0; k < ny; k++) {
+ convol2(&input[0*ny*nx*nt+k*nx*nt], &f1p[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, fmin, fmax, 1);
+ convol2(&input[1*ny*nx*nt+k*nx*nt], &f1m[l*ny*nx*nt+k*nx*nt], tmp1, nx, nt, dt, fmin, fmax, 1);
+ for (i=0; i<nx; i++) {
+ for (j=0; j<nt/2; j++) {
+ HomG[(j+nt/2)*Nfoc+l] -= 2.0*scl*(conv[i*nt+j] + tmp1[i*nt+j])/rho;
+ HomG[j*Nfoc+l] -= 2.0*scl*(conv[i*nt+(j+nt/2)] + tmp1[i*nt+(j+nt/2)])/rho;
+ }
+ }
+ }
+ }
+ else if (scheme==9) { // f1- redatuming 0=f1p 1=f1m
+ for (k = 0; k < ny; k++) {
+ convol2(&input[0*ny*nx*nt+k*nx*nt], &f1m[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, fmin, fmax, 1);
+ convol2(&input[1*ny*nx*nt+k*nx*nt], &f1p[l*ny*nx*nt+k*nx*nt], tmp1, nx, nt, dt, fmin, fmax, 1);
+ for (i=0; i<nx; i++) {
+ for (j=0; j<nt/2; j++) {
+ HomG[(j+nt/2)*Nfoc+l] += 2.0*scl*(conv[i*nt+j] + tmp1[i*nt+j])/rho;
+ HomG[j*Nfoc+l] += 2.0*scl*(conv[i*nt+(j+nt/2)] + tmp1[i*nt+(j+nt/2)])/rho;
+ }
+ }
+ }
+ }
+ else if (scheme==10) { // 2i IM(f1) redatuming
+ for (k = 0; k < ny; k++) {
+ convol2(&input[k*nx*nt] , &f1m[l*ny*nx*nt+k*nx*nt], conv, nx, nt, dt, fmin, fmax, 2);
+ convol2(&inputjkz[k*nx*nt], &f1p[l*ny*nx*nt+k*nx*nt], tmp1, nx, nt, dt, fmin, fmax, 2);
+ for (i=0; i<nx; i++) {
+ for (j=0; j<nt/2; j++) {
+ HomG[(j+nt/2)*Nfoc+l] += 4.0*scl*(conv[i*nt+j] - tmp1[i*nt+j])/rho;
+ HomG[j*Nfoc+l] += 4.0*scl*(conv[i*nt+(j+nt/2)] - tmp1[i*nt+(j+nt/2)])/rho;
+ }
+ }
+ }
+ }
}
free(conv);
- if (scheme==1) {
+ if (scheme==5) {
free(tmp1);
free(tmp2);
}
- if (scheme==3) free(tmp1);
- if (scheme==4) free(tmp1);
+ if (scheme==6) free(tmp1);
+ if (scheme==7) free(tmp1);
}
- if (scheme==1) {
+ if (scheme==5) {
free(input);
free(inputjkz);
free(greenjkz);
diff --git a/marchenko3D/imaging3D.c b/marchenko3D/imaging3D.c
index 6245d6e07f7da9f3c88332cfef0ec63629941acd..f7bd20feef88067e6ac9840b6485b66a577daba0 100644
--- a/marchenko3D/imaging3D.c
+++ b/marchenko3D/imaging3D.c
@@ -7,6 +7,14 @@
#include <assert.h>
#include <genfft.h>
+/*
+The imaging is computed using the double-focusing method
+For more information about the double-focusing method see:
+Myrna Staring, Roberto Pereira, Huub Douma, Joost van der Neut, and Kees Wapenaar, (2018),
+"Source-receiver Marchenko redatuming on field data using an adaptive double-focusing method," GEOPHYSICS 83: S579-S590.
+https://doi.org/10.1190/geo2017-0796.1
+*/
+
double wallclock_time(void);
void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout);
diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index fdfb59058a57b16ceec72fe586553d81078c0ab4..6581addd64a8d94ea68e0b29b24826fdf39d7841 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -95,9 +95,13 @@ char *sdoc[] = {
" ",
" Required parameters: ",
" ",
+" First arrival input options:",
" file_tinv= ............... direct arrival from focal point: G_d",
" file_ray= ................ direct arrival from raytimes",
-" file_shot= ............... Reflection response: R",
+" Shot data input options:",
+" file_shot= ............... Reflection response (time data): R(t)",
+" file_shotw= .............. Reflection response (frequency data): R(w)",
+" file_shotzfp= ............ Reflection response (frequency compressed data): zfp[R(w)]",
" ",
" Optional parameters: ",
" ",
@@ -112,32 +116,35 @@ char *sdoc[] = {
" MUTE-WINDOW ",
" file_amp= ................ amplitudes for the raytime estimation",
" file_wav= ................ Wavelet applied to the raytime data",
-" above=0 .................. mute above(1), around(0) or below(-1) the first travel times of file_tinv",
+" above=0 .................. mute above(1), around(0) or below(-1) the travel times of the first arrival",
" shift=12 ................. number of points above(positive) / below(negative) travel time for mute",
" hw=8 ..................... window in time samples to look for maximum in next trace",
" smooth=5 ................. number of points to smooth mute with cosine window",
" plane_wave=0 ............. enable plane-wave illumination function"
-// angles to be implemented currently only angle=0 works.
-//" src_angle=0 .............. angle of plane source array",
-//" src_velo=1500 ............ velocity to use in src_angle definition",
" REFLECTION RESPONSE CORRECTION ",
-" scale=2 .................. scale factor of R for summation of Ni with G_d",
+" scale=2 .................. scale factor of R for summation of Ni with G_d (only for time shot data)",
" 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",
" HOMOGENEOUS GREEN'S FUNCTION RETRIEVAL OPTIONS ",
+" file_homg= ............... output file with homogeneous Green's function ",
+" The homogeneous Green's function is computed if a filename is given",
" file_inp= ................ Input source function for the retrieval",
-" scheme=0 ................. Scheme for homogeneous Green's function retrieval",
-" .......................... scheme=0 Marchenko HomG retrieval with time-reversal",
-" .......................... scheme=1 Classical retrieval scheme",
-" .......................... scheme=2 Marchenko HomG retrieval without time-reversal",
-" .......................... scheme=3 Back propagation with multiple sources",
-" .......................... scheme=4 Marchenko HomG retrieval with multiple sources",
-" kxwfilt=0 ................ Apply a dip filter before integration",
-" alpha=65.0 ............... Alpha filter for the kxwfilter",
-" perc=0.15 ................ Percentage for the kxwfilter",
+" scheme=0 ................. Scheme for the retrieval",
+" .......................... scheme=0 Marchenko homogeneous Green's function retrieval with G source",
+" .......................... scheme=1 Marchenko homogeneous Green's function retrieval with f2 source",
+" .......................... scheme=2 Marchenko Green's function retrieval with source depending on virtual receiver location",
+" .......................... scheme=3 Marchenko Green's function retrieval with G source",
+" .......................... scheme=4 Marchenko Green's function retrieval with f2 source",
+" .......................... scheme=5 Classical homogeneous Green's function retrieval",
+" .......................... scheme=6 Marchenko homogeneous Green's function retrieval with multiple G sources",
+" .......................... scheme=7 Marchenko Green's function retrieval with multiple G sources",
+" .......................... scheme=8 f1+ redatuming",
+" .......................... scheme=9 f1- redatuming",
+" .......................... scheme=10 2i IM(f1) redatuming",
" cp=1000.0 ................ Velocity of upper layer for certain operations",
" rho=1000.0 ............... Density of upper layer for certain operations",
+" IMAGING",
+" file_imag= ............... output file with image ",
+" The image is computed if a filename is given",
" OUTPUT DEFINITION ",
" file_green= .............. output file with full Green function(s)",
" file_gplus= .............. output file with G+ ",
@@ -147,12 +154,10 @@ char *sdoc[] = {
" file_f2= ................. output file with f2 (=p+) ",
" file_pplus= .............. output file with p+ ",
" file_pmin= ............... output file with p- ",
-" file_imag= ............... output file with image ",
-" file_homg= ............... output file with homogeneous Green's function ",
" file_ampscl= ............. output file with estimated amplitudes ",
" file_iter= ............... output file with -Ni(-t) for each iteration",
" compact=0 ................ Write out homg and imag in compact format",
-" .......................... WARNING! This write-out cannot be displayed with SU"
+" .......................... WARNING! This write-out cannot be displayed with SU",
" verbose=0 ................ silent option; >0 displays info",
" ",
" ",
@@ -169,16 +174,16 @@ int main (int argc, char **argv)
long i, j, l, k, ret, nshots, nxshot, nyshot, Nfoc, nt, nx, ny, nts, nxs, nys, ngath;
long size, n1, n2, n3, ntap, tap, dxi, dyi, ntraces, pad, *sx, *sy, *sz;
long nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
- long reci, countmin, mode, n2out, n3out, verbose, ntfft;
+ long reci, mode, n2out, n3out, verbose, ntfft;
long iter, niter, tracf, *muteW, *tsynW, ampest, plane_wave, t0shift;
long hw, smooth, above, shift, *ixpos, *iypos, npos, ix, iy, nzim, nxim, nyim, compact;
- long nshots_r, *isxcount, *reci_xsrc, *reci_xrcv;
+ long *isxcount, *reci_xsrc, *reci_xrcv;
float fmin, fmax, *tapersh, *tapersy, fxf, fyf, dxf, dyf, *xsrc, *ysrc, *xrcv, *yrcv, *zsyn, *zsrc, *xrcvsyn, *yrcvsyn;
double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi, *energyN0;
float d1, d2, d3, f1, f2, f3, fxsb, fxse, fysb, fyse, ft, fx, fy, *xsyn, *ysyn, dxsrc, dysrc;
float *green, *f2p, *pmin, *G_d, dt, dx, dy, dxs, dys, scl, mem;
float *f1plus, *f1min, *iRN, *Ni, *Nig, *trace, *Gmin, *Gplus, *HomG;
- float xmin, xmax, ymin, ymax, scale, tsq, Q, f0, *tmpdata;
+ float scale, *tmpdata;
float *ixmask, *iymask, *ampscl, *Gd, *Image, dzim;
float grad2rad, p, src_angle, src_velo, *mutetest;
complex *Refl, *Fop;
@@ -225,16 +230,13 @@ int main (int argc, char **argv)
if (!getparfloat("fmax", &fmax)) fmax = 70.0;
if (!getparlong("reci", &reci)) reci = 0;
if (!getparfloat("scale", &scale)) scale = 2.0;
- if (!getparfloat("tsq", &tsq)) tsq = 0.0;
- if (!getparfloat("Q", &Q)) Q = 0.0;
- if (!getparfloat("f0", &f0)) f0 = 0.0;
if (!getparlong("tap", &tap)) tap = 0;
if (!getparlong("ntap", &ntap)) ntap = 0;
if (!getparlong("pad", &pad)) pad = 0;
if (!getparlong("ampest", &est)) ampest = 0;
if(!getparlong("niter", &niter)) niter = 10;
- if(!getparlong("hw", &hw)) hw = 15;
+ if(!getparlong("hw", &hw)) hw = 8;
if(!getparlong("smooth", &smooth)) smooth = 5;
if(!getparlong("above", &above)) above = 0;
if(!getparlong("shift", &shift)) shift=12;
@@ -305,7 +307,6 @@ int main (int argc, char **argv)
nw_high = MIN((long)(fmax*ntfft*dt), nfreq-1);
nw = nw_high - nw_low + 1;
scl = 1.0/((float)ntfft);
- if (!getparlong("countmin", &countmin)) countmin = 0.3*nx*ny;
/*================ Allocating all data arrays ================*/
@@ -393,19 +394,6 @@ int main (int argc, char **argv)
}
}
- // mutetest = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
- // for (i=0; i<Nfoc*nys*nxs*ntfft; i++) {
- // mutetest[i] = 1.0;
- // }
- // applyMute3D(mutetest, muteW, smooth, above, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift);
-
- // fp_out = fopen( "mute.bin", "w+" );
- // for (i=0; i<nx*ny; i++) {
- // fprintf(fp_out,"%.7f\n",mutetest[i]);
- // }
- // fclose(fp_out);
-
-
/* define tapers to taper edges of acquisition */
if (tap == 1 || tap == 3) {
for (j = 0; j < ntap; j++)
@@ -768,6 +756,7 @@ int main (int argc, char **argv)
}
}
}
+ if (above==-2) applyMute3D(f1min, muteW, smooth, 0, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift);
}
else { /* plane wave scheme */
for (l = 0; l < Nfoc; l++) {
@@ -781,6 +770,7 @@ int main (int argc, char **argv)
}
}
}
+ if (above==-2) applyMute3D(f1min, muteW, smooth, 0, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift);
}
}
else {/* odd iterations update: => f_1^+(t) */
@@ -1063,7 +1053,7 @@ int main (int argc, char **argv)
hdrs_Nfoc[0].sx = sx[0];
hdrs_Nfoc[0].sy = sy[0];
hdrs_Nfoc[0].sdepth = sz[0];
- hdrs_Nfoc[0].f1 = roundf(zsyn[0]*100.0)/1000.0;
+ hdrs_Nfoc[0].f1 = roundf(zsyn[0]*1000.0)/1000.0;
hdrs_Nfoc[0].f2 = roundf(xsyn[0]*1000.0)/1000.0;
hdrs_Nfoc[0].ungpow = roundf(ysyn[0]*1000.0)/1000.0;
hdrs_Nfoc[0].d1 = roundf(dzim*1000.0)/1000.0;
diff --git a/marchenko3D/readTinvData3D.c b/marchenko3D/readTinvData3D.c
index 572ed3dcd1c41dd62c2ad923e84681932acbfe4d..47a46b96742b6d55d5d5234c61ea081bb9664243 100644
--- a/marchenko3D/readTinvData3D.c
+++ b/marchenko3D/readTinvData3D.c
@@ -194,9 +194,11 @@ long readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
/* search forward in y-trace direction from maximum in file */
for (i = iymax+1; i < ny1; i++) {
+ tstart = MAX(0, (maxval[isyn*nxy+(i-1)*nx+ixmax]-hw));
+ tend = MIN(nt-1, (maxval[isyn*nxy+(i-1)*nx+ixmax]+hw));
tmax=0.0;
- jmax = 0;
- for(j = 0; j <= nt; j++) {
+ jmax = tstart;
+ for(j = tstart; j <= tend; j++) {
lmax = fabs(tinv[ig+i*nx*ntfft+ixmax*ntfft+j]);
if (lmax > tmax) {
jmax = j;
@@ -239,9 +241,11 @@ long readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
/* search backward in y-trace direction from maximum in file */
for (i = iymax-1; i >= 0; i--) {
+ tstart = MAX(0, (maxval[isyn*nxy+(i+1)*nx+ixmax]-hw));
+ tend = MIN(nt-1, (maxval[isyn*nxy+(i+1)*nx+ixmax]+hw));
tmax=0.0;
- jmax = 0;
- for(j = 0; j <= nt; j++) {
+ jmax = tstart;
+ for(j = tstart; j <= tend; j++) {
lmax = fabs(tinv[ig+i*nx*ntfft+ixmax*ntfft+j]);
if (lmax > tmax) {
jmax = j;
diff --git a/marchenko3D/synthesis3D.c b/marchenko3D/synthesis3D.c
index 41a66b0dfc787a2f1e24ec1044e6b5a5851507e0..70dab54944648a14c90ee6077f1c1828ba6f9d0c 100644
--- a/marchenko3D/synthesis3D.c
+++ b/marchenko3D/synthesis3D.c
@@ -143,10 +143,12 @@ long *ixpos, long *iypos, long npos, double *tfft, long *isxcount, long *reci_xs
long i, j, l, m, iw, ix, iy, k, isrc, il, ik, nxy, nxys;
float *rtrace, idxs, idys, fxb, fyb, fxe, fye;
complex *sum, *ctrace;
- long npe;
+ long npe, norm;
static long first=1, *ixrcv, *iyrcv;
static double t0, t1, t;
+ if (!getparlong("norm", &norm)) norm = 0;
+
if (fxsb < 0) fxb = 1.001*fxsb;
else fxb = 0.999*fxsb;
if (fysb < 0) fyb = 1.001*fysb;
@@ -164,7 +166,12 @@ long *ixpos, long *iypos, long npos, double *tfft, long *isxcount, long *reci_xs
/* scale factor 1/N for backward FFT,
* scale dt for correlation/convolution along time,
* scale dx*dy (or dxsrc*dysrc) for integration over receiver (or shot) coordinates */
- scl = 1.0*dt/((float)ntfft);
+ if (norm==0) { //pressure normalization
+ scl = (1.0*dt*dx*dy)/((float)ntfft);
+ }
+ else { // flux normalization
+ scl = 1.0/((float)ntfft);
+ }
#ifdef _OPENMP
npe = (long)omp_get_max_threads();
@@ -268,7 +275,7 @@ long *ixpos, long *iypos, long npos, double *tfft, long *isxcount, long *reci_xs
/* place result at source position ixsrc; dx = receiver distance */
for (j = 0; j < nts; j++)
- iRN[l*size+isrc*nts+j] += rtrace[j]*scl*dx*dy;
+ iRN[l*size+isrc*nts+j] += rtrace[j]*scl;
} /* end of Nfoc loop */
diff --git a/utils/getFileInfo3D.c b/utils/getFileInfo3D.c
index d914fc776ef8bba6003480715f7ad4893014f943..eef408c21a9af9892c1cceb9550e78fdd382144a 100644
--- a/utils/getFileInfo3D.c
+++ b/utils/getFileInfo3D.c
@@ -48,6 +48,10 @@ long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, fl
if (ret<0) perror("fseeko");
bytes = ftello( fp );
+ if (hdr.scalco < 0) scl = 1.0/fabs((float)hdr.scalco);
+ else if (hdr.scalco == 0) scl = 1.0;
+ else scl = hdr.scalco;
+
*n1 = hdr.ns;
if ( (hdr.trid == 1) && (hdr.dt != 0) ) {
*d1 = ((float) hdr.dt)*1.e-6;
@@ -58,16 +62,12 @@ long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, fl
*f1 = hdr.f1;
}
*f2 = hdr.f2;
- *f3 = hdr.gy;
+ *f3 = hdr.gy*scl;
data_sz = sizeof(float)*(*n1);
trace_sz = sizeof(float)*(*n1)+TRCBYTES;
ntraces = (long) (bytes/trace_sz);
- if (hdr.scalco < 0) scl = 1.0/fabs((float)hdr.scalco);
- else if (hdr.scalco == 0) scl = 1.0;
- else scl = hdr.scalco;
-
*sclsxgxsygy = scl;
/* check to find out number of traces in shot gather */
diff --git a/utils/makeR1D.c b/utils/makeR1D.c
index acfcaba5965794cc80322820594c7e284c64f23a..53cb1b4ed9a9a8aa517d13955fd918fd878f4869 100644
--- a/utils/makeR1D.c
+++ b/utils/makeR1D.c
@@ -147,7 +147,7 @@ int main (int argc, char **argv)
nysrc = ((nyin+1)/2-1)/dysf+1;
if (verbose) vmess("No nysrc given, setting it to maximum possibility (%li)",((nyin+1)/2-1)/dysf+1);
}
- else if (nysrc>(((nyin+1)/2-1)/dysf+1)) verr("nxsrc given (%li) is higher than the maximum (%li)",nysrc,((nyin+1)/2-1)/dysf+1);
+ else if (nysrc>(((nyin+1)/2-1)/dysf+1)) verr("nysrc given (%li) is higher than the maximum (%li)",nysrc,((nyin+1)/2-1)/dysf+1);
if (nxrcv==-1) {
nxrcv = ((nxin+1)/2-1)/dxrf+1;
if (verbose) vmess("No nxrcv given, setting it to maximum possibility (%li)",((nxin+1)/2-1)/dxrf+1);
@@ -155,12 +155,12 @@ int main (int argc, char **argv)
else if (nxrcv>(((nxin+1)/2-1)/dxrf+1)) verr("nxrcv given (%li) is higher than the maximum (%li)",nxrcv,((nxin+1)/2-1)/dxrf+1);
if (nyrcv==-1) {
nyrcv = ((nyin+1)/2-1)/dyrf+1;
- if (verbose) vmess("No nysrc given, setting it to maximum possibility (%li)",((nyin+1)/2-1)/dyrf+1);
+ if (verbose) vmess("No nyrcv given, setting it to maximum possibility (%li)",((nyin+1)/2-1)/dyrf+1);
}
- else if (nyrcv>(((nyin+1)/2-1)/dyrf+1)) verr("nrcv given (%li) is higher than the maximum (%li)",nyrcv,((nyin+1)/2-1)/dyrf+1);
+ else if (nyrcv>(((nyin+1)/2-1)/dyrf+1)) verr("nyrcv given (%li) is higher than the maximum (%li)",nyrcv,((nyin+1)/2-1)/dyrf+1);
t0out = t0in;
- x0rcv = x0in+((nxrcv-1)/2)*dxrcv; x0src = x0in+((nxsrc-1)/2)*dxsrc;
- y0rcv = y0in+((nyrcv-1)/2)*dyrcv; y0src = y0in+((nysrc-1)/2)*dysrc;
+ x0rcv = x0in+((nxin-1)/2)*dxin-((nxrcv-1)/2)*dxrcv; x0src = x0in+((nxin-1)/2)*dxin-((nxsrc-1)/2)*dxsrc;
+ y0rcv = y0in+((nyin-1)/2)*dyin-((nyrcv-1)/2)*dyrcv; y0src = y0in+((nyin-1)/2)*dyin-((nysrc-1)/2)*dysrc;
x0out = MINL(x0rcv,x0src);
y0out = MINL(y0rcv,y0src);
if (verbose) {