diff --git a/.gitignore b/.gitignore
index e99b1a90490e08865aad83e059ecdd17e03bf3e3..8049b95aea6d5a25f23e8f92fbcff25aee2d62f6 100644
--- a/.gitignore
+++ b/.gitignore
@@ -23,6 +23,12 @@ utils/makewave
utils/mat2su
utils/syn2d
utils/green3D
+utils/HomG
+utils/MuteSnap
+utils/Snap2Shot
+utils/combine
+utils/combine_induced
+utils/reshape_su
fdelmodc/fdelmodc
include/genfft.h
Make_include
@@ -48,4 +54,7 @@ marchenko_applications/iba
marchenko_applications/marchenko_app
marchenko_applications/reshape_su
raytime3d/raytime3d
+raytime/Raytime
MDD/mdd
+marchenko3D/fmute3D
+marchenko3D/marchenko3D
diff --git a/marchenko/demo/WS15/README.3 b/marchenko/demo/WS15/README.3
index 7cbb22a48af0f894044bdb300197189d287fc852..68509bd14999e8e2100dcbf3be0860de84f1e731 100644
--- a/marchenko/demo/WS15/README.3
+++ b/marchenko/demo/WS15/README.3
@@ -1,9 +1,13 @@
Plane waves
-Marchenko emthod worls as well with plane-waves, in stead of modeling a foccus point response for the initial operaror and plane wave is modeled.
+Marchenko method works as well with plane-waves, in stead of modeling a foccus point response for the initial operaror and plane wave is modeled.
Submit the initialPlane.scr to the queue to compute the intial operator. The marchenkoPlane.scr script computes the Macrehnko solution for the plane-wave response.
-
-
The script epsPlane.scr generates .eps and .png files that can be copied back to your local laptop for display.
+
+DISCUSSION:
+In epsPlane.scr there is a convolution between f1minPlane and p0lpus;
+ - why is this 'needed'?
+ - does the result have the correct amplitude?
+ - would this also work with focal points?
diff --git a/marchenko/demo/WS15/README.4 b/marchenko/demo/WS15/README.4
index 071e902dc11e406408ab6c18cc9c437f6cc3456b..5c9553936f58f0dfd245be000fe63bdbe18f2232 100644
--- a/marchenko/demo/WS15/README.4
+++ b/marchenko/demo/WS15/README.4
@@ -3,3 +3,9 @@ To generate a shot record with Primaries reflections only submit the demo/twoD/p
The script epsPrimaries.scr generates .eps and .png files that can be copied back to your local laptop for display.
+The directory demo/invisible contains the Reflection data to analyse:
+
+shotsdx4_rp.su
+
+Try to find out, by ONLY using the program marchenko_primaries, how many reflectors there are.
+
diff --git a/marchenko/demo/WS15/README.5 b/marchenko/demo/WS15/README.5
index ff4310971b52a9e4fde1a1fb2127886f84c2f780..f15516b746814f8cbd28c0ae16dd0ff2a9c08927 100644
--- a/marchenko/demo/WS15/README.5
+++ b/marchenko/demo/WS15/README.5
@@ -1,4 +1,4 @@
-To generate the raytimes and the estimated amplitudes for the virtual source and receivers, run:
+To generate the raytimes and the estimated amplitudes for the virtual source and receivers, run in demo/twoD
rayvsp.scr
diff --git a/marchenko/demo/twoD/epsPlane.scr b/marchenko/demo/twoD/epsPlane.scr
index e5abc5efff2f269ae36c4181a595c0528c0cf8f6..7f2e0be40b19e9a80cbd5e287460485dd641b94d 100755
--- a/marchenko/demo/twoD/epsPlane.scr
+++ b/marchenko/demo/twoD/epsPlane.scr
@@ -1,34 +1,53 @@
#!/bin/bash
#intial plane operator
-file=iniPlane_z1100_rp.su
+file=iniPlane_z1800_rp.su
file_base=${file%.su}
sumax < $file mode=abs outpar=nep
clipref=`cat nep | awk '{print $1/2}'`
supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
label1="time (s)" label2="lateral distance (m)" \
- n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 \
+ n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 title="${file_base}" \
f2=-3000 f2num=-3000 d2num=1000 clip=$clipref > ${file_base}.eps
convert -quality 90 -antialias ${file_base}.eps ${file_base}.jpg
supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < p0plus.su \
label1="time (s)" label2="lateral distance (m)" \
- n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 \
+ n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 title="${file_base}_muted" \
f2=-3000 f2num=-3000 d2num=1000 clip=$clipref > ${file_base}_muted.eps
-convert ${file_base}_muted.eps ${file_base}_muted.png
+convert ${file_base}_muted.eps ${file_base}_muted.jpg
-#Marchencko computed Greens function
-file=f2Plane.su
+#Zero-offset section
+suwind key=offset min=0 max=0 < shots/refl_rp.su > zo.su
+makewave fp=30 dt=0.004 file_out=wave4.su nt=1024 t0=0.0 scale=1
+fconv file_in1=zo.su file_in2=wave4.su file_out=zow.su
+
+file=zow.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/4}'`
+
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
+ label1="time (s)" label2="lateral distance (m)" \
+ n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 title="${file_base}" \
+ f2=-3000 f2num=-3000 d2num=1000 clip=$clipref > ${file_base}.eps
+convert ${file_base}.eps ${file_base}.jpg
+
+#Marchencko computed result based on f1min
+fconv file_in2=p0plus.su file_in1=f1minPlane.su file_out=convPlane.su verbose=3 mode=conv
+fconv file_in1=convPlane.su file_in2=wave4.su file_out=convPlanew.su
+
+file=convPlanew.su
file_base=${file%.su}
sumax < $file mode=abs outpar=nep
clipref=`cat nep | awk '{print $1/4}'`
supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
label1="time (s)" label2="lateral distance (m)" \
- n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 \
+ n1tic=2 x1beg=0 x1end=2.0 d1num=0.4 title="${file_base}" \
f2=-3000 f2num=-3000 d2num=1000 clip=$clipref > ${file_base}.eps
-convert ${file_base}.eps ${file_base}.png
+convert ${file_base}.eps ${file_base}.jpg
diff --git a/marchenko/demo/twoD/initialPlane.scr b/marchenko/demo/twoD/initialPlane.scr
index fbdd55e1fe68ffa5a729e4eff4ddbc8178d852d4..84fbd947bf812049d4ea4205ca88a599c5379c6e 100755
--- a/marchenko/demo/twoD/initialPlane.scr
+++ b/marchenko/demo/twoD/initialPlane.scr
@@ -13,7 +13,8 @@ makemod sizex=6000 sizez=2000 dx=$dx dz=$dx cp0=1900 ro0=1200 \
intt=def x=-3000,3000 z=780,780 poly=0 cp=2230 ro=1700 \
intt=def x=-3000,-2200,-1500,0,1300,2100,3000 z=520,580,680,840,680,600,500 poly=2 cp=2400 ro=2800 \
-makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+#makewave fp=30 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3 scale=0
dxshot=10
ishot=300
@@ -22,12 +23,12 @@ nshots=301
export KMP_AFFINITY=disabled
export OMP_NUM_THREADS=16
-file_rcv=iniPlane_z1100.su
+file_rcv=iniPlane_z1800.su
fdelmodc \
file_cp=synclDown_cp.su ischeme=1 iorder=4 \
file_den=synclDown_ro.su \
- file_src=wave.su \
+ file_src=wavefw.su \
file_rcv=$file_rcv \
src_type=1 \
src_orient=1 \
@@ -42,7 +43,7 @@ fdelmodc \
dxrcv=10.0 \
xrcv1=-3000 xrcv2=3000 \
zrcv1=0 zrcv2=0 \
- xsrc=0 zsrc=1100 plane_wave=1 nsrc=2401 \
+ xsrc=0 zsrc=1800 plane_wave=1 nsrc=2401 \
ntaper=200 \
left=2 right=2 top=2 bottom=2
diff --git a/marchenko/demo/twoD/marchenkoPlane.scr b/marchenko/demo/twoD/marchenkoPlane.scr
index 379837415576fac7d036956b8697017c1c4f56dc..2a16df8cb97758a733c44af6d32636e965a11f46 100755
--- a/marchenko/demo/twoD/marchenkoPlane.scr
+++ b/marchenko/demo/twoD/marchenkoPlane.scr
@@ -5,11 +5,10 @@ export PATH=$HOME/src/OpenSource/bin:$PATH:
export OMP_NUM_THREADS=40
#mute all events below the first arrival to get the intial focusing field
-fmute file_shot=iniPlane_z1100_rp.su file_out=p0plus.su above=-1 shift=-10 verbose=1 check=1 hw=4
+fmute file_shot=iniPlane_z1800_rp.su file_out=p0plus.su above=-1 shift=-10 verbose=1 check=1 hw=4
#apply the Marchenko algorithm
marchenko file_shot=shots/refl_rp.su file_tinv=p0plus.su nshots=601 verbose=1 \
- tap=0 niter=15 hw=8 shift=7 smooth=3 \
- file_green=pPlane.su \
- file_f1plus=f1plusPlane.su file_f1min=f1minPlane.su file_f2=f2Plane.su
+ tap=0 niter=15 hw=8 shift=7 smooth=3 rotate=0 \
+ file_f1plus=f1plusPlane.su file_f1min=f1minPlane.su file_green=greenPlane.su
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index d7d95c0f51b698d42421ebf983ac1604e6e6e8f5..75031101de7ecc97285f06535a2407c00eab9137 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -105,6 +105,7 @@ char *sdoc[] = {
" file_pplus= .............. output file with p+ ",
" file_pmin= ............... output file with p- ",
" file_iter= ............... output file with -Ni(-t) for each iteration",
+" rotate=1 ................. 1: t=0 at nt/2 (middle) 0: t=0 at sample 0 for f1+,-",
" verbose=0 ................ silent option; >0 displays info",
" ",
" ",
@@ -118,7 +119,7 @@ int main (int argc, char **argv)
FILE *fp_out, *fp_f1plus, *fp_f1min;
FILE *fp_gmin, *fp_gplus, *fp_f2, *fp_pmin;
int i, j, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath;
- int size, n1, n2, ntap, tap, di, ntraces, pad;
+ int size, n1, n2, ntap, tap, di, ntraces, pad, rotate;
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
int reci, countmin, mode, n2out, verbose, ntfft;
int iter, niter, tracf, *muteW, *tsynW;
@@ -170,6 +171,7 @@ int main (int argc, char **argv)
if (!getparint("tap", &tap)) tap = 0;
if (!getparint("ntap", &ntap)) ntap = 0;
if (!getparint("pad", &pad)) pad = 0;
+ if (!getparint("rotate", &rotate)) rotate = 1;
if(!getparint("niter", &niter)) niter = 10;
if(!getparint("hw", &hw)) hw = 15;
@@ -699,14 +701,16 @@ sqrt(energyNi/energyN0));
}
if (file_f1plus != NULL) {
/* rotate to get t=0 in the middle */
- for (i = 0; i < n2; i++) {
- hdrs_out[i].f1 = -n1*0.5*dt;
- memcpy(&trace[0],&f1plus[l*size+i*nts],nts*sizeof(float));
- for (j = 0; j < n1/2; j++) {
- f1plus[l*size+i*nts+n1/2+j] = trace[j];
- }
- for (j = n1/2; j < n1; j++) {
- f1plus[l*size+i*nts+j-n1/2] = trace[j];
+ if (rotate==1) {
+ for (i = 0; i < n2; i++) {
+ hdrs_out[i].f1 = -n1*0.5*dt;
+ memcpy(&trace[0],&f1plus[l*size+i*nts],nts*sizeof(float));
+ for (j = 0; j < n1/2; j++) {
+ f1plus[l*size+i*nts+n1/2+j] = trace[j];
+ }
+ for (j = n1/2; j < n1; j++) {
+ f1plus[l*size+i*nts+j-n1/2] = trace[j];
+ }
}
}
ret = writeData(fp_f1plus, (float *)&f1plus[l*size], hdrs_out, n1, n2);
@@ -714,14 +718,16 @@ sqrt(energyNi/energyN0));
}
if (file_f1min != NULL) {
/* rotate to get t=0 in the middle */
- for (i = 0; i < n2; i++) {
- hdrs_out[i].f1 = -n1*0.5*dt;
- memcpy(&trace[0],&f1min[l*size+i*nts],nts*sizeof(float));
- for (j = 0; j < n1/2; j++) {
- f1min[l*size+i*nts+n1/2+j] = trace[j];
- }
- for (j = n1/2; j < n1; j++) {
- f1min[l*size+i*nts+j-n1/2] = trace[j];
+ if (rotate==1) {
+ for (i = 0; i < n2; i++) {
+ hdrs_out[i].f1 = -n1*0.5*dt;
+ memcpy(&trace[0],&f1min[l*size+i*nts],nts*sizeof(float));
+ for (j = 0; j < n1/2; j++) {
+ f1min[l*size+i*nts+n1/2+j] = trace[j];
+ }
+ for (j = n1/2; j < n1; j++) {
+ f1min[l*size+i*nts+j-n1/2] = trace[j];
+ }
}
}
ret = writeData(fp_f1min, (float *)&f1min[l*size], hdrs_out, n1, n2);
diff --git a/marchenko/synthesis.c b/marchenko/synthesis.c
index 2118510d70aac94a0bd077cbab3b499152ce1ddb..a1197a2a8022aff389aa6f6b2e7e14abfb423b99 100644
--- a/marchenko/synthesis.c
+++ b/marchenko/synthesis.c
@@ -173,11 +173,11 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
for (j = 0; j < nts; j++)
iRN[l*size+ixsrc*nts+j] += rtrace[j]*scl*dx;
- } /* end of parallel Nfoc loop */
+ } /* end of Nfoc loop */
if (verbose>4) vmess("*** Shot gather %d processed ***", k);
- } /* end of nshots (k) loop */
+ } /* end of nparallel shots (k) loop */
free(sum);
free(rtrace);