adding script to reproduce results in papers

scripts/SolidEarth2018/ScatterModel/ClassicBacpProp.scr

+#!/bin/bash
+#
+#SBATCH -J ClassicTimeReverse
+#SBATCH --cpus-per-task=8
+#SBATCH --ntasks=1
+#SBATCH --time=1:00:00
+
+export PATH=:$HOME/src/OpenSource/bin:$HOME/bin:$PATH:
+which makewave
+which makemod
+which fdelmodc
+
+cd /vardim/home/thorbcke/data/Kees/Marchenko/SolidEarth/ScatterModel
+
+#makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
+#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
+depth=700
+
+ix1a=1
+ix1b=$(echo "scale=0; ${ix1a}+6000/${dx}" | bc -l)
+base=$(echo "scale=0; ${depth}/${dx}" | bc -l)
+
+makewave fp=25 dt=$dt file_out=wave.su nt=1024 t0=0.1 scale=1
+
+file_mod=scat
+
+export OMP_NUM_THREADS=8
+
+app=2250
+fileshom=shom_${app}.su
+suwind < shom_rp.su key=gx min=-${app}000 max=${app}000 > $fileshom
+
+fdelmodc \
+    file_cp=hom_cp.su ischeme=1 iorder=4 \
+    file_den=scatter_ro.su \
+    file_src=$fileshom \
+    grid_dir=1 \
+    src_type=1 \
+    src_injectionrate=1 \
+	src_orient=1 \
+	rec_delay=0.0 \
+    verbose=2 \
+    tmod=4.51 \
+	file_snap=${file_mod}_BackProp.su \
+	tsnap1=3.4950 dtsnap=0.05 tsnap2=4.5005 sna_type_vz=0 \
+	dxsnap=$dx dzsnap=$dx zsnap1=0 zsnap2=1200 xsnap1=-2250 xsnap2=2250 \
+    npml=250 \
+    left=2 right=2 top=2 bottom=2
+
+#	tsnap1=4.0905 dtsnap=0.0005 tsnap2=4.1005 sna_type_vz=0 \
+#	dxsnap=$dx dzsnap=$dx zsnap1=0 zsnap2=1000 xsnap1=-500 xsnap2=500 \
+
+
+#    curve=curve1,curve2,curve3 npair=25,25,25 curvecolor=black,black,black curvedash=3,3,3 \
+
+for file in ${file_mod}_BackProp 
+do
+	sumax < ${file}_sp.su mode=abs outpar=nep
+	clip=`cat nep | awk '{print $1/10}'`
+	echo $clip
+
+	for fldr in 10 13 16
+	do
+		times=$(echo "scale=2; 0.05*(13-${fldr})" | bc -l)
+		atime=`printf "%4.2f" $times`
+		suwind key=fldr min=$fldr max=$fldr < ${file}_sp.su | \
+    		supsimage hbox=4 wbox=4.7 labelsize=10 \
+    		x1beg=0 x1end=1200 clip=$clip \
+    		n1tic=4 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file}_$atime.eps
+	done
+	
+	suwind itmin=$base itmax=$base key=fldr min=13 max=13 < ${file}_sp.su | sumax mode=abs outpar=nep
+	scl=`cat nep | awk '{print 1.0/$1}'`
+	echo scale for trace = $scl 
+	suwind itmin=$base itmax=$base key=fldr min=13 max=13 < ${file}_sp.su | sugain scale=$scl | sustrip > trace.bin
+    suaddhead < trace.bin n1=1801 dt=$dx | supsgraph hbox=2 wbox=6 labelsize=10 \
+    f1=-2250 d1=$dx x1beg=-500 x1end=500 f1num=-500 d1num=500 style=normal > ${file}_z${depth}_t0.eps
+
+	(( imin = base - 50 ))
+	(( imax = base + 50 ))
+	echo $base $imin $imax
+	suwind key=fldr min=13 max=13 < ${file}_sp.su | \
+	suwind itmin=$imin itmax=$imax key=gx min=-125000 max=125000 | \
+	sustrip > ${file}_t0.bin
+
+	python3 readbin.py ${file}_t0.bin
+done
+
+rm nep trace.bin

scripts/SolidEarth2018/ScatterModel/ClassicTimeReverse.scr

+#!/bin/bash
+#
+#SBATCH -J ClassicTimeReverse
+#SBATCH --cpus-per-task=12
+#SBATCH --ntasks=1
+#SBATCH --time=1:00:00
+#SBATCH --hint=nomultithread
+
+export PATH=:$HOME/src/OpenSource/bin:$HOME/bin:$PATH:
+which makewave
+which makemod
+which fdelmodc
+
+cd /vardim/home/thorbcke/data/Kees/Marchenko/SolidEarth/ScatterModel
+
+#makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
+#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
+depth=850
+
+ix1a=1
+ix1b=$(echo "scale=0; ${ix1a}+6000/${dx}" | bc -l)
+base=$(echo "scale=0; ${depth}/${dx}" | bc -l)
+
+makewave fp=25 dt=$dt file_out=wave.su nt=1024 t0=0.1 scale=1
+
+file_mod=scat
+
+export OMP_NUM_THREADS=12
+
+#forward model of scattered response of source at depth 
+fdelmodc \
+    file_cp=hom_cp.su ischeme=1 iorder=4 \
+    file_den=scatter_ro.su \
+    file_src=wave.su \
+    file_rcv=ctr.su \
+    src_type=1 \
+    src_injectionrate=1 \
+	src_orient=1 \
+    rec_type_vz=0 \
+    rec_type_p=1 \
+    rec_int_vz=2 \
+    dtrcv=$dt \
+	rec_delay=0.1 \
+    verbose=2 \
+    tmod=4.1955 \
+    dxrcv=$dx \
+    xrcv1=-2250 xrcv2=2250 \
+    zrcv1=0 zrcv2=0 \
+    xsrc=0 zsrc=$depth \
+    npml=250 \
+    left=2 right=2 top=2 bottom=2
+
+#impulse response through scattered medium
+fdelmodc \
+    file_cp=hom_cp.su ischeme=1 iorder=4 \
+    file_den=scatter_ro.su \
+    file_src=wave.su \
+    file_rcv=ctr_impulse_response.su \
+    src_type=1 \
+    src_injectionrate=1 \
+	src_orient=1 \
+    rec_type_vz=0 \
+    rec_type_p=1 \
+    rec_int_vz=2 \
+	rec_delay=0.1 \
+    tmod=4.1940 \
+   	dtrcv=0.004 \
+   	verbose=2 \
+   	dxrcv=10.0 \
+   	xrcv1=-2250 xrcv2=2250 \
+   	zrcv1=$depth zrcv2=$depth \
+   	xsrc=$xsrc zsrc=0 \
+   	ntaper=250 \
+   	left=2 right=2 top=2 bottom=2
+
+for offset in -1000 -500 0 500 1000
+do
+suwind key=offset min=$offset max=$offset < ctr_impulse_response_rp.su | \
+	supswigp wbox=1 hbox=5 titlesize=-1 labelsize=-1 axescolor=white > impulse_response_off${offset}.eps
+done
+
+
+#Forward model of homogenoeus response of source at depth 
+fdelmodc \
+    file_cp=hom_cp.su ischeme=1 iorder=4 \
+    file_den=hom_ro.su \
+    file_src=wave.su \
+    file_rcv=shom.su \
+    src_type=1 \
+    src_injectionrate=1 \
+	src_orient=1 \
+    rec_type_vz=0 \
+    rec_type_p=1 \
+    rec_int_vz=2 \
+    dtrcv=$dt \
+	rec_delay=0.1 \
+    verbose=2 \
+    tmod=4.1955 \
+    dxrcv=$dx \
+    xrcv1=-2250 xrcv2=2250 \
+    zrcv1=0 zrcv2=0 \
+    xsrc=0 zsrc=$depth \
+    npml=250 \
+    left=2 right=2 top=2 bottom=2
+
+app=300
+filectr=ctr_${app}.su
+fileshom=shom_${app}.su
+suwind < ctr_rp.su key=gx min=-${app}000 max=${app}000 > $filectr
+suwind < shom_rp.su key=gx min=-${app}000 max=${app}000 > $fileshom
+
+#Time reverse of scattered field through scattered medium
+fdelmodc \
+    file_cp=hom_cp.su ischeme=1 iorder=4 \
+    file_den=scatter_ro.su \
+    file_src=$filectr \
+    grid_dir=1 \
+    src_type=1 \
+    src_injectionrate=1 \
+	src_orient=1 \
+	rec_delay=0.0 \
+    verbose=2 \
+    tmod=4.51 \
+	file_snap=${file_mod}_timerev_scat_scat${app}.su \
+	tsnap1=3.4950 dtsnap=0.05 tsnap2=4.5005 sna_type_vz=0 \
+	dxsnap=$dx dzsnap=$dx zsnap1=0 zsnap2=1200 xsnap1=-2250 xsnap2=2250 \
+    npml=250 \
+    left=2 right=2 top=2 bottom=2
+
+#Time reverse of homogenoeus field through homogenoeus medium
+fdelmodc \
+    file_cp=hom_cp.su ischeme=1 iorder=4 \
+    file_den=hom_ro.su \
+    file_src=$fileshom \
+    grid_dir=1 \
+    src_type=1 \
+    src_injectionrate=1 \
+	src_orient=1 \
+	rec_delay=0.0 \
+    verbose=2 \
+    tmod=4.51 \
+	file_snap=${file_mod}_timerev_hom_hom${app}.su \
+	tsnap1=3.4950 dtsnap=0.05 tsnap2=4.5005 sna_type_vz=0 \
+	dxsnap=$dx dzsnap=$dx zsnap1=0 zsnap2=1200 xsnap1=-2250 xsnap2=2250 \
+    npml=250 \
+    left=2 right=2 top=2 bottom=2
+
+#Time reverse of scattered field through homogenoeus medium
+fdelmodc \
+    file_cp=hom_cp.su ischeme=1 iorder=4 \
+    file_den=hom_ro.su \
+    file_src=$filectr \
+    grid_dir=1 \
+    src_type=1 \
+    src_injectionrate=1 \
+	src_orient=1 \
+	rec_delay=0.0 \
+    verbose=2 \
+    tmod=4.51 \
+	file_snap=${file_mod}_timerev_hom_scat${app}.su \
+	tsnap1=3.4950 dtsnap=0.05 tsnap2=4.5005 sna_type_vz=0 \
+	dxsnap=$dx dzsnap=$dx zsnap1=0 zsnap2=1200 xsnap1=-2250 xsnap2=2250 \
+    npml=250 \
+    left=2 right=2 top=2 bottom=2
+
+#    curve=curve1,curve2,curve3 npair=25,25,25 curvecolor=black,black,black curvedash=3,3,3 \
+
+#sumax < ${file_mod}_back_hom_ctr${app}_sp.su mode=abs outpar=nep
+#x2end=`cat nep | awk '{print $1}'`
+#echo $x2end
+
+for file in ${file_mod}_timerev_hom_scat${app} ${file_mod}_timerev_hom_hom${app} ${file_mod}_timerev_scat_scat${app}
+do
+	sumax < ${file}_sp.su mode=abs outpar=nep
+	clip=`cat nep | awk '{print $1/10}'`
+	echo $file has clip $clip
+
+	for fldr in 10 13 16
+	do
+		times=$(echo "scale=2; 0.05*(13-${fldr})" | bc -l)
+		atime=`printf "%4.2f" $times`
+		suwind key=fldr min=$fldr max=$fldr < ${file}_sp.su | \
+    		supsimage hbox=4 wbox=4.7 labelsize=10 \
+    		x1beg=0 x1end=1200 clip=$clip \
+    		n1tic=4 x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file}_$atime.eps
+	done
+	
+	suwind itmin=$base itmax=$base key=fldr min=13 max=13 < ${file}_sp.su | sumax mode=abs outpar=nep
+	scl=`cat nep | awk '{print 1.0/$1}'`
+	echo scale for trace = $scl 
+	suwind itmin=$base itmax=$base key=fldr min=13 max=13 < ${file}_sp.su | sugain scale=$scl | sustrip > trace.bin
+    suaddhead < trace.bin n1=1801 dt=$dx | supsgraph hbox=2 wbox=6 labelsize=10 \
+    f1=-2250 d1=$dx x1beg=-500 x1end=500 f1num=-500 d1num=500 style=normal > ${file}_z${depth}_t0.eps
+
+    suaddhead < trace.bin n1=1801 dt=$dx > ${file}_z${depth}_t0.su
+
+	(( imin = base - 50 ))
+	(( imax = base + 50 ))
+	echo $base $imin $imax
+	suwind key=fldr min=13 max=13 < ${file}_sp.su | \
+	suwind itmin=$imin itmax=$imax key=gx min=-125000 max=125000 | \
+	sustrip > ${file}_t0.bin
+
+	python3 readbin.py ${file}_t0.bin
+done
+
+
+(cat ${file_mod}_timerev_hom_hom${app}_z${depth}_t0.su; cat ${file_mod}_timerev_scat_scat${app}_z${depth}_t0.su;cat ${file_mod}_timerev_hom_scat${app}_z${depth}_t0.su ) | \
+	supsgraph hbox=2 wbox=6 labelsize=10 \
+    f1=-2250 d1=$dx x1beg=-500 x1end=500 f1num=-500 d1num=500 x2beg=-0.1 \
+	style=normal linecolor=red,blue,green > ${file_mod}_timerev_z${depth}_t0.eps
+
+
+rm nep trace.bin
+
+exit;
+
+xgraph < trace.bin n=451 pairs=2 d1=10 title=hom
+suwind itmin=75 itmax=75 key=fldr min=13 max=13 < snap_back_ctr_sp.su | sustrip > trace.bin
+xgraph < trace.bin n=451 pairs=2 d1=10 title=scatter

scripts/SolidEarth2018/ScatterModel/addNoiseShots.scr

+#!/bin/bash
+
+itmax=16383 
+suwind < shotsNoise/shots_0_rp.su key=tracl min=452 max=902 itmax=$itmax | \
+	suzero itmax=$itmax | \
+	suaddnoise f=0,5,80,100 amps=0,1,1,0 | \
+	sushw key=fldr a=1 c=1 j=1 > noiseTraces.su
+
+zsrc=0
+dxshot=10
+ishot=0
+nshots=451
+
+rm /tmp/all*.su
+
+while (( ishot < nshots ))
+do
+	(( xsrc = -2250 + ${ishot}*${dxshot} ))
+	(( ishot = $ishot + 1))
+
+	file_rcv=shotsNoise/shots_${xsrc}_rp.su
+	suwind < $file_rcv key=tracl min=452 max=902 > /tmp/sh.su
+	suwind < noiseTraces.su key=fldr min=$ishot max=$ishot > /tmp/tr.su
+	fconv file_in1=/tmp/sh.su file_in2=/tmp/tr.su ntfft=$itmax >> /tmp/allA.su
+
+	suwind < $file_rcv key=tracl min=903 max=1353 > /tmp/sh.su
+	suwind < noiseTraces.su key=fldr min=$ishot max=$ishot > /tmp/tr.su
+	fconv file_in1=/tmp/sh.su file_in2=/tmp/tr.su ntfft=$itmax >> /tmp/allB.su
+
+done
+
+susort gx fldr < /tmp/allA.su | sustack key=gx > NoiseSourcesA.su
+susort gx fldr < /tmp/allB.su | sustack key=gx > NoiseSourcesB.su
+
+rm /tmp/sh.su /tmp/tr.su /tmp/all*.su
+

scripts/SolidEarth2018/ScatterModel/corr.scr

+#!/bin/bash
+
+suwind key=gx min=0 max=0 < NoiseSourcesA.su > traceA.su
+
+fconv mode=cor1 file_in1=NoiseSourcesB.su file_in2=traceA.su file_out=BwithTraceA.su 
+
+#correlate with wavelet 
+makewave fp=25 dt=0.004 file_out=wavedt4.su nt=4096 t0=0.0 scale=1
+cp BwithTraceA.su nep.su
+fconv mode=cor1 file_in1=nep.su file_in2=wavedt4.su  file_out=BwithTraceA.su 
+
+sumax < BwithTraceA.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/3}'`
+
+supsimage hbox=2 wbox=4 labelsize=10 < BwithTraceA.su \
+	x1beg=0 x1end=1.0 clip=$clip f1=0 \
+	x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > SI_BwithTraceA.eps
+
+sudipfilt < BwithTraceA.su slopes=-0.0007,-0.0001,0,0.0001,0.0007 amps=0,1,1,1,0 | \
+	supsimage hbox=2 wbox=4 labelsize=10 \
+	x1beg=0 x1end=1.0 clip=$clip f1=0 \
+	x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > SI_BwithTraceA_dipfilt.eps
+
+sudipfilt < BwithTraceA.su slopes=-0.0007,-0.0001,0,0.0001,0.0007 amps=0,1,1,1,0 | \
+	suwind j=5 s=1 | \
+	supswigp hbox=2 wbox=4.0 labelsize=10 linewidth=0.0 \
+	x1beg=0 x1end=1.0 clip=$clip f1=0 n1tic=2 \
+	x2beg=-1000 f2num=-1000 d2=50 d2num=500 x2end=1000 > SI_BwithTraceA_dipfilt_wiggle.eps
+
+
+#optimal results obtained by placing recievers at level of TraceA (virtual shot position) in time-reversal computation 
+#fdelmodc \
+    file_cp=hom_cp.su ischeme=1 iorder=4 \
+    file_den=scatter_ro.su \
+    file_src=ctr_rp.su \
+	file_rcv=SI_reference_timerev.su \
+    grid_dir=1 \
+    src_type=1 \
+    src_injectionrate=1 \
+    src_orient=1 \
+	rec_type_vz=1 \
+    rec_delay=4.0955 \
+    verbose=2 \
+    tmod=5.1 \
+    xrcv1=-2250 xrcv2=2250 zrcv1=1000 zrcv2=1000 \
+    npml=250 \
+    left=2 right=2 top=2 bottom=2
+
+sumax < SI_reference_timerev_rp.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/5}'`
+
+supsimage hbox=2 wbox=4 labelsize=10 < SI_reference_timerev_rp.su \
+	x1beg=0 x1end=1.0 clip=$clip f1=0 \
+	x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > SI_reference_timerev_rp.eps
+
+sudipfilt < SI_reference_timerev_rp.su slopes=-0.0007,-0.0001,0,0.0001,0.0007 amps=0,1,1,1,0 | \
+	supsimage hbox=2 wbox=4 labelsize=10 \
+	x1beg=0 x1end=1.0 clip=$clip f1=0 \
+	x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > SI_reference_timerev_rp_dipfilt.eps
+
+#plot noisetraces
+
+for offset in -1000 -500 0 500 1000
+do
+suwind key=offset min=$offset max=$offset < noiseTraces.su | \
+    supswigp wbox=1 hbox=5 titlesize=-1 labelsize=-1 axescolor=white > noiseTraces_off${offset}.eps
+done
+
+

scripts/SolidEarth2018/ScatterModel/model.scr

+#!/bin/bash
+#
+#SBATCH -J model
+#SBATCH --cpus-per-task=8
+#SBATCH --ntasks=1
+#SBATCH --time=1:00:00
+
+export PATH=:$HOME/src/OpenSource/utils:$HOME/bin:$PATH:
+
+cd /vardim/home/thorbcke/data/Kees/Marchenko/SolidEarth/ScatterModel
+
+dx=2.5
+dt=0.0005
+
+#shots3=var=3000,5
+#shots=var=6000,5
+
+#high contrast model
+makemod sizex=5000 sizez=1200 dx=$dx dz=$dx cp0=1900 ro0=1200 \
+        orig=-2500,0 file_base=scatter.su \
+        intt=randdf x=-2500 z=200 cp=1900 ro=4700 var=4001,11 \
+        intt=def x=-2500,0,2500 z=700,700,700 cp=1900 ro=1200 \
+        verbose=4
+
+#low contrast model
+#makemod sizex=5000 sizez=1200 dx=$dx dz=$dx cp0=1900 ro0=1200 \
+#        orig=-2500,0 file_base=scatter.su \
+#        intt=randdf x=-2500 z=200 cp=1900 ro=2400 var=4001,11 \
+#        intt=def x=-2500,0,2500 z=700,700,700 cp=1900 ro=1200 \
+#        verbose=4
+
+makemod sizex=5000 sizez=1200 dx=$dx dz=$dx cp0=1900 ro0=1200 \
+        orig=-2500,0 file_base=hom verbose=2 
+
+makewave w=fw fmin=0 flef=5 frig=80 fmax=100  dt=$dt file_out=wavefw.su nt=4096 t0=0.3 scale=0
+
+supsimage hbox=4 wbox=6 labelsize=10 < scatter_ro.su\
+        x1beg=0 x1end=1000.0 legend=0 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=2400 wclip=1500 blockinterp=1 legend=1 \
+        n1tic=5 x2beg=-2500 f2num=-2000 d2num=1000 x2end=2500 > modelscatter_ro.eps
+
+exit;
+
+export OMP_NUM_THREADS=8
+
+makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
+
+~/src/OpenSource/bin/fdelmodc \
+    file_cp=scatter_cp.su ischeme=1 iorder=4 \
+    file_den=scatter_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=2.5 \
+    xrcv1=-2250 xrcv2=2250 \
+    zrcv1=0 zrcv2=0 \
+    xsrc=0 zsrc=0 \
+    ntaper=200 \
+    left=2 right=2 top=2 bottom=2 
+
+
+#    tsnap1=3.1 tsnap2=2.5 dtsnap=0.1 \
+
+makemod sizex=5000 sizez=1000 dx=$dx dz=$dx cp0=1900  ro0=1200 \
+        orig=-2500,0 file_base=hom.su 
+
+~/src/OpenSource/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=200 \
+    left=2 right=2 top=2 bottom=2 
+
+sudiff shot_fd_rp.su shot_hom_fd_rp.su > shot_rp.su
+

scripts/SolidEarth2018/ScatterModel/shotsNoise.scr

+#!/bin/bash
+#
+
+export PATH=:$HOME/src/OpenSource/bin:$HOME/bin:$PATH:
+
+cd /vardim/home/thorbcke/data/Kees/Marchenko/SolidEarth/ScatterModel
+
+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
+
+mkdir -p shotsNoise
+mkdir -p jobs
+
+zsrc=0
+dxshot=10
+ishot=0
+nshots=451
+
+while (( ishot < nshots ))
+do
+
+		(( xsrc = -2250 + ${ishot}*${dxshot} ))
+
+		echo ishot=$ishot xsrc=$xsrc zsrc=$zsrc
+
+  cat << EOF > jobs/job_$ishot.job 
+#!/bin/bash
+#
+#SBATCH -J scat_${xsrc}
+#SBATCH --cpus-per-task=8
+#SBATCH --ntasks=1
+#SBATCH --time=0:10:00
+
+cd \$SLURM_SUBMIT_DIR
+
+export PATH=:\$HOME/src/OpenSource/bin:\$HOME/bin:\$PATH:
+
+export OMP_NUM_THREADS=8
+file_rcv=shotsNoise/shots_${xsrc}.su
+
+fdelmodc \
+   		file_cp=scatter_cp.su ischeme=1 iorder=4 \
+   		file_den=scatter_ro.su \
+   		file_src=wavefw.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.3 \
+   		dtrcv=0.004 \
+   		verbose=2 \
+   		tmod=4.394 \
+   		dxrcv=10.0 \
+   		xrcv1=-2250,-2250,-2250 xrcv2=2250,2250,2250 \
+   		zrcv1=0,850,1000 zrcv2=0,850,1000 \
+   		xsrc=$xsrc zsrc=0 \
+   		ntaper=250 \
+   		left=2 right=2 top=2 bottom=2
+EOF
+
+sbatch jobs/job_$ishot.job
+
+   		(( ishot = $ishot + 1))
+
+done
+

scripts/SolidEarth2018/noFault/ClassicTimeReverse.scr

+#!/bin/bash
+#
+#SBATCH -J ClassicTimeReverse
+#SBATCH --cpus-per-task=12
+#SBATCH --ntasks=1
+#SBATCH --time=2:00:00
+
+export PATH=:$HOME/src/OpenSource/bin:$HOME/bin:$PATH:
+which makewave
+which makemod
+which fdelmodc
+
+cd /vardim/home/thorbcke/data/Kees/Marchenko/SolidEarth/noFault
+
+#makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
+#makewave w=fw fmin=0 flef=5 frig=80 fmax=100  dt=$dt file_out=wavefw.su nt=4096 t0=0.3
+
+dx=1.0
+dt=0.0002
+depth=1450
+
+base=$(echo "scale=0; ${depth}/${dx}" | bc -l)
+
+makewave fp=25 dt=$dt file_out=wave.su nt=1024 t0=0.1 scale=1
+
+file_mod=nofault
+
+export OMP_NUM_THREADS=12
+
+#fdelmodc \
+    file_cp=${file_mod}_cp.su ischeme=1 iorder=4 \
+    file_den=${file_mod}_ro.su \
+    file_src=wave.su \
+    file_rcv=ctr.su \
+    src_type=1 \
+    src_injectionrate=1 \
+	src_orient=1 \
+    rec_type_vz=1 \
+    rec_type_p=1 \
+    rec_int_vz=2 \
+    dtrcv=$dt \
+	rec_delay=0.1 \
+    verbose=2 \
+    tmod=4.1958 \
+    dxrcv=$dx \
+    xrcv1=-5000 xrcv2=5000 \
+    zrcv1=0 zrcv2=0 \
+    xsrc=0 zsrc=$depth \
+    npml=250 \
+    left=2 right=2 top=2 bottom=2
+
+
+scale=-1
+filectr=ctr_dn.su
+sugain scale=$scale < ctr_rvz.su > $filectr
+
+#fdelmodc \
+    file_cp=${file_mod}_cp.su ischeme=1 iorder=4 \
+    file_den=${file_mod}_ro.su \
+    file_src=$filectr \
+    grid_dir=1 \
+    src_type=1 \
+    src_injectionrate=1 \
+	src_orient=1 \
+	rec_delay=0.0 \
+    verbose=2 \
+    tmod=4.51 \
+	file_snap=${file_mod}_timerev.su \
+    tsnap1=3.4956 dtsnap=0.05 tsnap2=4.5005 sna_type_vz=0 \
+	dxsnap=$dx dzsnap=$dx zsnap1=0 zsnap2=2750 xsnap1=-1650 xsnap2=1650 \
+    npml=250 \
+    left=2 right=2 top=2 bottom=2
+
+#    curve=curve1,curve2,curve3 npair=25,25,25 curvecolor=black,black,black curvedash=3,3,3 \
+
+filectr=ctr_rp.su
+skip=250 
+sumax < $filectr mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/2}'`
+echo $clip
+suwind j=$skip s=1 < $filectr | \
+	supswigp hbox=4 wbox=6.6 d2=$skip labelsize=10 \
+	x1end=2.5 clip=$clip > ${file_mod}_P_wiggle.eps
+
+suwind j=$skip s=1 < $filectr tmax=2.5 |  \
+	supswigp hbox=4 wbox=6.6 d2=$skip labelsize=10 \
+	f1=2.5 x1beg=2.5 x1end=0 d1=-0.0002 clip=$clip > ${file_mod}_P_wiggle_flip.eps
+
+filectr=ctr_dn.su
+sumax < $filectr mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/2}'`
+echo $clip
+suwind j=$skip s=1 < $filectr | \
+	supswigp hbox=4 wbox=6.6 d2=$skip labelsize=10 \
+	x1end=2.5 clip=$clip > ${file_mod}_Vn_wiggle.eps
+
+suwind j=$skip s=1 < $filectr tmax=2.5 |  \
+	supswigp hbox=4 wbox=6.6 d2=$skip labelsize=10 \
+	f1=2.5 x1beg=2.5 x1end=0 d1=-0.0002 clip=$clip > ${file_mod}_Vn_wiggle_flip.eps
+
+
+for file in ${file_mod}_timerev
+do
+	sumax < ${file}_sp.su mode=abs outpar=nep
+	clip=`cat nep | awk '{print $1/10}'`
+	echo $clip
+
+	for fldr in 10 13 16
+	do
+		times=$(echo "scale=2; -0.05*(13-${fldr})" | bc -l)
+		atime=`printf "%4.2f" $times`
+		suwind key=fldr min=$fldr max=$fldr < ${file}_sp.su | \
+    		supsimage hbox=4.9 wbox=6.6 labelsize=10 \
+    		x1beg=300 x1end=2750 clip=$clip \
+    		x2beg=-1650 f2num=-1500 d2num=500 x2end=1650 > ${file}_$atime.eps
+	done
+	
+    suwind itmin=$base itmax=$base key=fldr min=13 max=13 < ${file}_sp.su | sumax mode=abs outpar=nep
+    scl=`cat nep | awk '{print 1.0/$1}'`
+    echo scale for trace = $scl 
+
+    suwind itmin=$base itmax=$base key=fldr min=13 max=13 < ${file}_sp.su | sugain scale=$scl | sustrip > trace.bin
+    suaddhead < trace.bin n1=3001 dt=$dx | supsgraph hbox=2 wbox=6 labelsize=10 \
+    f1=-1650 d1=$dx x1beg=-500 x1end=500 f1num=-500 d1num=500 style=normal > ${file}_z${depth}_t0.eps
+
+    (( imin = base - 125 ))
+    (( imax = base + 125 ))
+    echo $base $imin $imax
+    suwind key=fldr min=13 max=13 < ${file}_sp.su | \
+    suwind itmin=$imin itmax=$imax key=gx min=-125000 max=125000 | \
+    sustrip > ${file}_t0.bin
+
+    python3 readbin.py ${file}_t0.bin
+
+done
+
+# model slightly above depth to get Gd(x,xB,t)
+
+(( depth = depth - 400 ))
+echo $depth
+#fdelmodc \
+    file_cp=${file_mod}_cp.su ischeme=1 iorder=4 \
+    file_den=${file_mod}_ro.su \
+    file_src=wave.su \
+    file_rcv=T${depth}.su \
+    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.1 \
+    verbose=2 \
+    tmod=4.1920 \
+    dxrcv=10 \
+    xrcv1=-5000 xrcv2=5000 \
+    zrcv1=0 zrcv2=0 \
+    xsrc=0 zsrc=$depth \
+    npml=250 \
+    left=2 right=2 top=2 bottom=2
+
+fmute file_shot=T${depth}_rp.su file_out=Gd${depth}.su above=-1 shift=-10 verbose=1 check=1 hw=4
+
+sumax < T${depth}_rp.su  mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/2}'`
+skip=25
+suwind j=$skip s=1 < Gd${depth}.su  tmax=2.5 |  \
+	supswigp hbox=4 wbox=6.6 d2=250 labelsize=10 \
+	f1=2.5 x1beg=2.5 x1end=0 d1=-0.0002 clip=$clip > ${file_mod}_Gd${depth}_flip.eps
+
+suwind j=$skip s=1 < Gd${depth}.su  |  \
+	supswigp hbox=4 wbox=6.6 d2=250 labelsize=10 \
+	x1beg=0 x1end=2.5 d1=0.0002 clip=$clip > ${file_mod}_Gd${depth}.eps
+

scripts/SolidEarth2018/noFault/addNoiseShots.scr

+#!/bin/bash
+
+itmax=64530 
+suwind < shotsNoise/shots_0_rp.su key=tracl min=1 max=901 itmax=$itmax | \
+	suzero itmax=$itmax | \
+	suaddnoise f=0,5,80,100 amps=0,1,1,0 | \
+	sushw key=fldr a=1 c=1 j=1 > noiseTraces.su
+
+zsrc=0
+dxshot=10
+ishot=0
+nshots=901
+
+rm /tmp/all*.su
+
+while (( ishot < nshots ))
+do
+	(( xsrc = -4500 + ${ishot}*${dxshot} ))
+	(( ishot = $ishot + 1))
+
+	file_rcv=shotsNoise/shots_${xsrc}_rp.su
+    echo $file_rcv
+	suwind < $file_rcv key=tracl min=1 max=901 > /tmp/sh.su
+	suwind < noiseTraces.su key=fldr min=$ishot max=$ishot > /tmp/tr.su
+	fconv file_in1=/tmp/sh.su file_in2=/tmp/tr.su ntfft=$itmax >> /tmp/allA.su
+
+done
+
+susort gx fldr < /tmp/allA.su | sustack key=gx > NoiseSources.su
+
+rm /tmp/sh.su /tmp/tr.su /tmp/all*.su
+

scripts/SolidEarth2018/noFault/backpropf1plus.scr

+#!/bin/bash
+#SBATCH -J nofault_${xsrc}
+#SBATCH --cpus-per-task=12
+#SBATCH --ntasks=1
+#SBATCH --time=0:55:00
+#SBATCH --hint=nomultithread
+
+export PATH=:$HOME/src/OpenSource/bin:$PATH:
+which makewave
+which makemod
+
+cd /vardim/home/thorbcke/data/Kees/Marchenko/SolidEarth/noFault
+
+#makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
+#makewave w=fw fmin=0 flef=5 frig=80 fmax=100  dt=$dt file_out=wavefw.su nt=4096 t0=0.3
+
+dx=1.0
+dt=0.0002
+#dt=0.000125
+depth=1450
+
+#model in truncated medium: homogeneous below focal point 
+makemod sizex=10000 sizez=2000 dx=$dx dz=$dx cp0=1500 ro0=1000 \
+        orig=-5000,0 file_base=nofaultTL.su verbose=2 \
+        intt=def x=-5000,-2000,-1000,0,1000,2000,5000 z=220,170,190,210,200,240,220 poly=2 cp=1950 ro=1800 grad=0 \
+        intt=def x=-5000,-2200,-1500,0,1300,2100,5000 z=520,540,590,600,540,600,630 poly=2 cp=2000 ro=1000 grad=0 \
+        intt=def x=-5000,-1800,0,2200,5000 z=1020,1000,900,1000,1010 poly=2 cp=2300 ro=1600 grad=0 \
+		intt=def x=-5000,-1900,-1300,-800,-200,300,500,1000,1800,5000 z=1500,1450,1400,1350,1340,1350,1370,1390,1400,1450 poly=2 cp=2400 ro=2700 grad=0 
+
+#make model a bit smaller
+suwind key=gx min=-2500000 max=2500000 < nofaultTL_cp.su | sushw key=f2 a=-2500 > nofaultT_cp.su
+suwind key=gx min=-2500000 max=2500000 < nofaultTL_ro.su | sushw key=f2 a=-2500 > nofaultT_ro.su
+suwind key=gx min=-2500000 max=2500000 < f1plus.su | sushw key=f2 a=-2500 > f1plusS.su
+
+file_mod=nofaultT
+
+export OMP_NUM_THREADS=12
+
+#fdelmodc \
+    file_cp=${file_mod}_cp.su ischeme=1 iorder=4 \
+    file_den=${file_mod}_ro.su \
+    file_src=f1plusS.su \
+	dt=$dt \
+    file_rcv=backprop_f1plusz${depth}.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.002 \
+	rec_delay=0.0 \
+    verbose=2 \
+    tmod=8.004 \
+    dxrcv=10.0 \
+	tsnap1=3.0 tsnap2=4.0 dtsnap=0.05 sna_type_vz=0 \
+    xrcv1=-2250 xrcv2=2250 \
+    zrcv1=$depth zrcv2=$depth \
+    xsrc=0 zsrc=0 \
+    ntaper=250 \
+    left=2 right=2 top=2 bottom=2
+
+
+file_mod=nofault
+
+suwind < backprop_f1plusz${depth}_rp.su itmax=4095 > nep.su 
+basop choice=8 file_in=nep.su file_out=backprop_f1plusz${depth}_kz.su c=2400 verbose=1 fmax=120 
+
+sumax < backprop_f1plusz${depth}_kz.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/5}'`
+echo $clip
+skip=12 
+(( d2 = skip*10 ))
+
+suwind s=1 j=$skip < backprop_f1plusz${depth}_kz.su | \
+        supswigp hbox=4 wbox=6.6 labelsize=10 linewidth=0.0 \
+        n1tic=2 d2=$d2 f1=-4.004 x1beg=-1.5 x1end=1.5 \
+        f2=-2250 f2num=-2000 d2num=1000 clip=$clip > ${file_mod}_f1plusz${depth}_kz.eps
+
+sumax < backprop_f1plusz${depth}_rp.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/5}'`
+echo $clip
+suwind s=1 j=$skip < backprop_f1plusz${depth}_rp.su | \
+        supswigp hbox=4 wbox=6.6 labelsize=10 linewidth=0.0 \
+        n1tic=2 d2=$d2 f1=-4.004 x1beg=-1.5 x1end=1.5 \
+        f2=-2250 f2num=-2000 d2num=1000 clip=$clip > ${file_mod}_f1plusz${depth}_rp.eps
+
+#use same clip as in f2
+sumax < f2.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/5}'`
+#sumax < f1plus.su mode=abs outpar=nep
+#clip=`cat nep | awk '{print $1/10}'`
+echo $clip
+skip=12 
+(( d2 = skip*10 ))
+suwind s=1 j=$skip  < f1plus.su | \
+	suwind key=gx min=-4500000 max=4500000 | \
+    supswigp hbox=4 wbox=6.6 labelsize=10 linewidth=0.0 \
+    n1tic=2 d2=$d2 f1=-4.004 x1beg=-2.2 x1end=2.2 \
+    f2=-2500 f2num=-2000 d2num=1000 clip=$clip > ${file_mod}_f1plus.eps
+
+#sumax < f1min.su mode=abs outpar=nep
+#clip=`cat nep | awk '{print $1/6}'`
+echo $clip
+suwind s=1 j=$skip  < f1min.su | \
+	suwind key=gx min=-4500000 max=4500000 | \
+    supswigp hbox=4 wbox=6.6 labelsize=10 linewidth=0.0 \
+    n1tic=2 d2=$d2 f1=-4.004 x1beg=-2.2 x1end=2.2 \
+    f2=-2500 f2num=-2000 d2num=1000 clip=$clip > ${file_mod}_f1min.eps
+

scripts/SolidEarth2018/noFault/backpropf2.scr

+#!/bin/bash
+#SBATCH -J nofault_${xsrc}
+#SBATCH --cpus-per-task=12
+#SBATCH --ntasks=1
+#SBATCH --time=0:55:00
+#SBATCH --hint=nomultithread
+
+export PATH=$HOME/OpenSource/bin/:$PATH:
+
+dx=1.0
+dt=0.0002
+#dt=0.000125
+
+file_cp=nofaultS_cp.su
+file_ro=nofaultS_ro.su
+depth=1450
+
+export OMP_NUM_THREADS=12
+
+# 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}'`
+shift=$(echo "scale=6; ($dtrcv*(($ns)/2.0-1)+$dt)" | bc -l)
+echo $shift
+suwind key=gx min=-2250000 max=2250000 < f2.su > nep.su
+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=4.004
+
+#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 \
+    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=4.10 \
+    dxrcv=10.0 \
+    npml=250 \
+    sna_type_vz=0 \
+	dxsnap=$dx dzsnap=$dx zsnap1=0 zsnap2=2750 xsnap1=-1650 xsnap2=1650 \
+	file_snap=backpropf2.su tsnap1=3.992 dtsnap=0.0002 tsnap2=4.040 \
+    sna_type_vz=0 \
+    sna_type_p=1 \
+    left=2 right=2 top=2 bottom=2
+
+file_mod=nofault
+
+sumax < f2.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/5}'`
+echo $clip
+skip=12 
+(( d2 = skip*10 ))
+
+suwind s=1 j=$skip < f2.su > nep.su
+ns=`surange <f2.su | grep ns | awk '{print $2}'` 
+dtrcv=`surange < f2.su | grep dt | awk '{print $2/1000000.0}'`
+shift=$(echo "scale=6; ($dtrcv*(($ns)/2.0-1))" | bc -l)
+echo $shift
+basop choice=shift shift=$shift file_in=nep.su verbose=1 > nep2.su
+
+supswigp hbox=4 wbox=6.6 labelsize=10 linewidth=0.0 < nep2.su \
+        n1tic=2 d2=$d2 f1=-4.004 x1beg=-1.5 x1end=1.5 \
+        f2=-2500 f2num=-2000 d2num=1000 clip=$clip > ${file_mod}_f2_d${depth}.eps
+
+sumax < Td_rp.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/3}'`
+echo $clip
+skip=12
+(( d2 = skip*10 ))
+ns=1024
+dtrcv=0.004
+suwind j=$skip s=1 < Td_rp.su | \
+suwind itmax=2048 | suwind key=gx min=-2500000 max=2500000 > nep.su 
+shift=$(echo "scale=6; ($dtrcv*(($ns)/2.0-1))" | bc -l)
+echo $shift
+basop choice=shift shift=$shift file_in=nep.su verbose=1 > nep2.su
+
+supswigp hbox=4 wbox=6.6 labelsize=10 linewidth=0.0 < nep2.su\
+        n1tic=2 d2=$d2 f1=-$shift x1beg=-1.5 x1end=1.5 \
+        f2=-2500 f2num=-2000 d2num=1000 clip=$clip > ${file_mod}_P_cent0.eps
+
+#select snapshots at t=0
+suwind key=fldr min=30 max=30 < backpropf2_sp.su > snapt0.su
+
+sumax < snapt0.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/4}'`
+echo $clip
+
+supsimage hbox=4.9 wbox=6.6 labelsize=10 < snapt0.su \
+	x1beg=300 x1end=2750 clip=$clip \
+	x2beg=-1650 f2num=-1500 d2num=500 x2end=1650 > ${file_mod}_f2snapt0.eps
+
+
+
+#f2 at 1050
+depth=1050
+sumax < f2$depth.su mode=abs outpar=nep
+
+clip=`cat nep | awk '{print $1/5}'`
+echo $clip
+skip=12 
+(( d2 = skip*10 ))
+
+suwind s=1 j=$skip < f2$depth.su > nep.su
+ns=`surange <f2.su | grep ns | awk '{print $2}'` 
+dtrcv=`surange < f2.su | grep dt | awk '{print $2/1000000.0}'`
+shift=$(echo "scale=6; ($dtrcv*(($ns)/2.0-1))" | bc -l)
+echo $shift
+basop choice=shift shift=$shift file_in=nep.su verbose=1 > nep2.su
+
+supswigp hbox=4 wbox=6.6 labelsize=10 linewidth=0.0 < nep2.su \
+        n1tic=2 d2=$d2 f1=-4.004 x1beg=-1.5 x1end=1.5 \
+        f2=-2500 f2num=-2000 d2num=1000 clip=$clip > ${file_mod}_f2_d${depth}.eps
+

scripts/SolidEarth2018/noFault/corr.scr

+#!/bin/bash
+
+suwind key=gx min=0 max=0 < NoiseSources.su > traceA.su
+
+fconv mode=cor1 file_in1=NoiseSources.su file_in2=traceA.su file_out=BwithTraceA.su 
+
+sumax < NoiseSources.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/1}'`
+
+suwind j=5 s=1 < NoiseSources.su | \
+	supswigp hbox=2 wbox=4.0 labelsize=10 linewidth=0.0 \
+	x1beg=0 x1end=5.0 clip=$clip f1=0 n1tic=2 \
+	x2beg=-1000 f2num=-1000 d2=50 d2num=500 x2end=1000 > nofault_SI_noiseSources.eps
+
+#convolve with wavelet
+dt=0.004
+makewave fp=25 dt=$dt file_out=wavedt.su nt=1024 t0=0.0 scale=1
+suwind < BwithTraceA.su itmax=1023 > nep1.su 
+fconv mode=cor1 file_in1=nep1.su file_in2=wavedt.su file_out=nep.su 
+mv nep.su BwithTraceA.su
+rm nep1.su 
+
+sumax < BwithTraceA.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/5}'`
+
+supsimage hbox=2 wbox=4 labelsize=10 < BwithTraceA.su \
+	x1beg=0 x1end=1.0 clip=$clip f1=0 \
+	x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > nofault_SI_BwithTraceA.eps
+
+sudipfilt < BwithTraceA.su slopes=-0.0007,-0.0001,0,0.0001,0.0007 amps=0,1,1,1,0 | \
+	supsimage hbox=2 wbox=4 labelsize=10 \
+	x1beg=0 x1end=1.0 clip=$clip f1=0 \
+	x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > nofault_SI_BwithTraceA_dipfilt.eps
+
+sudipfilt < BwithTraceA.su slopes=-0.0007,-0.0001,0,0.0001,0.0007 amps=0,1,1,1,0 | \
+	suwind j=5 s=1 | \
+	supswigp hbox=2 wbox=4.0 labelsize=10 linewidth=0.0 \
+	x1beg=0 x1end=1.0 clip=$clip f1=0 n1tic=2 \
+	x2beg=-1000 f2num=-1000 d2=50 d2num=500 x2end=1000 > nofault_SI_BwithTraceA_dipfilt_wiggle.eps
+
+suwind j=5 s=1 < BwithTraceA.su | \
+	supswigp hbox=2 wbox=4.0 labelsize=10 linewidth=0.0 \
+	x1beg=0 x1end=1.0 clip=$clip f1=0 n1tic=2 \
+	x2beg=-1000 f2num=-1000 d2=50 d2num=500 x2end=1000 > nofault_SI_BwithTraceA_dipfilt_wiggle.eps
+
+export OMP_NUM_THREADS=8
+file_rcv=shots_SIref.su
+
+#fdelmodc \
+   		file_cp=nofaultS_cp.su ischeme=1 iorder=4 \
+   		file_den=nofaultS_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=1.320 \
+   		dxrcv=10.0 \
+   		xrcv1=-2500 xrcv2=2500 \
+   		zrcv1=100 zrcv2=100 \
+   		xsrc=0 zsrc=100 \
+   		ntaper=250 \
+   		left=2 right=2 top=1 bottom=2
+
+sumax < shots_SIref_rp.su mode=abs outpar=nep
+clip=`cat nep | awk '{print $1/5}'`
+
+suwind j=5 s=1 < shots_SIref_rp.su | \
+	supswigp hbox=2 wbox=4.0 labelsize=10 linewidth=0.0 \
+	x1beg=0 x1end=1.0 clip=$clip f1=0 n1tic=2 \
+	x2beg=-1000 f2num=-1000 d2=50 d2num=500 x2end=1000 > nofault_SI_reference_wiggle.eps
+
+supsimage hbox=2 wbox=4 labelsize=10 < shots_SIref_rp.su \
+	x1beg=0 x1end=1.0 clip=$clip f1=0 \
+	x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > nofault_SI_reference.eps
+

scripts/SolidEarth2018/noFault/marchenko.scr

+#!/bin/bash
+#SBATCH -J nofault_${xsrc}
+#SBATCH --cpus-per-task=12
+#SBATCH --ntasks=1
+#SBATCH --time=0:55:00
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+
+export OMP_NUM_THREADS=1
+
+cd /vardim/home/thorbcke/data/Kees/Marchenko/SolidEarth/noFault
+
+#compute Td
+dx=1.0
+dt=0.0002
+depth=1450
+
+makewave fp=25 dt=$dt file_out=wave.su nt=1024 t0=0.1 scale=1
+
+file_mod=nofault
+
+export OMP_NUM_THREADS=12
+
+#fdelmodc \
+    file_cp=${file_mod}_cp.su ischeme=1 iorder=4 \
+    file_den=${file_mod}_ro.su \
+    file_src=wave.su \
+    file_rcv=Td.su \
+    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.1 \
+    verbose=2 \
+    tmod=4.1920 \
+    dxrcv=10 \
+    xrcv1=-5000 xrcv2=5000 \
+    zrcv1=0 zrcv2=0 \
+    xsrc=0 zsrc=$depth \
+    npml=250 \
+    left=2 right=2 top=2 bottom=2
+
+fmute file_shot=Td_rp.su file_out=p0plus.su above=-1 shift=-10 verbose=1 check=1 hw=4
+
+#apply the Marchenko algorithm
+marchenko file_shot=shotsRefl/refl_rp.su file_tinv=p0plus.su nshots=901 verbose=2 \
+	tap=0 niter=15 hw=8 shift=10 smooth=3 pad=0 \
+    file_green=pgreen.su file_gplus=Gplus.su file_gmin=Gmin.su  \
+    file_f1plus=f1plus.su file_f1min=f1min.su file_f2=f2.su
+
+#apply the Marchenko algorithm for depth -= 400 => 1050
+(( depth = depth - 400 ))
+set -x 
+marchenko file_shot=shotsRefl/refl_rp.su file_tinv=Gd${depth}.su nshots=901 verbose=2 \
+	tap=0 niter=15 hw=4 shift=10 smooth=3 pad=46 \
+    file_green=pgreen${depth}.su \
+    file_f1plus=f1plus${depth}.su file_f1min=f1min${depth}.su file_f2=f2${depth}.su
+
+
+

scripts/SolidEarth2018/noFault/model_nofault.scr

+#!/bin/bash
+
+dx=1
+
+makemod sizex=10000 sizez=3000 dx=$dx dz=$dx cp0=1500 ro0=1000 supersmooth=1 \
+        orig=-5000,0 file_base=nofault.su verbose=2 \
+        intt=def x=-5000,-2000,-1000,0,1000,2000,5000 z=220,170,190,210,200,240,220 poly=2 cp=1950 ro=1800 grad=0 \
+        intt=def x=-5000,-2200,-1500,0,1300,2100,5000 z=520,540,590,600,540,600,630 poly=2 cp=2000 ro=1000 grad=0 \
+        intt=def x=-5000,-1800,0,2200,5000 z=1020,1000,900,1000,1010 poly=2 cp=2300 ro=1600 grad=0 \
+		intt=def x=-5000,-1900,-1300,-800,-200,300,500,1000,1800,5000 z=1500,1450,1400,1350,1340,1350,1370,1390,1400,1450 poly=2 cp=2400 ro=2700 grad=0 \
+		intt=def x=-5000,-2350,-1750,-1250,-650,-150,150,650,1450,5000 z=1850,1800,1750,1700,1690,1700,1720,1740,1750,1800 poly=2 cp=2600 ro=4400 grad=0 \
+		intt=def x=-5000,-2450,-1850,-1350,-750,-250,50,550,1350,5000 z=1950,1900,1850,1800,1790,1800,1820,1840,1850,1900 poly=2 cp=2700 ro=2200 grad=0 \
+		intt=def x=-5000,-2000,-500,0,500,1000,2000,5000 z=2650,2400,2520,2420,2540,2530,2720,2550 poly=2 cp=2800 ro=3000 grad=0 \
+
+
+supsimage hbox=4 wbox=6 labelsize=10 < nofault_ro.su \
+        x1beg=0 x1end=2750.0 legend=0 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 legend=1 \
+        n1tic=5 x2beg=-5000 f2num=-5000 d2num=1000 x2end=5000 > nofault_ro.eps
+
+supsimage hbox=4 wbox=6 labelsize=10 < nofault_cp.su \
+        x1beg=0 x1end=2750.0 legend=0 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 legend=1 \
+        n1tic=5 x2beg=-5000 f2num=-5000 d2num=1000 x2end=5000 > nofault_cp.eps
+

scripts/SolidEarth2018/noFault/model_nofault2.scr

+#!/bin/bash
+
+dx=1
+
+makemod sizex=10000 sizez=3000 dx=$dx dz=$dx cp0=1500 ro0=1000 supersmooth=1 \
+        orig=-5000,0 file_base=nofault.su verbose=2 \
+        intt=def x=-5000,-2000,-1000,0,1000,2000,5000 z=220,170,190,210,200,240,220 poly=2 cp=1950 ro=1800 grad=0 \
+        intt=def x=-5000,-2200,-1500,0,1300,2100,5000 z=520,540,590,600,540,600,630 poly=2 cp=2000 ro=1000 grad=0 \
+        intt=def x=-5000,-1800,0,2200,5000 z=1020,1000,900,1000,1010 poly=2 cp=2300 ro=1600 grad=0 \
+		intt=def x=-5000,-1900,-1300,-800,-200,300,500,1000,1800,5000 z=1500,1450,1400,1350,1340,1350,1370,1390,1400,1450 poly=2 cp=2400 ro=2700 grad=0 \
+		intt=def x=-5000,-2350,-1750,-1250,-650,-150,150,650,1450,5000 z=1850,1800,1750,1700,1690,1700,1720,1740,1750,1800 poly=2 cp=2600 ro=4400 grad=0 \
+		intt=def x=-5000,-2450,-1850,-1350,-750,-250,50,550,1350,5000 z=1950,1900,1850,1800,1790,1800,1820,1840,1850,1900 poly=2 cp=2700 ro=2200 grad=0 \
+		intt=def x=-5000,-2000,-500,0,500,1000,2000,5000 z=2650,2400,2520,2420,2540,2530,2720,2550 poly=2 cp=2800 ro=3000 grad=0 \
+
+
+supsimage hbox=4 wbox=6 labelsize=10 < nofault_ro.su \
+        x1beg=0 x1end=2750.0 legend=0 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 legend=1 \
+        n1tic=5 x2beg=-5000 f2num=-5000 d2num=1000 x2end=5000 > nofault_ro.eps
+
+supsimage hbox=4 wbox=6 labelsize=10 < nofault_cp.su \
+        x1beg=0 x1end=2750.0 legend=0 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 legend=1 \
+        n1tic=5 x2beg=-5000 f2num=-5000 d2num=1000 x2end=5000 > nofault_cp.eps
+

scripts/SolidEarth2018/noFault/shotsQueue.scr

+#!/bin/bash
+#
+
+export PATH=:$HOME/src/OpenSource/bin:$HOME/bin:$PATH:
+
+cd /vardim/home/thorbcke/data/Kees/Marchenko/SolidEarth/noFault
+#make model a bit smaller
+suwind key=gx min=-2500000 max=2500000 < nofault_cp.su | sushw key=f2 a=-2500 > nofaultS_cp.su
+suwind key=gx min=-2500000 max=2500000 < nofault_ro.su | sushw key=f2 a=-2500 > nofaultS_ro.su
+
+dt=0.0002
+makewave w=fw fmin=0 flef=5 frig=80 fmax=100  dt=$dt file_out=wavefw.su nt=4096 t0=0.3 scale=0
+
+mkdir -p shotsNoise
+mkdir -p jobs
+
+zsrc=2100
+dxshot=10
+ishot=0
+nshots=901
+
+while (( ishot < nshots ))
+do
+
+		(( xsrc = -4500 + ${ishot}*${dxshot} ))
+
+		echo ishot=$ishot xsrc=$xsrc zsrc=$zsrc
+
+  cat << EOF > jobs/job_$ishot.job 
+#!/bin/bash
+#
+#SBATCH -J nofault_${xsrc}
+#SBATCH --cpus-per-task=8
+#SBATCH --ntasks=1
+#SBATCH --time=1:50:00
+#SBATCH --hint=nomultithread
+
+cd \$SLURM_SUBMIT_DIR
+
+export PATH=:\$HOME/src/OpenSource/bin:\$HOME/bin:\$PATH:
+
+export OMP_NUM_THREADS=8
+file_rcv=shotsNoise/shots_${xsrc}.su
+
+fdelmodc \
+   		file_cp=nofault_cp.su ischeme=1 iorder=4 \
+   		file_den=nofault_ro.su \
+   		file_src=wavefw.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.3 \
+   		dtrcv=0.004 \
+   		verbose=2 \
+   		tmod=4.392 \
+   		dxrcv=10.0 \
+   		xrcv1=-4500 xrcv2=4500 \
+   		zrcv1=100 zrcv2=100 \
+   		xsrc=$xsrc zsrc=$zsrc \
+   		ntaper=250 \
+   		left=2 right=2 top=1 bottom=2
+EOF
+
+sbatch jobs/job_$ishot.job
+
+   		(( ishot = $ishot + 1))
+
+done
+