diff --git a/marchenko/demo/oneD/epsBackprop.scr b/marchenko/demo/oneD/epsBackprop.scr
index a4b4f17fa8a32c90bf368c8dd6d867f2e0485582..dff9ca0800cf58fe694cee59db162327891de6e6 100755
--- a/marchenko/demo/oneD/epsBackprop.scr
+++ b/marchenko/demo/oneD/epsBackprop.scr
@@ -2,6 +2,22 @@
export PATH=$HOME/src/OpenSource/bin/:$PATH:
+# Add interface line to postscript file of model
+cat << EOF1 > line1
+400 -2500
+400 2500
+EOF1
+
+cat << EOF2 > line2
+700 -2500
+700 2500
+EOF2
+
+cat << EOF3 > line3
+1100 -2500
+1100 2500
+EOF3
+
dx=2.5
dt=0.0005
diff --git a/marchenko/demo/oneD/epsMarchenkoIter.scr b/marchenko/demo/oneD/epsMarchenkoIter.scr
index 0c795ede15d626bfa83df06867c64e3f15ea7524..919d883d31df9cd27322b5046223ec4e7609a35d 100755
--- a/marchenko/demo/oneD/epsMarchenkoIter.scr
+++ b/marchenko/demo/oneD/epsMarchenkoIter.scr
@@ -34,7 +34,13 @@ fmute file_shot=iniFocus_rp.su file_out=nep.su above=0 shift=8 verbose=1 check=1
#set same clip factor for iteration updates
file=iter_001.su
sumax < $file mode=abs outpar=nep
-clipiter=`cat nep | awk '{print $1/6}'`
+clipiter=`cat nep | awk '{print $1/8}'`
+
+#set same clip factor for Green;s function updates
+file=pgreen_004.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+clipgreen=`cat nep | awk '{print $1/4}'`
#iterations
for (( iter=1; iter<=4; iter+=1 ))
@@ -57,7 +63,7 @@ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
file=f1min_$piter.su
file_base=${file%.su}
sumax < $file mode=abs outpar=nep
-clipref=`cat nep | awk '{print $1/4}'`
+clipref=`cat nep | awk '{print $1/5}'`
supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file\
n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
@@ -65,7 +71,7 @@ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file\
file=f1plus_$piter.su
file_base=${file%.su}
sumax < $file mode=abs outpar=nep
-clipref=`cat nep | awk '{print $1/4}'`
+clipref=`cat nep | awk '{print $1/5}'`
supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file\
n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
@@ -73,28 +79,29 @@ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file\
file=pgreen_$piter.su
file_base=${file%.su}
sumax < $file mode=abs outpar=nep
-clipref=`cat nep | awk '{print $1/3}'`
+clipref=`cat nep | awk '{print $1/4}'`
suwind key=gx min=-2250000 max=2250000 < $file | \
supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
- f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipgreen > $file_base.eps
#compare Green's funtions on Marhcenko and reference result
suwind key=gx min=0 max=0 itmax=511 < pgreen_$piter.su | sumax mode=abs outpar=nepmg
suwind key=gx min=0 max=0 itmax=511 < referenceP_rp.su | sumax mode=abs outpar=neprf
-mg=`cat nepmg | awk '{print $1}'`
-rf=`cat neprf | awk '{print $1}'`
+mg1=`cat nepmg | awk '{print $1}'`
+rf1=`cat neprf | awk '{print $1}'`
value=${value/[eE][+][0]/*10^}
-mg=${mg/[eE][+][0]/*10^}
-rf=${rf/[eE][+][0]/*10^}
+mg=${mg1/[eE][+][0]/*10^}
+rf=${rf1/[eE][+][0]/*10^}
rm nep*
scale=$(echo "scale=3; ($rf)/($mg)" | bc -l)
-#echo $scale
+scale=2.0
+echo $scale
(suwind key=gx min=0 max=0 itmax=511 < pgreen_$piter.su | sugain scale=$scale; \
suwind key=gx min=0 max=0 < referenceP_rp.su) | \
- supsgraph hbox=6 wbox=1 labelsize=10 linegray=0.3,0.0 style=seismic \
- lineon=3.5,1.0 lineoff=1.3,0.0 > compare_$piter.eps
+ supsgraph hbox=6 wbox=2 labelsize=10 linegray=0.3,0.0 style=seismic \
+ lineon=3.5,1.0 x2beg=-$rf1 x2end=$rf1 lineoff=1.3,0.0 > compare_$piter.eps
done
diff --git a/marchenko/demo/oneD/epsModel.scr b/marchenko/demo/oneD/epsModel.scr
index 6a31b126a03c9db8ec78547b1fe2abff5d9a419d..dbe3040a94c1d6742e0c56b6018e666c7fea49c4 100755
--- a/marchenko/demo/oneD/epsModel.scr
+++ b/marchenko/demo/oneD/epsModel.scr
@@ -2,7 +2,7 @@
export PATH=$HOME/src/OpenSource/bin/:$PATH:
-#Postscript file of model and shot record at xsrc=0
+# Add interface line to postscript file of model
cat << EOF1 > line1
400 -2500
400 2500
diff --git a/marchenko/demo/oneD/marchenko.scr b/marchenko/demo/oneD/marchenko.scr
index 16502d5cfb9fb64cf5154af5cfd51ca9c27d64f8..d17b5f4dcc3e8d1071635734cfb1c03f9bfeb723 100755
--- a/marchenko/demo/oneD/marchenko.scr
+++ b/marchenko/demo/oneD/marchenko.scr
@@ -4,11 +4,11 @@ export PATH=$HOME/src/OpenSource/bin:$PATH:
export OMP_NUM_THREADS=1
#mute all events below the first arrival to get the intial focusing field
-fmute file_shot=iniFocus_rp.su file_out=p0plus.su above=-1 shift=-8 verbose=1 check=0 hw=4
+fmute file_shot=iniFocus_rp.su file_out=p0plus.su above=-1 shift=-8 verbose=1 check=0 hw=8
#apply the Marchenko algorithm
marchenko file_shot=shotsdx5_rp.su file_tinv=p0plus.su nshots=901 verbose=1 \
- tap=0 ntap=41 niter=8 hw=12 shift=8 smooth=5 \
+ tap=0 niter=8 hw=8 shift=12 smooth=3 \
file_green=pgreen.su file_gplus=Gplus0.su file_gmin=Gmin0.su \
file_f1plus=f1plus0.su file_f1min=f1min0.su file_f2=f2.su
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index dc31009bcdc5494e9885111662e29dd84313bfec..bc36df906ffe04198df49401dda991ae19480cf6 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -100,7 +100,7 @@ int main (int argc, char **argv)
int iter, niter, tracf, *muteW;
int hw, smooth, above, shift, *ixpossyn, npossyn, ix;
float fmin, fmax, *tapersh, *tapersy, fxf, dxf, fxs2, *xsrc, *xrcv, *zsyn, *zsrc, *xrcvsyn;
- double t0, t1, t2, t3, tsyn, tread, tfft, tcopy;
+ double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi;
float d1, d2, f1, f2, fxs, ft, fx, *xsyn, dxsrc;
float *green, *f2p, *pmin, *G_d, dt, dx, dxs, scl, mem;
float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus;
@@ -279,7 +279,7 @@ int main (int argc, char **argv)
/* check consistency of header values */
fxf = xsrc[0];
- if (nx > 1) dxf = (xrcv[0] - xrcv[nx-1])/(float)(nx-1);
+ if (nx > 1) dxf = (xrcv[nx-1] - xrcv[0])/(float)(nx-1);
else dxf = d2;
if (NINT(dx*1e3) != NINT(fabs(dxf)*1e3)) {
vmess("dx in hdr.d1 (%.3f) and hdr.gx (%.3f) not equal",dx, dxf);
@@ -405,11 +405,15 @@ int main (int argc, char **argv)
j = 0;
Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+j];
pmin[l*nxs*nts+i*nts+j] += iRN[l*nxs*nts+i*nts+j];
+ energyNi = sqrt(iRN[l*nxs*nts+i*nts+j]*iRN[l*nxs*nts+i*nts+j]);
for (j = 1; j < nts; j++) {
Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+nts-j];
pmin[l*nxs*nts+i*nts+j] += iRN[l*nxs*nts+i*nts+j];
+ energyNi = sqrt(iRN[l*nxs*nts+i*nts+j]*iRN[l*nxs*nts+i*nts+j]);
}
}
+ if (verbose >=2) vmess(" - operator %d at iteration %d has energy %e", l, iter, energyNi);
+
}
/* apply mute window based on times of direct arrival (in muteW) */