diff --git a/marchenko3D/Makefile b/marchenko3D/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..237c134e904d7c4b7311e63065c4c09282ace67f
--- /dev/null
+++ b/marchenko3D/Makefile
@@ -0,0 +1,67 @@
+# Makefile
+
+include ../Make_include
+
+LIBS += -L$L -lgenfft -lm $(LIBSM)
+#OPTC += -g -O0 -Wall
+
+#ALL: fmute marchenko marchenko2
+
+ALL: fmute marchenko
+
+SRCJ = fmute.c \
+ getFileInfo.c \
+ readData.c \
+ applyMute.c \
+ writeData.c \
+ wallclock_time.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ getpars.c
+
+SRCH = marchenko.c \
+ getFileInfo.c \
+ readData.c \
+ readShotData.c \
+ readTinvData.c \
+ applyMute.c \
+ writeData.c \
+ writeDataIter.c \
+ wallclock_time.c \
+ name_ext.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ getpars.c
+
+OBJJ = $(SRCJ:%.c=%.o)
+
+fmute: $(OBJJ)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o fmute $(OBJJ) $(LIBS)
+
+OBJH = $(SRCH:%.c=%.o)
+
+marchenko: $(OBJH)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko $(OBJH) $(LIBS)
+
+OBJH2 = $(SRCH2:%.c=%.o)
+
+marchenko2: $(OBJH2)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko2 $(OBJH2) $(LIBS)
+
+install: fmute marchenko
+ cp fmute $B
+ cp marchenko $B
+
+# cp marchenko2 $B
+
+clean:
+ rm -f core fmute $(OBJJ) marchenko $(OBJH) marchenko2 $(OBJH2)
+
+realclean: clean
+ rm -f $B/fmute $B/marchenko $B/marchenko2
+
+
+
+
diff --git a/marchenko3D/applyMute.c b/marchenko3D/applyMute.c
new file mode 100644
index 0000000000000000000000000000000000000000..df98e8bc39f8b132c10b237f7e3e01f3167cdf7d
--- /dev/null
+++ b/marchenko3D/applyMute.c
@@ -0,0 +1,115 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <assert.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *ixpos, int npos, int shift)
+{
+ int i, j, l, isyn;
+ float *costaper, scl;
+ int imute, tmute;
+
+ if (smooth) {
+ costaper = (float *)malloc(smooth*sizeof(float));
+ scl = M_PI/((float)smooth);
+ for (i=0; i<smooth; i++) {
+ costaper[i] = 0.5*(1.0+cos((i+1)*scl));
+ }
+ }
+
+ for (isyn = 0; isyn < Nfoc; isyn++) {
+ if (above==1) {
+ for (i = 0; i < npos; i++) {
+ imute = ixpos[i];
+ tmute = mute[isyn*nxs+imute];
+ for (j = 0; j < MAX(0,tmute-shift-smooth); j++) {
+ data[isyn*nxs*nt+i*nt+j] = 0.0;
+ }
+ for (j = MAX(0,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
+ }
+ }
+ }
+ else if (above==0){
+ for (i = 0; i < npos; i++) {
+ imute = ixpos[i];
+ tmute = mute[isyn*nxs+imute];
+ if (tmute >= nt/2) {
+ memset(&data[isyn*nxs*nt+i*nt],0, sizeof(float)*nt);
+ continue;
+ }
+ for (j = MAX(0,tmute-shift),l=0; j < MAX(0,tmute-shift+smooth); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[l];
+ }
+ for (j = MAX(0,tmute-shift+smooth)+1; j < MIN(nt,nt+1-tmute+shift-smooth); j++) {
+ data[isyn*nxs*nt+i*nt+j] = 0.0;
+ }
+ for (j = MIN(nt,nt-tmute+shift-smooth),l=0; j < MIN(nt,nt-tmute+shift); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
+ }
+ }
+ }
+ else if (above==-1){
+ for (i = 0; i < npos; i++) {
+ imute = ixpos[i];
+ tmute = mute[isyn*nxs+imute];
+ for (j = MAX(0,tmute-shift),l=0; j < MAX(0,tmute-shift+smooth); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[l];
+ }
+ for (j = MAX(0,tmute-shift+smooth); j < nt; j++) {
+ data[isyn*nxs*nt+i*nt+j] = 0.0;
+ }
+ }
+ }
+ else if (above==4) { //Psi gate which is the inverse of the Theta gate (above=0)
+ for (i = 0; i < npos; i++) {
+ imute = ixpos[i];
+ tmute = mute[isyn*nxs+imute];
+ for (j = MAX(0,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
+ }
+ for (j = 0; j < MAX(0,tmute-shift-smooth-1); j++) {
+ data[isyn*nxs*nt+i*nt+j] = 0.0;
+ }
+ for (j = MIN(nt,nt+1-tmute+shift+smooth); j < nt; j++) {
+ data[isyn*nxs*nt+i*nt+j] = 0.0;
+ }
+ for (j = MIN(nt,nt-tmute+shift),l=0; j < MIN(nt,nt-tmute+shift+smooth); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[l];
+ }
+ }
+ }
+ else if (above==2){//Separates the direct part of the wavefield from the coda
+ for (i = 0; i < npos; i++) {
+ imute = ixpos[i];
+ tmute = mute[isyn*nxs+imute];
+ for (j = 0; j < MAX(0,tmute-shift-smooth); j++) {
+ data[isyn*nxs*nt+i*nt+j] = 0.0;
+ }
+ for (j = MAX(0,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
+ }
+ for (j = MAX(0,tmute+shift),l=0; j < MAX(0,tmute+shift+smooth); j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[l];
+ }
+ for (j = MAX(0,tmute+shift+smooth); j < nt; j++) {
+ data[isyn*nxs*nt+i*nt+j] = 0.0;
+ }
+ }
+ }
+ }
+
+ if (smooth) free(costaper);
+
+ return;
+}
+
diff --git a/marchenko3D/atopkge.c b/marchenko3D/atopkge.c
new file mode 120000
index 0000000000000000000000000000000000000000..5107e2b2ccd382ede29d397838d8fad88126a516
--- /dev/null
+++ b/marchenko3D/atopkge.c
@@ -0,0 +1 @@
+../utils/atopkge.c
\ No newline at end of file
diff --git a/marchenko3D/demo/README b/marchenko3D/demo/README
new file mode 100644
index 0000000000000000000000000000000000000000..f5a7c129e2a168d2325d162c9f886a2cc27ae43b
--- /dev/null
+++ b/marchenko3D/demo/README
@@ -0,0 +1,4 @@
+The scripts to reproduce the Figures in the manuscript can be found in the directory oneD. The oneD/README explains how to run the scripts.
+
+A more complicated model can be found in the directory twoD and will takes several hours to model the reflection data.
+
diff --git a/marchenko3D/demo/ScientificReports/NatureSnapshots.tex b/marchenko3D/demo/ScientificReports/NatureSnapshots.tex
new file mode 100644
index 0000000000000000000000000000000000000000..dc80d8cbb6bf25a0446c6f1be7e2cb4eeec777dd
--- /dev/null
+++ b/marchenko3D/demo/ScientificReports/NatureSnapshots.tex
@@ -0,0 +1,194 @@
+\documentclass[10pt]{article}
+\usepackage{graphpap, epsfig, amssymb, rotating, graphpap, pifont}
+\usepackage{amsmath}
+\usepackage{pstricks,placeins,pstricks-add}
+\usepackage{hyperref}
+%\usepackage{pst-node, pst-plot}
+\usepackage{afterpage}
+\usepackage{psfrag}
+\usepackage[text={16cm,24cm},centering]{geometry}
+
+\pagenumbering{gobble}% Remove page numbers (and reset to 1)
+
+% don't complain about overfull vboxes less than 5 points
+\vfuzz = 5pt
+\sloppy
+
+% now define realbold and undertilde
+\def\realbold#1{{\mbox{\boldmath $\bf #1$}}}
+\def\realboldcal#1{{\mbox{\boldmath $\bf \cal #1$}}}
+\def\@undertilde#1{\oalign{{\realbold{#1}}\crcr\hidewidth
+\vbox to .2ex{\hbox{\bf\char126}\vss}\hidewidth}}
+
+%Define helvetica bold for discrete vectors and matrices
+\DeclareMathAlphabet{\mathhelb}{OT1}{phv}{b}{n}
+\DeclareMathAlphabet{\mathhelbo}{OT1}{phv}{b}{sl}
+
+\def\vector#1{{\mbox{\boldmath $#1$}}}%
+\def\tensor#1{\realbold{#1}}%
+\def\op#1{\mathcal{#1}}%
+
+\def\dvector#1{\mathhelbo{#1}}%
+\def\dtensor#1{\mathhelb{#1}}%
+\def\dop#1{\mathhelb{#1}}%
+
+\let\BM\realbold
+
+\definecolor{LightBlue}{rgb}{0.10,0.55,1.0}
+\definecolor{LightBlue}{rgb}{0.0,1.0,1.0}
+
+% set unitlength to 1cm for picture environments
+\setlength{\unitlength}{1cm}
+
+% spacing between lines
+\renewcommand{\baselinestretch}{1.08}
+
+%Jan's definitions
+\newcommand{\dd}{\mathrm{d}}
+\newcommand{\x}{\vector x}
+\newcommand{\dA}{\dd^2\x}
+\newcommand{\dV}{\dd^3\x}
+\newcommand{\pD}{\partial D}
+\newcommand{\bx}{\mathbf{x}}
+
+
+\begin{document}
+
+\title{Note: Snapshots for article}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%% Figure with different time axis annotation
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\begin{figure}[ht]
+\begin{pspicture}(16,22)(0,0)
+%time axis top row
+%\put(-1.2, 17.5){$+t$}
+%\put(-1.2, 16.1){$-t$}
+\put(-1.5, 16.88){$t=0$}
+\psline{->}(-1.7,15.8)(-1.7,17.8) \rput{L}(-1.9, 16.6){ time}
+% time axis gray-scales
+\put(-1.2, 14.0){$-t_2$}
+\put(-1.2, 11.5){$-t_1$}
+\put(-1.5, 9.0){$t=0$}
+\put(-1.2, 6.5){$+t_1$}
+\put(-1.2, 4.0){$+t_2$}
+\psline{->}(-1.7,14.5)(-1.7,3.5) \rput{L}(-1.9, 8.8){ time}
+%
+%frame wiggle traces
+\psframe(-0.5,15.60)(3.45,18.2)
+\psframe(3.6,15.60)(7.55,18.2)
+\psframe(7.7,15.60)(11.65,18.2)
+%\psframe(11.8,15.60)(15.75,18.2)
+%
+% blue background on wiggles
+\pscurve[linecolor=LightBlue,linewidth=4pt]{c-c}(-0.4,16.3)(1.5,16.55)(3.4,16.3)
+\pscurve[linecolor=LightBlue,linewidth=4pt]{c-c}(3.7,16.3)(5.6,16.55)(7.5,16.3)
+\pscurve[linecolor=LightBlue,linewidth=4pt]{c-c}(7.8,16.3)(9.7,16.55)(11.6,16.3)
+%\pscurve[linecolor=LightBlue,linewidth=4pt]{c-c}(11.9,16.3)(13.8,16.55)(15.7,16.3)
+%
+\put(-0.57, 14.8){\epsfig{file={inj_rate_surf_dx0.5_rvz_wiggle}.eps, width=4.09cm }}
+\put(3.53, 14.8){\epsfig{file={inj_rate_surf_dx0.5_rvz_wiggle}.eps, width=4.09cm }}
+\put(7.63, 14.8){\epsfig{file={pplus_wiggle}.eps, width=4.09cm }}
+%\put(11.73, 14.8){\epsfig{file={pplus_wiggle}.eps, width=4.09cm }}
+%
+\put(-0.5, 12.0){\epsfig{file={snapinj_planevzvxsum_-0.60}.eps, width=3.9cm }}
+\put(3.6, 12.0){\epsfig{file={snapinj_surf_-0.60}.eps, width=3.9cm }}
+\put(7.7, 12.0){\epsfig{file={snapinj_f2_-0.60}.eps, width=3.9cm }}
+\put(11.8, 12.0){\epsfig{file={snapinj_f2sum_0.60}.eps, width=3.9cm }}
+%
+\put(-0.5, 9.5){\epsfig{file={snapinj_planevzvxsum_-0.30}.eps, width=3.9cm }}
+\put(3.6, 9.5){\epsfig{file={snapinj_surf_-0.30}.eps, width=3.9cm }}
+\put(7.7, 9.5){\epsfig{file={snapinj_f2_-0.30}.eps, width=3.9cm }}
+\put(11.8, 9.5){\epsfig{file={snapinj_f2sum_0.30}.eps, width=3.9cm }}
+%
+\put(-0.5, 7.0){\epsfig{file={snapinj_planevzvxsum_0.00}.eps, width=3.9cm }}
+\put(3.6, 7.0){\epsfig{file={snapinj_surf_0.00}.eps, width=3.9cm }}
+\put(7.7, 7.0){\epsfig{file={snapinj_f2_0.00}.eps, width=3.9cm }}
+\put(11.8, 7.0){\epsfig{file={snapinj_f2_0.00}.eps, width=3.9cm }}
+%
+\put(-0.5, 4.5){\epsfig{file={snapinj_planevzvxsum_0.30}.eps, width=3.9cm }}
+\put(3.6, 4.5){\epsfig{file={snapinj_surf_0.30}.eps, width=3.9cm }}
+\put(7.7, 4.5){\epsfig{file={snapinj_f2_0.30}.eps, width=3.9cm }}
+\put(11.8, 4.5){\epsfig{file={snapinj_f2sum_0.30}.eps, width=3.9cm }}
+%
+\put(-0.5, 2.0){\epsfig{file={snapinj_planevzvxsum_0.60}.eps, width=3.9cm }}
+\put(3.6, 2.0){\epsfig{file={snapinj_surf_0.60}.eps, width=3.9cm }}
+\put(7.7, 2.0){\epsfig{file={snapinj_f2_0.60}.eps, width=3.9cm }}
+\put(11.8, 2.0){\epsfig{file={snapinj_f2sum_0.60}.eps, width=3.9cm }}
+%
+% source position
+\multirput(1.44,13.75)(4.1,0){4}{\psdot[linecolor=yellow,fillcolor=yellow,dotstyle=o,dotsize=4pt](0,0)\put(0.1,0.0){${\bf s}$}}
+\multirput(1.44,11.25)(4.1,0){4}{\psdot[linecolor=yellow,fillcolor=yellow,dotstyle=o,dotsize=4pt](0,0)\put(0.1,0.0){${\bf s}$}}
+\multirput(1.44,8.75)(4.1,0){4}{\put(0.1,0.0){${\bf s}$}}
+\multirput(1.44,6.25)(4.1,0){4}{\psdot[linecolor=yellow,fillcolor=yellow,dotstyle=o,dotsize=4pt](0,0)\put(0.1,0.0){${\bf s}$}}
+\multirput(1.44,3.75)(4.1,0){4}{\psdot[linecolor=yellow,fillcolor=yellow,dotstyle=o,dotsize=4pt](0,0)\put(0.1,0.0){${\bf s}$}}
+%
+%source function
+\rput(1.44, 17.9){\psframebox[fillstyle=solid]{$V({\bf x},{\bf s},-t)$}}
+\rput(5.54, 17.9){\psframebox[fillstyle=solid]{$V({\bf x},{\bf s},-t)$}}
+\rput(9.64, 17.9){\psframebox[fillstyle=solid]{$F({\bf x},{\bf s},t)$}}
+%\rput(13.74, 17.9){\psframebox[fillstyle=solid]{$F({\bf x},{\bf s},t)$}}
+%
+%times in snapshots
+%\multirput(2.78,13.03)(4.1,0){4}{\put(0.0,0.0){$-t_2$}}
+%\multirput(2.78,10.53)(4.1,0){4}{\put(0.0,0.0){$-t_1$}}
+%\multirput(2.78,8.03)(4.1,0){4}{\put(-0.3,0.0){$t=0$}}
+%\multirput(2.78,5.53)(4.1,0){4}{\put(0.2,0.0){$t_1$}}
+%\multirput(2.78,3.03)(4.1,0){4}{\put(0.2,0.0){$t_2$}}
+%
+% Source postions (red dots) at top snapshots
+\psset{arrowscale=1.1,ArrowFill=true,linecolor=red,linestyle=none}
+\psline[ArrowInside=-*,ArrowInsidePos=0.0,ArrowInsideNo=13]{-}(3.4,15.4)(3.4,12.7)
+\psline[ArrowInside=-*,ArrowInsidePos=0.0,ArrowInsideNo=13]{-}(-0.5,15.4)(-0.5,12.7)
+\psline[ArrowInside=-*,ArrowInsidePos=0.0,ArrowInsideNo=21]{-}(-0.67,15.22)(3.5,15.22)
+\psline[ArrowInside=-*,ArrowInsidePos=0.0,ArrowInsideNo=21]{-}(-0.67,12.80)(3.5,12.80)
+%
+\psline[ArrowInside=-*,ArrowInsidePos=0.0,ArrowInsideNo=21]{-}(3.43,15.22)(7.6,15.22)
+\psline[ArrowInside=-*,ArrowInsidePos=0.0,ArrowInsideNo=21]{-}(7.53,15.22)(11.7,15.22)
+%\psline[ArrowInside=-*,ArrowInsidePos=0.0,ArrowInsideNo=21]{-}(11.63,15.22)(15.8,15.22)
+%
+% source arrows
+\psset{arrowscale=1.0,ArrowFill=true,linecolor=black,linestyle=solid}
+\multirput(-0.37,15.6)(0.184,0){21}{\psline{->}(0.0,0.0)(0.0,-0.36)}
+\multirput(3.73,15.6)(0.184,0){21}{\psline{->}(0.0,0.0)(0.0,-0.36)}
+\multirput(7.83,15.6)(0.184,0){21}{\psline{->}(0.0,0.0)(0.0,-0.36)}
+%\multirput(11.93,15.6)(0.184,0){21}{\psline{->}(0.0,0.0)(0.0,-0.36)}
+%
+%bottom labels
+\put(0.3, 2.0){omnidirectional}
+\put(0.4, 1.6){time reversal}
+\put(1.3, 0.8){(a)}
+\put(4.6, 2.0){single-sided}
+\put(4.5, 1.6){time reversal}
+\put(5.4, 0.8){(b)}
+\put(8.6, 2.0){single-sided}
+\put(8.8, 1.6){focusing}
+\put(9.5, 0.8){(c)}
+\put(12.8, 2.0){symmetrized}
+\put(12.9, 1.6){single-sided}
+\put(13.1, 1.2){focusing}
+\put(13.6, 0.8){(d)}
+%
+% x-axis top row
+\multirput(2.6,18.35)(4.1,0.0){3}{\psline{->}(0.0,0.0)(0.5,0.0) \put(0.5,-0.05){$x$}}
+% x-axis snapshots
+\multirput(2.6,15.05)(0,-2.5){5}{\psline{->}(0.0,0.0)(0.5,0.0) \put(0.5,-0.05){$x$}}
+\multirput(6.7,15.05)(0,-2.5){5}{\psline{->}(0.0,0.0)(0.5,0.0) \put(0.5,-0.05){$x$}}
+\multirput(10.8,15.05)(0,-2.5){5}{\psline{->}(0.0,0.0)(0.5,0.0) \put(0.5,-0.05){$x$}}
+\multirput(14.9,15.05)(0,-2.5){5}{\psline{->}(0.0,0.0)(0.5,0.0) \put(0.5,-0.05){$x$}}
+%
+% z-axis snapshots
+\multirput(-0.3,13.2)(0,-2.5){5}{\psline{->}(0.0,0.0)(0.0,-0.5) \put(-0.1,-0.7){$z$}}
+\multirput(3.8,13.2)(0,-2.5){5}{\psline{->}(0.0,0.0)(0.0,-0.5) \put(-0.1,-0.7){$z$}}
+\multirput(7.9,13.2)(0,-2.5){5}{\psline{->}(0.0,0.0)(0.0,-0.5) \put(-0.1,-0.7){$z$}}
+\multirput(12.0,13.2)(0,-2.5){5}{\psline{->}(0.0,0.0)(0.0,-0.5) \put(-0.1,-0.7){$z$}}
+%
+%\psgrid
+\end{pspicture}
+%\caption{Different... }\label{model}
+\end{figure}
+
+
+
+\end{document}
+
diff --git a/marchenko3D/demo/ScientificReports/README b/marchenko3D/demo/ScientificReports/README
new file mode 100644
index 0000000000000000000000000000000000000000..cbecad73b225273369f72c33d01774720013e790
--- /dev/null
+++ b/marchenko3D/demo/ScientificReports/README
@@ -0,0 +1,61 @@
+The scripts in this directory create Figure 1 of the paper in Scientific Reports
+
+Forward model the data
+1a) model.scr computes the gridded model
+ => tutodel_dx0.5_ro.su, tutodel_dx0.5_cp.su
+1b) shots_slurm.scr creates jobs to model the shots and submit jobs to slurm workload manager
+ => shotsnew/shots_x*_rp.su : ranging from -3000 to 3000 with dxshot=10
+ => to model one shots takes ~25 minutes, 601 shots are modeled
+1c) check.scr after the jobs on shots_*.scr are finished checks if all shots are there
+2) direct.scr creates the direct arrival to be removed from the shots
+ => direct_rp.su
+3) remove_direct.scr remove the direct wave from the shots
+ => shotsnew/refl_rp.su
+4) initialFocus.scr model G_d the intitial focusing function
+ => iniFocusz800x0_rp.su
+
+Apply the marchenko method to compute f2
+5) marchenko.scr perform the Marchenko scheme
+ => f2.su, f1plus.su, f1min.su, Gmin.su, Gplus.su, pgreen.su
+
+Backpropagation of the recorded data from one side and all-sides: snaphots of the wavefield are recorded
+6) back_injrate_planes.scr
+ => snapinj_planevzvxsum_sp.su : backpropagated from all 4 sides of the model
+ => snapinj_surf_sp.su : backpropagated from surface only
+
+Backpropagation of the Marchenko computed result
+7) backpropf2.scr
+ => snapinj_f2_sp.su : backpropagated f2
+
+Generate the postscript files of Figure 1
+8) epsBack.scr
+ => results of backpropagation from all 4 sides: first column of Figure 1
+snapinj_planevzvxsum_-0.60.eps
+snapinj_planevzvxsum_-0.30.eps
+snapinj_planevzvxsum_-0.03.eps
+snapinj_planevzvxsum_0.00.eps
+snapinj_planevzvxsum_0.03.eps
+snapinj_planevzvxsum_0.30.eps
+snapinj_planevzvxsum_0.60.eps
+ => results of backpropagation from surface only: second column of Figure 1
+snapinj_surf_-0.60.eps
+snapinj_surf_-0.30.eps
+snapinj_surf_-0.03.eps
+snapinj_surf_0.00.eps
+snapinj_surf_0.03.eps
+snapinj_surf_0.30.eps
+snapinj_surf_0.60.eps
+ => results of backpropagation of f2 from surface only: third column of Figure 1
+snapinj_f2_-0.60.eps
+snapinj_f2_-0.30.eps
+snapinj_f2_-0.03.eps
+snapinj_f2_0.00.eps
+snapinj_f2_0.03.eps
+snapinj_f2_0.30.eps
+snapinj_f2_0.60.eps
+ => results of symmetrized backpropagation of f2 from surface only: fourth column of Figure 1
+snapinj_f2sum_0.00.eps
+snapinj_f2sum_0.03.eps
+snapinj_f2sum_0.30.eps
+snapinj_f2sum_0.60.eps
+
diff --git a/marchenko3D/demo/ScientificReports/back_injrate_planes.scr b/marchenko3D/demo/ScientificReports/back_injrate_planes.scr
new file mode 100755
index 0000000000000000000000000000000000000000..e5e0ebdb86d52a9334a3f9f36e886aa3caa9e106
--- /dev/null
+++ b/marchenko3D/demo/ScientificReports/back_injrate_planes.scr
@@ -0,0 +1,141 @@
+#!/bin/bash
+
+export PATH=:$HOME/src/OpenSource/bin:$HOME/bin:$PATH:
+
+export OMP_NUM_THREADS=8
+dt=0.00050
+dx=2.5
+
+dt=0.00010
+dx=0.5
+
+makewave fp=30 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+./model.scr $dx
+ntr=$(echo "scale=0; 1+6000/${dx}" | bc -l)
+echo $ntr
+
+ix1a=1
+ix1b=$(echo "scale=0; ${ix1a}+6000/${dx}" | bc -l)
+ix2a=$(echo "scale=0; ${ix1b}+1" | bc -l)
+ix2b=$(echo "scale=0; ${ix2a}+6000/${dx}" | bc -l)
+ix3a=$(echo "scale=0; ${ix2b}+1" | bc -l)
+ix3b=$(echo "scale=0; ${ix3a}+1200/${dx}" | bc -l)
+ix4a=$(echo "scale=0; ${ix3b}+1" | bc -l)
+ix4b=$(echo "scale=0; ${ix4a}+1200/${dx}" | bc -l)
+
+file_cp=tutodel_dx${dx}_cp.su
+file_ro=tutodel_dx${dx}_ro.su
+
+#model data to be propagated back into medium
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=wave.su \
+ file_rcv=inj_rate_plane_dx${dx}.su \
+ src_type=1 \
+ src_injectionrate=1 \
+ rec_type_vz=1 \
+ rec_type_vx=1 \
+ rec_type_p=0 \
+ rec_int_vz=2 \
+ rec_int_vx=2 \
+ dtrcv=$dt \
+ rec_delay=0.1 \
+ verbose=1 \
+ tmod=4.4000 \
+ xrcv1=-3000,-3000,-3000,3000 xrcv2=3000,3000,-3000,3000 zrcv1=0,1200,0,0 zrcv2=0,1200,1200,1200 \
+ dxrcv=$dx,$dx,0,0 dzrcv=0,0,$dx,$dx \
+ xsrc=0 zsrc=800 \
+ ntaper=400 \
+ left=2 right=2 top=2 bottom=2
+
+
+##################################
+# back propagate from all sides
+# scale with -1 for outward pointing vector
+scale=-1.0
+suwind key=tracl min=$ix1a max=$ix1b < inj_rate_plane_dx${dx}_rvz.su | sugain scale=$scale > inj_rate_plane_dx${dx}vz.su
+suwind key=tracl min=$ix2a max=$ix2b < inj_rate_plane_dx${dx}_rvz.su >> inj_rate_plane_dx${dx}vz.su
+
+suwind < inj_rate_plane_dx${dx}_rvx.su key=tracl min=$ix3a max=$ix3b | sugain scale=$scale > inj_rate_plane_dx${dx}vx.su
+suwind < inj_rate_plane_dx${dx}_rvx.su key=tracl min=$ix4a max=$ix4b >> inj_rate_plane_dx${dx}vx.su
+
+# at 4.3000 seconds (tmod - rec_delay) the focus is at t=0
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=inj_rate_plane_dx${dx}vz.su \
+ file_snap=snapinj_planevz.su \
+ grid_dir=1 \
+ src_type=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_p=0 \
+ dtrcv=$dt \
+ rec_delay=0.0 \
+ verbose=3 \
+ tmod=5.004 \
+ tsnap1=3.000 tsnap2=5.00 dtsnap=0.01 \
+ xsnap1=-1000 xsnap2=1000 dxsnap=2.5 dzsnap=2.5 \
+ sna_type_vz=0 \
+ sna_type_p=1 \
+ ntaper=400 \
+ left=2 right=2 top=2 bottom=2
+
+# tsnap1=4.200 tsnap2=4.50 dtsnap=0.004 \
+
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=inj_rate_plane_dx${dx}vx.su \
+ file_snap=snapinj_planevx.su \
+ grid_dir=1 \
+ src_type=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_p=0 \
+ dtrcv=$dt \
+ rec_delay=0.0 \
+ verbose=1 \
+ tmod=5.004 \
+ tsnap1=3.000 tsnap2=5.00 dtsnap=0.01 \
+ xsnap1=-1000 xsnap2=1000 dxsnap=2.5 dzsnap=2.5 \
+ sna_type_vz=0 \
+ sna_type_p=1 \
+ ntaper=400 \
+ left=2 right=2 top=2 bottom=2
+
+# tsnap1=4.200 tsnap2=4.50 dtsnap=0.004 \
+
+suop2 snapinj_planevz_sp.su snapinj_planevx_sp.su op=sum w1=1 w2=1 > snapinj_planevzvxsum_sp.su
+
+##################################
+# back propagate from surface only
+scale=-1.0
+suwind key=tracl min=1 max=$ntr < inj_rate_plane_dx${dx}_rvz.su | sutaper tr1=100 tr2=100 ntr=$ntr | sugain scale=$scale > inj_rate_surf_dx${dx}_rvz.su
+
+# at 4.3000 seconds the focus is at t=0
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=inj_rate_surf_dx${dx}_rvz.su \
+ file_snap=snapinj_surf.su \
+ grid_dir=1 \
+ src_type=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_p=0 \
+ dtrcv=$dt \
+ rec_delay=0.0 \
+ verbose=1 \
+ tmod=5.004 \
+ tsnap1=3.000 tsnap2=5.00 dtsnap=0.01 \
+ xsnap1=-1000 xsnap2=1000 dxsnap=2.5 dzsnap=2.5 \
+ sna_type_vz=0 \
+ sna_type_p=1 \
+ ntaper=400 \
+ left=2 right=2 top=2 bottom=2
+
diff --git a/marchenko3D/demo/ScientificReports/backpropf2.scr b/marchenko3D/demo/ScientificReports/backpropf2.scr
new file mode 100755
index 0000000000000000000000000000000000000000..a19fca25cc61b0157ecdaf86b0cc9c52d3ce726e
--- /dev/null
+++ b/marchenko3D/demo/ScientificReports/backpropf2.scr
@@ -0,0 +1,57 @@
+#!/bin/bash
+
+export PATH=:$HOME/src/OpenSource/bin:$HOME/bin:$PATH:
+
+export OMP_NUM_THREADS=8
+dx=2.5
+dt=0.00050
+dx=0.5
+dt=0.000125
+
+ns=`surange < f2.su | grep ns | awk '{print $2}'`
+dtrcv=`surange < f2.su | grep dt | awk '{print $2/1000000.0}'`
+#suwind key=gx min=-2250000 max=2250000 itmax=1023 < f2.su | sushw key=f2 a=-2250 > nep.su
+#shift=$(echo "scale=6; ($dtrcv*($ns/2.0-1)+0.5*$dt-0.000250)" | bc -l)
+shift=$(echo "scale=6; ($dtrcv*($ns/2.0-1))" | bc -l)
+echo $shift
+basop choice=shift shift=$shift file_in=f2.su verbose=1 > pplus.su
+
+# the f2.su is sampled with 4ms the FD program need 0.5ms
+# time axis is interpolated by making use of FFT's: sinc interpolation
+# this is now done in fdelmodc
+
+#ftr1d file_in=pplus.su file_out=freq.su
+#sushw < freq.su key=nhs,dt a=8192,500 > fr.su
+#ftr1d file_in=fr.su n1=8194 file_out=pplusdt.su verbose=1
+
+file_cp=tutodel_dx${dx}_cp.su
+file_ro=tutodel_dx${dx}_ro.su
+
+# back propagate with f2 from marchenko.scr
+
+set -x
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=pplus.su \
+ file_snap=snapinj_f2.su \
+ grid_dir=0 \
+ src_type=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=0 \
+ rec_int_p=0 \
+ verbose=4 \
+ dt=$dt \
+ tmod=3.0000 \
+ tsnap1=0.744125 tsnap2=2.744125 dtsnap=0.01 \
+ xsnap1=-1000 xsnap2=1000 dxsnap=2.5 dzsnap=2.5 \
+ sna_type_vz=0 \
+ sna_type_p=1 \
+ fmax=120 \
+ ntaper=400 \
+ left=2 right=2 top=2 bottom=2
+
+#2.044125
+# tmod=2.1000 \
+# tsnap1=2.0400 tsnap2=2.0500 dtsnap=0.000125 \
diff --git a/marchenko3D/demo/ScientificReports/check.scr b/marchenko3D/demo/ScientificReports/check.scr
new file mode 100755
index 0000000000000000000000000000000000000000..b369721f40779a714ff6a020a404bacfc846e21c
--- /dev/null
+++ b/marchenko3D/demo/ScientificReports/check.scr
@@ -0,0 +1,29 @@
+#!/bin/bash
+#PBS -N fdelmod
+#PBS -q long
+#PBS -l nodes=1
+#PBS -k eo
+#PBS -j eo
+
+dxshot=10
+ishot=0
+nshots=601
+zsrc=0
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+
+ file_rcv=shotsnew/shots_${xsrc}_rp.su
+
+ if [ ! -e "$file_rcv" ]
+ then
+ echo $xsrc is missing
+ sbatch jobs/slurm_$ishot.job
+ fi
+
+ (( ishot = $ishot + 1))
+
+done
+
diff --git a/marchenko3D/demo/ScientificReports/clean b/marchenko3D/demo/ScientificReports/clean
new file mode 100755
index 0000000000000000000000000000000000000000..26987b43f2604488b03655488d6daaad171eb7ad
--- /dev/null
+++ b/marchenko3D/demo/ScientificReports/clean
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+rm *.su *.bin *.eps nep line* *.asci curve* *.mod
+
diff --git a/marchenko3D/demo/ScientificReports/direct.scr b/marchenko3D/demo/ScientificReports/direct.scr
new file mode 100755
index 0000000000000000000000000000000000000000..fa06b4215d8a8908a7bcbe37ff16177f265d146a
--- /dev/null
+++ b/marchenko3D/demo/ScientificReports/direct.scr
@@ -0,0 +1,40 @@
+#!/bin/bash
+
+cd /vardim/home/thorbcke/data/Kees/Marchenko/Tutorial
+
+dx=2.5
+dt=0.00050
+fast="fast"
+
+dx=0.5
+dt=0.0001
+fast=""
+
+makemod sizex=12000 sizez=250 dx=$dx dz=$dx cp0=1500 ro0=1000 \
+ orig=-6000,0 file_base=noContrast.su
+
+export OMP_NUM_THREADS=8
+makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw${fast}.su nt=8192 t0=0.4 scale=0 scfft=1
+
+fdelmodc \
+ file_cp=noContrast_cp.su ischeme=1 iorder=4 \
+ file_den=noContrast_ro.su \
+ file_src=wavefw${fast}.su \
+ file_rcv=direct${fast}.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ rec_delay=0.4 \
+ dtrcv=0.004 \
+ verbose=2 \
+ tmod=4.492 \
+ dxrcv=10.0 \
+ xrcv1=-6000 xrcv2=6000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=0 \
+ ntaper=200 \
+ left=2 right=2 top=2 bottom=2
+
diff --git a/marchenko3D/demo/ScientificReports/epsBack.scr b/marchenko3D/demo/ScientificReports/epsBack.scr
new file mode 100755
index 0000000000000000000000000000000000000000..fad63523dfa74f535cbaf7dd71847af4d4b3349c
--- /dev/null
+++ b/marchenko3D/demo/ScientificReports/epsBack.scr
@@ -0,0 +1,131 @@
+#!/bin/bash
+
+export PATH=:$HOME/src/OpenSource/bin:$HOME/bin:$PATH:
+
+dx=0.5
+file_snap_all=snapinj_planevzvxsum_sp.su
+file_snap_surf=snapinj_surf_sp.su
+file_snap_f2=snapinj_f2_sp.su
+
+
+sumax < ${file_snap_all} mode=abs outpar=nep
+clip_all=`cat nep | awk '{print $1/25}'`
+clip_t0all=`cat nep | awk '{print $1/10}'`
+sumax < ${file_snap_surf} mode=abs outpar=nep
+clip_surf=`cat nep | awk '{print $1/10}'`
+clip_t0surf=`cat nep | awk '{print $1/4}'`
+sumax < ${file_snap_f2} mode=abs outpar=nep
+clip_f2=`cat nep | awk '{print $1/10}'`
+clip_t0f2=`cat nep | awk '{print $1/4}'`
+echo $clip_all $clip_surf $clip_f2
+
+#131 => t=0
+for file_snap in $file_snap_all $file_snap_surf $file_snap_f2
+do
+ if [ $file_snap == $file_snap_all ]; then clip=$clip_all; fi;
+ if [ $file_snap == $file_snap_surf ]; then clip=$clip_surf; fi;
+ if [ $file_snap == $file_snap_f2 ]; then clip=$clip_f2; fi;
+ echo $clip
+ file_base=${file_snap%_sp.su}
+ suwind key=fldr min=61 max=201 < $file_snap | sugain scale=-1 | sustrip > $file_base.bin
+ for fldr in 71 101 128 131 134 161 191;
+ do
+ if [ $fldr == 131 ];
+ then
+ clipt=$(echo "scale=2; ${clip}*4.0" | bc -l);
+ else
+ clipt=$clip;
+ fi;
+ echo clip=$clip clipt=$clipt
+ times=$(echo "scale=2; 0.01*(${fldr}-131)" | bc -l)
+ atime=`printf "%4.2f" $times`
+ suwind key=fldr min=$fldr max=$fldr < ${file_snap} | \
+ supsimage hbox=2.5 wbox=4 labelsize=-1 \
+ x1beg=0 x1end=1250.0 clip=$clipt \
+ curve=curve1,curve2,curve3 npair=25,25,25 curvecolor=black,black,black curvedash=3,3,3 \
+ n1tic=5 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_base}_$atime.eps
+ done
+done
+
+file_base=${file_snap_f2%_sp.su}
+rm ${file_base}sum.su
+for shot in `seq 61 1 131`;
+do
+ suwind key=fldr min=$shot max=$shot < ${file_snap_f2} > neg.su
+ (( a = 131+(131-$shot) ))
+ echo $shot $a
+ suwind key=fldr min=$a max=$a < ${file_snap_f2} > pos.su
+ susum neg.su pos.su | sushw key=fldr a=$shot >> ${file_base}sum.su
+done
+
+for shot in `seq 132 1 201`;
+do
+ suwind key=fldr min=$shot max=$shot < ${file_snap_f2}> pos.su
+ (( a = 131+(131-$shot) ))
+ echo $shot $a
+ suwind key=fldr min=$a max=$a < ${file_snap_f2} > neg.su
+ susum neg.su pos.su | sushw key=fldr a=$fldr >> ${file_base}sum.su
+done
+
+sugain scale=-1 < ${file_base}sum.su | sustrip > ${file_base}sum.bin
+
+exit
+
+#select files for snapshot between -0.7 => 0 <= +0.07 (fldr 31-101-171)
+#add pos and negative times to get response of homogenoeus Green's function
+
+file_base=${file_snap%_sp.su}
+for fldr in 71 101 128 131;
+do
+ times=$(echo "scale=2; 0.01*(131-${fldr})" | bc -l)
+ atime=`printf "%4.2f" $times`
+ suwind key=fldr min=$fldr max=$fldr < ${file_snap_f2} > neg.su
+ (( fldr = 131+(131-$fldr) ))
+ suwind key=fldr min=$fldr max=$fldr < ${file_snap_f2}> pos.su
+ susum neg.su pos.su | \
+ supsimage hbox=2.5 wbox=4 labelsize=-1 \
+ x1beg=0 x1end=1250.0 clip=$clip_f2 \
+ curve=curve1,curve2,curve3 npair=25,25,25 curvecolor=black,black,black curvedash=3,3,3 \
+ n1tic=5 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_base}sum_$atime.eps
+done
+
+#supsgraph < wave.su hbox=6 wbox=2 labelsize=10 f1=0 x1end=0.4 > wave.eps
+#basop file_in=wave.su choice=15 | supsgraph hbox=6 wbox=2 labelsize=10 f1=0 x1end=0.4 > wavesjom.eps
+
+file=inj_rate_surf_dx${dx}_rvz.su
+file_base=${file%.su}
+sumax < ${file} mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/5}'`
+
+ns=`surange < $file | grep ns | awk '{print $2}'`
+dtrcv=`surange < $file | grep dt | awk '{print $2/1000000.0}'`
+shift=$(echo "scale=6; ($dtrcv*($ns/2.0-1))" | bc -l)
+skip=$(echo "scale=0; (10/${dx})" | bc -l)
+echo $shift $skip
+basop choice=shift shift=$shift file_in=$file verbose=1 | suwind s=1 j=$skip | sushw key=d2,tracl a=10,1 b=0,1 > nep.su
+supsimage hbox=2.5 wbox=4 labelsize=-1 < nep.su \
+ x1beg=-0.6 x1end=1.5 clip=$clip grid1=dot \
+ n1tic=5 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_base}.eps
+
+suwind s=1 j=10 < nep.su | \
+ supswigp hbox=2.5 wbox=4 labelsize=-1 d2=100 \
+ x1beg=-1.5 x1end=1.5 xcur=1.5 grid1=dot n1tic=1 d1num=2.145 \
+ f2=-3000 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_base}_wiggle.eps
+
+file=pplus.su
+file_base=${file%.su}
+sumax < ${file} mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/5}'`
+
+suflip flip=3 < $file | \
+ supsimage hbox=2.5 wbox=4 labelsize=-1 x2beg=-1000 x2end=1000 \
+ x1beg=-0.6 x1end=1.5 clip=$clip grid1=dot \
+ n1tic=5 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_base}.eps
+
+suflip flip=3 < $file | suwind s=1 j=10 | \
+ supswigp hbox=2.5 wbox=4 labelsize=-1 d2=100 \
+ x1beg=-1.5 x1end=1.5 xcur=1.5 grid1=dot n1tic=1 d1num=2.044 \
+ f2=-3000 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_base}_wiggle.eps
+
+
+
diff --git a/marchenko3D/demo/ScientificReports/initialFocus.scr b/marchenko3D/demo/ScientificReports/initialFocus.scr
new file mode 100755
index 0000000000000000000000000000000000000000..4353514c42aafbaaa8a4799b4a31c68a43d3be36
--- /dev/null
+++ b/marchenko3D/demo/ScientificReports/initialFocus.scr
@@ -0,0 +1,86 @@
+#!/bin/bash
+#PBS -N fdelmod
+#PBS -q verylong
+#PBS -l nodes=1
+#PBS -k eo
+#PBS -j eo
+
+export PATH=$HOME/bin64:$HOME/src/OpenSource/utils:$PATH:
+which makewave
+which makemod
+
+dx=2.5
+dt=0.0005
+dx=0.5
+dt=0.0001
+
+makemod sizex=6000 sizez=810 dx=$dx dz=$dx cp0=1500 ro0=1000 \
+ orig=-3000,0 file_base=tutodeldown.su verbose=2 \
+ intt=def x=-3000,-2000,-1000,0,1000,2000,3000 z=240,130,250,300,350,380,320 poly=2 cp=1950 ro=4500 grad=0 \
+ intt=def x=-3000,-2200,-1500,0,1300,2100,3000 z=620,640,590,600,740,700,600 poly=2 cp=2000 ro=1400 grad=0
+
+file_cp=tutodeldown_cp.su
+file_ro=tutodeldown_cp.su
+
+makewave fp=30 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+dxshot=10
+ishot=300
+nshots=301
+
+export OMP_NUM_THREADS=1
+mkdir -p shots${fast}
+mkdir -p jobs
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+ echo xsrc=$xsrc
+ file_rcv=iniFocusz800x${xsrc}.su
+
+ cat << EOF > jobs/job_$ishot.slurm
+#!/bin/bash
+#
+#SBATCH -J marchenko46
+#SBATCH --cpus-per-task=8
+#SBATCH --ntasks=1
+#SBATCH --time=0:45:00
+#SBATCH -o iniFocus-%A.out
+
+cd \$SLURM_SUBMIT_DIR
+
+export OMP_NUM_THREADS=8
+
+ fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=wave.su \
+ file_rcv=$file_rcv \
+ src_type=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ rec_delay=0.1 \
+ dtrcv=0.004 \
+ verbose=2 \
+ tmod=4.192 \
+ dxrcv=10.0 \
+ xrcv1=-3000 xrcv2=3000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=$xsrc zsrc=800 \
+ ntaper=200 \
+ left=2 right=2 top=2 bottom=2
+EOF
+
+ chmod +x jobs/job_$ishot.slurm
+ jobs/job_$ishot.slurm
+ #sbatch jobs/job_$ishot.slurm
+
+ (( ishot = $ishot + 1))
+done
+
+
+
diff --git a/marchenko3D/demo/ScientificReports/marchenko.scr b/marchenko3D/demo/ScientificReports/marchenko.scr
new file mode 100755
index 0000000000000000000000000000000000000000..427f2c096b3da9633d1f0a4bee3c86007f2d17a4
--- /dev/null
+++ b/marchenko3D/demo/ScientificReports/marchenko.scr
@@ -0,0 +1,22 @@
+#!/bin/bash -x
+#PBS -N fdelmod
+#PBS -q verylong
+#PBS -l nodes=1
+#PBS -k eo
+#PBS -j eo
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+export OMP_NUM_THREADS=8
+
+#for backprop.scr use fast model (with stair-case interfaces)
+#for backprop.scr use pml model (with 0.5 fine-grids)
+
+smooth=3
+
+fmute file_shot=iniFocusz800x0_rp.su file_out=p0plus.su above=-1 shift=-10 verbose=1 check=1 hw=4
+
+marchenko file_shot=shotsnew/refl_rp.su file_tinv=p0plus.su nshots=601 \
+ file_green=pgreen.su verbose=1 tap=0 ntap=10 niter=15 hw=4 shift=10 smooth=$smooth \
+ file_gplus=Gplus.su file_gmin=Gmin.su file_f1plus=f1plus.su file_f1min=f1min.su \
+ file_f2=f2.su fmax=90 file_f1plus=f1plus.su
+
diff --git a/marchenko3D/demo/ScientificReports/model.scr b/marchenko3D/demo/ScientificReports/model.scr
new file mode 100755
index 0000000000000000000000000000000000000000..7a37739be2a78183a3f2582b7f7409d7d8001327
--- /dev/null
+++ b/marchenko3D/demo/ScientificReports/model.scr
@@ -0,0 +1,98 @@
+#!/bin/bash -x
+
+#cd /vardim/home/thorbcke/data/Kees/Marchenko/Nature_SAGA
+
+if [ -n "$1" ]
+then
+ dx=$1
+else
+ dx=2.5
+fi
+
+makemod sizex=6000 sizez=1250 dx=$dx dz=$dx cp0=1500 ro0=1000 \
+ orig=-3000,0 file_base=tutodel_dx$dx.su rayfile=1 skip=100 verbose=2 \
+ intt=def x=-3000,-2000,-1000,0,1000,2000,3000 z=240,130,250,300,350,380,320 poly=2 cp=1950 ro=4500 \
+ intt=def x=-3000,-2200,-1500,0,1300,2100,3000 z=620,640,590,600,740,700,600 poly=2 cp=2000 ro=1400 \
+ intt=def x=-3000,-1800,0,2200,3000 z=920,1000,900,1000,1010 poly=2 cp=2300 ro=1600
+
+suwind key=gx min=0 max=0 < tutodel_dx${dx}_cp.su > tracemodel.su
+
+supsgraph < tracemodel.su hbox=6 wbox=3 labelsize=10 style=seismic \
+ x1beg=0 d1=$dx f1=0 d2num=100 label1="depth (m)" label2="velocity (m/s)"> tracemodel_x_0.eps
+
+supsimage hbox=4 wbox=6 labelsize=10 < tutodel_dx2.5_cp.su \
+ x1beg=0 x1end=1200.0 legend=1 threecolor=1 \
+ d1s=1.0 d2s=1.0 \
+ wrgb=1.0,.5,0 grgb=0,.7,1.0 brgb=0,1.0,0 \
+ bps=24 bclip=2300 wclip=1800 \
+ curve=curve1,curve2,curve3 npair=25,25,25 curvecolor=black,black,black curvedash=3,3,3 \
+ n1tic=5 x2beg=-3000 f2num=-3000 d2num=1000 x2end=3000 > tutodel_cp.eps
+
+#itmax=$(echo "scale=0; (1250/${dx}+1)" | bc -l)
+suwind key=gx min=-1000000 max=1000000 < tutodel_dx${dx}_cp.su | sugain scale=-1 dt=1 | sustrip > tutodelsnap_cp.bin
+
+exit;
+#model reference for SI results
+makewave fp=30 dt=$dt file_out=wave${dx}.su nt=4096 t0=0.1 scale=1
+
+dx=0.5
+file_cp=tutodelbelow_dx${dx}_cp.su
+file_ro=tutodelbelow_dx${dx}_ro.su
+
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=wave${dx}.su \
+ file_rcv=reference_${dx}Z800.su \
+ src_type=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_p=0 \
+ dtrcv=0.004 \
+ dxrcv=2.5 \
+ rec_delay=0.1 \
+ verbose=1 \
+ tmod=1.104 \
+ xrcv1=-3000 xrcv2=3000 zrcv1=800 zrcv2=800 \
+ xsrc=0 zsrc=800 \
+ ntaper=400 \
+ left=4 right=4 top=4 bottom=4
+
+sumax < reference_${dx}Z800_rp.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/30}'`
+sugain epow=0.8 < reference_${dx}Z800_rp.su | \
+ supsimage hbox=4 wbox=8 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=2.5 f1=0.0 x1beg=0.0 x1end=1.004 \
+ f2=-3000 f2num=-3000 d2num=1000 clip=$clip > reference_${dx}Z800_rp.eps
+
+file_cp=tutodel_dx${dx}_cp.su
+file_ro=tutodel_dx${dx}_ro.su
+
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=wave${dx}.su \
+ file_rcv=green_${dx}Z800.su \
+ src_type=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_p=0 \
+ dtrcv=0.004 \
+ dxrcv=2.5 \
+ rec_delay=0.1 \
+ verbose=1 \
+ tmod=1.104 \
+ xrcv1=-3000 xrcv2=3000 zrcv1=800 zrcv2=800 \
+ xsrc=0 zsrc=800 \
+ ntaper=400 \
+ left=4 right=4 top=4 bottom=4
+
+sumax < green_${dx}Z800_rp.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/30}'`
+sugain epow=0.8 < green_${dx}Z800_rp.su | \
+ supsimage hbox=4 wbox=8 labelsize=10 linewidth=0.0 \
+ n1tic=2 d2=2.5 f1=0.0 x1beg=0.0 x1end=1.004 \
+ f2=-3000 f2num=-3000 d2num=1000 clip=$clip > green_${dx}Z800_rp.eps
+
diff --git a/marchenko3D/demo/ScientificReports/remove_direct.scr b/marchenko3D/demo/ScientificReports/remove_direct.scr
new file mode 100755
index 0000000000000000000000000000000000000000..71a99c127a193cd2d9c3cc79ef4e72b7ee95f368
--- /dev/null
+++ b/marchenko3D/demo/ScientificReports/remove_direct.scr
@@ -0,0 +1,42 @@
+#!/bin/bash
+#PBS -N fdelmod
+#PBS -q verylong
+#PBS -l nodes=1
+#PBS -k eo
+#PBS -j eo
+
+export PATH=$HOME/bin:$HOME/src/OpenSource/bin:$PATH:
+
+dxshot=10
+ishot=0
+nshots=601
+
+file_R=shotsnew/refl_rp.su
+rm $file_R
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+ (( sx = ${xsrc}*1000 ))
+ (( iishot = ${ishot}*${dxshot}/10 ))
+ (( tr1 = 601 - ${iishot} ))
+ (( tr2 = ${tr1} + 600 ))
+ echo xsrc=$xsrc tr1=$tr1 tr2=$tr2
+
+# direct wave
+ suwind < direct_rp.su key=tracl min=$tr1 max=$tr2 > direct.su
+
+# 2D shot
+ file_rcv=shotsnew/shots_${xsrc}_rp.su
+ #suwind key=tracl min=1 max=601 < $file_rcv > shotz0.su
+ sudiff $file_rcv direct.su > refl.su
+
+ (( ishot = $ishot + 1))
+
+ sushw < refl.su key=fldr a=$ishot | sugain scale=0.5 >> $file_R
+
+done
+
+rm direct.su refl.su
+
diff --git a/marchenko3D/demo/ScientificReports/shotslurm.scr b/marchenko3D/demo/ScientificReports/shotslurm.scr
new file mode 100755
index 0000000000000000000000000000000000000000..fb5feee5c5172108e37479485121bc23c6a9c2ba
--- /dev/null
+++ b/marchenko3D/demo/ScientificReports/shotslurm.scr
@@ -0,0 +1,67 @@
+#!/bin/bash
+
+export PATH=:$HOME/src/OpenSource/bin:$PATH:
+
+dt=0.0001
+makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=8192 t0=0.4 scale=0 scfft=1
+
+./model.scr 0.5
+
+mkdir -p shotsnew
+mkdir -p jobs
+
+dxshot=10
+ishot=0
+nshots=601
+zsrc=0
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+
+ echo ishot=$ishot xsrc=$xsrc zsrc=$zsrc
+
+ cat << EOF > jobs/slurm_$ishot.job
+#!/bin/bash
+#SBATCH -J shot_${xsrc}
+#SBATCH --cpus-per-task=6
+#SBATCH --ntasks=1
+#SBATCH --time=03:00:00
+#SBATCH -o shot_${xsrc}.out
+
+cd \$SLURM_SUBMIT_DIR
+
+ export PATH=:\$HOME/src/OpenSource/bin:\$PATH:
+ export OMP_NUM_THREADS=6
+ file_rcv=shotsnew/shots_${xsrc}.su
+
+ fdelmodc \
+ file_cp=tutodel_cp.su ischeme=1 iorder=4 \
+ file_den=tutodel_ro.su \
+ file_src=wavefw.su \
+ file_rcv=\$file_rcv \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ rec_delay=0.4 \
+ dtrcv=0.004 \
+ verbose=2 \
+ tmod=4.492 \
+ dxrcv=10.0 \
+ xrcv1=-3000 xrcv2=3000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=$xsrc zsrc=$zsrc \
+ ntaper=200 \
+ left=2 right=2 top=2 bottom=2
+EOF
+
+ sbatch jobs/slurm_$ishot.job
+
+ (( ishot = $ishot + 1))
+
+done
+
diff --git a/marchenko3D/demo/old/README b/marchenko3D/demo/old/README
new file mode 100644
index 0000000000000000000000000000000000000000..5fc50362f6ee58da1104668356b8ac28380d6013
--- /dev/null
+++ b/marchenko3D/demo/old/README
@@ -0,0 +1,9 @@
+Description of files:
+1) shots.scr create the shots
+2) model.scr computes the model
+3) direct_wave.scr crate the direct wave to be removed from the shots
+4) remove_direct.scr remove the direct wave from the shots and scale them
+5) first_arrival.scr computes the first arrival
+6) marchenko.scr perform the Marchenko scheme
+7) referenceShot.scr creates the reference Green's function
+
diff --git a/marchenko3D/demo/old/direct.scr b/marchenko3D/demo/old/direct.scr
new file mode 100755
index 0000000000000000000000000000000000000000..a27121e6f21e354e0d16670ffee3e40e321b7bf9
--- /dev/null
+++ b/marchenko3D/demo/old/direct.scr
@@ -0,0 +1,37 @@
+#!/bin/bash
+
+export PATH=$HOME/bin:$HOME/src/OpenSource/bin:$PATH:
+
+cd /home/thorbcke/data/Kees/MultElim/ModelImageBackprop
+
+dx=2.5
+dt=0.0005
+
+makemod sizex=12000 sizez=4000 dx=$dx dz=$dx cp0=1900 ro0=1200 \
+ orig=-6000,-1000 file_base=noContrast.su
+
+export OMP_NUM_THREADS=8
+makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3 scale=0
+
+$HOME/bin/fdelmodc \
+ file_cp=noContrast_cp.su ischeme=1 iorder=4 \
+ file_den=noContrast_ro.su \
+ file_src=wavefw.su \
+ file_rcv=direct.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=1 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ rec_delay=0.3 \
+ dtrcv=0.004 \
+ verbose=2 \
+ tmod=4.394 \
+ dxrcv=10.0 \
+ xrcv1=-6000 xrcv2=6000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=0 \
+ ntaper=400 \
+ left=4 right=4 top=4 bottom=4
+
diff --git a/marchenko3D/demo/old/first_arrival.scr b/marchenko3D/demo/old/first_arrival.scr
new file mode 100755
index 0000000000000000000000000000000000000000..721721ddfae96deb0f734c57e6d10eb55e37dc04
--- /dev/null
+++ b/marchenko3D/demo/old/first_arrival.scr
@@ -0,0 +1,92 @@
+#!/bin/bash
+#PBS -N fdelmod
+#PBS -q verylong
+#PBS -l nodes=1
+#PBS -k eo
+#PBS -j eo
+
+export PATH=$HOME/bin:$HOME/src/OpenSource/bin:$PATH:
+
+cd /home/thorbcke/data/Kees/MultElim/ModelImageBackprop/Redatum
+
+dx=2.5
+dt=0.0005
+
+makemod sizex=6000 sizez=2000 dx=$dx dz=$dx cp0=1900 ro0=1200 \
+ orig=-3000,0 file_base=synclDown.su verbose=2 \
+ intt=def x=-3000,500,3000 z=195,195,195 poly=1 cp=1950 ro=3700 \
+ intt=def x=-3000,3000 z=600,600 poly=0 cp=2050 ro=1750 \
+ intt=def x=-3000,3000 z=680,680 poly=0 cp=2150 ro=2220 \
+ 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 \
+
+# intt=def x=-3000,0,3000 z=1110,1110,1110 poly=0 cp=2300 ro=1950 \
+# intt=def x=-3000,3000 z=1180,1180 poly=0 cp=2480 ro=1820 \
+# intt=def x=-3000,0,3000 z=1290,1290,1370 poly=0 cp=2600 ro=2000 \
+# intt=def x=-3000,3000 z=1380,1380 poly=0 cp=2720 ro=2050 \
+# intt=def x=-3000,3000 z=1480,1480 poly=0 cp=2800 ro=1850
+
+makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+#smooth file_in=synclDown_cp.su power=-1.0 ntsm=29 nxsm=29 niter=15 file_out=syncls_cp.su
+#smooth file_in=synclDown_ro.su power=-1.0 ntsm=29 nxsm=29 niter=15 file_out=syncls_ro.su
+
+dxshot=10
+ishot=300
+nshots=301
+
+export OMP_NUM_THREADS=1
+mkdir -p shots
+mkdir -p jobs
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+# (( xsrc = -1100 + ${ishot}*${dxshot} ))
+ echo xsrc=$xsrc
+ file_rcv=shots/shotsmonPz1100_${xsrc}.su
+
+ cat << EOF > jobs/pbs_$ishot.job
+#!/bin/bash
+#
+#PBS -q medium
+#PBS -N mod_${xsrc}
+#PBS -j eo
+#PBS -m n
+#PBS -l nodes=1
+#PBS -V
+
+cd /home/thorbcke/data/Kees/MultElim/ModelImageBackprop/Redatum
+export OMP_NUM_THREADS=1
+
+ $HOME/bin64/fdelmodc \
+ file_cp=synclDown_cp.su ischeme=1 iorder=4 \
+ file_den=synclDown_ro.su \
+ file_src=wave.su \
+ file_rcv=$file_rcv \
+ src_type=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ rec_delay=0.1 \
+ dtrcv=0.004 \
+ verbose=2 \
+ tmod=2.100 \
+ dxrcv=10.0 \
+ xrcv1=-3000 xrcv2=3000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=$xsrc zsrc=1100 \
+ ntaper=300 \
+ left=4 right=4 top=4 bottom=4
+EOF
+
+ qsub jobs/pbs_$ishot.job
+
+ (( ishot = $ishot + 1))
+done
+
+
+
diff --git a/marchenko3D/demo/old/marchenko.scr b/marchenko3D/demo/old/marchenko.scr
new file mode 100755
index 0000000000000000000000000000000000000000..bee3690f0b9482997ab64f4a652c3244cdda6c9a
--- /dev/null
+++ b/marchenko3D/demo/old/marchenko.scr
@@ -0,0 +1,39 @@
+#!/bin/bash -x
+#PBS -N fdelmod
+#PBS -q verylong
+#PBS -l nodes=1
+#PBS -k eo
+#PBS -j eo
+
+export PATH=$HOME/bin:$HOME/src/OpenSource/bin:$PATH:
+which makewave
+which makemod
+which fmute
+which syn2d
+export OMP_NUM_THREADS=8
+
+cd /home/thorbcke/data/Kees/MultElim/ModelImageBackprop/Redatum
+
+tmpdir=/tmp/shotI
+mkdir -p $tmpdir
+#for dt=0.004 with modeling at 0.0005
+scale=1.0
+w1=1
+smooth=3
+#smooth=0
+
+fmute file_shot=shots/shotsmonPz1100_0_rp.su file_out=p0plus.su above=-1 shift=-10 verbose=1 check=1 hw=4
+
+suwind itmax=1023 < p0plus.su | \
+suwind key=gx min=-3000000 max=3000000 | \
+sushw key=fldr a=1 > p0plussx.su
+
+#~/bin/marchenko file_shot=../shots/refl_rp.su file_tinv=p0plussx.su nshots=601 file_green=pgreen_1.su verbose=1 tap=0 ntap=0 reci=0 niter=1 hw=8 shift=7 smooth=$smooth w=1 file_gplus=Gplus0.su file_gmin=Gmin0.su file_f1plus=f1plus0_1.su file_f1min=f1min0_1.su file_pplus=Pplus0_1.su
+#
+#~/bin/marchenko file_shot=../shots/refl_rp.su file_tinv=p0plussx.su nshots=601 file_green=pgreen_4.su verbose=1 tap=0 ntap=0 reci=0 niter=4 hw=8 shift=7 smooth=$smooth w=1 file_gplus=Gplus0.su file_gmin=Gmin0.su file_f1plus=f1plus0_4.su file_f1min=f1min0_4.su file_pplus=Pplus0_4.su
+
+#for backpropagating pplus in marchenko scheme must be written to file
+~/bin/marchenko file_shot=../shots/refl_rp.su file_tinv=p0plussx.su nshots=601 file_green=pgreen.su verbose=1 tap=0 ntap=10 niter=15 hw=8 shift=7 smooth=$smooth file_gplus=Gplus0.su file_gmin=Gmin0.su file_f1plus=f1plus0.su file_f1min=f1min0.su file_pplus=Pplus0.su
+
+exit;
+
diff --git a/marchenko3D/demo/old/model.scr b/marchenko3D/demo/old/model.scr
new file mode 100755
index 0000000000000000000000000000000000000000..777a556cff71898882ebca61466f17fb7c09aaba
--- /dev/null
+++ b/marchenko3D/demo/old/model.scr
@@ -0,0 +1,212 @@
+#!/bin/bash
+
+export PATH=$HOME/bin:$HOME/src/OpenSource/bin:$PATH:
+
+dx=2.5
+dt=0.0005
+
+makemod sizex=6000 sizez=2000 dx=$dx dz=$dx cp0=1900 ro0=1200 \
+ orig=-3000,0 file_base=syncl.su verbose=2 \
+ intt=def x=-3000,500,3000 z=195,195,195 poly=1 cp=1950 ro=3700 \
+ intt=def x=-3000,3000 z=600,600 poly=0 cp=2050 ro=1750 \
+ intt=def x=-3000,3000 z=680,680 poly=0 cp=2150 ro=2220 \
+ 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 \
+ intt=def x=-3000,0,3000 z=1110,1110,1110 poly=0 cp=2300 ro=1950 \
+ intt=def x=-3000,3000 z=1180,1180 poly=0 cp=2480 ro=1820 \
+ intt=def x=-3000,0,3000 z=1290,1290,1370 poly=0 cp=2600 ro=2000 \
+ intt=def x=-3000,3000 z=1380,1380 poly=0 cp=2720 ro=2050 \
+ intt=def x=-3000,3000 z=1480,1480 poly=0 cp=2800 ro=1850
+
+ #intt=diffr x=-2000,-1000,0,1000,2000 z=1800,1800,1800,1800,1800 cp=0,0,0,0,0 ro=5000,5000,5000,5000,5000
+
+
+makemod sizex=6000 sizez=2000 dx=10 dz=5 cp0=1900 ro0=1200 \
+ orig=-3000,0 file_base=syncl_migr.su verbose=2 \
+ intt=def x=-3000,500,3000 z=195,195,195 poly=1 cp=1950 ro=3700 \
+ intt=def x=-3000,3000 z=600,600 poly=0 cp=2050 ro=1750 \
+ intt=def x=-3000,3000 z=680,680 poly=0 cp=2150 ro=2220 \
+ 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 \
+ intt=def x=-3000,0,3000 z=1110,1110,1110 poly=0 cp=2300 ro=1950 \
+ intt=def x=-3000,3000 z=1180,1180 poly=0 cp=2480 ro=1820 \
+ intt=def x=-3000,0,3000 z=1290,1290,1370 poly=0 cp=2600 ro=2000 \
+ intt=def x=-3000,3000 z=1380,1380 poly=0 cp=2720 ro=2050 \
+ intt=def x=-3000,3000 z=1480,1480 poly=0 cp=2800 ro=1850
+
+exit
+
+#example FD modeling with model defined above
+makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3 scale=0
+
+export OMP_NUM_THREADS=1
+
+makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+~/bin/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=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.0010 \
+ 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=400 \
+ tsnap1=3.1 tsnap2=2.5 dtsnap=0.1 \
+ left=4 right=4 top=4 bottom=4
+
+
+exit
+
+makemod sizex=6000 sizez=2000 dx=$dx dz=$dx cp0=1900 ro0=1200 \
+ orig=-3000,-1000 file_base=hom.su
+
+~/bin/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=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.0010 \
+ 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=400 \
+ left=4 right=4 top=4 bottom=4
+
+sudiff shot_fd_rp.su shot_hom_fd_rp.su > shot_rp.su
+
+migr file_shot=shot_rp.su file_vel=syncl_migr_cp.su imc=0 file_image=image.su verbose=1
+
+exit
+
+makemod sizex=6000 sizez=4000 dx=$dx dz=$dx cp0=$cp cs0=$cs ro0=$rho \
+ orig=-3000,-1000 file_base=synclTop.su \
+ intt=def x=-3000,0,3000 z=200,200,200 poly=0 cp=1800 ro=5000 \
+ intt=def x=-3000,-2000,-1000,-400,0,200,900,1800,3000 z=520,520,560,670,950,790,600,520,500 poly=2 cp=2300 ro=1800 \
+
+~/bin/fdelmodc \
+ file_cp=synclTop_cp.su ischeme=1 iorder=4 \
+ file_den=synclTop_ro.su \
+ file_src=wave.su \
+ file_rcv=p0.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.0010 \
+ rec_delay=0.1 \
+ verbose=2 \
+ tmod=2.100 \
+ dxrcv=10.0 \
+ xrcv1=-2250 xrcv2=2250 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=1100 \
+ ntaper=300 \
+ left=4 right=4 top=4 bottom=4 &
+
+~/bin/fdelmodc \
+ file_cp=synclTop_cp.su ischeme=1 iorder=4 \
+ file_den=synclTop_ro.su \
+ file_src=wave.su \
+ file_rcv=shot_top.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.0010 \
+ 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=300 \
+ left=4 right=4 top=4 bottom=4 &
+
+makemod sizex=6000 sizez=4000 dx=$dx dz=$dx cp0=2300 ro0=1800 \
+ orig=-3000,-1000 file_base=synclBot.su \
+ intt=def x=-3000,0,3000 z=1310,1310,1310 poly=0 cp=2450 ro=1950 \
+ intt=def x=-3000,3000 z=1380,1380 poly=0 cp=2460 ro=1820 \
+ intt=def x=-3000,0,3000 z=1490,1490,1570 poly=0 cp=2470 ro=2100 \
+ intt=def x=-3000,3000 z=1580,1580 poly=0 cp=2480 ro=2000 \
+ intt=def x=-3000,3000 z=1680,1680 poly=0 cp=2490 ro=1850
+
+~/bin/fdelmodc \
+ file_cp=synclBot_cp.su ischeme=1 iorder=4 \
+ file_den=synclBot_ro.su \
+ file_src=wave.su \
+ file_rcv=pRef.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.0010 \
+ rec_delay=0.1 \
+ verbose=2 \
+ tmod=2.100 \
+ dxrcv=10.0 \
+ xrcv1=-2250 xrcv2=2250 \
+ zrcv1=1100 zrcv2=1100 \
+ xsrc=0 zsrc=1100 \
+ ntaper=300 \
+ left=4 right=4 top=4 bottom=4 &
+
+makemod sizex=6000 sizez=4000 dx=$dx dz=$dx cp0=2300 ro0=1800 \
+ orig=-3000,-1000 file_base=synclBotHom.su
+
+~/bin/fdelmodc \
+ file_cp=synclBotHom_cp.su ischeme=1 iorder=4 \
+ file_den=synclBotHom_ro.su \
+ file_src=wave.su \
+ file_rcv=pRefHom.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.0010 \
+ rec_delay=0.1 \
+ verbose=2 \
+ tmod=2.100 \
+ dxrcv=10.0 \
+ xrcv1=-2250 xrcv2=2250 \
+ zrcv1=1100 zrcv2=1100 \
+ xsrc=0 zsrc=1100 \
+ ntaper=300 \
+ left=4 right=4 top=4 bottom=4 &
+
+wait
+
+sudiff shot_fd_rp.su shot_hom_fd_rp.su > shot_rp.su
+sudiff shot_top_rp.su shot_hom_fd_rp.su > shotTop_rp.su
+sudiff pRef_rp.su pRefHom_rp.su > pref_rp.su
+
diff --git a/marchenko3D/demo/old/referenceShot.scr b/marchenko3D/demo/old/referenceShot.scr
new file mode 100755
index 0000000000000000000000000000000000000000..b1ea70e1a0824cc1e60be117ff15110469f57942
--- /dev/null
+++ b/marchenko3D/demo/old/referenceShot.scr
@@ -0,0 +1,44 @@
+#!/bin/bash
+#PBS -N fdelmod
+#PBS -q verylong
+#PBS -l nodes=1
+#PBS -k eo
+#PBS -j eo
+
+export PATH=$HOME/bin:$HOME/src/OpenSource/bin:$PATH:
+
+cd /home/thorbcke/data/Kees/MultElim/ModelImageBackprop/Redatum
+
+#makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3
+
+dx=2.5
+dt=0.0005
+
+makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+export OMP_NUM_THREADS=8
+
+$HOME/bin/fdelmodc \
+ file_cp=../syncl_cp.su ischeme=1 iorder=4 \
+ file_den=../syncl_ro.su \
+ file_src=wave.su \
+ file_rcv=virtual_shot_fd_P_zsrc1100.su \
+ src_type=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_ud=1 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.1 \
+ verbose=2 \
+ tmod=2.144 \
+ dxrcv=10.0 \
+ xrcv1=-3000 xrcv2=3000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=1100 \
+ file_snap=backpropref.su tsnap1=0.1 dtsnap=0.010 tsnap2=2.100 dxsnap=10 dzsnap=10 zsnap1=0 zsnap2=2000 xsnap1=-2250 xsnap2=2250 \
+ ntaper=400 \
+ left=4 right=4 top=4 bottom=4
+
+
diff --git a/marchenko3D/demo/old/remove_direct.scr b/marchenko3D/demo/old/remove_direct.scr
new file mode 100755
index 0000000000000000000000000000000000000000..24a3f41e2c686f981ab12a38ed6ba4fc9eca4a1e
--- /dev/null
+++ b/marchenko3D/demo/old/remove_direct.scr
@@ -0,0 +1,38 @@
+#!/bin/bash
+#PBS -N fdelmod
+#PBS -q verylong
+#PBS -l nodes=1
+#PBS -k eo
+#PBS -j eo
+
+export PATH=$HOME/bin:$HOME/src/OpenSource/bin:$PATH:
+
+cd /home/thorbcke/data/Kees/MultElim/ModelImageBackprop
+
+dxshot=10
+ishot=0
+nshots=601
+
+rm shots/refl_rp.su
+
+while (( ishot < nshots ))
+do
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+ (( iishot = ${ishot}*${dxshot}/10 ))
+ (( tr1 = 601 - ${iishot} ))
+ (( tr2 = ${tr1} + 600 ))
+ echo xsrc=$xsrc tr1=$tr1 tr2=$tr2
+ suwind < direct_rp.su key=tracl min=$tr1 max=$tr2 > direct.su
+
+ file_rcv=shots/shots_${xsrc}_rp.su
+ suwind key=tracl min=1 max=601 < $file_rcv > shotz0.su
+
+ sudiff shotz0.su direct.su > refl.su
+
+ (( ishot = $ishot + 1))
+
+ sushw < refl.su key=fldr a=$ishot | \
+ suwind itmax=1023 >> shots/refl_rp.su
+
+done
+
diff --git a/marchenko3D/demo/old/shots.scr b/marchenko3D/demo/old/shots.scr
new file mode 100755
index 0000000000000000000000000000000000000000..907744a972c55c118a56b15237c55a028f269ad4
--- /dev/null
+++ b/marchenko3D/demo/old/shots.scr
@@ -0,0 +1,77 @@
+#!/bin/bash
+#PBS -N fdelmod
+#PBS -q long
+#PBS -l nodes=1
+#PBS -k eo
+#PBS -j eo
+
+export PATH=$HOME/bin:$HOME/src/OpenSource/bin:$PATH:
+which makewave
+which makemod
+
+cd /home/thorbcke/data/Kees/MultElim/ModelImageBackprop
+
+dt=0.0005
+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
+
+./model.scr
+
+mkdir -p shots
+mkdir -p jobs
+
+dxshot=10
+ishot=0
+nshots=601
+zsrc=0
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+
+ echo ishot=$ishot xsrc=$xsrc zsrc=$zsrc
+
+ cat << EOF > jobs/pbs_$ishot.job
+#!/bin/bash
+#
+#PBS -q medium
+#PBS -N mod_${xsrc}
+#PBS -j eo
+#PBS -m n
+#PBS -l nodes=1
+#PBS -V
+
+cd /home/thorbcke/data/Kees/MultElim/ModelImageBackprop
+
+ export OMP_NUM_THREADS=1
+ file_rcv=shots/shots_${xsrc}.su
+
+ $HOME/bin/fdelmodc \
+ file_cp=syncl_cp.su ischeme=1 iorder=4 \
+ file_den=syncl_ro.su \
+ file_src=wavefw.su \
+ file_rcv=\$file_rcv \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=1 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ rec_delay=0.3 \
+ dtrcv=0.004 \
+ verbose=2 \
+ tmod=4.394 \
+ dxrcv=10.0 \
+ xrcv1=-3000,-3000,-3000 xrcv2=3000,3000,3000 \
+ zrcv1=0,1000,1600 zrcv2=0,1000,1600 \
+ xsrc=$xsrc zsrc=$zsrc \
+ ntaper=400 \
+ left=4 right=4 top=4 bottom=4
+EOF
+
+qsub jobs/pbs_$ishot.job
+
+ (( ishot = $ishot + 1))
+
+done
+
diff --git a/marchenko3D/demo/oneD/README b/marchenko3D/demo/oneD/README
new file mode 100644
index 0000000000000000000000000000000000000000..fdc4653ad6e4c112a1c9c8d54806480859efa139
--- /dev/null
+++ b/marchenko3D/demo/oneD/README
@@ -0,0 +1,199 @@
+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) initialFocus.scr model G_d the initial focusing function => iniFocus_rp.su
+4) referenceShot.scr creates the reference Green's function at focal point => referenceP_rp.su
+5) marchenko.scr perform the Marchenko scheme => pgreen.su, f1plus0.su, f1min0.su, f2.su
+
+extra scripts
++) marchenkoIter.scr : to make the figure with "Four iterations of the Marchenko method."
++) backpropf2.scr : to make Figure "Snapshots of back-propagation of f_2."
++) eps*.scr : reproduce the postscript files of the manuscript using SU postscript plotting programs.
++) backProp_f2sum_movie.scr : produces a snapshot move of f2(-t) + f2(t) ; Figure 10 3'rd column
++) clean : remove all produced files and start with a clean directory
+
+
+To reproduce the Figures in the Manuscript:
+
+--------------------------
+* Figure 2: Wavelet
+* Figure 3: Model + Initial wavefield
+
+==> run model.scr to generate the data .su files: this will take 3-4 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 initialFocus.scr to compute the direct arrival of the transmission response G_d. This will take 1-2 minutes.
+ - modelup_cp.su
+ - modelup_ro.su
+ - iniFocus_rp.su
+Note if you model the initial Focusing operator also with a w=fw wavelet the length of the wavelet becomes very long. The
+mute-windows applied in Marchenko will then also mute a big part of this very long fw wavelet and will not converge anymore.
+
+
+==> run epsModel.scr to generate the postscript files of Figure 2 and 3
+
+wavefw.eps => Figure 2a
+wavefw_freq.eps => Figure 2b
+
+model_cp_line.eps => Figure 3a
+model_ro_line.eps => Figure 3b
+shotx0_rp.eps => Figure 3c
+iniFocus_rp.eps => Figure 3d
+
+
+--------------------------
+* Figure 4: Initialisation
+* Figure 5: first update
+* Figure 6: first 4 iterations
+
+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, together with iniFocus_rp.su, is the input of the Marchenko algorithm
+
+==> run marchenkoIter.scr to compute the first 4 iteration of the Marchenko algorithm. This will take 1-2 minutes. The generated files are:
+ - p0plus.su
+ - pgreen_001.su
+ - f1plus_001.su
+ - f1min_001.su
+ - Gplus_001.su
+ - Gmin_001.su
+ - pgreen_002.su
+ - f1plus_002.su
+ - f1min_002.su
+ - Gplus_002.su
+ - Gmin_002.su
+ - pgreen_003.su
+ - f1plus_003.su
+ - f1min_003.su
+ - Gplus_003.su
+ - Gmin_003.su
+ - pgreen_004.su
+ - f1plus_004.su
+ - f1min_004.su
+ - Gplus_004.su
+ - Gmin_004.su
+
+To Compute the reference Green's function at x=0 z=900 m in the actual model
+==> run referenceShot.scr This will take 1 minute and generates the file;
+ - referenceP_rp.su
+
+To generate all postscript files for Figure 4, 5 and 6
+
+==> run epsMarchenkoIter.scr
+
+shotx0_rp.eps => Figure 4 R == Figure 3c
+p0plus.eps => Figure 4 G_d
+iter_001.eps => Figure 4 N_0
+
+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.
+
+--------------------------
+* Figure 7: Comparison of Marchenko result with reference
+
+To compute the marchenko results for 8 iterations.
+
+==> 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
+
+
+At the end of the run the script will display in X11 a comparison of the middle trace.
+
+To make the postscript figure
+
+==> run epsCompare.scr
+
+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
+
diff --git a/marchenko3D/demo/oneD/backProp_f2sum_movie.scr b/marchenko3D/demo/oneD/backProp_f2sum_movie.scr
new file mode 100755
index 0000000000000000000000000000000000000000..1c5168b5c82d5fa6d7d100588a102858861619aa
--- /dev/null
+++ b/marchenko3D/demo/oneD/backProp_f2sum_movie.scr
@@ -0,0 +1,34 @@
+#!/bin/bash
+
+# $1 input-file of backpropf2_sp.su format
+
+export PATH=$HOME/OpenSource/bin/:$PATH:
+
+#rm prop_movie.su
+
+for fldr in $(seq 1 101);
+do
+ times=$(echo "scale=2; -0.01*(${fldr}-101)" | bc -l)
+ echo " Adding ${fldr}-th frame"
+ suwind key=fldr min=$fldr max=$fldr < backpropf2_sp.su > neg.su
+ (( fldr = 101+(101-$fldr) ))
+ suwind key=fldr min=$fldr max=$fldr < backpropf2_sp.su > pos.su
+ susum neg.su pos.su >prop_frame.su
+
+ if [ "$fldr" != "201" ]; then
+ cat prop_movie.su prop_frame.su >temp_movie.su
+ mv temp_movie.su prop_movie.su
+ else
+ cp prop_frame.su prop_movie.su
+ fi
+done
+
+suchw <prop_movie.su key1=ntr key2=fldr b=1001 >temp_movie.su
+
+susort <temp_movie.su -fldr tracf >prop_movie.su
+
+rm neg.su pos.su temp_movie.su
+
+n2=`surange <backpropf2_sp.su | grep tracf | awk '{print $3}'`
+
+suxmovie <prop_movie.su loop=1 height=500 width=1000 title=%g n2=$n2
diff --git a/marchenko3D/demo/oneD/backpropf2.scr b/marchenko3D/demo/oneD/backpropf2.scr
new file mode 100755
index 0000000000000000000000000000000000000000..2041570cce37829889067b1a50a72277402c7c5f
--- /dev/null
+++ b/marchenko3D/demo/oneD/backpropf2.scr
@@ -0,0 +1,59 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin/:$PATH:
+
+dx=2.5
+dt=0.0005
+
+file_cp=model10_cp.su
+file_ro=model10_ro.su
+
+export OMP_NUM_THREADS=4
+
+# t=0 focal time is at 2.0445 seconds back=propagating
+# shift f2.su such that t=0 is positioned in the middle of the time axis
+# the extra shift of 0.000250 is needed because of the staggered time implementation of the Finite Difference program.
+ns=1024
+dtrcv=`surange < f2.su | grep dt | awk '{print $2/1000000.0}'`
+suwind key=gx min=-2250000 max=2250000 itmax=1023 < f2.su > nep.su
+shift=$(echo "scale=6; ($dtrcv*($ns/2.0-1)+$dt)" | bc -l)
+echo $shift
+basop choice=shift shift=$shift file_in=nep.su verbose=1 > pplus.su
+
+
+# the f2.su is sampled with 4ms the FD program need 0.5ms
+# time axis is interpolated by making use of FFT's: sinc interpolation
+#ftr1d file_in=pplus.su file_out=freq.su
+#sushw < freq.su key=nhs,dt a=8192,500 > fr.su
+#ftr1d file_in=fr.su n1=8194 file_out=pplusdt.su verbose=1
+#this is now done within the fdelmodc code:
+
+#backpropagate f2.su and collect snapshots
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=pplus.su \
+ dt=$dt \
+ file_rcv=backprop_f2_z900.su \
+ grid_dir=0 \
+ src_type=1 \
+ src_injectionrate=1 \
+ src_orient=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.0 \
+ verbose=2 \
+ tmod=3.10 \
+ dxrcv=5.0 \
+ xrcv1=-2250 xrcv2=2250 \
+ zrcv1=900 zrcv2=900 \
+ zsrc=0 xsrc=0 \
+ npml=101 \
+ file_snap=backpropf2.su tsnap1=1.040 dtsnap=0.01 tsnap2=3.040 dxsnap=5 dzsnap=5 zsnap1=0 zsnap2=1250 xsnap1=-1000 xsnap2=1000 \
+ sna_type_vz=0 \
+ sna_type_p=1 \
+ left=2 right=2 top=2 bottom=2
+
+
diff --git a/marchenko3D/demo/oneD/clean b/marchenko3D/demo/oneD/clean
new file mode 100755
index 0000000000000000000000000000000000000000..3890128152ba3f4b11471dfdb5ddd1399840bc08
--- /dev/null
+++ b/marchenko3D/demo/oneD/clean
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+rm *.su *.bin *.eps nep line* *.asci
+
diff --git a/marchenko3D/demo/oneD/conv.gnp b/marchenko3D/demo/oneD/conv.gnp
new file mode 100644
index 0000000000000000000000000000000000000000..119341bef971d8e8dc3e7e4123c32a64f674f5d3
--- /dev/null
+++ b/marchenko3D/demo/oneD/conv.gnp
@@ -0,0 +1,15 @@
+set style data linespoints
+set mytics 10
+set xlabel 'number of iterations'
+set ylabel 'convergence rate'
+set size 2.0,2.0
+set size ratio 0.6
+set grid
+
+set log y
+set nolog x
+
+set term postscript eps font 'Helvetica,12' linewidth 4 fontscale 3
+set output 'convergence.eps'
+plot 'conv.txt' using 1:($2) lw 3 notitle
+
diff --git a/marchenko3D/demo/oneD/conv.txt b/marchenko3D/demo/oneD/conv.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f2e99f36e05c42e47cbbae161bfad8327bae2530
--- /dev/null
+++ b/marchenko3D/demo/oneD/conv.txt
@@ -0,0 +1,16 @@
+0 1.000000e+00
+1 8.104102e-01
+2 2.776407e-01
+3 1.775258e-01
+4 1.278046e-01
+5 8.376110e-02
+6 6.221900e-02
+7 4.089906e-02
+8 3.275844e-02
+9 2.070254e-02
+10 1.920658e-02
+11 1.091778e-02
+12 1.282995e-02
+13 6.060715e-03
+14 9.706275e-03
+15 3.603180e-03
diff --git a/marchenko3D/demo/oneD/epsBackprop.scr b/marchenko3D/demo/oneD/epsBackprop.scr
new file mode 100755
index 0000000000000000000000000000000000000000..5c2ecbcc92358b2cb0fe58a5914f54aabf1b4dc6
--- /dev/null
+++ b/marchenko3D/demo/oneD/epsBackprop.scr
@@ -0,0 +1,68 @@
+#!/bin/bash
+
+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=5
+file_snap="backpropf2"
+dtsnap=0.01
+nsnap=101
+
+sumax < ${file_snap}_sp.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/2}'`
+
+#first snap-shot with labels
+# fldr=71
+# times=$(echo "scale=2; $dtsnap*(${fldr}-$nsnap)" | bc -l)
+# atime=`printf "%4.2f" $times`
+# suwind key=fldr min=$fldr max=$fldr < ${file_snap}_sp.su | \
+# supsimage hbox=4 wbox=6 labelsize=10 \
+# label1="depth (m)" label2="lateral distance (m)" \
+# x1beg=0 x1end=1250.0 clip=${clip} \
+# curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black curvedash=3,3,3 \
+# n1tic=4 f2=-1000 d2=$dx x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_snap}_${atime}_labels.eps
+
+for fldr in 71 86 98 99 101 103 104 116 131;
+do
+ times=$(echo "scale=2; $dtsnap*(${fldr}-$nsnap)" | bc -l)
+ atime=`printf "%4.2f" $times`
+ suwind key=fldr min=$fldr max=$fldr < ${file_snap}_sp.su | \
+ supsimage hbox=4 wbox=6 labelsize=10 \
+ x1beg=0 x1end=1250.0 clip=${clip} \
+ curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black curvedash=3,3,3 \
+ n1tic=4 f2=-1000 d2=$dx x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_snap}_$atime.eps
+done
+
+#select files for snapshot between -0.7 => 0 <= +0.07 (fldr 31-101-171)
+#add pos and negative times to get response of homogenoeus Green's function
+
+file_snap="backpropf2"
+for fldr in 71 86 98 99 101;
+do
+ times=$(echo "scale=2; -0.01*(${fldr}-101)" | bc -l)
+ atime=`printf "%4.2f" $times`
+ suwind key=fldr min=$fldr max=$fldr < ${file_snap}_sp.su > neg.su
+ (( fldr = 101+(101-$fldr) ))
+ suwind key=fldr min=$fldr max=$fldr < ${file_snap}_sp.su > pos.su
+ susum neg.su pos.su | \
+ supsimage hbox=4 wbox=6 labelsize=10 \
+ x1beg=0 x1end=1250.0 clip=${clip} \
+ curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black curvedash=3,3,3 \
+ n1tic=4 x2beg=-1000 d2=$dx f2num=-1000 d2num=500 x2end=1000 > ${file_snap}sum_$atime.eps
+done
+
diff --git a/marchenko3D/demo/oneD/epsCompare.scr b/marchenko3D/demo/oneD/epsCompare.scr
new file mode 100755
index 0000000000000000000000000000000000000000..e9dae68ee730f2cceb5b36018040ceca5b1f49fa
--- /dev/null
+++ b/marchenko3D/demo/oneD/epsCompare.scr
@@ -0,0 +1,37 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin/:$PATH:
+
+#mke figures for reference and Marchenko result an merge into one file
+
+file=diffref.su
+file_base=${file%.su}
+sumax < referenceP_rp.su mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1}'`
+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
+suwind key=gx min=-2250000 max=2250000 < referenceP_rp.su | \
+ 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 > referenceP_rp.eps
+
+suwind < pgreen512.su j=50 s=1 | \
+ supswigp n2=19 fill=0 \
+ hbox=4 wbox=8 labelsize=10 linewidth=1.0 \
+ label1="time (s)" label2="lateral distance (m)" \
+ n1tic=2 d2=250 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=500 > green.eps
+suwind < referenceP_rp.su j=50 s=1 | \
+ supswigp n2=19 fill=0 tracecolor=#F \
+ hbox=4 wbox=8 labelsize=10 linewidth=2.0 \
+ label1="time (s)" label2="lateral distance (m)" \
+ n1tic=2 d2=250 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=500 > ref.eps
+
+sed -i.old -e "s/%%EndProlog/[ 1 1 ] 0 setdash %%EndProlog/" green.eps
+sed -i.old -e "s/0.5 0.5 0.5 setrgbcolor/0.65 0.65 0.65 setrgbcolor /" ref.eps
+
+psmerge in=ref.eps in=green.eps > mergeGreenRef.eps
+
diff --git a/marchenko3D/demo/oneD/epsIterwithLabels.scr b/marchenko3D/demo/oneD/epsIterwithLabels.scr
new file mode 100755
index 0000000000000000000000000000000000000000..cfb5a6a0e6e0a2515986007b94d8ee53faa4a74d
--- /dev/null
+++ b/marchenko3D/demo/oneD/epsIterwithLabels.scr
@@ -0,0 +1,76 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin/:$PATH:
+
+#mute to get pslinepos.asci files used in plotting only
+fmute file_shot=iniFocus_rp.su file_out=nep.su above=0 shift=8 verbose=1 check=1 hw=4
+
+#set same clip factor for iteration updates
+file=iter_001.su
+sumax < $file mode=abs outpar=nep
+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 ))
+do
+piter=$(printf %03d $iter)
+echo $piter
+
+file=iter_$piter.su
+#ns=`surange < iter_001.su | grep ns | awk '{print $2}'`
+#dtrcv=`surange < iter_001.su | grep dt | awk '{print $2/1000000.0}'`
+#shift=$(echo "scale=4; ($dtrcv*($ns/2.0-1))" | bc -l)
+#basop choice=shift shift=$shift file_in=$file | \
+file_base=${file%.su}
+clipref=$clipiter
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
+ label1="time (s)" label2="lateral distance (m)" \
+ n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
+ curve=pslinepos.asci,pslineneg.asci npair=901,901 curvewidth=2,2 curvecolor=black,black curvedash=3,3 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_labels.eps
+
+file=f1min_$piter.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/5}'`
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file\
+ label1="time (s)" label2="lateral distance (m)" \
+ n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_labels.eps
+
+file=f1plus_$piter.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/5}'`
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file\
+ label1="time (s)" label2="lateral distance (m)" \
+ n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_labels.eps
+
+file=pgreen_$piter.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+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 \
+ label1="time (s)" label2="lateral distance (m)" \
+ n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipgreen > ${file_base}_labels.eps
+
+done
+
+
+#special treatment of f1+ zero-iteration: which is zero, to make a nice gray plot (and not black)
+file=f1plus_001.su
+file_base=${file%.su}
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file\
+ label1="time (s)" label2="lateral distance (m)" \
+ n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 bclip=1 wclip=-1> ${file_base}_labels.eps
+
diff --git a/marchenko3D/demo/oneD/epsMarchenkoIter.scr b/marchenko3D/demo/oneD/epsMarchenkoIter.scr
new file mode 100755
index 0000000000000000000000000000000000000000..b2a417474810933105a76d01eb2c37168367ffda
--- /dev/null
+++ b/marchenko3D/demo/oneD/epsMarchenkoIter.scr
@@ -0,0 +1,121 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin/:$PATH:
+
+
+#Direct field of transmission repsponse
+file=p0plus.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/2}'`
+
+ns=1024
+dtrcv=`surange < p0plus.su | grep dt | awk '{print $2/1000000.0}'`
+suwind key=gx min=-2250000 max=2250000 itmax=1023 < $file > nep.su
+shift=$(echo "scale=4; ($dtrcv*($ns/2.0-1))" | bc -l)
+basop choice=shift shift=$shift file_in=nep.su | \
+ suflip flip=3 | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ label1="time (s)" label2="lateral distance (m)" \
+ n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_flip.eps
+rm nep.su
+
+file=p0plus.su
+file_base=${file%.su}
+suwind key=gx min=-2250000 max=2250000 < $file | \
+ suflip flip=3 | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ label1="time (s)" label2="lateral distance (m)" \
+ n1tic=2 d2=5 f1=-2.044 f1num=-2.000 x1beg=-2.004 x1end=0.0 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+#convolution of G_d with middle shot from R - not used -
+#suwind key=gx min=-2250000 max=2250000 < shot5_rp.su > shot0.su
+#fconv file_in1=iniFocus_rp.su file_in2=shot0.su file_out=GdRconv.su
+
+#mute to get pslinepos.asci files used in plotting only
+fmute file_shot=iniFocus_rp.su file_out=nep.su above=0 shift=8 verbose=1 check=1 hw=4
+
+#set same clip factor for iteration updates
+file=iter_001.su
+sumax < $file mode=abs outpar=nep
+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 ))
+do
+piter=$(printf %03d $iter)
+echo $piter
+
+file=iter_$piter.su
+#ns=`surange < iter_001.su | grep ns | awk '{print $2}'`
+#dtrcv=`surange < iter_001.su | grep dt | awk '{print $2/1000000.0}'`
+#shift=$(echo "scale=4; ($dtrcv*($ns/2.0-1))" | bc -l)
+#basop choice=shift shift=$shift file_in=$file | \
+file_base=${file%.su}
+clipref=$clipiter
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
+ n1tic=2 d2=5 x1beg=-1.504 x1end=1.5 d1num=0.4 \
+ curve=pslinepos.asci,pslineneg.asci npair=901,901 curvewidth=2,2 curvecolor=black,black curvedash=3,3 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+file=f1min_$piter.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+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
+
+file=f1plus_$piter.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+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
+
+file=pgreen_$piter.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+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=$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
+mg1=`cat nepmg | awk '{print $1}'`
+rf1=`cat neprf | awk '{print $1}'`
+value=${value/[eE][+][0]/*10^}
+mg=${mg1/[eE][+][0]/*10^}
+rf=${rf1/[eE][+][0]/*10^}
+rm nep*
+scale=$(echo "scale=3; ($rf)/($mg)" | bc -l)
+scale=2.0
+echo $scale
+
+(suwind key=gx min=0 max=0 < referenceP_rp.su; \
+ suwind key=gx min=0 max=0 itmax=511 < pgreen_$piter.su | sugain scale=$scale ) | \
+ supsgraph hbox=6 wbox=2 labelsize=10 linegray=0.5,0.0 style=seismic \
+ lineon=1.0,1.0 lineoff=0.0,1.0 linewidth=1.0,1.0 x2beg=-$rf1 x2end=$rf1 > compare_$piter.eps
+
+done
+
+
+#special treatment of f1+ zero-iteration: which is zero, to make a nice gray plot (and not black)
+file=f1plus_001.su
+file_base=${file%.su}
+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 bclip=1 wclip=-1> $file_base.eps
+
diff --git a/marchenko3D/demo/oneD/epsModel.scr b/marchenko3D/demo/oneD/epsModel.scr
new file mode 100755
index 0000000000000000000000000000000000000000..5ae0b460f468bf00cb8804d8882d6fa35a4f7885
--- /dev/null
+++ b/marchenko3D/demo/oneD/epsModel.scr
@@ -0,0 +1,68 @@
+#!/bin/bash
+
+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
+
+#model
+supsimage hbox=4 wbox=6 labelsize=12 < model10_cp.su \
+ x1beg=0 x1end=1400.0 d1num=200 lstyle=vertright legend=1 threecolor=0 \
+ label1="depth (m)" label2="lateral distance (m)" \
+ curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black \
+ n1tic=5 x2beg=-2250 f2num=-2000 d2num=1000 x2end=2250 > model_cp_line.eps
+
+supsimage hbox=4 wbox=6 labelsize=12 < model10_ro.su \
+ x1beg=0 x1end=1400.0 d1num=200 lstyle=vertright legend=1 threecolor=0 \
+ label1="depth (m)" label2="lateral distance (m)" \
+ curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black \
+ n1tic=5 x2beg=-2250 f2num=-2000 d2num=1000 x2end=2250 > model_ro_line.eps
+
+#wavelet
+dt=0.0005
+supsgraph < wavefw.su \
+ labelsize=12 d1=$dt style=normal \
+ label1="time (s)" label2="amplitude" \
+ d1num=0.15 wbox=6 hbox=3 x1end=0.9 > wavefw.eps
+
+sufft < wavefw.su | suamp | sugain scale=$dt | supsgraph \
+ labelsize=12 style=normal \
+ label1="frequency (1/s)" label2="amplitude" \
+ d1num=10 wbox=6 hbox=3 x1end=125 x2end=1.1 > wavefw_freq.eps
+
+
+#shot record
+file=shot5_rp.su
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/3}'`
+suwind key=gx min=-2250000 max=2250000 < $file | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ label1="time (s)" label2="lateral distance (m)" \
+ n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > shotx0_rp.eps
+
+#Initial focusing operator
+file=iniFocus_rp.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/3}'`
+suwind key=gx min=-2250000 max=2250000 < $file | \
+ supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 \
+ label1="time (s)" label2="lateral distance (m)" \
+ n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+rm nep
diff --git a/marchenko3D/demo/oneD/figAppendix.scr b/marchenko3D/demo/oneD/figAppendix.scr
new file mode 100755
index 0000000000000000000000000000000000000000..295d2cf53c10926bd222f8bc05b310763bf2f7d1
--- /dev/null
+++ b/marchenko3D/demo/oneD/figAppendix.scr
@@ -0,0 +1,47 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin/:$PATH:
+
+file=iter_002.su
+file_base=${file%.su}
+
+ns=`surange < $file | grep ns | awk '{print $2}'`
+dtrcv=`surange < $file | grep dt | awk '{print $2/1000000.0}'`
+shift=$(echo "scale=4; ($dtrcv*($ns/2.0-1))" | bc -l)
+suzero < $file itmax=$ns | suaddnoise | sushw key=f1 a=0 > noise.su
+file_base=noise
+sumax < ${file_base}.su mode=abs outpar=nep
+clipiter=`cat nep | awk '{print $1/6}'`
+clipref=$clipiter
+
+#basop choice=shift shift=$shift file_in=$file file_out=${file_base}_t0.su
+
+for above in 0 1 -1 2 4
+do
+fmute file_mute=iniFocus_rp.su file_shot=${file_base}.su file_out=nep.su above=${above} shift=8 verbose=1 check=1 hw=4
+
+basop choice=shift shift=-$shift file_in=nep.su file_out=nep_t0.su
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < nep.su \
+ label1="time (s)" label2="lateral distance (m)" \
+ n1tic=2 d2=5 x1beg=0 d1num=0.5 \
+ curve=pslinepos.asci,pslineneg.asci npair=901,901 curvewidth=2,2 curvecolor=black,black curvedash=3,3 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_above${above}.eps
+done
+
+for shift in 0 20 -20
+do
+fmute file_mute=iniFocus_rp.su file_shot=${file_base}.su file_out=nep.su above=${above} shift=$shift verbose=1 check=1 hw=4
+mv pslinepos.asci pslinepos${shift}.asci
+done
+
+suzero < $file itmax=$ns | sushw key=f1 a=0 > zero.su
+sumax < iniFocus_rp.su mode=abs outpar=nep
+clipiter=`cat nep | awk '{print $1/6}'`
+clipref=$clipiter
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < iniFocus_rp.su \
+ label1="time (s)" label2="lateral distance (m)" \
+ n1tic=2 d2=5 x1beg=0 d1num=0.5 \
+ curve=pslinepos0.asci,pslinepos20.asci,pslinepos-20.asci npair=901,901,901 \
+ curvewidth=1,1,1 curvecolor=white,black,black curvedash=3,3,3 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > iniFocus_shifts.eps
+
diff --git a/marchenko3D/demo/oneD/initialFocus.scr b/marchenko3D/demo/oneD/initialFocus.scr
new file mode 100755
index 0000000000000000000000000000000000000000..4d4fd68aee89b203d0976c9bfa2a4ac18d2f4731
--- /dev/null
+++ b/marchenko3D/demo/oneD/initialFocus.scr
@@ -0,0 +1,39 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+dx=2.5
+dt=0.0005
+
+#the model upto 900 m depth, deeper reflections are not needed to model the direct transmission response
+makemod sizex=10000 sizez=1400 dx=$dx dz=$dx cp0=1800 ro0=1000 \
+ orig=-5000,0 file_base=modelup.su verbose=2 \
+ intt=def x=-5000,5000 z=400,400 poly=0 cp=2300 ro=3000 \
+ intt=def x=-5000,5000 z=700,700 poly=0 cp=2000 ro=1100
+
+makewave fp=25 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+export OMP_NUM_THREADS=1
+
+fdelmodc \
+ file_cp=modelup_cp.su ischeme=1 iorder=4 \
+ file_den=modelup_ro.su \
+ file_src=wave.su \
+ file_rcv=iniFocus.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.1 \
+ verbose=2 \
+ tmod=2.144 \
+ dxrcv=5 \
+ xrcv1=-2250 xrcv2=2250 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=900 \
+ ntaper=101 \
+ left=2 right=2 top=2 bottom=2
+
diff --git a/marchenko3D/demo/oneD/marchenko.scr b/marchenko3D/demo/oneD/marchenko.scr
new file mode 100755
index 0000000000000000000000000000000000000000..04e495d53dcb305750132b3fa55dc4606d444b46
--- /dev/null
+++ b/marchenko3D/demo/oneD/marchenko.scr
@@ -0,0 +1,42 @@
+#!/bin/bash -x
+
+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=8
+
+#apply the Marchenko algorithm
+marchenko file_shot=shotsdx5_rp.su file_tinv=p0plus.su nshots=901 verbose=2 \
+ 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
+
+exit
+
+#compare Green's funtions on Marhcenko and reference result
+suwind key=gx min=0 max=0 itmax=511 < pgreen.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}'`
+value=${value/[eE][+][0]/*10^}
+mg=${mg/[eE][+][0]/*10^}
+rf=${rf/[eE][+][0]/*10^}
+rm nep*
+scale=$(echo "scale=3; ($rf)/($mg)" | bc -l)
+echo $scale
+
+(suwind key=gx min=0 max=0 itmax=511 < pgreen.su | sugain scale=$scale; \
+ suwind key=gx min=0 max=0 < referenceP_rp.su) | suxgraph
+
+#suwind itmax=511 < pgreen.su > pgreen512.su
+#suop2 pgreen512.su referenceP_rp.su op=diff w2=1 w1=$scale > diffref.su
+
+# plot for convergence rate, the values in conv.txt are collected from the output of the marhenko program with verbose=2
+# marchenko: - iSyn 0: Ni at iteration 0 has energy 6.234892e+02; relative to N0 1.000000e+00
+#a2b < conv.txt | \
+#psgraph n=16 style=normal hbox=2 wbox=6 labelsize=10 \
+#label2='convergence rate' label1='iteration number' > convergence.eps
+
+# If guplot is installed: the same plot can also be produced by gnuplot this figure is used in the paper
+#gnuplot conv.gnp
diff --git a/marchenko3D/demo/oneD/marchenkoIter.scr b/marchenko3D/demo/oneD/marchenkoIter.scr
new file mode 100755
index 0000000000000000000000000000000000000000..401f97f7c2108e92e0ff5ca813d9fdfd2b4d183a
--- /dev/null
+++ b/marchenko3D/demo/oneD/marchenkoIter.scr
@@ -0,0 +1,21 @@
+#!/bin/bash -x
+
+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
+
+for (( iter=1; iter<=4; iter+=1 ))
+do
+echo "doing iteration $iter"
+piter=$(printf %03d $iter)
+
+#apply the Marchenko algorithm
+marchenko file_shot=shotsdx5_rp.su file_tinv=p0plus.su nshots=901 verbose=1 \
+ tap=0 ntap=41 niter=$iter hw=12 shift=8 smooth=5 \
+ file_green=pgreen_$piter.su file_iter=iter.su \
+ file_f1plus=f1plus_$piter.su file_f1min=f1min_$piter.su
+
+done
+
diff --git a/marchenko3D/demo/oneD/model.scr b/marchenko3D/demo/oneD/model.scr
new file mode 100755
index 0000000000000000000000000000000000000000..9165e8758862d4cef0b497ca585bea46336e9d13
--- /dev/null
+++ b/marchenko3D/demo/oneD/model.scr
@@ -0,0 +1,77 @@
+#!/bin/bash
+
+#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=1000 \
+ orig=-5000,0 file_base=model10.su verbose=2 \
+ intt=def x=-5000,5000 z=400,400 poly=0 cp=2300 ro=3000 \
+ intt=def x=-5000,5000 z=700,700 poly=0 cp=2000 ro=1100 \
+ intt=def x=-5000,5000 z=1100,1100 poly=0 cp=2500 ro=4000
+
+#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=4
+
+#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
+
+
diff --git a/marchenko3D/demo/oneD/p5all.scr b/marchenko3D/demo/oneD/p5all.scr
new file mode 100755
index 0000000000000000000000000000000000000000..333be5510ec6a203c098595abfdabe5cdba2466b
--- /dev/null
+++ b/marchenko3D/demo/oneD/p5all.scr
@@ -0,0 +1,31 @@
+#!/bin/bash -x
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+# Generate the full R matrix for a fixed spread geometry.
+
+dxshot=5000 # with scalco factor of 1000
+ishot=0
+nshots=901
+
+echo $1
+
+rm shotsdx5_rp.su
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -2250000 + ${ishot}*${dxshot} ))
+ (( tr1 = 901 - ${ishot} ))
+ (( tr2 = ${tr1} + 900 ))
+ 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,901 b=0,5000,0,0 j=0,901,0,0 | \
+ suchw key1=offset key2=gx key3=sx c=-1 d=1000 >> shotsdx5_rp.su
+
+done
+
diff --git a/marchenko3D/demo/oneD/referenceShot.scr b/marchenko3D/demo/oneD/referenceShot.scr
new file mode 100755
index 0000000000000000000000000000000000000000..b7a2b771341b3115d71bfebe2ec06e308846cbc6
--- /dev/null
+++ b/marchenko3D/demo/oneD/referenceShot.scr
@@ -0,0 +1,34 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+#Compute the reference Green's fucntion at x=0 z=900 m in the actual model
+dx=2.5
+dt=0.0005
+
+makewave fp=25 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+export OMP_NUM_THREADS=2
+
+fdelmodc \
+ file_cp=model10_cp.su ischeme=1 iorder=4 \
+ file_den=model10_ro.su \
+ file_src=wave.su \
+ file_rcv=referenceP.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.1 \
+ verbose=2 \
+ tmod=2.144 \
+ dxrcv=5.0 \
+ xrcv1=-2250 xrcv2=2250 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=900 \
+ ntaper=101 \
+ left=2 right=2 top=2 bottom=2
+
diff --git a/marchenko3D/demo/twoD/README b/marchenko3D/demo/twoD/README
new file mode 100644
index 0000000000000000000000000000000000000000..a4c7852f088a25f6f53418e9043f5cf5565d6bb4
--- /dev/null
+++ b/marchenko3D/demo/twoD/README
@@ -0,0 +1,10 @@
+Description of files:
+1a) model.scr computes the model
+1b) shots_slurm/pbs.scr creates the shots and submit jobs to slurm or PBS
+1c) check.scr after the jobs on shots_*.scr are finished checks if all shots are there
+2) direct.scr creates the direct arrival to be removed from the shots
+3) remove_direct.scr remove the direct wave from the shots
+4) initialFocus.scr model G_d the intitial focusing function => iniFocus_z1100_x0_rp.su
+5) referenceShot.scr creates the reference Green's function at focal point => referenceP_rp.su
+6) marchenko.scr perform the Marchenko scheme
+
diff --git a/marchenko3D/demo/twoD/backProp_makemovie.scr b/marchenko3D/demo/twoD/backProp_makemovie.scr
new file mode 100755
index 0000000000000000000000000000000000000000..d039c054b2c7fc73c3694c6e075bbbb495005b01
--- /dev/null
+++ b/marchenko3D/demo/twoD/backProp_makemovie.scr
@@ -0,0 +1,34 @@
+#!/bin/bash
+
+# $1 input-file of backpropf2_sp.su format
+
+export PATH=$HOME/OpenSource/bin/:$PATH:
+
+#rm prop_movie.su
+
+for fldr in $(seq 1 101);
+do
+ times=$(echo "scale=2; -0.01*(${fldr}-101)" | bc -l)
+ echo " Adding ${fldr}-th frame"
+ suwind key=fldr min=$fldr max=$fldr < backpropmar_sp.su > neg.su
+ (( fldr = 101+(101-$fldr) ))
+ suwind key=fldr min=$fldr max=$fldr < backpropmar_sp.su > pos.su
+ susum neg.su pos.su >prop_frame.su
+
+ if [ "$fldr" != "201" ]; then
+ cat prop_movie.su prop_frame.su >temp_movie.su
+ mv temp_movie.su prop_movie.su
+ else
+ cp prop_frame.su prop_movie.su
+ fi
+done
+
+suchw <prop_movie.su key1=ntr key2=fldr b=1001 >temp_movie.su
+
+susort <temp_movie.su -fldr tracf >prop_movie.su
+
+rm neg.su pos.su temp_movie.su
+
+n2=`surange <backpropmar_sp.su | grep tracf | awk '{print $3}'`
+
+suxmovie <prop_movie.su loop=1 height=500 width=1000 title=%g n2=$n2
diff --git a/marchenko3D/demo/twoD/backpropf2.scr b/marchenko3D/demo/twoD/backpropf2.scr
new file mode 100755
index 0000000000000000000000000000000000000000..0b7c281ba11e4fdfb8059547c0bdd3d15bc80861
--- /dev/null
+++ b/marchenko3D/demo/twoD/backpropf2.scr
@@ -0,0 +1,59 @@
+#!/bin/bash
+
+export PATH=$HOME/OpenSource/bin/:$PATH:
+
+dx=2.5
+dt=0.0005
+
+file_cp=syncl_cp.su
+file_ro=syncl_ro.su
+
+export OMP_NUM_THREADS=8
+
+# t=0 focal time is at 2.0445 seconds back=propagating (dtrcv*(ns/2-1)+dt)
+# shift f2.su such that t=0 is positioned in the middle of the time axis
+# the extra shift of 0.000250 is needed because of the staggered time implementation of the Finite Difference program.
+ns=`surange <f2.su | grep ns | awk '{print $2}'`
+dtrcv=`surange < f2.su | grep dt | awk '{print $2/1000000.0}'`
+suwind key=gx min=-2250000 max=2250000 itmax=1023 < f2.su > nep.su
+shift=$(echo "scale=6; ($dtrcv*(($ns)/2.0-1)+$dt)" | bc -l)
+echo $shift
+basop choice=shift shift=$shift file_in=nep.su verbose=1 > pplus.su
+
+# the f2.su is sampled with 4ms the FD program needs 0.5ms
+# time axis is interpolated by making use of FFT's: sinc interpolation
+#ftr1d file_in=pplus.su file_out=freq.su
+#sushw <freq.su key=nhs,dt a=8192,500 >fr.su
+#ftr1d file_in=fr.su n1=8194 file_out=pplusdt.su verbose=1
+
+midsnap=2.04
+
+#backpropagate f2.su and collect snapshots
+fdelmodc \
+ file_cp=$file_cp ischeme=1 iorder=4 \
+ file_den=$file_ro \
+ file_src=pplus.su \
+ file_rcv=backprop_f2_z1100.su \
+ dt=$dt \
+ grid_dir=0 \
+ src_type=1 \
+ src_injectionrate=1 \
+ src_orient=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.0 \
+ verbose=2 \
+ tmod=3.10 \
+ dxrcv=10.0 \
+ xrcv1=-3000 xrcv2=3000 \
+ zrcv1=1100 zrcv2=1100 \
+ zsrc=0 xsrc=0 \
+ npml=101 \
+ file_snap=backpropmar.su tsnap1=`echo "$midsnap-1" | bc -l` dtsnap=0.01 tsnap2=`echo "$midsnap+1" | bc -l` dxsnap=5 dzsnap=5 zsnap1=0 zsnap2=2000 xsnap1=-2250 xsnap2=2250 \
+ sna_type_vz=0 \
+ sna_type_p=1 \
+ left=2 right=2 top=2 bottom=2
+
+
diff --git a/marchenko3D/demo/twoD/check.scr b/marchenko3D/demo/twoD/check.scr
new file mode 100755
index 0000000000000000000000000000000000000000..f1a40feb5fd574e1f48f5c7a4acd5981e03e8aca
--- /dev/null
+++ b/marchenko3D/demo/twoD/check.scr
@@ -0,0 +1,24 @@
+#!/bin/bash
+
+dxshot=10
+ishot=0
+nshots=601
+zsrc=0
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+
+ file_rcv=shots/shots_${xsrc}_rp.su
+
+ if [ ! -e "$file_rcv" ]
+ then
+ echo $xsrc is missing
+ sbatch jobs/slurm_$ishot.job
+ fi
+
+ (( ishot = $ishot + 1))
+
+done
+
diff --git a/marchenko3D/demo/twoD/clean b/marchenko3D/demo/twoD/clean
new file mode 100755
index 0000000000000000000000000000000000000000..0d2611c04c3751b3c7d36314cb089b3fcdc6b864
--- /dev/null
+++ b/marchenko3D/demo/twoD/clean
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+rm *.su *.bin *.txt *.eps nep *.asci
+
diff --git a/marchenko3D/demo/twoD/direct.scr b/marchenko3D/demo/twoD/direct.scr
new file mode 100755
index 0000000000000000000000000000000000000000..48ef53e850d9ee80f8b91027b36fbcbd7d825037
--- /dev/null
+++ b/marchenko3D/demo/twoD/direct.scr
@@ -0,0 +1,35 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+dx=2.5
+dt=0.0005
+
+makemod sizex=12000 sizez=4000 dx=$dx dz=$dx cp0=1900 ro0=1200 \
+ orig=-6000,-1000 file_base=noContrast.su
+
+export OMP_NUM_THREADS=8
+makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3 scale=0
+
+fdelmodc \
+ file_cp=noContrast_cp.su ischeme=1 iorder=4 \
+ file_den=noContrast_ro.su \
+ file_src=wavefw.su \
+ file_rcv=direct.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ rec_delay=0.3 \
+ dtrcv=0.004 \
+ verbose=2 \
+ tmod=4.392 \
+ dxrcv=10.0 \
+ xrcv1=-6000 xrcv2=6000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=0 \
+ ntaper=200 \
+ left=2 right=2 top=2 bottom=2
+
diff --git a/marchenko3D/demo/twoD/initialFocus_pbs.scr b/marchenko3D/demo/twoD/initialFocus_pbs.scr
new file mode 100755
index 0000000000000000000000000000000000000000..eb4e0c1d48c1a8533905a0ba3a0f6092dc48897d
--- /dev/null
+++ b/marchenko3D/demo/twoD/initialFocus_pbs.scr
@@ -0,0 +1,77 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+dx=2.5
+dt=0.0005
+
+makemod sizex=6000 sizez=2000 dx=$dx dz=$dx cp0=1900 ro0=1200 \
+ orig=-3000,0 file_base=synclDown.su verbose=2 \
+ intt=def x=-3000,500,3000 z=195,195,195 poly=1 cp=1950 ro=3700 \
+ intt=def x=-3000,3000 z=600,600 poly=0 cp=2050 ro=1750 \
+ intt=def x=-3000,3000 z=680,680 poly=0 cp=2150 ro=2220 \
+ 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
+
+dxshot=10
+ishot=300
+nshots=301
+
+export OMP_NUM_THREADS=1
+mkdir -p shots
+mkdir -p jobs
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+ echo xsrc=$xsrc
+ file_rcv=shots/iniFocus_z1100_x${xsrc}.su
+
+cat << EOF > jobs/pbs_$ishot.job
+#!/bin/bash
+#
+#PBS -q medium
+#PBS -N mod_${xsrc}
+#PBS -j eo
+#PBS -m n
+#PBS -l nodes=1
+#PBS -V
+
+export PATH=\$HOME/src/OpenSource/bin:\$PATH:
+cd \$PBS_O_WORKDIR
+
+export OMP_NUM_THREADS=4
+
+fdelmodc \
+ file_cp=synclDown_cp.su ischeme=1 iorder=4 \
+ file_den=synclDown_ro.su \
+ file_src=wave.su \
+ file_rcv=$file_rcv \
+ src_type=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ rec_delay=0.1 \
+ dtrcv=0.004 \
+ verbose=2 \
+ tmod=2.100 \
+ dxrcv=10.0 \
+ xrcv1=-3000 xrcv2=3000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=$xsrc zsrc=1100 \
+ ntaper=200 \
+ left=2 right=2 top=2 bottom=2
+EOF
+
+ qsub jobs/pbs_$ishot.job
+
+ (( ishot = $ishot + 1))
+done
+
+
+
diff --git a/marchenko3D/demo/twoD/initialFocus_slurm.scr b/marchenko3D/demo/twoD/initialFocus_slurm.scr
new file mode 100755
index 0000000000000000000000000000000000000000..a94a11d1a159d2f739a004a7ce3c86b627088223
--- /dev/null
+++ b/marchenko3D/demo/twoD/initialFocus_slurm.scr
@@ -0,0 +1,75 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+dx=2.5
+dt=0.0005
+
+makemod sizex=6000 sizez=2000 dx=$dx dz=$dx cp0=1900 ro0=1200 \
+ orig=-3000,0 file_base=synclDown.su verbose=2 \
+ intt=def x=-3000,500,3000 z=195,195,195 poly=1 cp=1950 ro=3700 \
+ intt=def x=-3000,3000 z=600,600 poly=0 cp=2050 ro=1750 \
+ intt=def x=-3000,3000 z=680,680 poly=0 cp=2150 ro=2220 \
+ 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
+
+dxshot=10
+ishot=300
+nshots=301
+
+export OMP_NUM_THREADS=1
+mkdir -p shots
+mkdir -p jobs
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+ echo xsrc=$xsrc
+ file_rcv=shots/iniFocus_z1100_x${xsrc}.su
+
+cat << EOF > jobs/slurm_$ishot.job
+#!/bin/bash
+#
+#SBATCH -J mod_${xsrc}
+#SBATCH --cpus-per-task=4
+#SBATCH --ntasks=1
+#SBATCH --time=0:20:00
+
+export PATH=\$HOME/src/OpenSource/bin:\$PATH:
+cd \$SLURM_SUBMIT_DIR
+
+export OMP_NUM_THREADS=4
+
+fdelmodc \
+ file_cp=synclDown_cp.su ischeme=1 iorder=4 \
+ file_den=synclDown_ro.su \
+ file_src=wave.su \
+ file_rcv=$file_rcv \
+ src_type=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ rec_delay=0.1 \
+ dtrcv=0.004 \
+ verbose=2 \
+ tmod=2.100 \
+ dxrcv=10.0 \
+ xrcv1=-3000 xrcv2=3000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=$xsrc zsrc=1100 \
+ ntaper=200 \
+ left=2 right=2 top=2 bottom=2
+EOF
+
+ sbatch jobs/slurm_$ishot.job
+
+ (( ishot = $ishot + 1))
+done
+
+
+
diff --git a/marchenko3D/demo/twoD/marchenko.scr b/marchenko3D/demo/twoD/marchenko.scr
new file mode 100755
index 0000000000000000000000000000000000000000..fcbb45ccb77f50f142d07a4824e889363e309a95
--- /dev/null
+++ b/marchenko3D/demo/twoD/marchenko.scr
@@ -0,0 +1,30 @@
+#!/bin/bash -x
+
+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=shots/iniFocus_z1100_x0_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=pgreen.su file_gplus=Gplus0.su file_gmin=Gmin0.su \
+ file_f1plus=f1plus0.su file_f1min=f1min0.su file_f2=f2.su
+
+#compare Green's funtions on Marhcenko and reference result
+suwind key=gx min=0 max=0 itmax=511 < pgreen.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}'`
+value=${value/[eE][+][0]/*10^}
+mg=${mg/[eE][+][0]/*10^}
+rf=${rf/[eE][+][0]/*10^}
+rm nep*
+scale=$(echo "scale=3; ($rf)/($mg)" | bc -l)
+echo $scale
+
+(suwind key=gx min=0 max=0 itmax=511 < pgreen.su | sugain scale=$scale; \
+ suwind key=gx min=0 max=0 < referenceP_rp.su) | suxgraph
+
diff --git a/marchenko3D/demo/twoD/model.scr b/marchenko3D/demo/twoD/model.scr
new file mode 100755
index 0000000000000000000000000000000000000000..f41606ab0cddeb19ab11c3f5af5c5c16ec8d9e93
--- /dev/null
+++ b/marchenko3D/demo/twoD/model.scr
@@ -0,0 +1,82 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+dx=2.5
+dt=0.0005
+
+makemod sizex=6000 sizez=2000 dx=$dx dz=$dx cp0=1900 ro0=1200 \
+ orig=-3000,0 file_base=syncl.su verbose=2 \
+ intt=def x=-3000,500,3000 z=195,195,195 poly=1 cp=1950 ro=3700 \
+ intt=def x=-3000,3000 z=600,600 poly=0 cp=2050 ro=1750 \
+ intt=def x=-3000,3000 z=680,680 poly=0 cp=2150 ro=2220 \
+ 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 \
+ intt=def x=-3000,0,3000 z=1110,1110,1110 poly=0 cp=2300 ro=1950 \
+ intt=def x=-3000,3000 z=1180,1180 poly=0 cp=2480 ro=1820 \
+ intt=def x=-3000,0,3000 z=1290,1290,1370 poly=0 cp=2600 ro=2000 \
+ intt=def x=-3000,3000 z=1380,1380 poly=0 cp=2720 ro=2050 \
+ intt=def x=-3000,3000 z=1480,1480 poly=0 cp=2800 ro=1850
+
+exit
+
+#example FD modeling with model defined above
+makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3 scale=0
+
+export OMP_NUM_THREADS=4
+
+makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+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=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.0010 \
+ 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=200 \
+ tsnap1=3.1 tsnap2=2.5 dtsnap=0.1 \
+ left=2 right=2 top=2 bottom=2
+
+
+
+makemod sizex=6000 sizez=2000 dx=$dx dz=$dx cp0=1900 ro0=1200 \
+ orig=-3000,-1000 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=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.0010 \
+ 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=400 \
+ left=4 right=4 top=4 bottom=4
+
+sudiff shot_fd_rp.su shot_hom_fd_rp.su > shot_rp.su
+
+
diff --git a/marchenko3D/demo/twoD/referenceShot.scr b/marchenko3D/demo/twoD/referenceShot.scr
new file mode 100755
index 0000000000000000000000000000000000000000..4c015f6baa98c092e3210801e3a6caff68fd34ad
--- /dev/null
+++ b/marchenko3D/demo/twoD/referenceShot.scr
@@ -0,0 +1,38 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+#makewave w=fw fmin=0 flef=5 frig=80 fmax=100 dt=$dt file_out=wavefw.su nt=4096 t0=0.3
+
+dx=2.5
+dt=0.0005
+
+makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+export OMP_NUM_THREADS=8
+
+fdelmodc \
+ file_cp=syncl_cp.su ischeme=1 iorder=4 \
+ file_den=syncl_ro.su \
+ file_src=wave.su \
+ file_rcv=referenceP.su \
+ src_type=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_ud=1 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.1 \
+ verbose=2 \
+ tmod=2.144 \
+ dxrcv=10.0 \
+ xrcv1=-3000 xrcv2=3000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=1100 \
+ file_snap=backpropref.su tsnap1=0.1 dtsnap=0.010 tsnap2=2.100 \
+ dxsnap=10 dzsnap=10 zsnap1=0 zsnap2=2000 xsnap1=-2250 xsnap2=2250 sna_type_vz=0 \
+ ntaper=200 \
+ left=2 right=2 top=2 bottom=2
+
+
diff --git a/marchenko3D/demo/twoD/remove_direct.scr b/marchenko3D/demo/twoD/remove_direct.scr
new file mode 100755
index 0000000000000000000000000000000000000000..0881615c381aba1f42314c7cc226ba18edb05496
--- /dev/null
+++ b/marchenko3D/demo/twoD/remove_direct.scr
@@ -0,0 +1,33 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+dxshot=10
+ishot=0
+nshots=601
+
+rm shots/refl_rp.su
+
+while (( ishot < nshots ))
+do
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+ (( iishot = ${ishot}*${dxshot}/10 ))
+ (( tr1 = 601 - ${iishot} ))
+ (( tr2 = ${tr1} + 600 ))
+ echo xsrc=$xsrc tr1=$tr1 tr2=$tr2
+ suwind < direct_rp.su key=tracl min=$tr1 max=$tr2 > direct.su
+
+ file_rcv=shots/shots_${xsrc}_rp.su
+ suwind key=tracl min=1 max=601 < $file_rcv > shotz0.su
+
+ sudiff shotz0.su direct.su > refl.su
+
+ (( ishot = $ishot + 1))
+
+ sushw < refl.su key=fldr a=$ishot | \
+ suwind itmax=1023 >> shots/refl_rp.su
+
+done
+
+rm refl.su shotz0.su direct.su
+
diff --git a/marchenko3D/demo/twoD/shots_pbs.scr b/marchenko3D/demo/twoD/shots_pbs.scr
new file mode 100755
index 0000000000000000000000000000000000000000..4d1f22967a9946331f9ad903af522d324e9f956d
--- /dev/null
+++ b/marchenko3D/demo/twoD/shots_pbs.scr
@@ -0,0 +1,69 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+dt=0.0005
+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
+
+./model.scr
+
+mkdir -p shots
+mkdir -p jobs
+
+dxshot=10
+ishot=0
+nshots=601
+zsrc=0
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+
+ echo ishot=$ishot xsrc=$xsrc zsrc=$zsrc
+
+cat << EOF > jobs/pbs_$ishot.job
+#!/bin/bash
+#
+#PBS -q medium
+#PBS -N mod_${xsrc}
+#PBS -j eo
+#PBS -m n
+#PBS -l nodes=1
+#PBS -V
+
+export PATH=\$HOME/src/OpenSource/bin:\$PATH:
+cd \$PBS_O_WORKDIR
+
+export OMP_NUM_THREADS=4
+file_rcv=shots/shots_${xsrc}.su
+
+fdelmodc \
+ file_cp=syncl_cp.su ischeme=1 iorder=4 \
+ file_den=syncl_ro.su \
+ file_src=wavefw.su \
+ file_rcv=\$file_rcv \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ rec_delay=0.3 \
+ dtrcv=0.004 \
+ verbose=2 \
+ tmod=4.392 \
+ dxrcv=10.0 \
+ xrcv1=-3000 xrcv2=3000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=$xsrc zsrc=$zsrc \
+ ntaper=200 \
+ left=2 right=2 top=2 bottom=2
+EOF
+
+qsub jobs/pbs_$ishot.job
+
+ (( ishot = $ishot + 1))
+
+done
+
diff --git a/marchenko3D/demo/twoD/shots_slurm.scr b/marchenko3D/demo/twoD/shots_slurm.scr
new file mode 100755
index 0000000000000000000000000000000000000000..0aac6ca0c021699f37216c5af8d08e7028cffc4a
--- /dev/null
+++ b/marchenko3D/demo/twoD/shots_slurm.scr
@@ -0,0 +1,67 @@
+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+dt=0.0005
+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
+
+./model.scr
+
+mkdir -p shots
+mkdir -p jobs
+
+dxshot=10
+ishot=0
+nshots=601
+zsrc=0
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -3000 + ${ishot}*${dxshot} ))
+
+ echo ishot=$ishot xsrc=$xsrc zsrc=$zsrc
+
+cat << EOF > jobs/slurm_$ishot.job
+#!/bin/bash
+#
+#SBATCH -J mod_${xsrc}
+#SBATCH --cpus-per-task=4
+#SBATCH --ntasks=1
+#SBATCH --time=0:20:00
+
+export PATH=\$HOME/src/OpenSource/bin:\$PATH:
+cd \$SLURM_SUBMIT_DIR
+
+export OMP_NUM_THREADS=4
+file_rcv=shots/shots_${xsrc}.su
+
+fdelmodc \
+ file_cp=syncl_cp.su ischeme=1 iorder=4 \
+ file_den=syncl_ro.su \
+ file_src=wavefw.su \
+ file_rcv=\$file_rcv \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ rec_delay=0.3 \
+ dtrcv=0.004 \
+ verbose=2 \
+ tmod=4.392 \
+ dxrcv=10.0 \
+ xrcv1=-3000 xrcv2=3000 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=$xsrc zsrc=$zsrc \
+ ntaper=200 \
+ left=2 right=2 top=2 bottom=2
+EOF
+
+ sbatch jobs/slurm_$ishot.job
+
+ (( ishot = $ishot + 1))
+
+done
+
diff --git a/marchenko3D/docpkge.c b/marchenko3D/docpkge.c
new file mode 120000
index 0000000000000000000000000000000000000000..5384bb3801703c3f0db8fcc032235ca6130fa08b
--- /dev/null
+++ b/marchenko3D/docpkge.c
@@ -0,0 +1 @@
+../utils/docpkge.c
\ No newline at end of file
diff --git a/marchenko3D/fmute.c b/marchenko3D/fmute.c
new file mode 100644
index 0000000000000000000000000000000000000000..ba4f39acb407d3dacf414096dafc0b3ab67a2c8d
--- /dev/null
+++ b/marchenko3D/fmute.c
@@ -0,0 +1,370 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
+int readData(FILE *fp, float *data, segy *hdrs, int n1);
+int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
+int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
+void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *ixpos, int npos, int shift);
+double wallclock_time(void);
+
+/*********************** self documentation **********************/
+char *sdoc[] = {
+" ",
+" fmute - mute in time domain file_shot along curve of maximum amplitude in file_mute ",
+" ",
+" fmute file_shot= {file_mute=} [optional parameters]",
+" ",
+" Required parameters: ",
+" ",
+" file_mute= ................ input file with event that defines the mute line",
+" file_shot= ................ input data that is muted",
+" ",
+" Optional parameters: ",
+" ",
+" file_out= ................ output file",
+" above=0 .................. mute after(0), before(1) or around(2) the maximum times of file_mute",
+" .......................... options 4 is the inverse of 0 and -1 the inverse of 1",
+" shift=0 .................. number of points above(positive) / below(negative) maximum time for mute",
+" check=0 .................. plots muting window on top of file_mute: output file check.su",
+" scale=0 .................. scale data by dividing through maximum",
+" hw=15 .................... number of time samples to look up and down in next trace for maximum",
+" smooth=0 ................. number of points to smooth mute with cosine window",
+//" nxmax=512 ................ maximum number of traces in input file",
+//" ntmax=1024 ............... maximum number of samples/trace in input file",
+" verbose=0 ................ silent option; >0 display info",
+" ",
+" author : Jan Thorbecke : 2012 (janth@xs4all.nl)",
+" ",
+NULL};
+/**************** end self doc ***********************************/
+
+int main (int argc, char **argv)
+{
+ FILE *fp_in1, *fp_in2, *fp_out, *fp_chk, *fp_psline1, *fp_psline2;
+ int verbose, shift, k, nx1, nt1, nx2, nt2;
+ int ntmax, nxmax, ret, i, j, jmax, imax, above, check;
+ int size, ntraces, ngath, *maxval, hw, smooth;
+ int tstart, tend, scale, *xrcv;
+ float dt, d2, f1, f2, t0, t1, f1b, f2b, d1, d1b, d2b;
+ float w1, w2, dxrcv;
+ float *tmpdata, *tmpdata2, *costaper;
+ char *file_mute, *file_shot, *file_out;
+ float scl, sclsxgx, sclshot, xmin, xmax, tmax, lmax;
+ segy *hdrs_in1, *hdrs_in2;
+
+ t0 = wallclock_time();
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if(!getparstring("file_mute", &file_mute)) file_mute=NULL;
+ if(!getparstring("file_shot", &file_shot)) file_shot=NULL;
+ if(!getparstring("file_out", &file_out)) file_out=NULL;
+ if(!getparint("ntmax", &ntmax)) ntmax = 1024;
+ if(!getparint("nxmax", &nxmax)) nxmax = 512;
+ if(!getparint("above", &above)) above = 0;
+ if(!getparint("check", &check)) check = 0;
+ if(!getparint("scale", &scale)) scale = 0;
+ if(!getparint("hw", &hw)) hw = 15;
+ if(!getparint("smooth", &smooth)) smooth = 0;
+ if(!getparfloat("w1", &w1)) w1=1.0;
+ if(!getparfloat("w2", &w2)) w2=1.0;
+ if(!getparint("shift", &shift)) shift=0;
+ if(!getparint("verbose", &verbose)) verbose=0;
+
+/* Reading input data for file_mute */
+
+ if (file_mute != NULL) {
+ ngath = 1;
+ getFileInfo(file_mute, &nt1, &nx1, &ngath, &d1, &d2, &f1, &f2, &xmin, &xmax, &sclsxgx, &ntraces);
+
+ if (!getparint("ntmax", &ntmax)) ntmax = nt1;
+ if (!getparint("nxmax", &nxmax)) nxmax = nx1;
+ if (verbose>=2 && (ntmax!=nt1 || nxmax!=nx1))
+ vmess("dimensions overruled: %d x %d",ntmax,nxmax);
+ if(!getparfloat("dt", &dt)) dt=d1;
+
+ fp_in1 = fopen(file_mute, "r");
+ if (fp_in1 == NULL) verr("error on opening input file_mute=%s", file_mute);
+
+ size = ntmax * nxmax;
+ tmpdata = (float *)malloc(size*sizeof(float));
+ hdrs_in1 = (segy *) calloc(nxmax,sizeof(segy));
+
+ nx1 = readData(fp_in1, tmpdata, hdrs_in1, nt1);
+ if (nx1 == 0) {
+ fclose(fp_in1);
+ if (verbose) vmess("end of file_mute data reached");
+ }
+
+ if (verbose) {
+ disp_fileinfo(file_mute, nt1, nx1, f1, f2, dt, d2, hdrs_in1);
+ }
+ }
+
+/* Reading input data for file_shot */
+
+ ngath = 1;
+ getFileInfo(file_shot, &nt2, &nx2, &ngath, &d1b, &d2b, &f1b, &f2b, &xmin, &xmax, &sclshot, &ntraces);
+
+ if (!getparint("ntmax", &ntmax)) ntmax = nt2;
+ if (!getparint("nxmax", &nxmax)) nxmax = nx2;
+
+ size = ntmax * nxmax;
+ tmpdata2 = (float *)malloc(size*sizeof(float));
+ hdrs_in2 = (segy *) calloc(nxmax,sizeof(segy));
+
+ if (file_shot != NULL) fp_in2 = fopen(file_shot, "r");
+ else fp_in2=stdin;
+ if (fp_in2 == NULL) verr("error on opening input file_shot=%s", file_shot);
+
+ nx2 = readData(fp_in2, tmpdata2, hdrs_in2, nt2);
+ if (nx2 == 0) {
+ fclose(fp_in2);
+ if (verbose) vmess("end of file_shot data reached");
+ }
+ nt2 = hdrs_in2[0].ns;
+ f1b = hdrs_in2[0].f1;
+ f2b = hdrs_in2[0].f2;
+ d1b = (float)hdrs_in2[0].dt*1e-6;
+
+ if (verbose) {
+ disp_fileinfo(file_shot, nt2, nx2, f1b, f2b, d1b, d2b, hdrs_in2);
+ }
+
+ /* file_shot will be used as well to define the mute window */
+ if (file_mute == NULL) {
+ nx1=nx2;
+ nt1=nt2;
+ dt=d1b;
+ f1=f1b;
+ f2=f2b;
+ tmpdata = tmpdata2;
+ hdrs_in1 = hdrs_in2;
+ sclsxgx = sclshot;
+ }
+
+ if (verbose) vmess("sampling file_mute=%d, file_shot=%d", nt1, nt2);
+
+/*================ initializations ================*/
+
+ maxval = (int *)calloc(nx1,sizeof(int));
+ xrcv = (int *)calloc(nx1,sizeof(int));
+
+ if (file_out==NULL) fp_out = stdout;
+ else {
+ fp_out = fopen(file_out, "w+");
+ if (fp_out==NULL) verr("error on ceating output file");
+ }
+ if (check!=0){
+ fp_chk = fopen("check.su", "w+");
+ if (fp_chk==NULL) verr("error on ceating output file");
+ fp_psline1 = fopen("pslinepos.asci", "w+");
+ if (fp_psline1==NULL) verr("error on ceating output file");
+ fp_psline2 = fopen("pslineneg.asci", "w+");
+ if (fp_psline2==NULL) verr("error on ceating output file");
+
+ }
+ if (smooth) {
+ costaper = (float *)malloc(smooth*sizeof(float));
+ scl = M_PI/((float)smooth);
+ for (i=0; i<smooth; i++) {
+ costaper[i] = 0.5*(1.0+cos((i+1)*scl));
+/* fprintf(stderr,"costaper[%d]=%f\n",i,costaper[i]);*/
+ }
+ }
+
+/*================ loop over all shot records ================*/
+
+ k=1;
+ while (nx1 > 0) {
+ if (verbose) vmess("processing input gather %d", k);
+
+/*================ loop over all shot records ================*/
+
+ /* find consistent (one event) maximum related to maximum value */
+
+ /* find global maximum
+ xmax=0.0;
+ for (i = 0; i < nx1; i++) {
+ tmax=0.0;
+ jmax = 0;
+ for (j = 0; j < nt1; j++) {
+ lmax = fabs(tmpdata[i*nt1+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ if (lmax > xmax) {
+ imax = i;
+ xmax=lmax;
+ }
+ }
+ }
+ maxval[i] = jmax;
+ }
+ */
+
+ /* alternative find maximum at source position */
+ dxrcv = (hdrs_in1[nx1-1].gx - hdrs_in1[0].gx)*sclsxgx/(float)(nx1-1);
+ imax = NINT(((hdrs_in1[0].sx-hdrs_in1[0].gx)*sclsxgx)/dxrcv);
+ tmax=0.0;
+ jmax = 0;
+ xmax=0.0;
+ for (j = 0; j < nt1; j++) {
+ lmax = fabs(tmpdata[imax*nt1+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ if (lmax > xmax) {
+ xmax=lmax;
+ }
+ }
+ }
+ maxval[imax] = jmax;
+ if (verbose >= 3) vmess("Mute max at src-trace %d is sample %d", imax, maxval[imax]);
+
+ /* search forward */
+ for (i = imax+1; i < nx1; i++) {
+ tstart = MAX(0, (maxval[i-1]-hw));
+ tend = MIN(nt1-1, (maxval[i-1]+hw));
+ jmax=tstart;
+ tmax=0.0;
+ for(j = tstart; j <= tend; j++) {
+ lmax = fabs(tmpdata[i*nt1+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ }
+ }
+ maxval[i] = jmax;
+ }
+ /* search backward */
+ for (i = imax-1; i >=0; i--) {
+ tstart = MAX(0, (maxval[i+1]-hw));
+ tend = MIN(nt1-1, (maxval[i+1]+hw));
+ jmax=tstart;
+ tmax=0.0;
+ for(j = tstart; j <= tend; j++) {
+ lmax = fabs(tmpdata[i*nt1+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ }
+ }
+ maxval[i] = jmax;
+ }
+
+/* scale with maximum ampltiude */
+
+ if (scale==1) {
+ for (i = 0; i < nx2; i++) {
+ lmax = fabs(tmpdata2[i*nt2+maxval[i]]);
+ for (j = 0; j < nt2; j++) {
+ tmpdata2[i*nt2+j] = tmpdata2[i*nt2+j]/lmax;
+ }
+ }
+ }
+
+ for (i = 0; i < nx2; i++) xrcv[i] = i;
+
+/*================ apply mute window ================*/
+
+ applyMute(tmpdata2, maxval, smooth, above, 1, nx2, nt2, xrcv, nx2, shift);
+
+/*================ write result to output file ================*/
+
+ ret = writeData(fp_out, tmpdata2, hdrs_in2, nt2, nx2);
+ if (ret < 0 ) verr("error on writing output file.");
+
+ /* put mute window in file to check correctness of mute */
+ if (check !=0) {
+ for (i = 0; i < nx1; i++) {
+ jmax = maxval[i]-shift;
+ tmpdata[i*nt1+jmax] = 2*xmax;
+ }
+ if (above==0){
+ for (i = 0; i < nx1; i++) {
+ jmax = nt2-maxval[i]+shift;
+ tmpdata[i*nt1+jmax] = 2*xmax;
+ }
+ }
+ ret = writeData(fp_chk, tmpdata, hdrs_in1, nt1, nx1);
+ if (ret < 0 ) verr("error on writing check file.");
+ for (i=0; i<nx1; i++) {
+ jmax = maxval[i]-shift;
+ ret = fprintf(fp_psline1, "%.5f %.5f \n",jmax*dt,hdrs_in1[i].gx*sclshot);
+ jmax =-maxval[i]+shift;
+ ret = fprintf(fp_psline2, "%.5f %.5f \n",jmax*dt,hdrs_in1[i].gx*sclshot);
+ }
+ }
+
+/*================ Read next record for muting ================*/
+
+ if (file_mute != NULL) {
+ nx1 = readData(fp_in1, tmpdata, hdrs_in1, nt1);
+ if (nx1 == 0) {
+ fclose(fp_in1);
+ if (verbose) vmess("end of file_mute data reached");
+ fclose(fp_in2);
+ if (fp_out!=stdout) fclose(fp_out);
+ if (check!=0) fclose(fp_chk);
+ if (check!=0) {
+ fclose(fp_psline1);
+ fclose(fp_psline2);
+ }
+ break;
+ }
+ nt1 = (int)hdrs_in1[0].ns;
+ if (nt1 > ntmax) verr("n_samples (%d) greater than ntmax", nt1);
+ if (nx1 > nxmax) verr("n_traces (%d) greater than nxmax", nx1);
+ }
+
+/*================ Read next shot record(s) ================*/
+
+ nx2 = readData(fp_in2, tmpdata2, hdrs_in2, nt2);
+ if (nx2 == 0) {
+ if (verbose) vmess("end of file_shot data reached");
+ fclose(fp_in2);
+ break;
+ }
+ nt2 = (int)hdrs_in2[0].ns;
+ if (nt2 > ntmax) verr("n_samples (%d) greater than ntmax", nt2);
+ if (nx2 > nxmax) verr("n_traces (%d) greater than nxmax", nx2);
+
+ if (file_mute == NULL) {
+ nx1=nx2;
+ nt1=nt2;
+ hdrs_in1 = hdrs_in2;
+ tmpdata = tmpdata2;
+ }
+
+ k++;
+ }
+
+ t1 = wallclock_time();
+ if (verbose) vmess("Total CPU-time = %f",t1-t0);
+
+
+ return 0;
+}
+
diff --git a/marchenko3D/getFileInfo.c b/marchenko3D/getFileInfo.c
new file mode 120000
index 0000000000000000000000000000000000000000..ae38ea27f17697d65d7248c8e89038b632314182
--- /dev/null
+++ b/marchenko3D/getFileInfo.c
@@ -0,0 +1 @@
+../utils/getFileInfo.c
\ No newline at end of file
diff --git a/marchenko3D/getpars.c b/marchenko3D/getpars.c
new file mode 120000
index 0000000000000000000000000000000000000000..fa7dc3355428e8ea9013fafad6e319dde3a48ebb
--- /dev/null
+++ b/marchenko3D/getpars.c
@@ -0,0 +1 @@
+../utils/getpars.c
\ No newline at end of file
diff --git a/marchenko3D/marchenko.c b/marchenko3D/marchenko.c
new file mode 100644
index 0000000000000000000000000000000000000000..9812335d5eeca1883de2be2810350e053687f58b
--- /dev/null
+++ b/marchenko3D/marchenko.c
@@ -0,0 +1,1030 @@
+/*
+ * 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.
+ */
+
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+int omp_get_max_threads(void);
+int omp_get_num_threads(void);
+void omp_set_num_threads(int num_threads);
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+int compareInt(const void *a, const void *b)
+{ return (*(int *)a-*(int *)b); }
+
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int nshots,
+int nx, int nxs, float fxsb, float dxs, int ntfft, int mode, float scale, float tsq, float Q, float f0, int reci, int *nshots_r, int *isxcount, int *reci_xsrc, int *reci_xrcv, float *ixmask, int verbose);
+int readTinvData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nfoc, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
+int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn,
+float *zsyn, int *ixpos, int npos, int iter);
+
+void name_ext(char *filename, char *extension);
+
+void applyMute(float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *xrcvsyn, int npos, int shift);
+
+int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *ntraces);
+int readData(FILE *fp, float *data, segy *hdrs, int n1);
+int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
+int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
+double wallclock_time(void);
+
+void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int
+Nfoc, float *xrcv, float *xsrc, int *xnx, float fxse, float fxsb, float dxs, float dxsrc, float dx, int ntfft, int
+nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int npos, double *tfft, int *isxcount, int
+*reci_xsrc, int *reci_xrcv, float *ixmask, int verbose);
+
+void synthesisPosistions(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);
+
+int linearsearch(int *array, size_t N, int value);
+
+/*********************** self documentation **********************/
+char *sdoc[] = {
+" ",
+" MARCHENKO - Iterative Green's function and focusing functions retrieval",
+" ",
+" marchenko file_tinv= file_shot= [optional parameters]",
+" ",
+" Required parameters: ",
+" ",
+" file_tinv= ............... direct arrival from focal point: G_d",
+" file_shot= ............... Reflection response: R",
+" ",
+" Optional parameters: ",
+" ",
+" INTEGRATION ",
+" tap=0 .................... lateral taper focusing(1), shot(2) or both(3)",
+" ntap=0 ................... number of taper points at boundaries",
+" fmin=0 ................... minimum frequency in the Fourier transform",
+" fmax=70 .................. maximum frequency in the Fourier transform",
+" MARCHENKO ITERATIONS ",
+" niter=10 ................. number of iterations",
+" MUTE-WINDOW ",
+" above=0 .................. mute above(1), around(0) or below(-1) the first travel times of file_tinv",
+" 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",
+" REFLECTION RESPONSE CORRECTION ",
+" tsq=0.0 .................. scale factor n for t^n for true amplitude recovery",
+" Q=0.0 .......,............ Q correction factor",
+" f0=0.0 ................... ... for Q correction factor",
+" scale=2 .................. scale factor of R for summation of Ni with G_d",
+" pad=0 .................... amount of samples to pad the reflection series",
+" reci=0 ................... 1; add receivers as shots 2; only use receivers as shot positions",
+" countmin=0 ............... 0.3*nxrcv; minumum number of reciprocal traces for a contribution",
+" OUTPUT DEFINITION ",
+" file_green= .............. output file with full Green function(s)",
+" file_gplus= .............. output file with G+ ",
+" file_gmin= ............... output file with G- ",
+" file_f1plus= ............. output file with f1+ ",
+" file_f1min= .............. output file with f1- ",
+" file_f2= ................. output file with f2 (=p+) ",
+" file_pplus= .............. output file with p+ ",
+" file_pmin= ............... output file with p- ",
+" file_iter= ............... output file with -Ni(-t) for each iteration",
+" verbose=0 ................ silent option; >0 displays info",
+" ",
+" ",
+" author : Jan Thorbecke : 2016 (j.w.thorbecke@tudelft.nl)",
+" ",
+NULL};
+/**************** end self doc ***********************************/
+
+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 nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
+ int reci, countmin, mode, n2out, verbose, ntfft;
+ int iter, niter, tracf, *muteW;
+ int hw, smooth, above, shift, *ixpos, npos, ix;
+ int nshots_r, *isxcount, *reci_xsrc, *reci_xrcv;
+ float fmin, fmax, *tapersh, *tapersy, fxf, dxf, *xsrc, *xrcv, *zsyn, *zsrc, *xrcvsyn;
+ double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi, energyN0;
+ float d1, d2, f1, f2, fxsb, fxse, ft, fx, *xsyn, dxsrc;
+ float *green, *f2p, *pmin, *G_d, dt, dx, dxs, scl, mem;
+ float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus;
+ float xmin, xmax, scale, tsq, Q, f0;
+ float *ixmask;
+ complex *Refl, *Fop;
+ char *file_tinv, *file_shot, *file_green, *file_iter;
+ char *file_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_pmin;
+ segy *hdrs_out;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ tsyn = tread = tfft = tcopy = 0.0;
+ t0 = wallclock_time();
+
+ if (!getparstring("file_shot", &file_shot)) file_shot = NULL;
+ if (!getparstring("file_tinv", &file_tinv)) file_tinv = NULL;
+ if (!getparstring("file_f1plus", &file_f1plus)) file_f1plus = NULL;
+ if (!getparstring("file_f1min", &file_f1min)) file_f1min = NULL;
+ if (!getparstring("file_gplus", &file_gplus)) file_gplus = NULL;
+ if (!getparstring("file_gmin", &file_gmin)) file_gmin = NULL;
+ if (!getparstring("file_pplus", &file_f2)) file_f2 = NULL;
+ if (!getparstring("file_f2", &file_f2)) file_f2 = NULL;
+ if (!getparstring("file_pmin", &file_pmin)) file_pmin = NULL;
+ if (!getparstring("file_iter", &file_iter)) file_iter = NULL;
+ if (!getparint("verbose", &verbose)) verbose = 0;
+ if (file_tinv == NULL && file_shot == NULL)
+ verr("file_tinv and file_shot cannot be both input pipe");
+ if (!getparstring("file_green", &file_green)) {
+ if (verbose) vwarn("parameter file_green not found, assume pipe");
+ file_green = NULL;
+ }
+ if (!getparfloat("fmin", &fmin)) fmin = 0.0;
+ if (!getparfloat("fmax", &fmax)) fmax = 70.0;
+ if (!getparint("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 (!getparint("tap", &tap)) tap = 0;
+ if (!getparint("ntap", &ntap)) ntap = 0;
+ if (!getparint("pad", &pad)) pad = 0;
+
+ if(!getparint("niter", &niter)) niter = 10;
+ if(!getparint("hw", &hw)) hw = 15;
+ if(!getparint("smooth", &smooth)) smooth = 5;
+ if(!getparint("above", &above)) above = 0;
+ if(!getparint("shift", &shift)) shift=12;
+
+ if (reci && ntap) vwarn("tapering influences the reciprocal result");
+
+/*================ Reading info about shot and initial operator sizes ================*/
+
+ ngath = 0; /* setting ngath=0 scans all traces; n2 contains maximum traces/gather */
+ ret = getFileInfo(file_tinv, &n1, &n2, &ngath, &d1, &d2, &f1, &f2, &xmin, &xmax, &scl, &ntraces);
+ Nfoc = ngath;
+ nxs = n2;
+ nts = n1;
+ dxs = d2;
+ fxsb = f2;
+
+ ngath = 0; /* setting ngath=0 scans all traces; nx contains maximum traces/gather */
+ ret = getFileInfo(file_shot, &nt, &nx, &ngath, &d1, &dx, &ft, &fx, &xmin, &xmax, &scl, &ntraces);
+ nshots = ngath;
+ assert (nxs >= nshots);
+
+ if (!getparfloat("dt", &dt)) dt = d1;
+
+ ntfft = optncr(MAX(nt+pad, nts+pad));
+ nfreq = ntfft/2+1;
+ nw_low = (int)MIN((fmin*ntfft*dt), nfreq-1);
+ nw_low = MAX(nw_low, 1);
+ nw_high = MIN((int)(fmax*ntfft*dt), nfreq-1);
+ nw = nw_high - nw_low + 1;
+ scl = 1.0/((float)ntfft);
+ if (!getparint("countmin", &countmin)) countmin = 0.3*nx;
+
+/*================ Allocating all data arrays ================*/
+
+ Fop = (complex *)calloc(nxs*nw*Nfoc,sizeof(complex));
+ green = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ f2p = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ pmin = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ f1plus = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ f1min = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ iRN = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ Ni = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ G_d = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ muteW = (int *)calloc(Nfoc*nxs,sizeof(int));
+ trace = (float *)malloc(ntfft*sizeof(float));
+ tapersy = (float *)malloc(nxs*sizeof(float));
+ xrcvsyn = (float *)calloc(Nfoc*nxs,sizeof(float)); // x-rcv postions of focal points
+ xsyn = (float *)malloc(Nfoc*sizeof(float)); // x-src position of focal points
+ zsyn = (float *)malloc(Nfoc*sizeof(float)); // z-src position of focal points
+ xnxsyn = (int *)calloc(Nfoc,sizeof(int)); // number of traces per focal point
+ ixpos = (int *)calloc(nxs,sizeof(int)); // x-position of source of shot in G_d domain (nxs with dxs)
+
+ Refl = (complex *)malloc(nw*nx*nshots*sizeof(complex));
+ tapersh = (float *)malloc(nx*sizeof(float));
+ xrcv = (float *)calloc(nshots*nx,sizeof(float)); // x-rcv postions of shots
+ xsrc = (float *)calloc(nshots,sizeof(float)); //x-src position of shots
+ zsrc = (float *)calloc(nshots,sizeof(float)); // z-src position of shots
+ xnx = (int *)calloc(nshots,sizeof(int)); // number of traces per shot
+
+ if (reci!=0) {
+ reci_xsrc = (int *)malloc((nxs*nxs)*sizeof(int));
+ reci_xrcv = (int *)malloc((nxs*nxs)*sizeof(int));
+ isxcount = (int *)calloc(nxs,sizeof(int));
+ ixmask = (float *)calloc(nxs,sizeof(float));
+ }
+
+/*================ Read and define mute window based on focusing operator(s) ================*/
+/* G_d = p_0^+ = G_d (-t) ~ Tinv */
+
+ mode=-1; /* apply complex conjugate to read in data */
+ readTinvData(file_tinv, xrcvsyn, xsyn, zsyn, xnxsyn, Nfoc, nxs, ntfft,
+ mode, muteW, G_d, hw, verbose);
+ /* reading data added zero's to the number of time samples to be the same as ntfft */
+ nts = ntfft;
+
+ /* define tapers to taper edges of acquisition */
+ if (tap == 1 || tap == 3) {
+ for (j = 0; j < ntap; j++)
+ tapersy[j] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
+ for (j = ntap; j < nxs-ntap; j++)
+ tapersy[j] = 1.0;
+ for (j = nxs-ntap; j < nxs; j++)
+ tapersy[j] =(cos(PI*(j-(nxs-ntap))/ntap)+1)/2.0;
+ }
+ else {
+ for (j = 0; j < nxs; j++) tapersy[j] = 1.0;
+ }
+ if (tap == 1 || tap == 3) {
+ if (verbose) vmess("Taper for operator applied ntap=%d", ntap);
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < nxs; i++) {
+ for (j = 0; j < nts; j++) {
+ G_d[l*nxs*nts+i*nts+j] *= tapersy[i];
+ }
+ }
+ }
+ }
+
+ /* check consistency of header values */
+ if (xrcvsyn[0] != 0 || xrcvsyn[1] != 0 ) fxsb = xrcvsyn[0];
+ fxse = fxsb + (float)(nxs-1)*dxs;
+ dxf = (xrcvsyn[nxs-1] - xrcvsyn[0])/(float)(nxs-1);
+ if (NINT(dxs*1e3) != NINT(fabs(dxf)*1e3)) {
+ vmess("dx in hdr.d1 (%.3f) and hdr.gx (%.3f) not equal",d2, dxf);
+ if (dxf != 0) dxs = fabs(dxf);
+ vmess("dx in operator => %f", dxs);
+ }
+
+/*================ Reading shot records ================*/
+
+ mode=1;
+ readShotData(file_shot, xrcv, xsrc, zsrc, xnx, Refl, nw, nw_low, nshots, nx, nxs, fxsb, dxs, ntfft,
+ mode, scale, tsq, Q, f0, reci, &nshots_r, isxcount, reci_xsrc, reci_xrcv, ixmask, verbose);
+
+ tapersh = (float *)malloc(nx*sizeof(float));
+ if (tap == 2 || tap == 3) {
+ for (j = 0; j < ntap; j++)
+ tapersh[j] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
+ for (j = ntap; j < nx-ntap; j++)
+ tapersh[j] = 1.0;
+ for (j = nx-ntap; j < nx; j++)
+ tapersh[j] =(cos(PI*(j-(nx-ntap))/ntap)+1)/2.0;
+ }
+ else {
+ for (j = 0; j < nx; j++) tapersh[j] = 1.0;
+ }
+ if (tap == 2 || tap == 3) {
+ if (verbose) vmess("Taper for shots applied ntap=%d", ntap);
+ for (l = 0; l < nshots; l++) {
+ for (j = 1; j < nw; j++) {
+ for (i = 0; i < nx; i++) {
+ Refl[l*nx*nw+j*nx+i].r *= tapersh[i];
+ Refl[l*nx*nw+j*nx+i].i *= tapersh[i];
+ }
+ }
+ }
+ }
+ free(tapersh);
+
+ /* check consistency of header values */
+ fxf = xsrc[0];
+ if (nx > 1) dxf = (xrcv[nx-1] - xrcv[0])/(float)(nx-1);
+ else dxf = d2;
+ if (NINT(dx*1e3) != NINT(fabs(dxf)*1e3)) {
+ vmess("dx in hdr.d1 (%.3f) and hdr.gx (%.3f) not equal",dx, dxf);
+ if (dxf != 0) dx = fabs(dxf);
+ else verr("gx hdrs not set");
+ vmess("dx used => %f", dx);
+ }
+
+ dxsrc = (float)xsrc[1] - xsrc[0];
+ if (dxsrc == 0) {
+ vwarn("sx hdrs are not filled in!!");
+ dxsrc = dx;
+ }
+
+/*================ Check the size of the files ================*/
+
+ if (NINT(dxsrc/dx)*dx != NINT(dxsrc)) {
+ vwarn("source (%.2f) and receiver step (%.2f) don't match",dxsrc,dx);
+ if (reci == 2) vwarn("step used from operator (%.2f) ",dxs);
+ }
+ di = NINT(dxf/dxs);
+ if ((NINT(di*dxs) != NINT(dxf)) && verbose)
+ vwarn("dx in receiver (%.2f) and operator (%.2f) don't match",dx,dxs);
+ if (nt != nts)
+ vmess("Time samples in shot (%d) and focusing operator (%d) are not equal",nt, nts);
+ if (verbose) {
+ vmess("Number of focusing operators = %d", Nfoc);
+ vmess("Number of receivers in focusop = %d", nxs);
+ vmess("number of shots = %d", nshots);
+ vmess("number of receiver/shot = %d", nx);
+ vmess("first model position = %.2f", fxsb);
+ vmess("last model position = %.2f", fxse);
+ vmess("first source position fxf = %.2f", fxf);
+ vmess("source distance dxsrc = %.2f", dxsrc);
+ vmess("last source position = %.2f", fxf+(nshots-1)*dxsrc);
+ vmess("receiver distance dxf = %.2f", dxf);
+ vmess("direction of increasing traces = %d", di);
+ vmess("number of time samples (nt,nts) = %d (%d,%d)", ntfft, nt, nts);
+ vmess("time sampling = %e ", dt);
+ if (file_green != NULL) vmess("Green output file = %s ", file_green);
+ if (file_gmin != NULL) vmess("Gmin output file = %s ", file_gmin);
+ if (file_gplus != NULL) vmess("Gplus output file = %s ", file_gplus);
+ if (file_pmin != NULL) vmess("Pmin output file = %s ", file_pmin);
+ if (file_f2 != NULL) vmess("f2 (=pplus) output file = %s ", file_f2);
+ if (file_f1min != NULL) vmess("f1min output file = %s ", file_f1min);
+ if (file_f1plus != NULL)vmess("f1plus output file = %s ", file_f1plus);
+ if (file_iter != NULL) vmess("Iterations output file = %s ", file_iter);
+ }
+
+/*================ initializations ================*/
+
+ if (reci) n2out = nxs;
+ else n2out = nshots;
+ mem = Nfoc*n2out*ntfft*sizeof(float)/1048576.0;
+ if (verbose) {
+ vmess("number of output traces = %d", n2out);
+ vmess("number of output samples = %d", ntfft);
+ vmess("Size of output data/file = %.1f MB", mem);
+ }
+
+
+ /* dry-run of synthesis to get all x-positions calcalated by the integration */
+ synthesisPosistions(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("synthesisPosistions: nshots=%d npos=%d", nshots, npos);
+ }
+
+/*================ set variables for output data ================*/
+
+ n1 = nts; n2 = n2out;
+ f1 = ft; f2 = fxsb+dxs*ixpos[0];
+ d1 = dt;
+ if (reci == 0) d2 = dxsrc; // distance between sources in R
+ else if (reci == 1) d2 = dxs; // distance between traces in G_d
+ else if (reci == 2) d2 = dx; // distance between receivers in R
+
+ hdrs_out = (segy *) calloc(n2,sizeof(segy));
+ if (hdrs_out == NULL) verr("allocation for hdrs_out");
+ size = nxs*nts;
+
+ for (i = 0; i < n2; i++) {
+ hdrs_out[i].ns = n1;
+ hdrs_out[i].trid = 1;
+ hdrs_out[i].dt = dt*1000000;
+ hdrs_out[i].f1 = f1;
+ hdrs_out[i].f2 = f2;
+ hdrs_out[i].d1 = d1;
+ hdrs_out[i].d2 = d2;
+ hdrs_out[i].trwf = n2out;
+ hdrs_out[i].scalco = -1000;
+ hdrs_out[i].gx = NINT(1000*(f2+i*d2));
+ hdrs_out[i].scalel = -1000;
+ hdrs_out[i].tracl = i+1;
+ }
+ t1 = wallclock_time();
+ tread = t1-t0;
+
+/*================ initialization ================*/
+
+ memcpy(Ni, G_d, Nfoc*nxs*ntfft*sizeof(float));
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ ix = ixpos[i]; /* select the traces that have an output trace after integration */
+ f2p[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts+j];
+ f1plus[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts+j];
+ for (j = 1; j < nts; j++) {
+ f2p[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts+j];
+ f1plus[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts+j];
+ }
+ }
+ }
+
+/*================ start Marchenko iterations ================*/
+
+ for (iter=0; iter<niter; iter++) {
+
+ t2 = wallclock_time();
+
+/*================ construction of Ni(-t) = - \int R(x,t) Ni(t) ================*/
+
+ synthesis(Refl, Fop, Ni, iRN, nx, nt, nxs, nts, dt, xsyn, Nfoc,
+ xrcv, xsrc, xnx, fxse, fxsb, dxs, dxsrc, dx, ntfft, nw, nw_low, nw_high, mode,
+ reci, nshots, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv, ixmask, verbose);
+
+ t3 = wallclock_time();
+ tsyn += t3 - t2;
+
+ if (file_iter != NULL) {
+ writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, iter);
+ }
+ /* N_k(x,t) = -N_(k-1)(x,-t) */
+ /* p0^-(x,t) += iRN = (R * T_d^inv)(t) */
+ for (l = 0; l < Nfoc; l++) {
+ energyNi = 0.0;
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ ix = ixpos[i];
+ Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+ix*nts+j];
+ pmin[l*nxs*nts+i*nts+j] += iRN[l*nxs*nts+ix*nts+j];
+ energyNi += iRN[l*nxs*nts+ix*nts+j]*iRN[l*nxs*nts+ix*nts+j];
+ for (j = 1; j < nts; j++) {
+ Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+ix*nts+nts-j];
+ pmin[l*nxs*nts+i*nts+j] += iRN[l*nxs*nts+ix*nts+j];
+ energyNi += iRN[l*nxs*nts+ix*nts+j]*iRN[l*nxs*nts+ix*nts+j];
+ }
+ }
+ if (iter==0) energyN0 = energyNi;
+ if (verbose >=2) vmess(" - iSyn %d: Ni at iteration %d has energy %e; relative to N0 %e", l, iter, sqrt(energyNi),
+sqrt(energyNi/energyN0));
+ }
+
+ /* apply mute window based on times of direct arrival (in muteW) */
+ applyMute(Ni, muteW, smooth, above, Nfoc, nxs, nts, ixpos, npos, shift);
+
+ /* update f2 */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ f2p[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ f2p[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
+ }
+ }
+ }
+
+ if (iter % 2 == 0) { /* even iterations update: => f_1^-(t) */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j];
+ }
+ }
+ }
+ }
+ else {/* odd iterations update: => f_1^+(t) */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
+ }
+ }
+ }
+ }
+
+ t2 = wallclock_time();
+ tcopy += t2 - t3;
+
+ if (verbose) vmess("*** Iteration %d finished ***", iter);
+
+ } /* end of iterations */
+
+ free(Ni);
+ free(G_d);
+
+ /* compute full Green's function G = int R * f2(t) + f2(-t) = Pplus + Pmin */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ /* set green to zero if mute-window exceeds nt/2 */
+ if (muteW[l*nxs+ixpos[i]] >= nts/2) {
+ memset(&green[l*nxs*nts+i*nts],0, sizeof(float)*nt);
+ continue;
+ }
+ green[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+j] + pmin[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ green[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+nts-j] + pmin[l*nxs*nts+i*nts+j];
+ }
+ }
+ }
+
+ /* compute upgoing Green's function G^+,- */
+ if (file_gmin != NULL) {
+ Gmin = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+
+ /* use f1+ as operator on R in frequency domain */
+ mode=1;
+ synthesis(Refl, Fop, f1plus, iRN, nx, nt, nxs, nts, dt, xsyn, Nfoc,
+ xrcv, xsrc, xnx, fxse, fxsb, dxs, dxsrc, dx, ntfft, nw, nw_low, nw_high, mode,
+ reci, nshots, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv, ixmask, verbose);
+
+ /* compute upgoing Green's G^-,+ */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j=0;
+ ix = ixpos[i];
+ Gmin[l*nxs*nts+i*nts+j] = iRN[l*nxs*nts+ix*nts+j] - f1min[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ Gmin[l*nxs*nts+i*nts+j] = iRN[l*nxs*nts+ix*nts+j] - f1min[l*nxs*nts+i*nts+j];
+ }
+ }
+ }
+ /* Apply mute with window for Gmin */
+ applyMute(Gmin, muteW, smooth, 1, Nfoc, nxs, nts, ixpos, npos, shift);
+ } /* end if Gmin */
+
+ /* compute downgoing Green's function G^+,+ */
+ if (file_gplus != NULL) {
+ Gplus = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+
+ /* use f1-(*) as operator on R in frequency domain */
+ mode=-1;
+ synthesis(Refl, Fop, f1min, iRN, nx, nt, nxs, nts, dt, xsyn, Nfoc,
+ xrcv, xsrc, xnx, fxse, fxsb, dxs, dxsrc, dx, ntfft, nw, nw_low, nw_high, mode,
+ reci, nshots, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv, ixmask, verbose);
+
+ /* compute downgoing Green's G^+,+ */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j=0;
+ ix = ixpos[i];
+ Gplus[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+ix*nts+j] + f1plus[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ Gplus[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+ix*nts+j] + f1plus[l*nxs*nts+i*nts+nts-j];
+ }
+ }
+ }
+ } /* end if Gplus */
+
+ t2 = wallclock_time();
+ if (verbose) {
+ vmess("Total CPU-time marchenko = %.3f", t2-t0);
+ vmess("with CPU-time synthesis = %.3f", tsyn);
+ vmess("with CPU-time copy array = %.3f", tcopy);
+ vmess(" CPU-time fft data = %.3f", tfft);
+ vmess("and CPU-time read data = %.3f", tread);
+ }
+
+/*================ write output files ================*/
+
+
+ fp_out = fopen(file_green, "w+");
+ if (fp_out==NULL) verr("error on creating output file %s", file_green);
+ if (file_gmin != NULL) {
+ fp_gmin = fopen(file_gmin, "w+");
+ if (fp_gmin==NULL) verr("error on creating output file %s", file_gmin);
+ }
+ if (file_gplus != NULL) {
+ fp_gplus = fopen(file_gplus, "w+");
+ if (fp_gplus==NULL) verr("error on creating output file %s", file_gplus);
+ }
+ if (file_f2 != NULL) {
+ fp_f2 = fopen(file_f2, "w+");
+ if (fp_f2==NULL) verr("error on creating output file %s", file_f2);
+ }
+ if (file_pmin != NULL) {
+ fp_pmin = fopen(file_pmin, "w+");
+ if (fp_pmin==NULL) verr("error on creating output file %s", file_pmin);
+ }
+ if (file_f1plus != NULL) {
+ fp_f1plus = fopen(file_f1plus, "w+");
+ if (fp_f1plus==NULL) verr("error on creating output file %s", file_f1plus);
+ }
+ if (file_f1min != NULL) {
+ fp_f1min = fopen(file_f1min, "w+");
+ if (fp_f1min==NULL) verr("error on creating output file %s", file_f1min);
+ }
+
+
+ tracf = 1;
+ for (l = 0; l < Nfoc; l++) {
+ if (reci) f2 = fxsb;
+ else f2 = fxf;
+
+ for (i = 0; i < n2; i++) {
+ hdrs_out[i].fldr = l+1;
+ hdrs_out[i].sx = NINT(xsyn[l]*1000);
+ hdrs_out[i].offset = (long)NINT((f2+i*d2) - xsyn[l]);
+ hdrs_out[i].tracf = tracf++;
+ hdrs_out[i].selev = NINT(zsyn[l]*1000);
+ hdrs_out[i].sdepth = NINT(-zsyn[l]*1000);
+ hdrs_out[i].f1 = f1;
+ }
+
+ ret = writeData(fp_out, (float *)&green[l*size], hdrs_out, n1, n2);
+ if (ret < 0 ) verr("error on writing output file.");
+
+ if (file_gmin != NULL) {
+ ret = writeData(fp_gmin, (float *)&Gmin[l*size], hdrs_out, n1, n2);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ if (file_gplus != NULL) {
+ ret = writeData(fp_gplus, (float *)&Gplus[l*size], hdrs_out, n1, n2);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ if (file_f2 != NULL) {
+ ret = writeData(fp_f2, (float *)&f2p[l*size], hdrs_out, n1, n2);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ if (file_pmin != NULL) {
+ ret = writeData(fp_pmin, (float *)&pmin[l*size], hdrs_out, n1, n2);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ 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];
+ }
+ }
+ ret = writeData(fp_f1plus, (float *)&f1plus[l*size], hdrs_out, n1, n2);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ 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];
+ }
+ }
+ ret = writeData(fp_f1min, (float *)&f1min[l*size], hdrs_out, n1, n2);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ }
+ ret = fclose(fp_out);
+ if (file_gplus != NULL) {ret += fclose(fp_gplus);}
+ if (file_gmin != NULL) {ret += fclose(fp_gmin);}
+ if (file_f2 != NULL) {ret += fclose(fp_f2);}
+ if (file_pmin != NULL) {ret += fclose(fp_pmin);}
+ if (file_f1plus != NULL) {ret += fclose(fp_f1plus);}
+ if (file_f1min != NULL) {ret += fclose(fp_f1min);}
+ if (ret < 0) verr("err %d on closing output file",ret);
+
+ if (verbose) {
+ t1 = wallclock_time();
+ vmess("and CPU-time write data = %.3f", t1-t2);
+ }
+
+/*================ free memory ================*/
+
+ free(hdrs_out);
+ free(tapersy);
+
+ exit(0);
+}
+
+
+/*================ Convolution and Integration ================*/
+
+void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int
+Nfoc, float *xrcv, float *xsrc, int *xnx, float fxse, float fxsb, float dxs, float dxsrc, float dx, int ntfft, int
+nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int npos, double *tfft, int *isxcount, int
+*reci_xsrc, int *reci_xrcv, float *ixmask, int verbose)
+{
+ int nfreq, size, inx;
+ float scl;
+ int i, j, l, m, iw, ix, k, ixsrc, il, ik;
+ float *rtrace, idxs;
+ complex *sum, *ctrace;
+ int npe;
+ static int first=1, *ixrcv;
+ static double t0, t1, t;
+
+ 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 virtual source positions (Nfoc) */
+ if (npe > Nfoc) {
+ vmess("Number of OpenMP threads set to %d (was %d)", Nfoc, npe);
+ omp_set_num_threads(Nfoc);
+ }
+#endif
+
+ t0 = wallclock_time();
+
+ /* reset output data to zero */
+ memset(&iRN[0], 0, Nfoc*nxs*nts*sizeof(float));
+ ctrace = (complex *)calloc(ntfft,sizeof(complex));
+
+ if (!first) {
+ /* 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].r, 0, nxs*nw*2*sizeof(float));
+ for (i = 0; i < npos; i++) {
+ rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
+ ix = ixpos[i];
+ 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;
+ }
+ }
+ }
+ }
+ else { /* only for first call to synthesis using all nxs traces in G_d */
+ /* 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 */
+ memset(&Fop[l*nxs*nw].r, 0, nxs*nw*2*sizeof(float));
+ for (i = 0; i < nxs; i++) {
+ rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
+ for (iw=0; iw<nw; iw++) {
+ Fop[l*nxs*nw+iw*nxs+i].r = ctrace[nw_low+iw].r;
+ Fop[l*nxs*nw+iw*nxs+i].i = mode*ctrace[nw_low+iw].i;
+ }
+ }
+ }
+ 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);
+ }
+ }
+ }
+ free(ctrace);
+ t1 = wallclock_time();
+ *tfft += t1 - t0;
+
+/* Loop over total number of shots */
+ if (reci == 0 || reci == 1) {
+ for (k=0; k<nshots; k++) {
+ if ((xsrc[k] < 0.999*fxsb) || (xsrc[k] > 1.001*fxse)) continue;
+ ixsrc = NINT((xsrc[k] - fxsb)/dxs);
+ inx = xnx[k]; /* number of traces per shot */
+
+/*================ SYNTHESIS ================*/
+
+#pragma omp parallel default(none) \
+ shared(iRN, dx, npe, nw, verbose) \
+ shared(Refl, Nfoc, reci, xrcv, xsrc, xsyn, fxsb, fxse, nxs, dxs) \
+ shared(nx, dxsrc, inx, k, nfreq, nw_low, nw_high) \
+ shared(Fop, size, nts, ntfft, scl, ixrcv, ixsrc) \
+ private(l, ix, j, m, i, sum, rtrace)
+{ /* start of parallel region */
+ sum = (complex *)malloc(nfreq*sizeof(complex));
+ rtrace = (float *)calloc(ntfft,sizeof(float));
+
+#pragma omp for schedule(guided,1)
+ 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++)
+ iRN[l*size+ixsrc*nts+j] += rtrace[j]*scl*dx;
+
+ } /* end of parallel Nfoc loop */
+ free(sum);
+ free(rtrace);
+
+#ifdef _OPENMP
+#pragma omp single
+ npe = omp_get_num_threads();
+#endif
+} /* end of parallel region */
+
+ if (verbose>4) vmess("*** Shot gather %d processed ***", k);
+
+ } /* end of nshots (k) loop */
+ } /* end of if reci */
+
+/* if reciprocal traces are enabled start a new loop over reciprocal shot positions */
+ if (reci != 0) {
+ for (k=0; k<nxs; k++) {
+ if (isxcount[k] == 0) continue;
+ ixsrc = k;
+ inx = isxcount[ixsrc]; /* number of traces per reciprocal shot */
+
+#pragma omp parallel default(none) \
+ shared(iRN, dx, nw, verbose) \
+ shared(Refl, Nfoc, reci, xrcv, xsrc, xsyn, fxsb, fxse, nxs, dxs) \
+ shared(nx, dxsrc, inx, k, nfreq, nw_low, nw_high) \
+ shared(reci_xrcv, reci_xsrc, ixmask) \
+ shared(Fop, size, nts, ntfft, scl, ixrcv, ixsrc) \
+ private(l, ix, j, m, i, sum, rtrace, ik, il)
+{ /* start of parallel region */
+ sum = (complex *)malloc(nfreq*sizeof(complex));
+ rtrace = (float *)calloc(ntfft,sizeof(float));
+
+#pragma omp for schedule(guided,1)
+ 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++)
+ iRN[l*size+ixsrc*nts+j] = ixmask[ixsrc]*(iRN[l*size+ixsrc*nts+j]+rtrace[j]*scl*dxsrc);
+
+ } /* end of parallel Nfoc loop */
+
+ free(sum);
+ free(rtrace);
+
+ } /* end of parallel region */
+
+ } /* end of reciprocal shots (k) loop */
+ } /* end of if reci */
+
+ t = wallclock_time() - t0;
+ if (verbose) {
+ vmess("OMP: parallel region = %f seconds (%d threads)", t, npe);
+ }
+
+ return;
+}
+
+void synthesisPosistions(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)
+{
+ int i, j, l, ixsrc, ixrcv, dosrc, k, *count;
+ float x0, x1;
+
+ count = (int *)calloc(nxs,sizeof(int)); // number of traces that contribute to the integration over x
+
+/*================ SYNTHESIS ================*/
+
+ for (l = 0; l < 1; l++) { /* assuming all focal operators cover the same lateral area */
+// for (l = 0; l < Nfoc; l++) {
+ *npos=0;
+
+ if (reci == 0 || reci == 1) {
+ for (k=0; k<nshots; k++) {
+
+ ixsrc = NINT((xsrc[k] - fxsb)/dxs);
+ if (verbose>=3) {
+ vmess("source position: %.2f in operator %d", xsrc[k], ixsrc);
+ vmess("receiver positions: %.2f <--> %.2f", xrcv[k*nx+0], xrcv[k*nx+nx-1]);
+ vmess("focal point positions: %.2f <--> %.2f", fxsb, fxse);
+ }
+
+ if ((NINT(xsrc[k]-fxse) > 0) || (NINT(xrcv[k*nx+nx-1]-fxse) > 0) ||
+ (NINT(xrcv[k*nx+nx-1]-fxsb) < 0) || (NINT(xsrc[k]-fxsb) < 0) ||
+ (NINT(xrcv[k*nx+0]-fxsb) < 0) || (NINT(xrcv[k*nx+0]-fxse) > 0) ) {
+ vwarn("source/receiver positions are outside synthesis aperture");
+ vmess("xsrc = %.2f xrcv_1 = %.2f xrvc_N = %.2f", xsrc[k], xrcv[k*nx+0], xrcv[k*nx+nx-1]);
+ vmess("source position: %.2f in operator %d", xsrc[k], ixsrc);
+ vmess("receiver positions: %.2f <--> %.2f", xrcv[k*nx+0], xrcv[k*nx+nx-1]);
+ vmess("focal point positions: %.2f <--> %.2f", fxsb, fxse);
+ }
+
+ if ( (xsrc[k] >= 0.999*fxsb) && (xsrc[k] <= 1.001*fxse) ) {
+ j = linearsearch(ixpos, *npos, ixsrc);
+ if (j < *npos) { /* the position (at j) is already included */
+ count[j] += xnx[k];
+ }
+ else { /* add new postion */
+ ixpos[*npos]=ixsrc;
+ count[*npos] += xnx[k];
+ *npos += 1;
+ }
+ // vmess("source position %d is inside synthesis model %f *npos=%d count=%d", k, xsrc[k], *npos, count[*npos]);
+ }
+
+ } /* end of nshots (k) loop */
+ } /* end of reci branch */
+
+ /* if reci=1 or reci=2 source-receive reciprocity is used and new (reciprocal-)sources are added */
+ if (reci != 0) {
+ for (k=0; k<nxs; k++) { /* check count in total number of shots added by reciprocity */
+ if (isxcount[k] >= countmin) {
+ j = linearsearch(ixpos, *npos, k);
+ if (j < *npos) { /* the position (at j) is already included */
+ count[j] += isxcount[k];
+ }
+ else { /* add new postion */
+ ixpos[*npos]=k;
+ count[*npos] += isxcount[k];
+ *npos += 1;
+ }
+ }
+ else {
+ isxcount[k] = 0;
+ }
+ }
+ } /* end of reci branch */
+ } /* end of Nfoc loop */
+
+ if (verbose>=4) {
+ for (j=0; j < *npos; j++) {
+ vmess("ixpos[%d] = %d count=%d", j, ixpos[j], count[j]);
+ }
+ }
+ free(count);
+
+/* sort ixpos into increasing values */
+ qsort(ixpos, *npos, sizeof(int), compareInt);
+
+
+ return;
+}
+
+int linearsearch(int *array, size_t N, int value)
+{
+ int j;
+/* Check is position is already in array */
+ j = 0;
+ while (j < N && value != array[j]) {
+ j++;
+ }
+ return j;
+}
+
+/*
+void update(float *field, float *term, int Nfoc, int nx, int nt, int reverse, int ixpos)
+{
+ int i, j, l, ix;
+
+ if (reverse) {
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ Ni[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];
+ }
+ }
+ }
+ }
+ else {
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ Ni[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];
+ }
+ }
+ }
+ }
+ return;
+}
+*/
diff --git a/marchenko3D/name_ext.c b/marchenko3D/name_ext.c
new file mode 120000
index 0000000000000000000000000000000000000000..83ac1f8ddf2ec6a316557877ae7db38720a5ca53
--- /dev/null
+++ b/marchenko3D/name_ext.c
@@ -0,0 +1 @@
+../utils/name_ext.c
\ No newline at end of file
diff --git a/marchenko3D/par.h b/marchenko3D/par.h
new file mode 120000
index 0000000000000000000000000000000000000000..0fa273cea748f9ead16e0e231201941174a3dd46
--- /dev/null
+++ b/marchenko3D/par.h
@@ -0,0 +1 @@
+../utils/par.h
\ No newline at end of file
diff --git a/marchenko3D/readData.c b/marchenko3D/readData.c
new file mode 120000
index 0000000000000000000000000000000000000000..af43798573495d45a669aacf2dfe5d1094834bf8
--- /dev/null
+++ b/marchenko3D/readData.c
@@ -0,0 +1 @@
+../utils/readData.c
\ No newline at end of file
diff --git a/marchenko3D/readShotData.c b/marchenko3D/readShotData.c
new file mode 100644
index 0000000000000000000000000000000000000000..a619799113de01c77af1a667cb060e3de055731b
--- /dev/null
+++ b/marchenko3D/readShotData.c
@@ -0,0 +1,207 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "segy.h"
+#include <assert.h>
+
+typedef struct { /* complex number */
+ float r,i;
+} complex;
+
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+
+int optncr(int n);
+void cc1fft(complex *data, int n, int sign);
+void rc1fft(float *rdata, complex *cdata, int n, int sign);
+
+int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int nshots,
+int nx, int nxs, float fxsb, float dxs, int ntfft, int mode, float scale, float tsq, float Q, float f0, int reci, int *nshots_r, int *isxcount, int *reci_xsrc, int *reci_xrcv, float *ixmask, int verbose)
+{
+ FILE *fp;
+ segy hdr;
+ size_t nread;
+ int fldr_shot, sx_shot, itrace, one_shot, igath, iw;
+ int end_of_file, nt;
+ int *isx, *igx, k, l, m, j, nreci;
+ int samercv, samesrc, nxrk, nxrm, maxtraces, ixsrc;
+ float scl, scel, *trace, dt;
+ complex *ctrace;
+
+ /* Reading first header */
+
+ if (filename == NULL) fp = stdin;
+ else fp = fopen( filename, "r" );
+ if ( fp == NULL ) {
+ fprintf(stderr,"input file %s has an error\n", filename);
+ perror("error in opening file: ");
+ fflush(stderr);
+ return -1;
+ }
+
+ fseek(fp, 0, SEEK_SET);
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+ if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
+ else if (hdr.scalco == 0) scl = 1.0;
+ else scl = hdr.scalco;
+ if (hdr.scalel < 0) scel = 1.0/fabs(hdr.scalel);
+ else if (hdr.scalel == 0) scel = 1.0;
+ else scel = hdr.scalel;
+
+ fseek(fp, 0, SEEK_SET);
+
+ nt = hdr.ns;
+ dt = hdr.dt/(1E6);
+
+ trace = (float *)calloc(ntfft,sizeof(float));
+ ctrace = (complex *)malloc(ntfft*sizeof(complex));
+ isx = (int *)malloc((nx*nshots)*sizeof(int));
+ igx = (int *)malloc((nx*nshots)*sizeof(int));
+
+ end_of_file = 0;
+ one_shot = 1;
+ igath = 0;
+
+ /* Read shots in file */
+
+ while (!end_of_file) {
+
+ /* start reading data (shot records) */
+ itrace = 0;
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) { /* no more data in file */
+ break;
+ }
+
+/* ToDo Don't store the traces that are not in the aperture */
+/*
+ if ( (NINT(sx_shot*scl-fxse) > 0) || (NINT(-fxsb) > 0) ) {
+ vwarn("source positions are outside synthesis aperture");
+ vmess("xsrc = %.2f", xsrc[k], xrcv[k*nx+0], xrcv[k*nx+nx-1]);
+ }
+*/
+
+ sx_shot = hdr.sx;
+ fldr_shot = hdr.fldr;
+ isx[igath] = sx_shot;
+ xsrc[igath] = sx_shot*scl;
+ zsrc[igath] = hdr.selev*scel;
+ xnx[igath]=0;
+ while (one_shot) {
+ igx[igath*nx+itrace] = hdr.gx;
+ xrcv[igath*nx+itrace] = hdr.gx*scl;
+
+ nread = fread( trace, sizeof(float), nt, fp );
+ assert (nread == hdr.ns);
+
+ /* True Amplitude Recovery */
+ if (tsq != 0.0) {
+ for (iw=0; iw<nt; iw++) {
+ trace[iw] *= powf(dt*iw,tsq);
+ }
+ }
+
+ /* Q-correction */
+ if (Q != 0.0 && f0 != 0.0) {
+ for (iw=0; iw<nt; iw++) {
+ trace[iw] *= expf(((dt*iw)*M_PI*f0)/Q);
+ }
+ }
+
+ /* transform to frequency domain */
+ if (ntfft > hdr.ns)
+ memset( &trace[nt-1], 0, sizeof(float)*(ntfft-nt) );
+
+ rc1fft(trace,ctrace,ntfft,-1);
+ for (iw=0; iw<nw; iw++) {
+ cdata[igath*nx*nw+iw*nx+itrace].r = scale*ctrace[nw_low+iw].r;
+ cdata[igath*nx*nw+iw*nx+itrace].i = scale*mode*ctrace[nw_low+iw].i;
+ }
+ itrace++;
+ xnx[igath]+=1;
+
+ /* read next hdr of next trace */
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) {
+ one_shot = 0;
+ end_of_file = 1;
+ break;
+ }
+ if ((sx_shot != hdr.sx) || (fldr_shot != hdr.fldr) ) break;
+ }
+ if (verbose>2) {
+ vmess("finished reading shot %d (%d) with %d traces",sx_shot,igath,itrace);
+ }
+
+ if (itrace != 0) { /* end of shot record */
+ fseek( fp, -TRCBYTES, SEEK_CUR );
+ igath++;
+ }
+ else {
+ end_of_file = 1;
+ }
+ }
+
+ free(ctrace);
+ free(trace);
+
+/* if reci=1 or reci=2 source-receive reciprocity is used and traces are added */
+
+ if (reci != 0) {
+ for (k=0; k<nxs; k++) ixmask[k] = 1.0;
+ for (k=0; k<nshots; k++) {
+ ixsrc = NINT((xsrc[k] - fxsb)/dxs);
+ nxrk = xnx[k];
+ for (l=0; l<nxrk; l++) {
+ samercv = 0;
+ samesrc = 0;
+ for (m=0; m<nshots; m++) {
+ if (igx[k*nx+l] == isx[m] && reci == 1) { // receiver position already present as source position m
+ nxrm = xnx[m];
+ for (j=0; j<nxrm; j++) { // check if receiver l with source k is also present in shot m
+ if (isx[k] == igx[m*nx+j]) { // shot k with receiver l already known as receiver j in shot m: same data
+ samercv = 1;
+ break;
+ }
+ }
+ if (samercv == 0) { // source k of receiver l -> accept trace as new receiver position for source l
+ ixsrc = NINT((xrcv[k*nx+l] - fxsb)/dxs);
+ if ((ixsrc >= 0) && (ixsrc < nxs)) {
+ reci_xrcv[ixsrc*nxs+isxcount[ixsrc]] = k;
+ reci_xsrc[ixsrc*nxs+isxcount[ixsrc]] = l;
+ isxcount[ixsrc] += 1;
+ if (reci==1) ixmask[ixsrc] = 0.5; // traces are added to already existing traces and must be scaled
+ }
+ }
+ samesrc = 1;
+ break;
+ }
+ }
+ if (samesrc == 0) { // receiver l with source k -> accept trace as new source position l with receiver k
+ //fprintf(stderr,"not a samesrc for receiver l=%d for source k=%d\n", l,k);
+ ixsrc = NINT((xrcv[k*nx+l] - fxsb)/dxs);
+ if ((ixsrc >= 0) && (ixsrc < nxs)) { // is this correct or should k and l be reversed: rcv=l src=k
+ reci_xrcv[ixsrc*nxs+isxcount[ixsrc]] = k;
+ reci_xsrc[ixsrc*nxs+isxcount[ixsrc]] = l;
+ isxcount[ixsrc] += 1;
+ }
+ }
+ }
+ }
+ nreci = 0;
+ for (k=0; k<nxs; k++) { // count total number of shots added by reciprocity
+ if (isxcount[k] != 0) {
+ maxtraces = MAX(maxtraces,isxcount[k]);
+ nreci++;
+ if (verbose>1) vmess("reciprocal receiver at %f (%d) has %d sources contributing", k, k*dxs+fxsb, isxcount[k]);
+ }
+ }
+ *nshots_r = nreci;
+ }
+
+ return 0;
+}
+
+
diff --git a/marchenko3D/readTinvData.c b/marchenko3D/readTinvData.c
new file mode 100644
index 0000000000000000000000000000000000000000..028d4c21b4ddef73fe08edc7fcec0f1400b6d518
--- /dev/null
+++ b/marchenko3D/readTinvData.c
@@ -0,0 +1,240 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "segy.h"
+#include <assert.h>
+
+typedef struct { /* complex number */
+ float r,i;
+} complex;
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+void findShotInMute(float *xrcvMute, float xrcvShot, int nxs, int *imute);
+
+int readTinvData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nfoc, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose)
+{
+ FILE *fp;
+ segy hdr;
+ size_t nread;
+ int fldr_shot, sx_shot, itrace, one_shot, ig, isyn, i, j;
+ int end_of_file, nt, gx0, gx1;
+ int nx1, jmax, imax, tstart, tend;
+ float xmax, tmax, lmax;
+ float scl, scel, *trace, dxrcv;
+ complex *ctrace;
+
+ /* Reading first header */
+
+ if (filename == NULL) fp = stdin;
+ else fp = fopen( filename, "r" );
+ if ( fp == NULL ) {
+ fprintf(stderr,"input file %s has an error\n", filename);
+ perror("error in opening file: ");
+ fflush(stderr);
+ return -1;
+ }
+
+ fseek(fp, 0, SEEK_SET);
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+ if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
+ else if (hdr.scalco == 0) scl = 1.0;
+ else scl = hdr.scalco;
+ if (hdr.scalel < 0) scel = 1.0/fabs(hdr.scalel);
+ else if (hdr.scalel == 0) scel = 1.0;
+ else scel = hdr.scalel;
+ fseek(fp, 0, SEEK_SET);
+
+ nt = hdr.ns;
+ trace = (float *)calloc(ntfft,sizeof(float));
+ ctrace = (complex *)malloc(ntfft*sizeof(complex));
+
+ end_of_file = 0;
+ one_shot = 1;
+ isyn = 0;
+
+ /* Read shots in file */
+
+ while (!end_of_file) {
+
+ /* start reading data (shot records) */
+ itrace = 0;
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) { /* no more data in file */
+ break;
+ }
+
+ sx_shot = hdr.sx;
+ fldr_shot = hdr.fldr;
+ gx0 = hdr.gx;
+ xsrc[isyn] = sx_shot*scl;
+ zsrc[isyn] = hdr.selev*scel;
+ xnx[isyn] = 0;
+ ig = isyn*nx*ntfft;
+ while (one_shot) {
+ xrcv[isyn*nx+itrace] = hdr.gx*scl;
+ nread = fread( trace, sizeof(float), nt, fp );
+ assert (nread == hdr.ns);
+
+ /* copy trace to data array */
+ memcpy( &tinv[ig+itrace*ntfft], trace, nt*sizeof(float));
+
+ gx1 = hdr.gx;
+ itrace++;
+
+ /* read next hdr of next trace */
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) {
+ one_shot = 0;
+ end_of_file = 1;
+ break;
+ }
+ if ((sx_shot != hdr.sx) || (fldr_shot != hdr.fldr) ) break;
+ }
+ if (verbose>2) {
+ fprintf(stderr,"finished reading shot %d (%d) with %d traces\n",sx_shot,isyn,itrace);
+ //disp_fileinfo(filename, nt, xnx[isyn], hdr.f1, xrcv[isyn*nxm], d1, d2, &hdr);
+ }
+
+ /* look for maximum in shot record to define mute window */
+ /* find consistent (one event) maximum related to maximum value */
+ nx1 = itrace;
+ xnx[isyn]=nx1;
+ /* find global maximum
+ xmax=0.0;
+ for (i = 0; i < nx1; i++) {
+ tmax=0.0;
+ jmax = 0;
+ for (j = 0; j < nt; j++) {
+ lmax = fabs(tinv[ig+i*ntfft+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ if (lmax > xmax) {
+ imax = i;
+ xmax=lmax;
+ }
+ }
+ }
+ maxval[isyn*nx+i] = jmax;
+ }
+ */
+
+ /* alternative find maximum at source position */
+ dxrcv = (gx1 - gx0)*scl/(float)(nx1-1);
+ imax = NINT(((sx_shot-gx0)*scl)/dxrcv);
+ tmax=0.0;
+ jmax = 0;
+ for (j = 0; j < nt; j++) {
+ lmax = fabs(tinv[ig+imax*ntfft+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ if (lmax > xmax) {
+ xmax=lmax;
+ }
+ }
+ }
+ maxval[isyn*nx+imax] = jmax;
+ if (verbose >= 3) vmess("Mute max at src-trace %d is sample %d", imax, maxval[imax]);
+
+ /* search forward in trace direction from maximum in file */
+ for (i = imax+1; i < nx1; i++) {
+ tstart = MAX(0, (maxval[isyn*nx+i-1]-hw));
+ tend = MIN(nt-1, (maxval[isyn*nx+i-1]+hw));
+ jmax=tstart;
+ tmax=0.0;
+ for(j = tstart; j <= tend; j++) {
+ lmax = fabs(tinv[ig+i*ntfft+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ }
+ }
+ maxval[isyn*nx+i] = jmax;
+ }
+ /* search backward in trace direction from maximum in file */
+ for (i = imax-1; i >=0; i--) {
+ tstart = MAX(0, (maxval[isyn*nx+i+1]-hw));
+ tend = MIN(nt-1, (maxval[isyn*nx+i+1]+hw));
+ jmax=tstart;
+ tmax=0.0;
+ for(j = tstart; j <= tend; j++) {
+ lmax = fabs(tinv[ig+i*ntfft+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ }
+ }
+ maxval[isyn*nx+i] = jmax;
+ }
+
+ if (itrace != 0) { /* end of shot record, but not end-of-file */
+ fseek( fp, -TRCBYTES, SEEK_CUR );
+ isyn++;
+ }
+ else {
+ end_of_file = 1;
+ }
+
+ /* copy trace to data array for mode=-1 */
+ /* time reverse trace */
+ if (mode==-1) {
+ for (i = 0; i < nx1; i++) {
+ memcpy( trace, &tinv[ig+i*ntfft], ntfft*sizeof(float));
+ j=0;
+ tinv[ig+i*ntfft+j] = trace[j];
+ for (j=1; j<ntfft; j++) tinv[ig+i*ntfft+ntfft-j] = trace[j];
+ }
+ }
+ }
+
+ free(ctrace);
+ free(trace);
+
+ return 0;
+}
+
+
+/* simple sort algorithm */
+void findShotInMute(float *xrcvMute, float xrcvShot, int nxs, int *imute)
+{
+ int i, sign;
+ float diff1, diff2;
+
+ *imute=0;
+
+ if (xrcvMute[0] < xrcvMute[1]) sign = 1;
+ else sign = -1;
+
+ if (sign == 1) {
+ i = 0;
+ while (xrcvMute[i] < xrcvShot && i < nxs) {
+ i++;
+ }
+ /* i is now position larger than xrcvShot */
+ }
+ else {
+ i = 0;
+ while (xrcvMute[i] > xrcvShot && i < nxs) {
+ i++;
+ }
+ /* i is now position smaller than xrcvShot */
+ }
+
+ diff1 = fabsf(xrcvMute[i]-xrcvShot);
+ diff2 = fabsf(xrcvMute[i-1]-xrcvShot);
+ if (diff1 < diff2) *imute = i;
+ else *imute = i-1;
+
+ return;
+}
+
diff --git a/marchenko3D/segy.h b/marchenko3D/segy.h
new file mode 120000
index 0000000000000000000000000000000000000000..8eaebbdccb4f6c015d1ed7d5d3d227bb22ca55c8
--- /dev/null
+++ b/marchenko3D/segy.h
@@ -0,0 +1 @@
+../utils/segy.h
\ No newline at end of file
diff --git a/marchenko3D/verbosepkg.c b/marchenko3D/verbosepkg.c
new file mode 120000
index 0000000000000000000000000000000000000000..248253edebc2c7b207e139ecf16b68b318f057df
--- /dev/null
+++ b/marchenko3D/verbosepkg.c
@@ -0,0 +1 @@
+../utils/verbosepkg.c
\ No newline at end of file
diff --git a/marchenko3D/wallclock_time.c b/marchenko3D/wallclock_time.c
new file mode 120000
index 0000000000000000000000000000000000000000..0bd00b4c2878f007a8dc398f0af7c7cb44f50717
--- /dev/null
+++ b/marchenko3D/wallclock_time.c
@@ -0,0 +1 @@
+../utils/wallclock_time.c
\ No newline at end of file
diff --git a/marchenko3D/writeData.c b/marchenko3D/writeData.c
new file mode 120000
index 0000000000000000000000000000000000000000..b761f28f24545fb2e550406a85b67afe0410db7e
--- /dev/null
+++ b/marchenko3D/writeData.c
@@ -0,0 +1 @@
+../utils/writeData.c
\ No newline at end of file
diff --git a/marchenko3D/writeDataIter.c b/marchenko3D/writeDataIter.c
new file mode 100644
index 0000000000000000000000000000000000000000..e705736dba0adfb711d2b542d8c6ab224618b09e
--- /dev/null
+++ b/marchenko3D/writeDataIter.c
@@ -0,0 +1,65 @@
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include "segy.h"
+#include "par.h"
+
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+
+/**
+* writes an 2D array to a SU file
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+void name_ext(char *filename, char *extension);
+int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
+
+int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int iter)
+{
+ FILE *fp_iter;
+ size_t nwrite;
+ int i, l, j, ret, tracf, size, ix;
+ char number[16], filename[1024];
+ float *trace;
+
+ trace = (float *)malloc(n1*sizeof(float));
+ strcpy(filename, file_iter);
+ sprintf(number,"_%03d",(iter+1));
+ name_ext(filename, number);
+ fp_iter = fopen(filename, "w+");
+ if (fp_iter==NULL) verr("error on creating output file %s", filename);
+ tracf=1;
+ size=n1*n2;
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ ix = ixpos[i]; /* select proper position */
+ hdrs[i].fldr = l+1;
+ hdrs[i].sx = NINT(xsyn[l]*1000);
+ hdrs[i].offset = (long)NINT((f2+i*d2) - xsyn[l]);
+ hdrs[i].tracf = tracf++;
+ hdrs[i].selev = NINT(zsyn[l]*1000);
+ hdrs[i].sdepth = NINT(-zsyn[l]*1000);
+ /* rotate to get t=0 in the middle */
+ hdrs[i].f1 = -n1*0.5*hdrs[i].d1;
+ memcpy(&trace[0],&data[l*size+ix*n1],n1*sizeof(float));
+ for (j = 0; j < n1/2; j++) {
+ trace[n1/2+j] = data[l*size+ix*n1+j];
+ }
+ for (j = n1/2; j < n1; j++) {
+ trace[j-n1/2] = data[l*size+ix*n1+j];
+ }
+ nwrite = fwrite(&hdrs[i], 1, TRCBYTES, fp_iter);
+ assert(nwrite == TRCBYTES);
+ nwrite = fwrite(trace, sizeof(float), n1, fp_iter);
+ assert (nwrite == n1);
+ }
+ }
+ ret = fclose(fp_iter);
+ if (ret < 0 ) verr("error on writing output file.");
+ free(trace);
+
+ return 0;
+}