From 631c4181be0a11d85b7f531d18ef40438c58caf8 Mon Sep 17 00:00:00 2001
From: Jan Thorbecke <janth@xs4all.nl>
Date: Thu, 3 Jan 2019 10:31:04 +0100
Subject: [PATCH] adding script to reproduce results in papers
---
.../ScatterModel/ClassicBacpProp.scr | 92 ++++++++
.../ScatterModel/ClassicTimeReverse.scr | 223 ++++++++++++++++++
.../ScatterModel/addNoiseShots.scr | 36 +++
scripts/SolidEarth2018/ScatterModel/corr.scr | 69 ++++++
scripts/SolidEarth2018/ScatterModel/model.scr | 101 ++++++++
.../ScatterModel/shotsNoise.scr | 69 ++++++
.../noFault/ClassicTimeReverse.scr | 175 ++++++++++++++
.../SolidEarth2018/noFault/addNoiseShots.scr | 32 +++
.../SolidEarth2018/noFault/backpropf1plus.scr | 111 +++++++++
scripts/SolidEarth2018/noFault/backpropf2.scr | 132 +++++++++++
scripts/SolidEarth2018/noFault/corr.scr | 82 +++++++
scripts/SolidEarth2018/noFault/marchenko.scr | 63 +++++
.../SolidEarth2018/noFault/model_nofault.scr | 29 +++
.../SolidEarth2018/noFault/model_nofault2.scr | 29 +++
scripts/SolidEarth2018/noFault/shotsQueue.scr | 73 ++++++
15 files changed, 1316 insertions(+)
create mode 100755 scripts/SolidEarth2018/ScatterModel/ClassicBacpProp.scr
create mode 100755 scripts/SolidEarth2018/ScatterModel/ClassicTimeReverse.scr
create mode 100755 scripts/SolidEarth2018/ScatterModel/addNoiseShots.scr
create mode 100755 scripts/SolidEarth2018/ScatterModel/corr.scr
create mode 100755 scripts/SolidEarth2018/ScatterModel/model.scr
create mode 100755 scripts/SolidEarth2018/ScatterModel/shotsNoise.scr
create mode 100755 scripts/SolidEarth2018/noFault/ClassicTimeReverse.scr
create mode 100755 scripts/SolidEarth2018/noFault/addNoiseShots.scr
create mode 100755 scripts/SolidEarth2018/noFault/backpropf1plus.scr
create mode 100755 scripts/SolidEarth2018/noFault/backpropf2.scr
create mode 100755 scripts/SolidEarth2018/noFault/corr.scr
create mode 100755 scripts/SolidEarth2018/noFault/marchenko.scr
create mode 100755 scripts/SolidEarth2018/noFault/model_nofault.scr
create mode 100755 scripts/SolidEarth2018/noFault/model_nofault2.scr
create mode 100755 scripts/SolidEarth2018/noFault/shotsQueue.scr
diff --git a/scripts/SolidEarth2018/ScatterModel/ClassicBacpProp.scr b/scripts/SolidEarth2018/ScatterModel/ClassicBacpProp.scr
new file mode 100755
index 0000000..7272fd1
--- /dev/null
+++ b/scripts/SolidEarth2018/ScatterModel/ClassicBacpProp.scr
@@ -0,0 +1,92 @@
+#!/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
diff --git a/scripts/SolidEarth2018/ScatterModel/ClassicTimeReverse.scr b/scripts/SolidEarth2018/ScatterModel/ClassicTimeReverse.scr
new file mode 100755
index 0000000..5bd28b1
--- /dev/null
+++ b/scripts/SolidEarth2018/ScatterModel/ClassicTimeReverse.scr
@@ -0,0 +1,223 @@
+#!/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
diff --git a/scripts/SolidEarth2018/ScatterModel/addNoiseShots.scr b/scripts/SolidEarth2018/ScatterModel/addNoiseShots.scr
new file mode 100755
index 0000000..d88b243
--- /dev/null
+++ b/scripts/SolidEarth2018/ScatterModel/addNoiseShots.scr
@@ -0,0 +1,36 @@
+#!/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
+
diff --git a/scripts/SolidEarth2018/ScatterModel/corr.scr b/scripts/SolidEarth2018/ScatterModel/corr.scr
new file mode 100755
index 0000000..dbb30f0
--- /dev/null
+++ b/scripts/SolidEarth2018/ScatterModel/corr.scr
@@ -0,0 +1,69 @@
+#!/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
+
+
diff --git a/scripts/SolidEarth2018/ScatterModel/model.scr b/scripts/SolidEarth2018/ScatterModel/model.scr
new file mode 100755
index 0000000..be8f960
--- /dev/null
+++ b/scripts/SolidEarth2018/ScatterModel/model.scr
@@ -0,0 +1,101 @@
+#!/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
+
diff --git a/scripts/SolidEarth2018/ScatterModel/shotsNoise.scr b/scripts/SolidEarth2018/ScatterModel/shotsNoise.scr
new file mode 100755
index 0000000..580562c
--- /dev/null
+++ b/scripts/SolidEarth2018/ScatterModel/shotsNoise.scr
@@ -0,0 +1,69 @@
+#!/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
+
diff --git a/scripts/SolidEarth2018/noFault/ClassicTimeReverse.scr b/scripts/SolidEarth2018/noFault/ClassicTimeReverse.scr
new file mode 100755
index 0000000..d91f180
--- /dev/null
+++ b/scripts/SolidEarth2018/noFault/ClassicTimeReverse.scr
@@ -0,0 +1,175 @@
+#!/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
+
diff --git a/scripts/SolidEarth2018/noFault/addNoiseShots.scr b/scripts/SolidEarth2018/noFault/addNoiseShots.scr
new file mode 100755
index 0000000..91c8975
--- /dev/null
+++ b/scripts/SolidEarth2018/noFault/addNoiseShots.scr
@@ -0,0 +1,32 @@
+#!/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
+
diff --git a/scripts/SolidEarth2018/noFault/backpropf1plus.scr b/scripts/SolidEarth2018/noFault/backpropf1plus.scr
new file mode 100755
index 0000000..a59593c
--- /dev/null
+++ b/scripts/SolidEarth2018/noFault/backpropf1plus.scr
@@ -0,0 +1,111 @@
+#!/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
+
diff --git a/scripts/SolidEarth2018/noFault/backpropf2.scr b/scripts/SolidEarth2018/noFault/backpropf2.scr
new file mode 100755
index 0000000..d5dc987
--- /dev/null
+++ b/scripts/SolidEarth2018/noFault/backpropf2.scr
@@ -0,0 +1,132 @@
+#!/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
+
diff --git a/scripts/SolidEarth2018/noFault/corr.scr b/scripts/SolidEarth2018/noFault/corr.scr
new file mode 100755
index 0000000..3b34477
--- /dev/null
+++ b/scripts/SolidEarth2018/noFault/corr.scr
@@ -0,0 +1,82 @@
+#!/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
+
diff --git a/scripts/SolidEarth2018/noFault/marchenko.scr b/scripts/SolidEarth2018/noFault/marchenko.scr
new file mode 100755
index 0000000..333dd10
--- /dev/null
+++ b/scripts/SolidEarth2018/noFault/marchenko.scr
@@ -0,0 +1,63 @@
+#!/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
+
+
+
diff --git a/scripts/SolidEarth2018/noFault/model_nofault.scr b/scripts/SolidEarth2018/noFault/model_nofault.scr
new file mode 100755
index 0000000..00d530d
--- /dev/null
+++ b/scripts/SolidEarth2018/noFault/model_nofault.scr
@@ -0,0 +1,29 @@
+#!/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
+
diff --git a/scripts/SolidEarth2018/noFault/model_nofault2.scr b/scripts/SolidEarth2018/noFault/model_nofault2.scr
new file mode 100755
index 0000000..00d530d
--- /dev/null
+++ b/scripts/SolidEarth2018/noFault/model_nofault2.scr
@@ -0,0 +1,29 @@
+#!/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
+
diff --git a/scripts/SolidEarth2018/noFault/shotsQueue.scr b/scripts/SolidEarth2018/noFault/shotsQueue.scr
new file mode 100755
index 0000000..38cc727
--- /dev/null
+++ b/scripts/SolidEarth2018/noFault/shotsQueue.scr
@@ -0,0 +1,73 @@
+#!/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
+
--
GitLab