adding invisible demo

marchenko/demo/invisible/README

+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
+

marchenko/demo/invisible/clean

+#!/bin/bash
+
+rm *.su *.bin *.eps nep line* *.asci
+

marchenko/demo/invisible/marchenko.scr

+#!/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

marchenko/demo/invisible/model.scr

+#!/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
+
+

marchenko/demo/invisible/p4all.scr

+#!/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
+

marchenko/demo/invisible/primaries.scr

+#!/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 
+
+

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