adding plane-wave option to marchenko_primaries; cleanig files

.gitignore

@@ -3,6 +3,7 @@
 bin/*
 *.su
 *.bin
+*.asci
 */.DS*
 **/.DS*
 .DS*
@@ -60,3 +61,6 @@ marchenko3D/fmute3D
 marchenko3D/marchenko3D
 fdelmodc3D/fdelmodc3D
 fdacrtmc/fdacrtmc
+marchenko3D/TWtransform
+utils/mutesnap
+utils/snap2shot

marchenko/demo/oneD/README_PRIMARIES

+Description of files:
+
+
+1) model.scr computes the model and the 'basis' shot of R => shot5_rp.su
+2) p5all.scr create from basis shot full Reflection response matrix => shotsdx5_rp.su (3.3 GB)
+3) primaries.scr computes the internal multiple attenuated (middle) shot record.
+4) itertions.scr computes the intermediate results of the multiple attenutation scheme and produces all output files that are used in the manuscript.
+5) epsPrimaries.scr selected output form 4) are converted to .eps pictures that are used in the Figures to explain the method.
+   reproduce the postscript files of the manuscript using SU postscript plotting programs.
+
+extra scripts
++) primariesPlane.scr: computes the internal moval scheme for plane-waves (see Meles 2020)
++) clean: remove all produced files and start with a clean directory
+
+To reproduce the Figures in the Manuscript:
+
+--------------------------
+* Figure 2: Model + Initial wavefield
+
+==> run model.scr to generate the data .su files: this will take 3-4 minutes. The files generate are:
+	- hom_cp.su, hom_ro.su
+	- model10_cp.su, model10_ro.su
+	- shot5_fd_rp.su
+	- shot5_hom_fd_rp.su
+	- shot5_rp.su
+	- wave.su
+	- wavefw.su
+
+==> run './epsPrimaries.scr Figure2' to generate the postscript files of Figure 2 
+
+model_cp_line.eps 	=> Figure 2a 
+model_ro_line.eps 	=> Figure 2b
+shotx0_rp.eps 		=> Figure 2c
+
+--------------------------
+* Figure 3: First Iteration
+
+The full R matrix is build up from the the shot record computed with model.scr
+
+==> run p5all.scr to generate the full R matrix for a fixed spread geometry. This will take less than one minute. The file generated is
+	- shotsdx5_rp.su this file has a size of 3.3 GB
+
+This R is the only input of the Marchenko Primaries algorithm.
+
+
+==> run marchenkoIter.scr to compute the first 4 iteration of the Marchenko algorithm. This will take 1-2 minutes. The generated files are:
+	- p0plus.su
+	- pgreen_001.su
+	- f1plus_001.su
+	- f1min_001.su
+	- Gplus_001.su
+	- Gmin_001.su
+	- pgreen_002.su
+	- f1plus_002.su
+	- f1min_002.su
+	- Gplus_002.su
+	- Gmin_002.su
+	- pgreen_003.su
+	- f1plus_003.su
+	- f1min_003.su
+	- Gplus_003.su
+	- Gmin_003.su
+	- pgreen_004.su
+	- f1plus_004.su
+	- f1min_004.su
+	- Gplus_004.su
+	- Gmin_004.su
+
+To generate all postscript files for Figure 4, 5 and 6
+
+==> run './epsPrimaries.scr Figure3' to generate the postscript files of Figure 3 
+
+==> run epsMarchenkoIter.scr
+
+shotx0_rp.eps 		=> Figure 4 R == Figure 3c
+p0plus.eps	 	=> Figure 4 G_d
+iter_001.eps	 	=> Figure 4 N_0
+
+shotx0_rp.eps 		=> Figure 5 R == Figure 3c
+f1min_001.eps	 	=> Figure 5 f^-_1,0
+iter_002.eps	 	=> Figure 5 -N_1
+f1plus_002.eps	 	=> Figure 5 f^+_1,0
+
+-- Figure 6 column 1
+iter_001.eps
+iter_002.eps
+iter_003.eps
+iter_004.eps
+-- Figure 6 column 2
+f1min_001.eps
+f1min_002.eps
+f1min_003.eps
+f1min_004.eps
+-- Figure 6 column 3
+p0plus_flip.eps
+f1plus_002.eps
+f1plus_003.eps
+f1plus_004.eps
+-- Figure 6 column 4
+pgreen_001.eps
+pgreen_002.eps
+pgreen_003.eps
+pgreen_004.eps
+-- Figure 6 column 5
+compare_001.eps
+compare_002.eps
+compare_003.eps
+compare_004.eps
+
+
+Note that the script epsIterwithLabels.scr produces the same figures, but with axis-labels. 
+
+--------------------------
+* Figure 7: Comparison of Marchenko result with reference
+
+To compute the marchenko results for 8 iterations.  
+
+==> run marchenko.scr This will take less than 1 minute. The generated files are:
+	- pgreen.su, pgreen512.su
+	- diffref.su
+	- Gplus0.su
+	- Gmin0.su
+	- f1plus0.su
+	- f1min0.su
+	- f2.su 
+
+
+At the end of the run the script will display in X11 a comparison of the middle trace. 
+
+To make the postscript figure 
+
+==> run epsCompare.scr
+
+mergeGreenRef.eps 	=> Figure 7
+
+--------------------------
+* Figure 8: snapshots of back propagating f2 in actual medium 
+
+To compute the snapshots 
+
+==> run backpropf2.scr This will take about 1 minute. The generated output file is
+	- backpropf2_sp.su
+
+The postscript files of Figure 8 are generated with 
+
+==> run epsBackprop.scr
+
+-- Figure 8 column 1
+backpropf2_-0.30.eps
+backpropf2_-0.15.eps
+backpropf2_-0.03.eps
+backpropf2_-0.02.eps
+backpropf2_0.00.eps
+-- Figure 8 column 2
+backpropf2_0.30.eps
+backpropf2_0.15.eps
+backpropf2_0.03.eps
+backpropf2_0.02.eps
+backpropf2_0.00.eps
+-- Figure 8 column 3
+backpropf2sum_0.30.eps
+backpropf2sum_0.15.eps
+backpropf2sum_0.03.eps
+backpropf2sum_0.02.eps
+backpropf2_0.00.eps
+
+
+The figures in the appendix, to explain the different options in the programs, are reproduced by
+
+==> run figAppendi.scr
+
+-- Figure A-1
+noise_above0.eps
+noise_above1.eps
+noise_above-1.eps
+noise_above2.eps
+noise_above4.eps
+
+-- Figure A-2
+iniFocus_shifts.eps
+

marchenko/demo/oneD/epsModel.scr

@@ -50,7 +50,7 @@ clipref=`cat nep | awk '{print $1/15}'`
 suwind key=gx min=-2250000 max=2250000 < $file | \
         supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
     	label1="time sample number" label2="lateral distance (m)" \
-        n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=500 d1=1 f1=0 d1num=100 \
+        n1tic=2 d2=5 f1=0.0 x1beg=0 x1end=400 d1=1 f1=0 d1num=100 \
         f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > shotx0_rp.eps
 
 rm nep 

marchenko/demo/oneD/epsPrimaries.scr

+#!/bin/bash
+
+export PATH=$HOME/src/OpenSource/bin/:$PATH:
+
+echo "Making eps files for $1"
+
+if [[ "$1" == "Figure2" ]];
+then
+./epsModel.scr
+exit 0
+fi
+
+istart=276
+#set same clip factor for iteration updates
+file=Mi_${istart}002.su
+sumax < $file mode=abs outpar=nep
+clipiter=`cat nep | awk '{print $1/15}'`
+ns=`surange < Mi_276014.su | grep ns | awk '{print $2}'`
+nsd=400 #number of samples to display
+(( nstart = ns - nsd ))
+
+file=f1min_${istart}002.su
+sumax < $file mode=abs outpar=nep
+clipf1=`cat nep | awk '{print $1/11}'`
+clipm1=`cat nep | awk '{print $1/22}'`
+
+#M0
+file=M0_276000.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/15}'`
+
+echo "clipiter="$clipiter "clipref="$clipref "clipf1="$clipf1
+clipiter=$clipref
+clipf1=$clipref
+
+supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+	n1tic=2 d2=5 x1beg=$nstart x1end=$ns d1=1 f1=0 f1num=700 d1num=100 \
+	label1="time sample number" label2="lateral distance" \
+    f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+suflip < $file flip=3 | \
+supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+	n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 f1=0 f1num=0 d1num=100 \
+	label1="time sample number" label2="lateral distance" \
+    f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_flip.eps
+
+# First iteration
+#convolve M0 with middle shot record of R
+select=451
+#original shot record from Reflection matrix
+suwind key=fldr min=$select max=$select < shotsdx5_rp.su > shotsx0.su
+#first iteration
+cp M0_276000.su N0.su
+#compute R*N0
+fconv file_in1=shotsx0.su file_in2=N0.su file_out=fconvN0.su verbose=1 fmax=90
+
+file=fconvN0.su
+file_base=${file%.su}
+supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+	n1tic=2 d2=5 x1beg=624 x1end=1024 d1=1 f1=0 f1num=700 d1num=100 \
+	label1="time sample number" label2="lateral distance" \
+    f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+suflip < fconvN0.su flip=3 | sugain scale=-1 > fconvN0flip.su
+file=fconvN0flip.su
+file_base=${file%.su}
+supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+	n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 f1=0 f1num=0 d1num=100 \
+	label1="time sample number" label2="lateral distance" \
+    f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+file=fconvN0.su
+file_base=${file%.su}
+supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+	n1tic=2 d2=5 x1beg=0 d1=1 f1=0 f1num=0 d1num=200 \
+	label1="time sample number" label2="lateral distance" \
+    f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}fulltime.eps
+
+# second iteration
+suwind itmax=1023 <  Mi_276001.su > N0.su
+#compute R*N0
+fconv file_in1=shotsx0.su file_in2=N0.su file_out=fconvN1.su verbose=1 fmax=90
+
+suflip < fconvN1.su flip=3 | sugain scale=-1 > fconvN1flip.su
+file=fconvN1flip.su
+file_base=${file%.su}
+supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+	n1tic=2 d2=5 x1beg=624 x1end=1024 d1=1 f1=0 f1num=600 d1num=100 \
+	label1="time sample number" label2="lateral distance" \
+    f2=-2250 f2num=-2000 d2num=1000 clip=$clipm1 > $file_base.eps
+
+file=fconvN1.su
+file_base=${file%.su}
+supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+	n1tic=2 d2=5 x1beg=0 d1=1 f1=0 f1num=0 d1num=200 \
+	label1="time sample number" label2="lateral distance" \
+    f2=-2250 f2num=-2000 d2num=1000 clip=$clipm1 > ${file_base}fulltime.eps
+
+#RNi
+file=iter_276000.su
+file_base=${file%.su}
+#clipref=$clipiter
+supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+	n1tic=2 d2=5 x1beg=0 d1=1 d1num=100 \
+	label1="time sample number" label2="lateral distance" \
+    f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
+
+#Iterations
+for (( iter=0; iter<=21; iter+=2 ))
+do
+piter=$(printf %03d $iter)
+echo $piter
+file=Mi_${istart}${piter}.su
+file_base=${file%.su}
+supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+	n1tic=2 d2=5 x1beg=$nstart x1end=$ns d1=1 d1num=100 \
+	label1="time sample number" label2="lateral distance" \
+    f2=-2250 f2num=-2000 d2num=1000 clip=$clipiter > $file_base.eps
+done
+
+for (( iter=1; iter<=21; iter+=2 ))
+do
+piter=$(printf %03d $iter)
+echo $piter
+file=Mi_${istart}${piter}.su
+file_base=${file%.su}
+supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+	n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
+	label1="time sample number" label2="lateral distance" \
+    f2=-2250 f2num=-2000 d2num=1000 clip=$clipiter > $file_base.eps
+done
+
+#iterations
+for (( iter=0; iter<=31; iter+=2 ))
+do
+piter=$(printf %03d $iter)
+echo $piter
+
+file=Mi_${istart}${piter}.su
+#ns=`surange < iter_001.su | grep ns | awk '{print $2}'`
+#dtrcv=`surange < iter_001.su | grep dt | awk '{print $2/1000000.0}'`
+#shift=$(echo "scale=4; ($dtrcv*($ns/2.0-1))" | bc -l)
+#basop choice=shift shift=$shift file_in=$file | \
+file_base=${file%.su}
+suflip flip=3 < $file | supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+	n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
+	label1="time sample number" label2="lateral distance" \
+    f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > ${file_base}flip.eps
+
+file=f1min_${istart}${piter}.su
+file_base=${file%.su}
+supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 < $file\
+	n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
+	label1="time sample number" label2="lateral distance" \
+    f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > $file_base.eps
+
+done
+
+iter=32
+piter=$(printf %03d $iter)
+echo $piter
+for (( istart=246; istart<=316; istart+=10 ))
+do
+file=f1min_${istart}${piter}.su
+file_base=${file%.su}
+supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 < $file\
+	n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
+	label1="time sample number" label2="lateral distance" \
+    f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > $file_base.eps
+done
+
+exit;
+
+#Windows for odd and even iterations
+file=WindowOdd276.su
+file_base=${file%.su}
+suwind key=tracl min=1 max=1 < $file | suwind itmin=1024 | \
+supsgraph d1=1 style=normal f1=0 \
+    labelsize=12 label1="time sample number" label2="amplitude" \
+    d1num=100 wbox=6 hbox=3 x2end=1.2 > $file_base.eps
+file=WindowEven276.su
+file_base=${file%.su}
+suwind key=tracl min=1 max=1 < $file | suwind itmax=1024 | \
+supsgraph d1=1 style=normal f1=0 \
+    labelsize=12 label1="time sample number" label2="amplitude" \
+    d1num=100 wbox=6 hbox=3 x2end=1.2 > $file_base.eps
+

marchenko/demo/oneD/iterations.scr

@@ -9,17 +9,23 @@
 
 select=451
 
-for istart in 276
-#for istart in 176 200 276 400
-#for istart in 246 256 266 276 286 296 306 316
+#for istart in 276
+for istart in 246 256 266 276 286 296 306 316
 do
 (( iend = istart + 1 ))
 
 ../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
 	nshots=901 verbose=2 istart=$istart iend=$iend fmax=90 \
-	pad=0 niter=4 smooth=10 niterskip=600 shift=20 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su 
+	pad=512 niter=33 smooth=8 niterskip=600 shift=16 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su 
 
 #suximage < f1min_${istart}043.su  clip=1 title="${istart}" &
 
 done
 
+exit;
+#T-MME
+istart=276
+(( iend = istart + 1 ))
+../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
+	nshots=901 verbose=2 istart=$istart iend=$iend fmax=90 \
+	pad=512 niter=33 smooth=8 niterskip=600 shift=16 file_rr=pred_rr_${istart}.su T=1 file_iter=iterT.su 

marchenko/demo/oneD/primaries.scr

@@ -11,13 +11,13 @@ export PATH=$HOME/src/OpenSource/bin:$PATH:
 export OMP_NUM_THREADS=40
 
 #shot record to remove internal multiples
-select=51
+select=451
 
 makewave fp=20 dt=0.004 file_out=wave.su nt=2048 t0=0.0 scale=0 scfft=1
 
 ../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
 	nshots=901 verbose=2 istart=40 iend=500 fmax=90 pad=1024 \
-	niter=20 smooth=10 niterskip=1 shift=20 file_rr=pred_rr.su T=0
+	niter=31 smooth=10 niterec=2 niterskip=1 shift=20 file_rr=pred_rr.su T=0 file_update=update.su
 
 exit;
 

marchenko/demo/oneD/primariesPlane.scr

+#!/bin/bash -x
+#SBATCH -J marchenko_primaries
+#SBATCH --cpus-per-task=40
+#SBATCH --ntasks=1
+#SBATCH --time=2:40:00
+
+
+#cd $SLURM_SUBMIT_DIR
+
+export PATH=$HOME/src/OpenSource/bin:$PATH:
+export OMP_NUM_THREADS=40
+
+#shot record to remove internal multiples
+select=451
+
+makewave fp=20 dt=0.004 file_out=wave.su nt=1024 t0=0.0 scale=0 scfft=1
+
+#../../marchenko_primaries file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=10 \
+#	nshots=901 verbose=2 istart=276 iend=277 fmax=90 pad=1024 \
+#	xorig=0 niter=20 smooth=10 niterskip=1 file_iter=iterplane.su shift=20 file_rr=plane2_10_rr.su T=0
+#exit
+../../marchenko_primaries file_shot=shotsdx5_rp.su plane_wave=1 file_src=wave.su src_angle=-15 \
+	nshots=901 verbose=2 istart=40 iend=1000 fmax=90 pad=1024 \
+	niter=20 smooth=10 niterskip=100 shift=20 file_rr=plane2_10_rr.su T=0
+
+exit
+
+#for reference original shot record from Reflection matrix
+suwind key=fldr min=$select max=$select < shotsdx5_rp.su itmax=2047 > shotsx0.su
+fconv file_in1=shotsx0.su file_in2=wave.su file_out=shotw.su
+
+# for displaying results
+
+(suwind key=offset min=0 max=0 < pred_rr.su ; suwind key=offset min=0 max=0 < shotw.su) | suxgraph &
+
+sudiff shotw.su pred_rr.su > diff.su
+suximage < shotw.su  x1end=2.5 clip=1 title="original shot"&
+suximage < pred_rr.su  x1end=2.5 clip=1 title="shot with multiples removed"&
+suximage < diff.su   x1end=2.5 clip=1 title="removed multiples"&
+