working on primaries demo

marchenko/demo/oneD/marchenko.scr

@@ -7,8 +7,6 @@ export MKL_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
 
-which marchenko
-ldd marchenko
 #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 \

marchenko/demo/oneD/primaries.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=20
+
+#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 ishot=$select file_src=wave.su \
+	nshots=601 verbose=1 istart=40 fmax=90 \
+	niter=15 niterskip=600 shift=20 file_rr=pred_rr.su T=0
+
+#for reference original shot record from Reflection matrix
+suwind key=fldr min=$select max=$select < shotsdx5_rp.su > 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"&
+

marchenko/marchenko_primaries.c

@@ -523,11 +523,6 @@ int main (int argc, char **argv)
                 for (i = 0; i < npos; i++) {
 					j=0;
                     ix = ixpos[i];
-                    G_d[l*nxs*nts+i*nts+j] = -f1min[l*nxs*nts+i*nts+j];
-                    for (j = 1; j < nts; j++) {
-                        G_d[l*nxs*nts+i*nts+j] = -f1min[l*nxs*nts+i*nts+nts-j];
-                    }
-/* org
                     G_d[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts] - f1min[l*nxs*nts+i*nts+j];
                     for (j = 1; j < nts; j++) {
                         G_d[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j] - f1min[l*nxs*nts+i*nts+nts-j];
@@ -535,7 +530,6 @@ int main (int argc, char **argv)
                     for (j = 0; j < nts-ii+T*shift; j++) {
                         G_d[l*nxs*nts+i*nts+j] = 0.0;
                     }
-*/
                 }
             }
         }