diff --git a/marchenko/demo/invisible/README b/marchenko/demo/invisible/README
new file mode 100644
index 0000000000000000000000000000000000000000..c07026ae3ed3a12c3b0ffd4c3757929c75076c8d
--- /dev/null
+++ b/marchenko/demo/invisible/README
@@ -0,0 +1,22 @@
+This is a special model with only (carefully selected) density contrasts.
+
+Description of files:
+1) model.scr computes the model and the 'basis' shot of R => shot_rp.su
+2) p4all.scr create from basis shot full Reflection response matrix => shotsdx4_rp.su (2.3 GB)
+3) primaries.scr perform the Marchenko scheme to remove multiples from a shot record => pred_inv.su
+
+==> 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
+ - invisible_cp.su invisible_ro.su
+ - shot_fd_rp.su
+ - shot_hom_fd_rp.su
+ - shot_rp.su
+ - wavefw.su
+
+- Compare the invisible_ro.su with the number of refectors that are visible in shot_fd_rp.su
+The multiples overlap with the primary reflections is such a way that most primaries dissapear.
+
+The primaries.scr computes the primaries only from a selected shot record. => pred_inv.su
+
+suximage < pred_inv.su clip=1
+
diff --git a/marchenko/demo/invisible/clean b/marchenko/demo/invisible/clean
new file mode 100755
index 0000000000000000000000000000000000000000..3890128152ba3f4b11471dfdb5ddd1399840bc08
--- /dev/null
+++ b/marchenko/demo/invisible/clean
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+rm *.su *.bin *.eps nep line* *.asci
+
diff --git a/marchenko/demo/invisible/marchenko.scr b/marchenko/demo/invisible/marchenko.scr
new file mode 100755
index 0000000000000000000000000000000000000000..04e495d53dcb305750132b3fa55dc4606d444b46
--- /dev/null
+++ b/marchenko/demo/invisible/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/marchenko/demo/invisible/model.scr b/marchenko/demo/invisible/model.scr
new file mode 100755
index 0000000000000000000000000000000000000000..cbd96ebdb49565f710ebc13bbfdf7ac1a0befbf2
--- /dev/null
+++ b/marchenko/demo/invisible/model.scr
@@ -0,0 +1,83 @@
+#!/bin/bash
+
+#adjust this PATH to where the code is installed
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+dx=2.0
+dt=0.0005
+
+#define gridded model for FD computations
+makemod sizex=6000 sizez=1500 dx=$dx dz=$dx cp0=1000 ro0=1000 \
+ orig=-3000,0 file_base=invisible.su verbose=2 \
+ intt=def x=-3000,3000 z=100,100 poly=0 cp=1000 ro=2000 \
+ intt=def x=-3000,3000 z=200,200 poly=0 cp=1000 ro=300 \
+ intt=def x=-3000,3000 z=300,300 poly=0 cp=1000 ro=702.34 \
+ intt=def x=-3000,3000 z=400,400 poly=0 cp=1000 ro=412.57 \
+ intt=def x=-3000,3000 z=500,500 poly=0 cp=1000 ro=594.32 \
+ intt=def x=-3000,3000 z=600,600 poly=0 cp=1000 ro=457.94 \
+ intt=def x=-3000,3000 z=700,700 poly=0 cp=1000 ro=553.64 \
+ intt=def x=-3000,3000 z=800,800 poly=0 cp=1000 ro=481.49 \
+ intt=def x=-3000,3000 z=900,900 poly=0 cp=1000 ro=533.88 \
+ intt=def x=-3000,3000 z=1000,1000 poly=0 cp=1000 ro=494.5 \
+ intt=def x=-3000,3000 z=1100,1100 poly=0 cp=1000 ro=523.51 \
+ intt=def x=-3000,3000 z=1200,1200 poly=0 cp=1000 ro=501.77 \
+
+#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
+
+export OMP_NUM_THREADS=8
+
+#Model shot record in middle of model
+fdelmodc \
+ file_cp=invisible_cp.su ischeme=1 iorder=4 \
+ file_den=invisible_ro.su \
+ file_src=wavefw.su \
+ file_rcv=shot_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=4.0 \
+ xrcv1=-3000 xrcv2=3000 \
+ 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=6000 sizez=1200 dx=$dx dz=$dx cp0=1000 ro0=1000 \
+ orig=-3000,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=shot_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=4.0 \
+ xrcv1=-3000 xrcv2=3000 \
+ 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 shot_fd_rp.su shot_hom_fd_rp.su > shot_rp.su
+
+
diff --git a/marchenko/demo/invisible/p4all.scr b/marchenko/demo/invisible/p4all.scr
new file mode 100755
index 0000000000000000000000000000000000000000..8b1ecc310548c429d399ca277aec65f6a0e1a34d
--- /dev/null
+++ b/marchenko/demo/invisible/p4all.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=4000 # with scalco factor of 1000
+ishot=0
+nshots=751
+
+echo $1
+
+rm shotsdx4_rp.su
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -1500000 + ${ishot}*${dxshot} ))
+ (( tr1 = 751 - ${ishot} ))
+ (( tr2 = ${tr1} + 750 ))
+ echo xsrc=$xsrc tr1=$tr1 tr2=$tr2
+
+ (( ishot = $ishot + 1))
+
+ suwind < shot_rp.su key=tracl min=$tr1 max=$tr2 | \
+ sushw key=sx,gx,fldr,trwf \
+ a=$xsrc,-1500000,$ishot,751 b=0,$dxshot,0,0 j=0,751,0,0 | \
+ suchw key1=offset key2=gx key3=sx c=-1 d=1000 >> shotsdx4_rp.su
+
+done
+
diff --git a/marchenko/demo/invisible/primaries.scr b/marchenko/demo/invisible/primaries.scr
new file mode 100755
index 0000000000000000000000000000000000000000..26d4483628eadaaa521a6691576c7bcdb0f29344
--- /dev/null
+++ b/marchenko/demo/invisible/primaries.scr
@@ -0,0 +1,24 @@
+#!/bin/bash -x
+#SBATCH -J marchenko_primaries
+#SBATCH --cpus-per-task=40
+#SBATCH --ntasks=1
+#SBATCH --time=0:40:00
+
+#cd $SLURM_SUBMIT_DIR
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+export OMP_NUM_THREADS=20
+
+#shot record to remove internal multiples
+select=376
+
+makewave fp=25 dt=0.004 file_out=wave.su nt=1024 t0=0.0 scale=0 scfft=1
+suwind key=fldr min=$select max=$select < shotsdx4_rp.su > shot${select}.su
+fconv file_in1=shot${select}.su file_in2=wave.su file_out=shotw.su
+
+marchenko_primaries file_shot=shotsdx4_rp.su ishot=$select \
+ file_src=wave.su file_rr=pred_inv.su \
+ verbose=2 istart=21 iend=751 fmax=90 \
+ niterskip=1024 niter=21 shift=20
+
+
diff --git a/marchenko/demo/twoD/primaries.scr b/marchenko/demo/twoD/primaries.scr
index f874ca3f566617f23927d4a8eac3c3cf5558118e..df8d79f795d9be175bbc633fa0477068958babe3 100755
--- a/marchenko/demo/twoD/primaries.scr
+++ b/marchenko/demo/twoD/primaries.scr
@@ -11,16 +11,20 @@ export PATH=$HOME/src/OpenSource/bin:$PATH:
export OMP_NUM_THREADS=40
#export OMP_NUM_THREADS=1
+#shot record to remove internal multiples
+select=301
+
makewave fp=20 dt=0.004 file_out=wave.su nt=1024 t0=0.0 scale=0 scfft=1
-suwind key=sx min=0 max=0 < shots/refl_rp.su > shotsx0.su
+suwind key=fldr min=$select max=$select < shots/refl_rp.su > shotsx0.su
fconv file_in1=shotsx0.su file_in2=wave.su file_out=shotw.su
-marchenko_primaries file_shot=shots/refl_rp.su ishot=301 file_src=wave.su \
+marchenko_primaries file_shot=shots/refl_rp.su ishot=$select file_src=wave.su \
nshots=601 verbose=2 istart=40 iend=1024 fmax=90 \
niter=22 niterskip=50 shift=20 file_rr=pred_rr.su T=0
#alternative use shotw as input, must first be multiplied by -1 (see theory)
#sugain scale=-1 < shotw.su > pplus.su
+#
#marchenko_primaries file_shot=shots/refl_rp.su file_tinv=pplus.su \
# nshots=601 verbose=2 istart=40 iend=1024 fmax=90 \
# niter=22 niterskip=50 shift=20 file_rr=pred_rr.su T=0