Cosmetic changes and ready for marchenko paper R1

FFTlib/fft_mayer.c

@@ -64,6 +64,7 @@ REAL *fz;
  REAL *fi,*fn,*gi;
  TRIG_VARS;
 
+#pragma omp threadprivate (halsec,costab,coswrk,sintab,sinwrk)
  for (k1=1,k2=0;k1<n;k1++)
     {
      REAL a;

fdelmodc/acoustic4.c

@@ -94,7 +94,7 @@ int acoustic4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixs
 	}
 
 	/* calculate vz for all grid points except on the virtual boundary */
-#pragma omp for private (ix, iz) schedule(guided,1)
+#pragma omp for private (ix, iz) schedule(guided,1) 
 #pragma ivdep
 	for (ix=mod.ioZx; ix<mod.ieZx; ix++) {
 #pragma ivdep
@@ -120,7 +120,7 @@ int acoustic4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixs
     if (bnd.rig==2) mod.iePx -= bnd.npml;
 
 	/* calculate p/tzz for all grid points except on the virtual boundary */
-#pragma omp for private (ix, iz) schedule(guided,1)
+#pragma omp for private (ix, iz) schedule(guided,1) 
 //#pragma omp for private (ix, iz) schedule(dynamic) 
 #pragma ivdep
 	for (ix=mod.ioPx; ix<mod.iePx; ix++) {

fdelmodc/demo/fdelmodc_visco.scr

@@ -50,9 +50,8 @@ export fileqs=relax_cs.su
 	zrcv1=300 zrcv2=300 \
 	dtrcv=0.004 xsrc=1000 zsrc=300 nshot=1 \
 	src_type=1 \
-	tapfact=0.9 \
-	ntaper=200 \
-	left=4 right=4 bottom=4 top=4 \
+	ntaper=100 \
+	left=2 right=2 bottom=2 top=2 \
 	tmod=1.5 dt=0.001 
 
 #    tsnap1=0 tsnap2=1.5 dtsnap=0.05 \
@@ -89,9 +88,8 @@ export fileqs=relax_cs.su
 	zrcv1=300 zrcv2=300 \
 	dtrcv=0.004 xsrc=1000 zsrc=300 nshot=1 ntaper=100 \
 	src_type=1 \
-	tapfact=0.9 \
-	ntaper=200 \
-	left=4 right=4 bottom=4 top=4 \
+	ntaper=100 \
+	left=2 right=2 bottom=2 top=2 \
 	tmod=1.5 dt=0.001 
 
 # substract mean and plot in eps

marchenko/applyMute.c

@@ -12,7 +12,7 @@
 #endif
 #define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
 
-void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs, int nt, int *xrcvsyn, int npossyn, int shift)
+void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *xrcvsyn, int npossyn, int shift)
 {
      int i, j, l, isyn;
     float *costaper, scl;
@@ -26,7 +26,7 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs
         }
     }
 
-    for (isyn = 0; isyn < Nsyn; isyn++) {
+    for (isyn = 0; isyn < Nfoc; isyn++) {
         if (above==1) {
             for (i = 0; i < npossyn; i++) {
                 imute = xrcvsyn[i];

marchenko/demo/oneD/README

@@ -125,7 +125,7 @@ compare_003.eps
 compare_004.eps
 
 
-Note that the script epsIterwithLabels.scr produces the same figures, but with axis-labels. In figures 6 and 8 only one picture has axis-labels to avoid a busy picture. From that single label the labels of the other pictures, with the same axis, can be derived. 
+Note that the script epsIterwithLabels.scr produces the same figures, but with axis-labels. 
 
 --------------------------
 * Figure 7: Comparison of Marchenko result with reference

marchenko/demo/oneD/epsBackprop.scr

@@ -27,15 +27,15 @@ sumax < ${file_snap}_sp.su mode=abs outpar=nep
 clip=`cat nep | awk '{print $1/2}'`
 
 #first snap-shot with labels
-  fldr=71
-  times=$(echo "scale=2; $dtsnap*(${fldr}-$nsnap)" | bc -l)
-  atime=`printf "%4.2f" $times`
-  suwind key=fldr min=$fldr max=$fldr < ${file_snap}_sp.su | \
-	supsimage hbox=4 wbox=6 labelsize=10 \
-   	label1="depth (m)" label2="lateral distance (m)" \
-	x1beg=0 x1end=1250.0 clip=${clip} \
-    curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black curvedash=3,3,3 \
-	n1tic=4 f2=-1000 d2=$dx x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_snap}_${atime}_labels.eps
+#  fldr=71
+#  times=$(echo "scale=2; $dtsnap*(${fldr}-$nsnap)" | bc -l)
+#  atime=`printf "%4.2f" $times`
+#  suwind key=fldr min=$fldr max=$fldr < ${file_snap}_sp.su | \
+#	supsimage hbox=4 wbox=6 labelsize=10 \
+#   	label1="depth (m)" label2="lateral distance (m)" \
+#	x1beg=0 x1end=1250.0 clip=${clip} \
+#    curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black curvedash=3,3,3 \
+#	n1tic=4 f2=-1000 d2=$dx x2beg=-1000 f2num=-1000 d2num=500 x2end=1000 > ${file_snap}_${atime}_labels.eps
 
 for fldr in 71 86 98 99 101 103 104 116 131;
 do

marchenko/demo/oneD/epsCompare.scr

@@ -24,13 +24,14 @@ suwind < pgreen512.su j=50 s=1 | \
     n1tic=2 d2=250 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
     f2=-2250 f2num=-2000 d2num=500 > green.eps
 suwind < referenceP_rp.su j=50 s=1 | \
-	supswigp n2=19 fill=0 tracecolor=gray \
+	supswigp n2=19 fill=0 tracecolor=#F \
 	hbox=4 wbox=8 labelsize=10 linewidth=2.0 \
     label1="time (s)" label2="lateral distance (m)" \
     n1tic=2 d2=250 f1=0.0 x1beg=0 x1end=2.004 d1num=0.4 \
     f2=-2250 f2num=-2000 d2num=500 > ref.eps
 
-sed -i -e "s/%%EndProlog/[ 1 1 ]   0 setdash %%EndProlog/" green.eps
+sed -i.old -e "s/%%EndProlog/[ 1 1 ]   0 setdash %%EndProlog/" green.eps
+sed -i.old -e "s/0.5 0.5 0.5 setrgbcolor/0.65 0.65 0.65 setrgbcolor /" ref.eps
 
 psmerge in=ref.eps in=green.eps > mergeGreenRef.eps
 

marchenko/demo/oneD/epsModel.scr

@@ -40,7 +40,7 @@ supsgraph < wavefw.su \
  
 sufft < wavefw.su | suamp | sugain scale=$dt | supsgraph \
     labelsize=12 style=normal \
-    label1="time (s)" label2="amplitude" \
+    label1="frequency (1/s)" label2="amplitude" \
     d1num=10 wbox=6 hbox=3 x1end=125 x2end=1.1 > wavefw_freq.eps
  
 

marchenko/demo/oneD/marchenko.scr

@@ -27,6 +27,14 @@ 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 
+#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/twoD/README

 Description of files:
-1) shots.scr create the shots
-2) model.scr computes the model
-3) direct_wave.scr crate the direct wave to be removed from the shots
-4) remove_direct.scr remove the direct wave from the shots and scale them
-5) first_arrival.scr computes the first arrival
-6) marchenko.scr perform the Marchenko scheme
-7) referenceShot.scr creates the reference Green's function
+1a) model.scr computes the model
+1b) shots.scr creates the shots
+2) direct.scr creates the direct arrival to be removed from the shots
+3) remove_direct.scr remove the direct wave from the shots 
+4) initialFocus.scr model G_d the intitial focusing function => iniFocus_z1100_x0_rp.su
+5) referenceShot.scr creates the reference Green's function at focal point => referenceP_z1100_x0_rp.su
+5) marchenko.scr perform the Marchenko scheme
 

marchenko/demo/twoD/direct.scr

@@ -25,7 +25,7 @@ fdelmodc \
 	rec_delay=0.3 \
     dtrcv=0.004 \
     verbose=2 \
-    tmod=4.394 \
+    tmod=4.392 \
     dxrcv=10.0 \
     xrcv1=-6000 xrcv2=6000 \
     zrcv1=0 zrcv2=0 \

marchenko/demo/twoD/first_arrival.scr

-#!/bin/bash
-#PBS -N fdelmod
-#PBS -q verylong
-#PBS -l nodes=1
-#PBS -k eo
-#PBS -j eo
-
-export PATH=$HOME/src/OpenSource/bin:$PATH:
-
-dx=2.5
-dt=0.0005
-
-makemod sizex=6000 sizez=2000 dx=$dx dz=$dx cp0=1900  ro0=1200 \
-        orig=-3000,0 file_base=synclDown.su verbose=2 \
-        intt=def x=-3000,500,3000 z=195,195,195 poly=1 cp=1950 ro=3700 \
-        intt=def x=-3000,3000 z=600,600 poly=0 cp=2050 ro=1750 \
-        intt=def x=-3000,3000 z=680,680 poly=0 cp=2150 ro=2220 \
-        intt=def x=-3000,3000 z=780,780 poly=0 cp=2230 ro=1700 \
-        intt=def x=-3000,-2200,-1500,0,1300,2100,3000 z=520,580,680,840,680,600,500 poly=2 cp=2400 ro=2800 \
-
-#        intt=def x=-3000,0,3000 z=1110,1110,1110 poly=0 cp=2300 ro=1950 \
-#        intt=def x=-3000,3000 z=1180,1180 poly=0 cp=2480 ro=1820 \
-#        intt=def x=-3000,0,3000 z=1290,1290,1370 poly=0 cp=2600 ro=2000 \
-#        intt=def x=-3000,3000 z=1380,1380 poly=0 cp=2720 ro=2050 \
-#        intt=def x=-3000,3000 z=1480,1480 poly=0 cp=2800 ro=1850
-
-makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1 scale=1
-
-#smooth file_in=synclDown_cp.su power=-1.0 ntsm=29 nxsm=29 niter=15 file_out=syncls_cp.su
-#smooth file_in=synclDown_ro.su power=-1.0 ntsm=29 nxsm=29 niter=15 file_out=syncls_ro.su
-
-dxshot=10
-ishot=300
-nshots=301
-
-export OMP_NUM_THREADS=1
-mkdir -p shots
-mkdir -p jobs
-
-while (( ishot < nshots ))
-do
-
-		(( xsrc = -3000 + ${ishot}*${dxshot} ))
-		echo xsrc=$xsrc
-		file_rcv=shots/shotsmonPz1100_${xsrc}.su
-
-  cat << EOF > jobs/pbs_$ishot.job 
-#!/bin/bash
-#
-#PBS -q medium
-#PBS -N mod_${xsrc}
-#PBS -j eo 
-#PBS -m n 
-#PBS -l nodes=1
-#PBS -V
-
-
-export PATH=\$HOME/src/OpenSource/bin:\$PATH:
-cd \$PBS_O_WORKDIR
-
-export OMP_NUM_THREADS=4
-	
-		fdelmodc \
-			file_cp=synclDown_cp.su ischeme=1 iorder=4 \
-   			file_den=synclDown_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=2.100 \
-   			dxrcv=10.0 \
-   			xrcv1=-3000 xrcv2=3000 \
-   			zrcv1=0 zrcv2=0 \
-   			xsrc=$xsrc zsrc=1100 \
-   			ntaper=200 \
-   			left=2 right=2 top=2 bottom=2
-EOF
-
-        qsub jobs/pbs_$ishot.job 
-
-			(( ishot = $ishot + 1))
-done
-
-
-

marchenko/demo/twoD/marchenko.scr

 #!/bin/bash -x
-#PBS -N fdelmod
-#PBS -q verylong
-#PBS -l nodes=1
-#PBS -k eo
-#PBS -j eo
 
 export PATH=$HOME/src/OpenSource/bin:$PATH:
 
 export OMP_NUM_THREADS=1
 
-cd Redatum
-
-#for dt=0.004 with modeling at 0.0005
-smooth=3
-
-fmute file_shot=shots/shotsmonPz1100_0_rp.su file_out=p0plus.su above=-1 shift=-10 verbose=1 check=1 hw=4
-
-suwind itmax=1023 < p0plus.su | \
-suwind key=gx min=-3000000 max=3000000 | \
-sushw key=fldr a=1 > p0plussx.su
-
-#
-marchenko file_shot=../shots/refl_rp.su file_tinv=p0plussx.su nshots=601 file_green=pgreen.su verbose=1 tap=0 ntap=10 niter=15 hw=8 shift=7 smooth=$smooth file_gplus=Gplus0.su file_gmin=Gmin0.su  file_f1plus=f1plus0.su file_f1min=f1min0.su file_pplus=Pplus0.su
+#mute all events below the first arrival to get the intial focusing field
+fmute file_shot=shots/iniFocus_z1100_x0_rp.su file_out=p0plus.su above=-1 shift=-10 verbose=1 check=1 hw=4
 
+#apply the Marchenko algorithm
+marchenko file_shot=../shots/refl_rp.su file_tinv=p0plus.su nshots=601 verbose=1 \
+	tap=0 niter=15 hw=8 shift=7 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
 

marchenko/demo/twoD/referenceShot.scr

@@ -20,7 +20,7 @@ fdelmodc \
     file_cp=../syncl_cp.su ischeme=1 iorder=4 \
     file_den=../syncl_ro.su \
     file_src=wave.su \
-    file_rcv=virtual_shot_fd_P_zsrc1100.su \
+    file_rcv=referenceP_z1100_x0.su \
     src_type=1 \
 	src_orient=1 \
 	src_injectionrate=1 \

marchenko/demo/twoD/remove_direct.scr

 #!/bin/bash
-#PBS -N fdelmod
-#PBS -q verylong
-#PBS -l nodes=1
-#PBS -k eo
-#PBS -j eo
 
 export PATH=$HOME/src/OpenSource/bin:$PATH:
 

marchenko/demo/twoD/shots.scr

@@ -57,7 +57,7 @@ cd \$PBS_O_WORKDIR
 		rec_delay=0.3 \
    		dtrcv=0.004 \
    		verbose=2 \
-   		tmod=4.394 \
+   		tmod=4.392 \
    		dxrcv=10.0 \
    		xrcv1=-3000,-3000,-3000 xrcv2=3000,3000,3000 \
    		zrcv1=0,1000,1600 zrcv2=0,1000,1600 \

marchenko/fmute.c

@@ -25,7 +25,7 @@ int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *
 int readData(FILE *fp, float *data, segy *hdrs, int n1);
 int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
 int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
-void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs, int nt, int *xrcvsyn, int npossyn, int shift);
+void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *xrcvsyn, int npossyn, int shift);
 double wallclock_time(void);
 
 /*********************** self documentation **********************/
@@ -43,11 +43,12 @@ char *sdoc[] = {
 " Optional parameters: ",
 " ",
 "   file_out= ................ output file",
-"   above=0 .................. mute above(1), around(0) or below(-1) the maximum times of file_mute",
+"   above=0 .................. mute after(0), before(1) or around(2) the maximum times of file_mute",
+"   .......................... options 4 is the inverse of 0 and -1 the inverse of 1",
 "   shift=0 .................. number of points above(positive) / below(negative) maximum time for mute",
 "   check=0 .................. plots muting window on top of file_mute: output file check.su",
 "   scale=0 .................. scale data by dividing through maximum",
-"   hw=15 .................... window in time samples to look for maximum in next trace of file_mute",
+"   hw=15 .................... number of time samples to look up and down in next trace for maximum",
 "   smooth=0 ................. number of points to smooth mute with cosine window",
 //"   nxmax=512 ................ maximum number of traces in input file",
 //"   ntmax=1024 ............... maximum number of samples/trace in input file",

marchenko/readShotData.c

@@ -18,7 +18,7 @@ void rc1fft(float *rdata, complex *cdata, int n, int sign);
 int compare(const void *a, const void *b) 
 { return (*(float *)b-*(float *)a); }
 
-int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int ngath, int nx, int nxm, int ntfft, int mode, float weight, float tsq, int verbose)
+int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int ngath, int nx, int nxm, int ntfft, int mode, float scale, float tsq, int verbose)
 {
     FILE *fp;
     segy hdr;
@@ -95,8 +95,8 @@ int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx
 
             rc1fft(trace,ctrace,ntfft,-1);
             for (iw=0; iw<nw; iw++) {
-                cdata[igath*nx*nw+iw*nx+itrace].r = weight*ctrace[nw_low+iw].r;
-                cdata[igath*nx*nw+iw*nx+itrace].i = weight*mode*ctrace[nw_low+iw].i;
+                cdata[igath*nx*nw+iw*nx+itrace].r = scale*ctrace[nw_low+iw].r;
+                cdata[igath*nx*nw+iw*nx+itrace].i = scale*mode*ctrace[nw_low+iw].i;
             }
             itrace++;
             xnx[igath]+=1;