amplitude correction

marchenko/demo/twoD/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=40
+#export OMP_NUM_THREADS=1
+
+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
+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 \
+	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
+
+exit
+
+# for displaying results
+
+(suwind key=offset min=0 max=0 < pred_rr.su ; suwind key=offset min=0 max=0 < shotw.su) | suxgraph
+
+sudiff pred_rr.su shotw.su > diff.su
+suximage < shotw.su  x1end=2 clip=1 title="original shot"&
+suximage < pred_rr.su  x1end=2 clip=1 title="shot with multiples removed"&
+suximage < diff.su   x1end=2 clip=1 title="removed multiples"&
+

marchenko/marchenko_primaries.c

@@ -59,7 +59,7 @@ char *sdoc[] = {
 " Optional parameters: ",
 " ",
 " INTEGRATION ",
-"   ishot=nshots/2 ........... shot position(s) to remove internal multiples ",
+"   ishot=nshots/2 ........... shot number(s) to remove internal multiples ",
 "   file_src=spike ........... convolve ishot(s) with source wavelet",
 "   file_tinv= ............... use file_tinv to remove internal multiples",
 " COMPUTATION",
@@ -181,7 +181,8 @@ int main (int argc, char **argv)
     if(!getparint("iend", &iend)) iend=nt;
 
     if (file_tinv == NULL) {/* 'G_d' is one of the shot records */
-        if(!getparint("ishot", &ishot)) ishot=(nshots)/2;
+        if(!getparint("ishot", &ishot)) ishot=1+(nshots-1)/2;
+		ishot -= 1; /* shot numbering starts at 0 */
         Nfoc = 1;
         nxs  = nx;
         nts  = nt;
@@ -348,8 +349,7 @@ int main (int argc, char **argv)
         if (verbose) vmess("Selecting G_d from Refl of %s", file_shot);
         nts   = ntfft;
 
-        //scl   = 1.0/((float)ntfft);
-        scl   = 1.0;
+        scl   = 1.0/((float)2.0*ntfft);
         rtrace = (float *)calloc(ntfft,sizeof(float));
         ctrace = (complex *)calloc(nfreq+1,sizeof(complex));
         for (i = 0; i < xnx[ishot]; i++) {