diff --git a/marchenko3D/applyMute3D.c b/marchenko3D/applyMute3D.c
index 985b918f84121413d956a5f29b58896240b85dc8..8e51ff2851f1a771c7778f5e8d227f387f80306b 100644
--- a/marchenko3D/applyMute3D.c
+++ b/marchenko3D/applyMute3D.c
@@ -202,7 +202,7 @@ void applyMute3D_tshift( float *data, long *mute, long smooth, long above, long
for (j = MAX(0,tmute-shift),l=0; j < MAX(0,tmute-shift+smooth); j++,l++) {
Nig[j] *= costaper[l];
}
- for (j = MAX(0,tmute-shift+smooth+1); j < MIN(nt,nt+1-tmute+2*ts+shift-smooth); j++) {
+ for (j = MAX(0,tmute-shift+smooth); j < MIN(nt,nt-tmute+2*ts+shift-smooth); j++) {
Nig[j] = 0.0;
}
for (j = MIN(nt-1,nt-tmute+2*ts+shift-smooth),l=0; j < MIN(nt,nt-tmute+shift); j++,l++) {
@@ -220,6 +220,24 @@ void applyMute3D_tshift( float *data, long *mute, long smooth, long above, long
Nig[j] = 0.0;
}
}
+ else if (above==-2){ //New Theta window keeps Gd
+ imute = iypos[i]*nxs+ixpos[i];
+ tmute = mute[isyn*nxys+imute];
+ ts = tsynW[isyn*nxys+imute];
+ if (tmute >= nt/2) {
+ memset(&Nig[0],0, sizeof(float)*nt);
+ continue;
+ }
+ for (j = MAX(0,-2*ts+tmute+shift),l=0; j < MAX(0,-2*ts+tmute+shift+smooth); j++,l++) {
+ Nig[j] *= costaper[l];
+ }
+ for (j = MAX(0,-2*ts+tmute+shift+smooth); j < MIN(nt,nt-tmute-shift-smooth); j++) {
+ Nig[j] = 0.0;
+ }
+ for (j = MIN(nt,nt-tmute-shift-smooth),l=0; j < MIN(nt,nt-tmute-shift); j++,l++) {
+ Nig[j] *= costaper[smooth-l-1];
+ }
+ }
else if (above==4) { //Psi gate which is the inverse of the Theta gate (above=0)
imute = iypos[i]*nxs+ixpos[i];
tmute = mute[isyn*nxys+imute];
@@ -230,7 +248,7 @@ void applyMute3D_tshift( float *data, long *mute, long smooth, long above, long
for (j = 0; j < MAX(0,-2*ts+tmute-shift-smooth-1); j++) {
Nig[j] = 0.0;
}
- for (j = nt+1-tmute+shift+smooth; j < nt; j++) {
+ for (j = nt-tmute+shift+smooth; j < nt; j++) {
Nig[j] = 0.0;
}
for (j = nt-tmute+shift,l=0; j < nt-tmute+shift+smooth; j++,l++) {
diff --git a/marchenko3D/applyMute3D_tshift.c b/marchenko3D/applyMute3D_tshift.c
index 1111c29b164a7958b516dead018a258c6fcf4a51..e58381eae411a49846af9dc08f516792bbfe9a43 100644
--- a/marchenko3D/applyMute3D_tshift.c
+++ b/marchenko3D/applyMute3D_tshift.c
@@ -90,14 +90,32 @@ void applyMute3D_tshift( float *data, long *mute, long smooth, long above, long
Nig[j] = 0.0;
}
}
+ else if (above==-2){ //New Theta window keeps Gd
+ imute = iypos[i]*nxs+ixpos[i];
+ tmute = mute[isyn*nxys+imute];
+ ts = tsynW[isyn*nxys+imute];
+ if (tmute >= nt/2) {
+ memset(&Nig[0],0, sizeof(float)*nt);
+ continue;
+ }
+ for (j = MAX(0,-2*ts+tmute+shift),l=0; j < MAX(0,-2*ts+tmute+shift+smooth); j++,l++) {
+ Nig[isyn*nxys*nt+i*nt+j] *= costaper[l];
+ }
+ for (j = MAX(0,-2*ts+tmute+shift+smooth)+1; j < MIN(nt,nt+1-tmute-shift-smooth); j++) {
+ Nig[isyn*nxys*nt+i*nt+j] = 0.0;
+ }
+ for (j = MIN(nt,nt-tmute-shift-smooth),l=0; j < MIN(nt,nt-tmute-shift); j++,l++) {
+ Nig[isyn*nxys*nt+i*nt+j] *= costaper[smooth-l-1];
+ }
+ }
else if (above==4) { //Psi gate which is the inverse of the Theta gate (above=0)
imute = iypos[i]*nxs+ixpos[i];
tmute = mute[isyn*nxys+imute];
ts = tsynW[isyn*nxys+imute];
- for (j = -2*ts+tmute-shift-smooth,l=0; j < -2*ts+tmute-shift; j++,l++) {
+ for (j = MAX(0,-2*ts+tmute-shift-smooth),l=0; j < MAX(0,-2*ts+tmute-shift); j++,l++) {
Nig[j] *= costaper[smooth-l-1];
}
- for (j = 0; j < -2*ts+tmute-shift-smooth-1; j++) {
+ for (j = 0; j < MAX(0,-2*ts+tmute-shift-smooth-1); j++) {
Nig[j] = 0.0;
}
for (j = nt+1-tmute+shift+smooth; j < nt; j++) {
diff --git a/marchenko3D/demo/marchenko3D/oneD/README b/marchenko3D/demo/marchenko3D/oneD/README
index 52d17c72cb350b8b428f3c3944492402f14d5a17..3f79ca853bbb8d2f0577bcc0a1a83048ad53851b 100644
--- a/marchenko3D/demo/marchenko3D/oneD/README
+++ b/marchenko3D/demo/marchenko3D/oneD/README
@@ -4,3 +4,33 @@ Due to the size of the data, one requires around 256 GB of RAM. Keep this limita
To avoid the need for large storage space on your directory, we reccomend using the ZFP compressed option.
This is the standard option included in the data. You can just execute execute.sh and everything will be done for you. It will take a considerable amount of time, depending on your machine.
+
+In order to run all of the code the folders fdelmodc3D/ and utils/ in the main directory should have been compiled.
+Additionally, some steps require the installation of Seismic Unix
+
+model.sh will create a layered 2D model and convert it into a 3D model (requires Seimsic Unix and utils/)
+
+shotmodel.sh will model a shot of reflection data in the model as well as its direct arrival. It will subtract this first arrival
+from the reflection data (requires fdelmodc3D)
+
+makeR.sh will create a full 3D reflection response from the single shot by shifting traces. The file contains a warning due to the
+large size of the dataset. (requires utils/)
+
+farrmod.sh models a first arrival in a truncated version of the subsurface model. (requires fdelmodc3D/)
+
+fmute3D.sh separates the first arrival from the coda (requires utils/)
+
+zfp.sh transforms the reflection data to the frequency data and applies ZFP compression in order to shrink the datasize of the
+reflection data. It will also delete the original reflection data once the compression has been applied. (requires marchenko3D/)
+
+marchenko3d.sh will perform the Marchenko method using the data created with the previous scripts. (requires marchenko3D)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+Additionally, there are two demo's for running plane-wave Marchenko in the directories plane_horizontal/ and plane_tilted/
+In order to run these examples, the scripts model.sh shotmodel.sh makeR.sh in the main directory should all have been run
+
+plane_horizontal/ will run the 3D Marchenko method on a horizontal plane-wave
+
+plane_tilted/ will run the 3D Marchenko method on a tilted plane-wave. In order to construct a full Green's function
+the method needs to be run twice at opposing dips. For more information see
+Meles, G. A., Wapenaar, K., & Thorbecke, J. (2018). Virtual plane-wave imaging via Marchenko redatuming. Geophysical Journal International, 214(1), 508-519.
diff --git a/marchenko3D/demo/marchenko3D/oneD/execute.sh b/marchenko3D/demo/marchenko3D/oneD/execute.sh
index f183f9e9dc334595fd842b125d6a07c7f710b8df..5e9cee4dd628287ddb6aa864cdf3128c166d31c4 100755
--- a/marchenko3D/demo/marchenko3D/oneD/execute.sh
+++ b/marchenko3D/demo/marchenko3D/oneD/execute.sh
@@ -1,15 +1,15 @@
#!/bin/bash
-./model.sh
+#./model.sh
-./shotmodel.sh
+#./shotmodel.sh
-./makeR.sh
+#./makeR.sh
+
+./zfp.sh
./farrmod.sh
./fmute3D.sh
-./zfp.sh
-
./marchenko3d.sh
diff --git a/marchenko3D/demo/marchenko3D/oneD/farrmod.sh b/marchenko3D/demo/marchenko3D/oneD/farrmod.sh
index 4cfbf5211d87618c33e6987dff124e4250fc669b..b14fe894ed44cd8aef1c1b21a9cb8d69dc024985 100755
--- a/marchenko3D/demo/marchenko3D/oneD/farrmod.sh
+++ b/marchenko3D/demo/marchenko3D/oneD/farrmod.sh
@@ -7,7 +7,7 @@ makewave fp=15 fmax=22 dt=${dt} file_out=wavefpmod.su nt=8192 t0=0.2 scale=1 scf
export OMP_NUM_THREADS=8
#Model shot record in middle of model
-fdelmodc3D \
+../../../../fdelmodc3D/fdelmodc3D \
file_cp=cp3d.su ischeme=1 iorder=4 \
file_den=ro3d.su \
file_src=wavefpmod.su \
diff --git a/marchenko3D/demo/marchenko3D/oneD/fmute3D.sh b/marchenko3D/demo/marchenko3D/oneD/fmute3D.sh
index 8b89c7354eaef93cff259c4df8e0240660626654..62a62acf0cc8f0755931df80d5b27bf1d049ab41 100755
--- a/marchenko3D/demo/marchenko3D/oneD/fmute3D.sh
+++ b/marchenko3D/demo/marchenko3D/oneD/fmute3D.sh
@@ -1,4 +1,4 @@
#!/bin/bash
#Mute the coda from the modeled first arrival
-fmute3D file_mute=farrmod_rp.su file_shot=farrmod_rp.su file_out=farr.su above=2 shift=18 smooth=10 hw=5
+../../../fmute3D file_mute=farrmod_rp.su file_shot=farrmod_rp.su file_out=farr.su above=2 shift=18 smooth=10 hw=5
diff --git a/marchenko3D/demo/marchenko3D/oneD/makeR.sh b/marchenko3D/demo/marchenko3D/oneD/makeR.sh
index 87d1a122ac47f57f99eae844fac113bf49e9364e..a8a54719e5e7c6672508cc9e3794d4bc69a6d63f 100755
--- a/marchenko3D/demo/marchenko3D/oneD/makeR.sh
+++ b/marchenko3D/demo/marchenko3D/oneD/makeR.sh
@@ -4,8 +4,8 @@ echo "WARNING! The data size of the 3D reflection data is very large (>200 GB)!!
echo "When you are positive your machine has enough memory,"
echo "delete the exit; statement from makeR.sh"
-exit;
+#exit;
#create the reflection data that was modeled in 1 1D medium
-makeR1D file_in=shotx10y10.su file_out=reflx10y10.su verbose=2 nxrcv=201 nyrcv=61 nxsrc=201 nysrc=61
+../../../../utils/makeR1D file_in=shotx10y10.su file_out=reflx10y10.su verbose=2 nxrcv=201 nyrcv=61 nxsrc=201 nysrc=61
diff --git a/marchenko3D/demo/marchenko3D/oneD/marchenko3d.sh b/marchenko3D/demo/marchenko3D/oneD/marchenko3d.sh
index 68439928f64b33cb290aa83fdbebce4712b8c16f..16efd2b679210b7158962f6d777d975fa007fc0b 100755
--- a/marchenko3D/demo/marchenko3D/oneD/marchenko3d.sh
+++ b/marchenko3D/demo/marchenko3D/oneD/marchenko3d.sh
@@ -1,19 +1,19 @@
#!/bin/bash
#Marchenko using time data
-#marchenko3D file_shot=reflx10y10.su verbose=2 \
+#../../../marchenko3D file_shot=reflx10y10.su verbose=2 \
# file_tinv=farr.su \
# niter=10 shift=15 smooth=10 hw=5 \
# file_green=green_timeshot.su file_file_f2=f2_timeshot.su
#Marchenko using frequency data
-#marchenko3D file_shotw=reflx10y10_W.bin verbose=2 \
+#../../../marchenko3D file_shotw=reflx10y10_W.bin verbose=2 \
# file_tinv=farr.su \
# niter=10 shift=15 smooth=10 hw=5 \
# file_green=green_freqshot.su file_file_f2=f2_freqshot.su
#Marchenko using zfp data
-marchenko3D file_shotzfp=reflx10y10_zfp.bin verbose=2 \
+../../../marchenko3D file_shotzfp=reflx10y10_zfp.bin verbose=2 \
file_tinv=farr.su \
niter=10 shift=15 smooth=10 hw=5 \
file_green=green_zfpshot.su file_file_f2=f2_zfpshot.su
diff --git a/marchenko3D/demo/marchenko3D/oneD/model.sh b/marchenko3D/demo/marchenko3D/oneD/model.sh
index 3725e0aff7c937946560dccc3f1a64d6720f4379..b5affbc7a72492b89d2f2d591cb6eba3b6cd4743 100755
--- a/marchenko3D/demo/marchenko3D/oneD/model.sh
+++ b/marchenko3D/demo/marchenko3D/oneD/model.sh
@@ -8,14 +8,14 @@ ny=181
dy=10000
#define gridded model for FD computations
-makemod sizex=5000 sizez=1400 dx=$dx dz=$dx cp0=1800 ro0=1000 \
+../../../../utils/makemod sizex=5000 sizez=1400 dx=$dx dz=$dx cp0=1800 ro0=1000 \
orig=-2500,0 file_base=model5.su verbose=2 \
intt=def x=-2500,2500 z=400,400 poly=0 cp=2300 ro=3000 \
intt=def x=-2500,2500 z=700,700 poly=0 cp=2000 ro=1100 \
intt=def x=-2500,2500 z=1100,1100 poly=0 cp=2500 ro=4000
#define homogenoeus model to compute direct wave only
-makemod sizex=5000 sizez=300 dx=$dx dz=$dx cp0=1800 ro0=1000 \
+../../../../utils/makemod sizex=5000 sizez=300 dx=$dx dz=$dx cp0=1800 ro0=1000 \
orig=-2500,0 file_base=hom.su verbose=2
rm cp3d.su
diff --git a/marchenko3D/demo/marchenko3D/oneD/plane_horizontal/execute_plane.sh b/marchenko3D/demo/marchenko3D/oneD/plane_horizontal/execute_plane.sh
new file mode 100755
index 0000000000000000000000000000000000000000..7a5d864cf1433b2f4aad5979a73fae77a13c46ba
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/plane_horizontal/execute_plane.sh
@@ -0,0 +1,9 @@
+#!/bin/bash
+
+echo "Before running this file, run model.sh, makeR.sh and zfp.sh in the folder above"
+
+./farrmod_plane.sh
+
+./fmute3D_plane.sh
+
+./marchenko3d_plane.sh
diff --git a/marchenko3D/demo/marchenko3D/oneD/plane_horizontal/farrmod_plane.sh b/marchenko3D/demo/marchenko3D/oneD/plane_horizontal/farrmod_plane.sh
new file mode 100755
index 0000000000000000000000000000000000000000..87e1ba2b5b5135abae9475cc7dd3855a1934a199
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/plane_horizontal/farrmod_plane.sh
@@ -0,0 +1,36 @@
+#!/bin/bash
+
+dt=0.001
+
+makewave fp=15 fmax=22 dt=${dt} file_out=wavefpmod.su nt=8192 t0=0.2 scale=1 scfft=0
+
+export OMP_NUM_THREADS=8
+
+#Model shot record in middle of model
+../../../../../fdelmodc3D/fdelmodc3D \
+ file_cp=../cp3d.su ischeme=1 iorder=4 \
+ file_den=../ro3d.su \
+ file_src=wavefpmod.su \
+ file_rcv=farrmod_plane.su \
+ plane_wave=1 npxsrc=201 npysrc=61 \
+ src_anglex=0 src_angley=0 \
+ src_velox=2170 src_veloy=2170 \
+ src_nxwindow=0 src_nywindow=0 \
+ src_type=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.2 \
+ verbose=2 \
+ tmod=1.2 \
+ dxrcv=10.0 \
+ xrcv1=-1000 xrcv2=1000 \
+ dyrcv=10.0 \
+ yrcv1=-300 yrcv2=300 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 ysrc=0 zsrc=900 \
+ ntaper=61 \
+ left=4 right=4 top=4 bottom=4 front=4 back=4 \
diff --git a/marchenko3D/demo/marchenko3D/oneD/plane_horizontal/fmute3D_plane.sh b/marchenko3D/demo/marchenko3D/oneD/plane_horizontal/fmute3D_plane.sh
new file mode 100755
index 0000000000000000000000000000000000000000..a6a56dc878fa615c35d52d36e09e48cb14082978
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/plane_horizontal/fmute3D_plane.sh
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+#Mute the coda from the modeled first arrival
+../../../../fmute3D file_mute=farrmod_plane_rp.su file_shot=farrmod_plane_rp.su file_out=farr_plane.su above=2 shift=15 smooth=10 hw=5
diff --git a/marchenko3D/demo/marchenko3D/oneD/plane_horizontal/marchenko3d_plane.sh b/marchenko3D/demo/marchenko3D/oneD/plane_horizontal/marchenko3d_plane.sh
new file mode 100755
index 0000000000000000000000000000000000000000..b7c2d3a7376a360610c6806cdd479d1e997c9cb4
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/plane_horizontal/marchenko3d_plane.sh
@@ -0,0 +1,22 @@
+#!/bin/bash
+
+#Marchenko using time data
+#../../../../marchenko3D file_shot=../reflx10y10.su verbose=2 \
+# file_tinv=farr_plane.su \
+# plane_wave=1 src_anglex=0 src_angley=0 src_velox=2170 src_veloy=2170 \
+# niter=10 shift=15 smooth=10 hw=5 \
+# file_green=green_plane_timeshot.su file_file_f2=f2_plane_timeshot.su
+
+#Marchenko using frequency data
+#../../../../marchenko3D file_shotw=../reflx10y10_W.bin verbose=2 \
+# file_tinv=farr_plane.su \
+# plane_wave=1 src_anglex=0 src_angley=0 src_velox=2170 src_veloy=2170 \
+# niter=10 shift=15 smooth=10 hw=5 \
+# file_green=green_plane_freqshot.su file_file_f2=f2_plane_freqshot.su
+
+#Marchenko using zfp data
+../../../../marchenko3D file_shotzfp=../reflx10y10_zfp.bin verbose=2 \
+ file_tinv=farr_plane.su \
+ plane_wave=1 src_anglex=0 src_angley=0 src_velox=2170 src_veloy=2170 \
+ niter=10 shift=15 smooth=10 hw=5 \
+ file_green=green_plane_zfpshot.su file_file_f2=f2_plane_zfpshot.su
diff --git a/marchenko3D/demo/marchenko3D/oneD/plane_tilted/execute_tilted.sh b/marchenko3D/demo/marchenko3D/oneD/plane_tilted/execute_tilted.sh
new file mode 100755
index 0000000000000000000000000000000000000000..a8762fc23185c36ec41ad7fcb726869c493b98d0
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/plane_tilted/execute_tilted.sh
@@ -0,0 +1,18 @@
+#!/bin/bash
+
+echo "Before running this file, run model.sh, makeR.sh and zfp.sh in the folder above"
+
+./farrmod_tilted_pos.sh
+
+./farrmod_tilted_neg.sh
+
+./fmute3D_tilted_pos.sh
+
+./fmute3D_tilted_neg.sh
+
+./marchenko3d_tilted_pos.sh
+
+./marchenko3d_tilted_neg.sh
+
+susum gplus_tilted_neg_zfpshot.su gmin_tilted_pos_zfpshot.su > green_total_zfpshot_neg.su
+susum gplus_tilted_pos_zfpshot.su gmin_tilted_neg_zfpshot.su > green_total_zfpshot_pos.su
diff --git a/marchenko3D/demo/marchenko3D/oneD/plane_tilted/farrmod_tilted_neg.sh b/marchenko3D/demo/marchenko3D/oneD/plane_tilted/farrmod_tilted_neg.sh
new file mode 100755
index 0000000000000000000000000000000000000000..2ade699a3a4f1442796f73e77d43ed44e62e225e
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/plane_tilted/farrmod_tilted_neg.sh
@@ -0,0 +1,36 @@
+#!/bin/bash
+
+dt=0.001
+
+makewave fp=15 fmax=22 dt=${dt} file_out=wavefpmod.su nt=8192 t0=0.2 scale=1 scfft=0
+
+export OMP_NUM_THREADS=8
+
+#Model shot record in middle of model
+../../../../../fdelmodc3D/fdelmodc3D \
+ file_cp=../cp3d.su ischeme=1 iorder=4 \
+ file_den=../ro3d.su \
+ file_src=wavefpmod.su \
+ file_rcv=farrmod_tilted_neg.su \
+ plane_wave=1 npxsrc=201 npysrc=61 \
+ src_anglex=-10 src_angley=-5 \
+ src_velox=2170 src_veloy=2170 \
+ src_nxwindow=0 src_nywindow=0 \
+ src_type=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.2 \
+ verbose=2 \
+ tmod=1.2 \
+ dxrcv=10.0 \
+ xrcv1=-1000 xrcv2=1000 \
+ dyrcv=10.0 \
+ yrcv1=-300 yrcv2=300 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 ysrc=0 zsrc=900 \
+ ntaper=61 \
+ left=4 right=4 top=4 bottom=4 front=4 back=4 \
diff --git a/marchenko3D/demo/marchenko3D/oneD/plane_tilted/farrmod_tilted_pos.sh b/marchenko3D/demo/marchenko3D/oneD/plane_tilted/farrmod_tilted_pos.sh
new file mode 100755
index 0000000000000000000000000000000000000000..0de3ae65b2705ab39cd53eaa062bbe2413769c36
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/plane_tilted/farrmod_tilted_pos.sh
@@ -0,0 +1,36 @@
+#!/bin/bash
+
+dt=0.001
+
+makewave fp=15 fmax=22 dt=${dt} file_out=wavefpmod.su nt=8192 t0=0.2 scale=1 scfft=0
+
+export OMP_NUM_THREADS=8
+
+#Model shot record in middle of model
+../../../../../fdelmodc3D/fdelmodc3D \
+ file_cp=../cp3d.su ischeme=1 iorder=4 \
+ file_den=../ro3d.su \
+ file_src=wavefpmod.su \
+ file_rcv=farrmod_tilted_pos.su \
+ plane_wave=1 npxsrc=201 npysrc=61 \
+ src_anglex=10 src_angley=5 \
+ src_velox=2170 src_veloy=2170 \
+ src_nxwindow=0 src_nywindow=0 \
+ src_type=1 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.2 \
+ verbose=2 \
+ tmod=1.2 \
+ dxrcv=10.0 \
+ xrcv1=-1000 xrcv2=1000 \
+ dyrcv=10.0 \
+ yrcv1=-300 yrcv2=300 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 ysrc=0 zsrc=900 \
+ ntaper=61 \
+ left=4 right=4 top=4 bottom=4 front=4 back=4 \
diff --git a/marchenko3D/demo/marchenko3D/oneD/plane_tilted/fmute3D_tilted_neg.sh b/marchenko3D/demo/marchenko3D/oneD/plane_tilted/fmute3D_tilted_neg.sh
new file mode 100755
index 0000000000000000000000000000000000000000..9d5a622e98591eaf969312441e480ecbdac856be
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/plane_tilted/fmute3D_tilted_neg.sh
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+#Mute the coda from the modeled first arrival
+../../../../fmute3D file_mute=farrmod_tilted_neg_rp.su file_shot=farrmod_tilted_neg_rp.su file_out=farr_tilted_neg.su above=2 shift=15 smooth=10 hw=5
diff --git a/marchenko3D/demo/marchenko3D/oneD/plane_tilted/fmute3D_tilted_pos.sh b/marchenko3D/demo/marchenko3D/oneD/plane_tilted/fmute3D_tilted_pos.sh
new file mode 100755
index 0000000000000000000000000000000000000000..a835de535bdbf82b04525710bb31c735ce584e29
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/plane_tilted/fmute3D_tilted_pos.sh
@@ -0,0 +1,4 @@
+#!/bin/bash
+
+#Mute the coda from the modeled first arrival
+../../../../fmute3D file_mute=farrmod_tilted_pos_rp.su file_shot=farrmod_tilted_pos_rp.su file_out=farr_tilted_pos.su above=2 shift=15 smooth=10 hw=5
diff --git a/marchenko3D/demo/marchenko3D/oneD/plane_tilted/marchenko3d_tilted_neg.sh b/marchenko3D/demo/marchenko3D/oneD/plane_tilted/marchenko3d_tilted_neg.sh
new file mode 100755
index 0000000000000000000000000000000000000000..bd7f41352a97e65718b8653d8f60e2a7946d3b8c
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/plane_tilted/marchenko3d_tilted_neg.sh
@@ -0,0 +1,25 @@
+#!/bin/bash
+
+#Marchenko using time data
+#../../../../marchenko3D file_shot=../reflx10y10.su verbose=2 \
+# file_tinv=farr_tilted_neg.su \
+# plane_wave=1 src_anglex=-10 src_angley=-5 src_velox=2170 src_veloy=2170 \
+# niter=10 shift=15 smooth=10 hw=5 \
+# file_green=green_tilted_neg_timeshot.su file_file_f2=f2_tilted_neg_timeshot.su \
+# file_gplus=gplus_tilted_neg_timeshot.su file_file_gmin=gmin_tilted_neg_timeshot.su
+
+#Marchenko using frequency data
+#../../../../marchenko3D file_shotw=../reflx10y10_W.bin verbose=2 \
+# file_tinv=farr_tilted_neg.su \
+# plane_wave=1 src_anglex=-10 src_angley=-5 src_velox=2170 src_veloy=2170 \
+# niter=10 shift=15 smooth=10 hw=5 \
+# file_green=green_tilted_neg_freqshot.su file_file_f2=f2_tilted_neg_freqshot.su \
+# file_gplus=gplus_tilted_neg_freqshot.su file_file_gmin=gmin_tilted_neg_freqshot.su
+
+#Marchenko using zfp data
+../../../../marchenko3D file_shotzfp=../reflx10y10_zfp.bin verbose=2 \
+ file_tinv=farr_tilted_neg.su \
+ plane_wave=1 src_anglex=-10 src_angley=-5 src_velox=2170 src_veloy=2170 \
+ niter=10 shift=15 smooth=10 hw=5 \
+ file_green=green_tilted_neg_zfpshot.su file_file_f2=f2_tilted_neg_zfpshot.su \
+ file_gplus=gplus_tilted_neg_zfpshot.su file_gmin=gmin_tilted_neg_zfpshot.su
diff --git a/marchenko3D/demo/marchenko3D/oneD/plane_tilted/marchenko3d_tilted_pos.sh b/marchenko3D/demo/marchenko3D/oneD/plane_tilted/marchenko3d_tilted_pos.sh
new file mode 100755
index 0000000000000000000000000000000000000000..3aa1a48dc34bf972f2ef8bb5065fd5265178f0bc
--- /dev/null
+++ b/marchenko3D/demo/marchenko3D/oneD/plane_tilted/marchenko3d_tilted_pos.sh
@@ -0,0 +1,25 @@
+#!/bin/bash
+
+#Marchenko using time data
+#../../../../marchenko3D file_shot=../reflx10y10.su verbose=2 \
+# file_tinv=farr_tilted_pos.su \
+# plane_wave=1 src_anglex=10 src_angley=5 src_velox=2170 src_veloy=2170 \
+# niter=10 shift=15 smooth=10 hw=5 \
+# file_green=green_tilted_pos_timeshot.su file_file_f2=f2_tilted_pos_timeshot.su \
+# file_gplus=gplus_tilted_pos_timeshot.su file_gmin=gmin_tilted_pos_timeshot.su
+
+#Marchenko using frequency data
+#../../../../marchenko3D file_shotw=../reflx10y10_W.bin verbose=2 \
+# file_tinv=farr_tilted_pos.su \
+# plane_wave=1 src_anglex=10 src_angley=5 src_velox=2170 src_veloy=2170 \
+# niter=10 shift=15 smooth=10 hw=5 \
+# file_green=green_tilted_pos_freqshot.su file_file_f2=f2_tilted_pos_freqshot.su
+# file_gplus=gplus_tilted_pos_freqshot.su file_gmin=gmin_tilted_pos_freqshot.su
+
+#Marchenko using zfp data
+../../../../marchenko3D file_shotzfp=../reflx10y10_zfp.bin verbose=2 \
+ file_tinv=farr_tilted_pos.su \
+ plane_wave=1 src_anglex=10 src_angley=5 src_velox=2170 src_veloy=2170 \
+ niter=10 shift=15 smooth=10 hw=5 \
+ file_green=green_tilted_pos_zfpshot.su file_file_f2=f2_tilted_pos_zfpshot.su \
+ file_gplus=gplus_tilted_pos_zfpshot.su file_gmin=gmin_tilted_pos_zfpshot.su
diff --git a/marchenko3D/demo/marchenko3D/oneD/shotmodel.sh b/marchenko3D/demo/marchenko3D/oneD/shotmodel.sh
index cd6c947964220e7a2bbfd5573a8d9adcae7eb9bc..4e8453ca11b0983e4430cddc3b9147eb26010836 100755
--- a/marchenko3D/demo/marchenko3D/oneD/shotmodel.sh
+++ b/marchenko3D/demo/marchenko3D/oneD/shotmodel.sh
@@ -9,7 +9,7 @@ makewave w=fw fmin=0 flef=5 frig=25 fmax=30 dt=$dt file_out=wavefw.su nt=4096 t0
export OMP_NUM_THREADS=8
#Model shot record in middle of model
-fdelmodc3D \
+../../../../fdelmodc3D/fdelmodc3D \
file_cp=cp3d.su ischeme=1 iorder=4 \
file_den=ro3d.su \
file_src=wavefw.su \
@@ -34,7 +34,7 @@ fdelmodc3D \
left=4 right=4 top=4 bottom=4 front=4 back=4
#Model direct wave only in middle of model
-fdelmodc3D \
+../../../../fdelmodc3D/fdelmodc3D \
file_cp=cphom3d.su ischeme=1 iorder=4 \
file_den=rohom3d.su \
file_src=wavefw.su \
diff --git a/marchenko3D/demo/marchenko3D/oneD/zfp.sh b/marchenko3D/demo/marchenko3D/oneD/zfp.sh
index 45f436827b4c75bfecd7b5c535053002b2930443..2fd8fb88d3b2651fac969afdd2487536b452ba0c 100755
--- a/marchenko3D/demo/marchenko3D/oneD/zfp.sh
+++ b/marchenko3D/demo/marchenko3D/oneD/zfp.sh
@@ -1,9 +1,15 @@
#!/bin/bash
# Pre-transform the data to the frequency domain
-#TWtransform file_in=reflx10y10.su file_out=reflx10y10_W.bin verbose=1 fmin=0 fmax=30 zfp=0 tolerance=1e-7
+#../../../TWtransform file_in=reflx10y10.su file_out=reflx10y10_W.bin verbose=1 fmin=0 fmax=30 zfp=0 tolerance=1e-7
# Pre-transform the data to the frequency domain and apply zfp compression
-TWtransform file_in=reflx10y10.su file_out=reflx10y10_zfp.bin verbose=1 fmin=0 fmax=30 zfp=1 tolerance=1e-7
+../../../TWtransform file_in=reflx10y10.su file_out=reflx10y10_zfp.bin verbose=1 fmin=0 fmax=30 zfp=1 tolerance=1e-7
+if [ -e "reflx10y10_zfp.bin" ]
+then
+echo "Data are compressed. Deleting uncompressed data"
rm reflx10y10.su
+else
+echo "Compressed data not found. Retaining uncompressed data"
+fi
diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index 4ed77c7e5d14e15ef267d351cdb5312deddfbbd5..f8fd501224a60f4f22d7b5e9b76a3c19a423ee4d 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -397,7 +397,7 @@ int main (int argc, char **argv)
px = sin(src_anglex*grad2rad)/src_velox;
py = sin(src_angley*grad2rad)/src_veloy;
- tshift = fabs((nys-1)*dys*py) + fabs((nxs-1)*dxs*px);
+ tshift = fabs((nys-1)*dys*py) + fabs((nxs-1)*dxs*px) - dt;
for (j = 0; j < Nfoc; j++) {
itmin[j] = nt;
@@ -602,7 +602,13 @@ int main (int argc, char **argv)
vmess("number of time samples (nt,nts) = %li (%li,%li)", ntfft, nt, nts);
vmess("frequency cutoffs = min:%.3f max:%.3f",fmin,fmax);
vmess("time sampling = %e ", dt);
- if (plane_wave) vmess("Plane wave focusing is applied");
+ if (plane_wave) {
+ vmess("Plane wave focusing is applied");
+ vmess("Plane wave angle = x:%.2f y%.2f",src_anglex,src_angley);
+ vmess("Plane wave velocity = x:%.2f y%.2f",src_velox,src_veloy);
+ vmess("Plane wave p value = x:%.6f y%.6f",px,py);
+ vmess("tshift = x:%.3f",tshift);
+ }
if (file_green != NULL) vmess("Green output file = %s ", file_green);
if (file_gmin != NULL) vmess("Gmin output file = %s ", file_gmin);
if (file_gplus != NULL) vmess("Gplus output file = %s ", file_gplus);
@@ -781,7 +787,14 @@ int main (int argc, char **argv)
}
}
}
- if (above==-2) applyMute3D(f1min, muteW, smooth, 0, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
+ if (above==-2) {
+ if (plane_wave==1) {
+ applyMute3D_tshift(f1min, muteW, smooth, 0, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, iter, tsynW);
+ }
+ else {
+ applyMute3D(f1min, muteW, smooth, 0, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
+ }
+ }
}
else {/* odd iterations update: => f_1^+(t) */
for (l = 0; l < Nfoc; l++) {
@@ -878,17 +891,6 @@ int main (int argc, char **argv)
if (plane_wave==1) {
applyMute3D_tshift(Gmin, muteW, smooth, 4, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, 0, tsynW);
/* for plane wave with angle shift itmin downward */
- for (l = 0; l < Nfoc; l++) {
- for (i = 0; i < npos; i++) {
- memcpy(&trace[0],&Gmin[l*nys*nxs*nts+i*nts],nts*sizeof(float));
- for (j = 0; j < itmin[l]; j++) {
- Gmin[l*nys*nxs*nts+i*nts+j] = 0.0;
- }
- for (j = 0; j < nts-itmin[l]; j++) {
- Gmin[l*nys*nxs*nts+i*nts+j+itmin[l]] = trace[j];
- }
- }
- }
timeShift(Gmin, nts, npos*Nfoc, dt, tshift, fmin, fmax);
}
else {
diff --git a/raytime3d/raytime3d.c b/raytime3d/raytime3d.c
index a689661c756702dbe87911eb696e6f78aed07060..801b4dfd73450291e5ab2a22c9a2048fd486cfb5 100644
--- a/raytime3d/raytime3d.c
+++ b/raytime3d/raytime3d.c
@@ -324,7 +324,7 @@ private (is,time0,ampl,nrx,nry,nrz,nr,cp_average,i,j,k,ix,iy,iz,hdrs,tmpdata,nwr
if (ret < 0 ) verr("error on writing output file.");
if (ray.geomspread==1) {
- ret = writeData3D(fpt, &l[0], hdrs, rec.nx, rec.nz*rec.ny);
+ ret = writeData3D(fpa, &l[0], hdrs, rec.nx, rec.nz*rec.ny);
if (ret < 0 ) verr("error on writing output file.");
}
diff --git a/utils/HomG.c b/utils/HomG.c
index 52cf1864bc34b440f49e4c60e3b92129a79d9a98..87712f6e181f442a5d6f9733ee3ca3aaef55839e 100644
--- a/utils/HomG.c
+++ b/utils/HomG.c
@@ -65,6 +65,8 @@ void scl_data3D(float *data, long nt, long nx, long ny, float scl, float *datout
void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout);
void pad_data(float *data, long nsam, long nrec, long nsamout, float *datout);
void convol(float *data1, float *data2, float *con, long nrec, long nsam, float dt, long shift);
+void deconv(float *data1, float *data2, float *decon, long nrec, long nsam,
+ float dt, float eps, float reps, long shift);
void corr(float *data1, float *data2, float *cov, long nrec, long nsam, float dt, long shift);
void timeDiff(float *data, long nsam, long nrec, float dt, float fmin, float fmax, long opt);
void depthDiff(float *data, long nsam, long nrec, float dt, float dx, float fmin, float fmax, float c, long opt);
@@ -97,6 +99,9 @@ char *sdoc[] = {
" zfpr=0 ................... virtual receiver data are in SU format (=0) or zfp compressed (=1)",
" cp=1000.0 ................ Velocity at the top of the medium in m/s",
" rho=1000.0 ............... Density at the top of the medium in kg/m^3",
+" file_wav= ................ Wavelet that is deconvolved from the virtual receiver data",
+" eps=0.0 .................. Absolute stabilization factor for deconvolution",
+" reps=0.01 ................ Relative stabilization factor for deconvolution (is multiplied by the maximum of the data)",
" scheme=0 ................. Scheme for the retrieval",
" .......................... scheme=0 Marchenko homogeneous Green's function retrieval with G source",
" .......................... scheme=1 Marchenko homogeneous Green's function retrieval with f2 source",
@@ -113,17 +118,17 @@ NULL};
int main (int argc, char **argv)
{
- FILE *fp_in, *fp_shot, *fp_out;
- char *fin, *fshot, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100], *direction;
+ FILE *fp_in, *fp_shot, *fp_out, *fp_wav;
+ char *fin, *fshot, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100], *direction, *file_wav;
float *rcvdata, *Ghom, *shotdata, *shotdata_jkz, rho, fmin, fmax;
float dt, dy, dx, t0, y0, x0, xmin, xmax1, sclsxgx, dxrcv, dyrcv, dzrcv;
- float *conv, *conv2, *tmp1, *tmp2, cp, shift;
+ float *conv, *conv2, *decon, *tmp1, *tmp2, cp, shift, eps, reps;
long nshots, ntvs, nyvs, nxvs, ntraces, ret, ix, iy, it, is, ir, ig, file_det, verbose;
long ntr, nxr, nyr, nsr, i, l, j, k, nxvr, nyvr, nzvr, count, num, isn;
- float dtr, dxr, dyr, ftr, fxr, fyr, sclr, scl;
+ float dtr, dxr, dyr, ftr, fxr, fyr, sclr, scl, *wavelet;
long pos1, npos, zmax, numb, dnumb, scheme, ntmax, ntshift, shift_num, zfps, zfpr, size;
long ixr, iyr, zsrc, zrcv, *xvr, *yvr, *zvr;
- segy *hdr_rcv, *hdr_out, *hdr_shot;
+ segy *hdr_rcv, *hdr_out, *hdr_shot, *hdr_wav;
initargs(argc, argv);
requestdoc(1);
@@ -144,6 +149,9 @@ int main (int argc, char **argv)
if (!getparlong("dnumb", &dnumb)) dnumb=1;
if (!getparlong("scheme", &scheme)) scheme = 0;
if (!getparlong("verbose", &verbose)) verbose = 0;
+ if (!getparstring("file_wav", &file_wav)) file_wav = NULL;
+ if(!getparfloat("eps", &eps)) eps=0.00;
+ if(!getparfloat("reps", &reps)) reps=0.01;
if (!getparlong("zfps", &zfps)) zfps = 0;
if (!getparlong("zfpr", &zfpr)) zfpr = 0;
if (!getparstring("direction", &direction)) direction = "z";
@@ -281,6 +289,24 @@ int main (int argc, char **argv)
yvr = (long *)malloc(nxvr*nyvr*nzvr*sizeof(long));
zvr = (long *)malloc(nxvr*nyvr*nzvr*sizeof(long));
+ /*----------------------------------------------------------------------------*
+ * Read in the wavelet and create the gather
+ *----------------------------------------------------------------------------*/
+
+ if (file_wav) {
+ if (verbose) vmess("Reading in wavelet for deconvolution");
+ wavelet = (float *)calloc(ntvs*nxvs*nyvs,sizeof(float));
+ hdr_wav = (segy *)calloc(nxvs*nyvs,sizeof(segy));
+ readSnapData3D(file_wav, &wavelet[0], &hdr_wav[0], 1, 1, 1, ntvs, 0, 1, 0, 1, 0, ntvs);
+ for (i = 1; i < nxvs*nyvs; i++){
+ for (j = 0; j < ntvs; j++) {
+ wavelet[i*ntvs+j] = wavelet[j];
+ hdr_wav[i] = hdr_wav[0];
+ }
+ }
+ free(hdr_wav);
+ }
+
/*----------------------------------------------------------------------------*
* Get the file info for the source position
*----------------------------------------------------------------------------*/
@@ -371,6 +397,7 @@ int main (int argc, char **argv)
tmp2 = (float *)calloc(nyr*nxr*ntr,sizeof(float));
}
if (scheme==8 || scheme==9 || scheme==10) tmp1 = (float *)calloc(nyr*nxr*ntr,sizeof(float));
+ if (file_wav!=NULL) decon = (float *)calloc(nyr*nxr*ntr,sizeof(float));
sprintf(fins,"%s%li",direction,ir*dnumb+numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
@@ -439,7 +466,13 @@ int main (int argc, char **argv)
else if (scheme==3) { //Marchenko representation without time-reversal G source
if (nyr>1) depthDiff3D(&rcvdata[l*nyr*nxr*ntr], ntr, nxr, nyr, dt, dx, dy, fmin, fmax, cp, 1);
else depthDiff(&rcvdata[l*nyr*nxr*ntr], ntr, nxr, dt, dx, fmin, fmax, cp, 1);
- convol2(&shotdata[0], &rcvdata[l*nyr*nxr*ntr], conv, nyr*nxr, ntr, dt, fmin, fmax, 1);
+ if (file_wav!=NULL){
+ deconv(&rcvdata[l*nyr*nxr*ntr], wavelet, decon, nyr*nxr, ntr, dt, eps, reps, 0);
+ convol2(&shotdata[0], decon, conv, nyr*nxr, ntr, dt, fmin, fmax, 1);
+ }
+ else {
+ convol2(&shotdata[0], &rcvdata[l*nyr*nxr*ntr], conv, nyr*nxr, ntr, dt, fmin, fmax, 1);
+ }
for (i=0; i<nyr*nxr; i++) {
for (j=0; j<ntr/2; j++) {
Ghom[(j+ntr/2)*nxvr*nyvr*nzvr+l*nzvr+ir] += 2.0*scl*conv[i*ntr+j]/rho;
@@ -542,9 +575,11 @@ int main (int argc, char **argv)
free(tmp1); free(tmp2);
}
if (scheme==8 || scheme==9 || scheme==10) free(tmp1);
+ if (file_wav!=NULL) free(decon);
}
free(shotdata);
+ if (file_wav!=NULL) free(wavelet);
if (strcmp(direction,"z") == 0) {
if (nxvr>1) dxrcv = (float)((xvr[nzvr] - xvr[0])/1000.0);
@@ -1203,7 +1238,7 @@ void getVirRec(char *filename, long *nxs, long *nys, long *nxr, long *nyr, long
fp = fopen( filename, "r" );
if ( fp == NULL ) verr("Could not open %s",filename);
nread = fread(&hdr, 1, TRCBYTES, fp);
- if (nread != TRCBYTES) verr("Could not read the header of the input file");
+ if (nread != TRCBYTES) verr("Could not read the header of the VR file");
*nxs = 1;
*nys = 1;
@@ -1655,4 +1690,105 @@ void scl_data3D(float *data, long nt, long nx, long ny, float scl, float *datout
}
}
}
+}
+
+void deconv(float *data1, float *data2, float *decon, long nrec, long nsam,
+ float dt, float eps, float reps, long shift)
+{
+ long i, j, n, optn, nfreq, sign;
+ float df, dw, om, tau, *den, scl;
+ float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1, *rdata2, maxden, leps;
+ complex *cdata1, *cdata2, *cdec, tmp;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+
+ cdata1 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata1 == NULL) verr("memory allocation error for cdata1");
+ cdata2 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata2 == NULL) verr("memory allocation error for cdata2");
+ cdec = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdec == NULL) verr("memory allocation error for ccov");
+
+ rdata1 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata1 == NULL) verr("memory allocation error for rdata1");
+ rdata2 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata2 == NULL) verr("memory allocation error for rdata2");
+ den = (float *)malloc(nfreq*nrec*sizeof(float));
+ if (den == NULL) verr("memory allocation error for rdata1");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data1, nsam, nrec, optn, rdata1);
+ pad_data(data2, nsam, nrec, optn, rdata2);
+
+ /* forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata1[0], &cdata1[0], optn, nrec, optn, nfreq, sign);
+ rcmfft(&rdata2[0], &cdata2[0], optn, nrec, optn, nfreq, sign);
+
+ /* apply deconvolution */
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ n = nrec*nfreq;
+ maxden=0.0;
+ for (j = 0; j < n; j++) {
+ den[j] = *p2r**p2r + *p2i**p2i;
+ maxden = MAX(den[j], maxden);
+ p2r += 2;
+ p2i += 2;
+ }
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ qr = (float *) &cdec[0].r;
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ qi = qr + 1;
+ leps = reps*maxden+eps;
+ for (j = 0; j < n; j++) {
+
+ if (fabs(*p2r)>=fabs(*p2i)) {
+ *qr = (*p2r**p1r+*p2i**p1i)/(den[j]+leps);
+ *qi = (*p2r**p1i-*p2i**p1r)/(den[j]+leps);
+ } else {
+ *qr = (*p1r**p2r+*p1i**p2i)/(den[j]+leps);
+ *qi = (*p1i**p2r-*p1r**p2i)/(den[j]+leps);
+ }
+ qr += 2;
+ qi += 2;
+ p1r += 2;
+ p1i += 2;
+ p2r += 2;
+ p2i += 2;
+ }
+ free(cdata1);
+ free(cdata2);
+ free(den);
+
+ if (shift) {
+ df = 1.0/(dt*optn);
+ dw = 2*PI*df;
+ tau = dt*(nsam/2);
+ for (j = 0; j < nrec; j++) {
+ om = 0.0;
+ for (i = 0; i < nfreq; i++) {
+ tmp.r = cdec[j*nfreq+i].r*cos(om*tau) + cdec[j*nfreq+i].i*sin(om*tau);
+ tmp.i = cdec[j*nfreq+i].i*cos(om*tau) - cdec[j*nfreq+i].r*sin(om*tau);
+ cdec[j*nfreq+i] = tmp;
+ om += dw;
+ }
+ }
+ }
+
+ /* inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = dt/(float)optn;
+ crmfft(&cdec[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata1,optn,nrec,scl,decon,nsam);
+
+ free(cdec);
+ free(rdata1);
+ free(rdata2);
+ return;
}
\ No newline at end of file
diff --git a/utils/combine_induced.c b/utils/combine_induced.c
index b1c63a8679c6c1bbb2a847fd255df9cd1ab9c43a..3a6229258614f9b9bcc43993520f6c878012843c 100644
--- a/utils/combine_induced.c
+++ b/utils/combine_induced.c
@@ -48,6 +48,7 @@ char *sdoc[] = {
" nshift=0 ................. Sample number shift",
" shift=0.0 ................ Base shift in seconds for the data",
" dtshift=0.0 .............. Shift per gather in seconds",
+" opt=0 .................... Option for muting acausal events (=0) or no muting (=1)",
" verbose=1 ................ Give information about the process (=1) or not (=0)",
NULL};
@@ -58,7 +59,7 @@ int main (int argc, char **argv)
float *indata, *outdata, *loopdata, shift, dtshift;
float dt, dz, dy, dx, t0, x0, y0, z0, scl, dxrcv, dyrcv, dzrcv;
long nt, nz, ny, nx, nxyz, ntr, ix, iy, it, is, iz, pos, file_det, nxs, nys, nzs;
- long numb, dnumb, ret, nzmax, verbose, nt_out, ishift, nshift, *sx, *sy, *sz;
+ long numb, dnumb, ret, nzmax, verbose, nt_out, ishift, nshift, *sx, *sy, *sz, opt;
segy *hdr_in, *hdr_loop, *hdr_out;
initargs(argc, argv);
@@ -74,6 +75,7 @@ int main (int argc, char **argv)
if (!getparlong("nzmax", &nzmax)) nzmax=0;
if (!getparlong("verbose", &verbose)) verbose=0;
if (!getparlong("nshift", &nshift)) nshift=0;
+ if (!getparlong("opt", &opt)) opt=0;
if (!getparfloat("shift", &shift)) shift=0.0;
if (!getparfloat("dtshift", &dtshift)) dtshift=0.0;
if (fin == NULL) verr("Incorrect downgoing input");
@@ -137,6 +139,8 @@ int main (int argc, char **argv)
vmess("Starting distance for x: %.3f, y: %.3f, z: %.3f",x0,y0,z0);
vmess("Sampling distance for x: %.3f, y: %.3f, z: %.3f",dx,dy,dz);
vmess("Number of virtual sources: %li",nzs);
+ if (opt==0) vmess("Muting is applied");
+ if (opt==1) vmess("Muting is not applied");
}
/*----------------------------------------------------------------------------*
@@ -165,6 +169,8 @@ int main (int argc, char **argv)
hdr_out = (segy *)calloc(nx*ny,sizeof(segy));
outdata = (float *)calloc(nt_out*nxyz,sizeof(float));
+ if (dtshift==0.0) dtshift=((float) nshift)*dt;
+
/*----------------------------------------------------------------------------*
* Parallel loop for reading in and combining the various shots
*----------------------------------------------------------------------------*/
@@ -189,12 +195,25 @@ int main (int argc, char **argv)
sz[is] = hdr_loop[0].sdepth;
ishift = nshift*is;
- if (verbose) vmess("Shifting %li timesteps for a total of %.3f seconds",ishift,shift+(dtshift*((float)iz)));
- for (it = ishift; it < nt_out; it++) {
- for (iy = 0; iy < ny; iy++) {
- for (ix = 0; ix < nx; ix++) {
- for (iz = 0; iz < nz; iz++) {
- outdata[it*nxyz+iy*nx*nz+ix*nz+iz] += loopdata[(it-ishift+(nt/2))*nxyz+iy*nx*nz+ix*nz+iz];
+ if (verbose) vmess("Shifting %li timesteps for a total of %.3f seconds",ishift,shift+(dtshift*((float)is)));
+ if (opt==0) {
+ for (it = ishift; it < nt_out; it++) {
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (iz = 0; iz < nz; iz++) {
+ outdata[it*nxyz+iy*nx*nz+ix*nz+iz] += loopdata[(it-ishift+(nt/2))*nxyz+iy*nx*nz+ix*nz+iz];
+ }
+ }
+ }
+ }
+ }
+ else {
+ for (it = 0; it < nt_out; it++) {
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (iz = 0; iz < nz; iz++) {
+ outdata[it*nxyz+iy*nx*nz+ix*nz+iz] += loopdata[(it-ishift+(nt/2))*nxyz+iy*nx*nz+ix*nz+iz];
+ }
}
}
}
diff --git a/utils/convhomg.c b/utils/convhomg.c
index 2bd57624d4b95dd4f68053bbaca5d862dfeb2cc7..e328f0ec666108871913735e2854ce70b4aef74f 100644
--- a/utils/convhomg.c
+++ b/utils/convhomg.c
@@ -31,6 +31,7 @@ long readSnapData3D(char *filename, float *data, segy *hdrs, long nsnaps, long n
long sx, long ex, long sy, long ey, long sz, long ez);
void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout);
void convol(float *data1, float *data2, float *con, long ntfft);
+void deconv(float *data1, float *data2, float *con, long ntfft, float eps, float reps);
char *sdoc[] = {
" ",
@@ -41,22 +42,25 @@ char *sdoc[] = {
" ",
" Required parameters: ",
"",
-" file_hom= ................. First file of the array of virtual receivers",
-" file_wav= ................. File containing the virtual source",
+" file_hom= ................. File containing the homogeneous Green's function",
+" file_wavcon= .............. File containing the wavelet that will be used for convolution",
" ",
" Optional parameters: ",
" ",
+" file_wavdecon= ........... File containing the wavelet that will be used for deconvolution",
+" eps=0.0 .................. Absolute stabilization factor for deconvolution",
+" reps=0.01 ................ Relative stabilization factor for deconvolution (is multiplied by the maximum of the data)",
" file_out= ................ Filename of the output",
NULL};
int main (int argc, char **argv)
{
FILE *fp_hom, *fp_wav, *fp_out;
- char *file_hom, *file_wav, *file_out;
+ char *file_hom, *file_wavcon, *file_wavdecon, *file_out;
long nt, nx, ny, nz, ntwav, ntfft, ntr, verbose;
long nt_wav, it, ipos;
float dt, dx, dy, dz, x0, y0, z0, scl, *Ghom, *wav, *wavelet, dt_wav;
- float *trace, *conv;
+ float *trace, *conv, *decon, *tmp, eps, reps;
segy *hdr_hom, hdr_wav;
size_t nread;
@@ -67,12 +71,15 @@ int main (int argc, char **argv)
* Get the parameters passed to the function
*----------------------------------------------------------------------------*/
if (!getparstring("file_hom", &file_hom)) file_hom = NULL;
- if (!getparstring("file_wav", &file_wav)) file_wav = NULL;
+ if (!getparstring("file_wavcon", &file_wavcon)) file_wavcon = NULL;
+ if (!getparstring("file_wavdecon", &file_wavdecon)) file_wavdecon = NULL;
if (!getparstring("file_out", &file_out)) file_out = "out.su";
+ if (!getparfloat("eps", &eps)) eps = 0.00;
+ if (!getparfloat("reps", &reps)) reps = 0.01;
if (!getparlong("verbose", &verbose)) verbose = 1;
if (file_hom==NULL) verr("Error file_hom is not given");
- if (file_wav==NULL) verr("Error file_wav is not given");
+ if (file_wavcon==NULL) verr("Error file_wav is not given");
/*----------------------------------------------------------------------------*
* Get the file info of the data and determine the indez of the truncation
@@ -91,8 +98,8 @@ int main (int argc, char **argv)
* Allocate and read in the data
*----------------------------------------------------------------------------*/
- fp_wav = fopen( file_wav, "r" );
- if (fp_wav == NULL) verr("File %s does not exist or cannot be opened", file_wav);
+ fp_wav = fopen( file_wavcon, "r" );
+ if (fp_wav == NULL) verr("File %s does not exist or cannot be opened", file_wavcon);
nread = fread( &hdr_wav, 1, TRCBYTES, fp_wav );
assert(nread == TRCBYTES);
nt_wav = hdr_wav.ns;
@@ -121,18 +128,51 @@ int main (int argc, char **argv)
wavelet = (float *)calloc(ntfft,sizeof(float));
trace = (float *)calloc(ntfft,sizeof(float));
+ tmp = (float *)calloc(ntfft,sizeof(float));
conv = (float *)calloc(ntfft,sizeof(float));
for (it = 0; it < nt_wav/2; it++) {
wavelet[it] = wav[it];
wavelet[ntfft-1-it] = wav[nt_wav-1-it];
}
- free(wav);
+
+ if (file_wavdecon!=NULL) {
+ if (verbose) vmess("Reading in deconvolution wavelet");
+ fp_wav = fopen( file_wavdecon, "r" );
+ if (fp_wav == NULL) verr("File %s does not exist or cannot be opened", file_wavdecon);
+ nread = fread( &hdr_wav, 1, TRCBYTES, fp_wav );
+ assert(nread == TRCBYTES);
+ nt_wav = hdr_wav.ns;
+ dt_wav = (float)(hdr_wav.dt/1E6);
+ wav = (float *)calloc(nt_wav,sizeof(float));
+ nread = fread(&wav[0], sizeof(float), nt_wav, fp_wav);
+ assert(nread==nt_wav);
+ fclose(fp_wav);
+
+ decon = (float *)calloc(ntfft,sizeof(float));
+
+ for (it = 0; it < nt_wav/2; it++) {
+ decon[it] = wav[it];
+ decon[ntfft-1-it] = wav[nt_wav-1-it];
+ }
+ vmess("Time sampling wavelet decon : %.3f",dt_wav);
+ if (dt_wav!=dt) vmess("WARNING! dt of HomG (%f) and deconvolution wavelet (%f) do not match.",dt,dt_wav);
+ }
+ free(wav);
+
for (ipos = 0; ipos<nx*ny*nz; ipos++) {
- for (it = 0; it < nt; it++) {
- trace[it] = Ghom[it*nx*ny*nz+ipos];
- }
- convol(trace,wavelet,conv,ntfft);
+ if (file_wavdecon!=NULL) {
+ for (it = 0; it < nt; it++) {
+ tmp[it] = Ghom[it*nx*ny*nz+ipos];
+ }
+ deconv(tmp,decon,trace,ntfft,eps,reps);
+ }
+ else{
+ for (it = 0; it < nt; it++) {
+ trace[it] = Ghom[it*nx*ny*nz+ipos];
+ }
+ }
+ convol(trace,wavelet,conv,ntfft);
for (it = 0; it < nt; it++) {
Ghom[it*nx*ny*nz+ipos] = conv[it];
}
@@ -143,6 +183,7 @@ int main (int argc, char **argv)
}
free(wavelet); free(trace); free(conv);
+ if (file_wavdecon!=NULL) free(decon);
fp_out = fopen(file_out, "w+");
@@ -158,7 +199,7 @@ int main (int argc, char **argv)
void convol(float *data1, float *data2, float *con, long ntfft)
{
- long i, j, n, nfreq, sign;
+ long i, j, nfreq, sign;
float scl;
float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1;
complex *cdata1, *cdata2, *ccon;
@@ -187,8 +228,7 @@ void convol(float *data1, float *data2, float *con, long ntfft)
p1i = p1r + 1;
p2i = p2r + 1;
qi = qr + 1;
- n = nfreq;
- for (j = 0; j < n; j++) {
+ for (j = 0; j < nfreq; j++) {
*qr = (*p2r**p1r-*p2i**p1i);
*qi = (*p2r**p1i+*p2i**p1r);
qr += 2;
@@ -212,6 +252,81 @@ void convol(float *data1, float *data2, float *con, long ntfft)
return;
}
+void deconv(float *data1, float *data2, float *con, long ntfft, float eps, float reps)
+{
+ long i, j, nfreq, sign;
+ float scl, maxden, *den, leps;
+ float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1;
+ complex *cdata1, *cdata2, *ccon;
+
+ nfreq = ntfft/2+1;
+
+ den = (float *)malloc(nfreq*sizeof(float));
+
+ cdata1 = (complex *)malloc(nfreq*sizeof(complex));
+ if (cdata1 == NULL) verr("memory allocation error for cdata1");
+ cdata2 = (complex *)malloc(nfreq*sizeof(complex));
+ if (cdata2 == NULL) verr("memory allocation error for cdata2");
+ ccon = (complex *)malloc(nfreq*sizeof(complex));
+ if (ccon == NULL) verr("memory allocation error for ccov");
+
+ rdata1 = (float *)malloc(ntfft*sizeof(float));
+ if (rdata1 == NULL) verr("memory allocation error for rdata1");
+
+ /* forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&data1[0], &cdata1[0], (int)ntfft, 1, (int)ntfft, (int)nfreq, (int)sign);
+ rcmfft(&data2[0], &cdata2[0], (int)ntfft, 1, (int)ntfft, (int)nfreq, (int)sign);
+
+ /* apply deconvolution */
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ maxden=0.0;
+ for (j = 0; j < nfreq; j++) {
+ den[j] = *p2r**p2r + *p2i**p2i;
+ maxden = MAX(den[j], maxden);
+ p2r += 2;
+ p2i += 2;
+ }
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ qr = (float *) &ccon[0].r;
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ qi = qr + 1;
+ leps = reps*maxden+eps;
+ for (j = 0; j < nfreq; j++) {
+
+ if (fabs(*p2r)>=fabs(*p2i)) {
+ *qr = (*p2r**p1r+*p2i**p1i)/(den[j]+leps);
+ *qi = (*p2r**p1i-*p2i**p1r)/(den[j]+leps);
+ } else {
+ *qr = (*p1r**p2r+*p1i**p2i)/(den[j]+leps);
+ *qi = (*p1i**p2r-*p1r**p2i)/(den[j]+leps);
+ }
+ qr += 2;
+ qi += 2;
+ p1r += 2;
+ p1i += 2;
+ p2r += 2;
+ p2i += 2;
+ }
+ free(cdata1);
+ free(cdata2);
+
+ /* inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/((float)(ntfft));
+ crmfft(&ccon[0], &rdata1[0], (int)ntfft, 1, (int)nfreq, (int)ntfft, (int)sign);
+ scl_data(rdata1,ntfft,1,scl,con,ntfft);
+
+ free(ccon);
+ free(rdata1);
+ return;
+}
+
void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout)
{
long it,ix;
diff --git a/utils/mutesnap.c b/utils/mutesnap.c
index 7afa393e19ddef9bd7f8c2354e8b6d742fff77b7..6102c1751490f45b362f9f7e615df4891f298242 100644
--- a/utils/mutesnap.c
+++ b/utils/mutesnap.c
@@ -61,7 +61,7 @@ int main (int argc, char **argv)
float *homdata, *snapdata, *outdata, *rtrace, *costaper, scl, tol, *timeval, dt;
float dxs, dys, dzs, scls, fzs, fxs, fys;
float dxh, dyh, dzh, sclh, fzh, fxh, fyh;
- float dxrcv, dyrcv, dzrcv, dxpos, offset;
+ float dxrcv, dyrcv, dzrcv, dxpos, offset, maxz;
long nts, nxs, nys, nzs, ntrs, nth, nxh, nyh, nzh, ntrh, opt;
long nxyz, nxy, ret, ix, iy, iz, it, ht, indrcv, shift, rmt, mode, smooth, verbose;
segy *hdr_hom, *hdr_snap, *hdrs_mute;
@@ -82,6 +82,7 @@ int main (int argc, char **argv)
if (!getparlong("verbose", &verbose)) verbose = 1;
if (!getparlong("opt", &opt)) opt = 0;
if (!getparfloat("tol", &tol)) tol = 5;
+ if (!getparfloat("maxz", &maxz)) maxz = 0;
if (fhom == NULL) verr("Incorrect G_hom input");
if (mode != 2) {
if (fsnap == NULL) verr("Incorrect snapshot input");
@@ -104,7 +105,7 @@ int main (int argc, char **argv)
/*----------------------------------------------------------------------------*
* Determine the parameters of the files and determine if they match
*----------------------------------------------------------------------------*/
- getFileInfo3D(fhom, &nzh, &nxh, &nyh, &nth, &dzh, &dxh, &dyh, &fzh, &fxh, &fyh, &sclh, &ntrh);
+ getFileInfo3D(fhom, &nzh, &nyh, &nxh, &nth, &dzh, &dxh, &dyh, &fzh, &fxh, &fyh, &sclh, &ntrh);
if (mode != 2) {
getFileInfo3D(fsnap, &nzs, &nxs, &nys, &nts, &dzs, &dxs, &dys, &fzs, &fxs, &fys, &scls, &ntrs);
@@ -122,6 +123,7 @@ int main (int argc, char **argv)
if (verbose) {
vmess("Number of virtual receivers: %li, nz: %li, nx: %li, ny: %li",nxyz,nzh,nxh,nyh);
vmess("Sampling distance in direction of z: %.3f, x: %.3f, y: %.3f",dzh,dxh,dyh);
+ vmess("Starting distance in direction of z: %.3f, x: %.3f, y: %.3f",fzh,fxh,fyh);
vmess("Number of time samples: %li",nth);
}
@@ -159,7 +161,7 @@ int main (int argc, char **argv)
fclose(fp_snap);
hdrs_mute = (segy *) calloc(nzh*nyh,sizeof(segy));
timeval = (float *)calloc(nxyz,sizeof(float));
- readSnapData3D(fray, timeval, hdrs_mute, nzh, 1, nyh, nxh, 0, 1, 0, nyh, 0, nxh);
+ readSnapData3D(fray, timeval, hdrs_mute, 1, nzh, nyh, nxh, 0, nzh, 0, nyh, 0, nxh);
dt = hdr_hom[0].dt/1E6;
}
diff --git a/utils/pwshift.c b/utils/pwshift.c
index 74ccaa1ba8f5590eeed6059a59862436de08251a..e72dd4a17ae2b59f30e800c9c092235ca8bf88f6 100644
--- a/utils/pwshift.c
+++ b/utils/pwshift.c
@@ -28,7 +28,7 @@ long writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2);
long readSnapData3D(char *filename, float *data, segy *hdrs, long nsnaps,
long nx, long ny, long nz, long sx, long ex, long sy, long ey, long sz, long ez);
-void timeShift(float *data, long nsam, long nrec, float dt, float *time, float *amp, float *delay, float fmin, float fmax);
+void timeShift(float *data, long nsam, long nrec, float dt, float *time, float *amp, float delay, float fmin, float fmax);
void pad_data(float *data, long nsam, long nrec, long nsamout, float *datout);
void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout);
@@ -64,9 +64,9 @@ int main (int argc, char **argv)
FILE *fp_gmin, *fp_ray, *fp_amp, *fp_out;
char *file_gmin, *file_ray, *file_amp, *file_out, fbr[100], fer[100], fba[100], fea[100], fins[100], fin2[100], numb1[100], *ptr;
float *gmin, *conv, *time, *amp, *image, fmin, fmax;
- float dt, dy, dx, dz, t0, y0, x0, scl, *shift, *block;
+ float dt, dy, dx, dz, t0, y0, x0, te, ye, xe, scl, *shift, *block, sy, sx;
float dt_ray, dy_ray, dx_ray, t0_ray, y0_ray, x0_ray, scl_ray, px, py, src_velox, src_veloy, src_anglex, src_angley, grad2rad;
- long nshots, nt, ny, nx, ntr, delay;
+ long nshots, nt, ny, nx, ntr, *delay, nx1, ny1, tmp;
long nray, nt_ray, ny_ray, nx_ray, ntr_ray;
long verbose, ix, iy, it, iz, is, *gx, *gy, *gz, numb, dnumb, pos, nzs, file_det;
size_t ret;
@@ -146,11 +146,16 @@ int main (int argc, char **argv)
*----------------------------------------------------------------------------*/
getFileInfo3D(file_gmin, &nt, &nx, &ny, &nshots, &dt, &dx, &dy, &t0, &x0, &y0, &scl, &ntr);
+ xe = x0 + ((float)(nx-1))*dx;
+ ye = y0 + ((float)(ny-1))*dy;
+ te = t0 + ((float)(nt-1))*dt;
+
if (verbose) {
vmess("************************ Gmin info ************************");
vmess("Number of depth levels : %li",nshots);
vmess("Number of samples x: %li, y: %li, t: %li",nx,ny,nt);
vmess("Starting distance for x: %.3f, y: %.3f, t: %.3f",x0,y0,t0);
+ vmess("Ending distance for x: %.3f, y: %.3f, t: %.3f",xe,ye,te);
vmess("Sampling distance for x: %.3f, y: %.3f, t: %.3f",dx,dy,dt);
vmess("***********************************************************");
}
@@ -186,6 +191,7 @@ int main (int argc, char **argv)
gz = (long *)calloc(nray*nshots,sizeof(long));
block = (float *)calloc(nray*nshots*nt,sizeof(float));
+ delay = (long *)calloc(nray,sizeof(long));
readSnapData3D(file_gmin, gmin, hdr_gmin, nshots, nx, ny, nt, 0, nx, 0, ny, 0, nt);
if (verbose) vmess("Read in Gmin data");
@@ -193,41 +199,67 @@ int main (int argc, char **argv)
/*----------------------------------------------------------------------------*
* Add the delay in case the plane wave is at an angle
*----------------------------------------------------------------------------*/
+
+ hdr_time = (segy *)calloc(ny*nray,sizeof(segy));
+ time = (float *)calloc(nx*ny*nray,sizeof(float));
+
+ sprintf(fins,"z%li",numb);
+ sprintf(fin2,"%s%s%s",fbr,fins,fer);
+ readSnapData3D(fin2, time, hdr_time, nray, 1, ny, nx, 0, 1, 0, ny, 0, nx);
+
+ ny1 = 1;
+ tmp = hdr_time[0].sy;
+ for (it=1; it<nray; it++) {
+ if (tmp!=hdr_time[it*ny].sy)
+ {
+ ny1++;
+ tmp = hdr_time[it*ny].sy;
+ }
+ }
+ nx1=nray/ny1;
+ vmess("ny1=%li nx1=%li nray=%li",ny1,nx1,nray);
- shift = (float *)calloc(nx*ny,sizeof(float));
+ shift = (float *)calloc(nray,sizeof(float));
grad2rad = 17.453292e-3;
- px = sin(src_anglex*grad2rad)/src_velox;
- py = sin(src_angley*grad2rad)/src_veloy;
- if (verbose) vmess("px value is %f and py value is %f",px,py);
-
- // if (py < 0.0) {
- // for (iy=0; iy<ny; iy++) {
- // if (px < 0.0) {
- // for (ix=0; ix<nx; ix++) {
- // shift[iy*nx+ix] = fabsf((nx-1-ix)*dx*px) + fabsf((ny-1-iy)*dy*py);
- // }
- // }
- // else {
- // for (ix=0; ix<nx; ix++) {
- // shift[iy*nx+ix] = ix*dx*px + fabsf((ny-1-iy)*dy*py);
- // }
- // }
- // }
- // }
- // else {
- // for (iy=0; iy<ny; iy++) {
- // if (px < 0.0) {
- // for (ix=0; ix<nx; ix++) {
- // shift[iy*nx+ix] = fabsf((nx-1-ix)*dx*px) + iy*dy*py;
- // }
- // }
- // else {
- // for (ix=0; ix<nx; ix++) {
- // shift[iy*nx+ix] = ix*dx*px + iy*dy*py;
- // }
- // }
- // }
- // }
+ if (src_anglex==0.0) px=0.0;
+ else px = sin(-1.0*src_anglex*grad2rad)/src_velox;
+ if (src_angley==0.0) py=0.0;
+ else py = sin(-1.0*src_angley*grad2rad)/src_veloy;
+ if (verbose) {
+ vmess("************************ plane wave info *************************");
+ vmess("Plane wave angle x: %.3f, y: %.3f",src_anglex,src_angley);
+ vmess("Plane wave velocity x: %.3f, y: %.3f",src_velox,src_veloy);
+ vmess("Plane wave slowness x: %.3f, y: %.3f",px,py);
+ vmess("***********************************************************");
+ }
+
+ for (iy=0; iy<ny1; iy++) {
+ for (ix=0; ix<nx1; ix++) {
+ sy = ((float)(hdr_time[(iy*nx1+ix)*ny].sy))/1000.0;
+ sx = ((float)(hdr_time[(iy*nx1+ix)*ny].sx))/1000.0;
+ //vmess("sy=%.3f sx=%.3f",sy,sx);
+ if (py < 0.0) {
+ if (px < 0.0) {
+ shift[iy*nx1+ix] = fabsf((xe-sx)*px) + fabsf((ye-sy)*py);
+ }
+ else {
+ shift[iy*nx1+ix] = fabsf((x0-sx)*px) + fabsf((ye-sy)*py);
+ }
+ }
+ else {
+ if (px < 0.0) {
+ shift[iy*nx1+ix] = fabsf((xe-sx)*px) + fabsf((y0-sy)*py);
+ }
+ else {
+ shift[iy*nx1+ix] = fabsf((x0-sx)*px) + fabsf((y0-sy)*py);
+ }
+ }
+ }
+ }
+
+
+ free(hdr_time); free(time);
+
/*----------------------------------------------------------------------------*
* Apply the imaging condition
@@ -247,7 +279,7 @@ int main (int argc, char **argv)
sprintf(fin2,"%s%s%s",fbr,fins,fer);
readSnapData3D(fin2, time, hdr_time, nray, 1, ny, nx, 0, 1, 0, ny, 0, nx);
sprintf(fin2,"%s%s%s",fba,fins,fea);
- readSnapData3D(file_amp, amp, hdr_amp, nray, 1, ny, nx, 0, 1, 0, ny, 0, nx);
+ readSnapData3D(file_amp, amp, hdr_amp, nray, 1, ny, nx, 0, 1, 0, ny, 0, nx);
for (is = 0; is < nray; is++) {
gx[is*nshots+iy] = hdr_time[is*ny].sx;
@@ -257,34 +289,35 @@ int main (int argc, char **argv)
x0_ray = ((float)gx[is*nshots+iy])/1000.0;
y0_ray = ((float)gy[is*nshots+iy])/1000.0;
- if (py < 0.0) {
- if (px < 0.0) {
- delay = NINT(fabsf((x0-x0_ray)*px)/dt) + NINT(fabsf((y0-y0_ray)*py)/dt);
- }
- else {
- delay = NINT((x0_ray-x0)*px/dt) + NINT(fabsf((y0-y0_ray)*py)/dt);
- }
- }
- else {
- if (px < 0.0) {
- delay = NINT(fabsf((x0-x0_ray)*px)/dt) + NINT((y0_ray-y0)*py/dt);
- }
- else {
- delay = NINT((x0_ray-x0)*px/dt) + NINT((y0_ray-y0)*py/dt);
- }
- }
-
- if (delay > nt-1) delay = nt-1;
+ // if (py < 0.0) {
+ // if (px < 0.0) {
+ // delay = NINT(fabsf((x0-x0_ray)*px)/dt) + NINT(fabsf((y0-y0_ray)*py)/dt);
+ // }
+ // else {
+ // delay = NINT((x0_ray-x0)*px/dt) + NINT(fabsf((y0-y0_ray)*py)/dt);
+ // }
+ // }
+ // else {
+ // if (px < 0.0) {
+ // delay = NINT(fabsf((x0-x0_ray)*px)/dt) + NINT((y0_ray-y0)*py/dt);
+ // }
+ // else {
+ // delay = NINT((x0_ray-x0)*px/dt) + NINT((y0_ray-y0)*py/dt);
+ // }
+ // }
+
+ // if (delay > nt-1) delay = nt-1;
for (it = 0; it < nx*ny*nt; it++) {
conv[it] = gmin[iy*nx*ny*nt+it];
}
- timeShift(&conv[0],nt,nx*ny,dt,&time[is*ny*nx],&[is*ny*nx],shift,fmin,fmax);
+ timeShift(&conv[0],nt,nx*ny,dt,&time[is*ny*nx],&[is*ny*nx],shift[is],fmin,fmax);
for (ix = 0; ix < ny*nx; ix++) {
- image[is*nshots+iy] += conv[ix*nt+delay+2]*dx*dy*dt;
+ image[is*nshots+iy] += conv[ix*nt]*dx*dy*dt;
for (it=0; it<nt; it++) {
block[is*nshots*nt+iy*nt+it] += conv[ix*nt+it];
}
+ block[is*nshots*nt+iy*nt] = 1e20;
}
}
@@ -324,37 +357,37 @@ int main (int argc, char **argv)
fclose(fp_out);
- fp_out = fopen("block.su", "w+");
-
- if (nshots>1) dz = ((float)(gz[1] - gz[0]))/1000.0;
- else dz = 1.0;
-
- for (is = 0; is < nshots; is++) {
- for (it = 0; it < nray; it++) {
- hdr_out[it*nshots+is].fldr = it+1;
- hdr_out[it*nshots+is].tracl = it+1;
- hdr_out[it*nshots+is].tracf = it+1;
- hdr_out[it*nshots+is].scalco = -1000;
- hdr_out[it*nshots+is].scalel = -1000;
- hdr_out[it*nshots+is].trid = 1;
- hdr_out[it*nshots+is].ns = nt;
- hdr_out[it*nshots+is].trwf = nray;
- hdr_out[it*nshots+is].ntr = nray;
- hdr_out[it*nshots+is].f1 = (((float)gz[0])/1000.0);
- hdr_out[it*nshots+is].f2 = (((float)gx[0])/1000.0);
- hdr_out[it*nshots+is].dt = ((int)(dt*1E6));
- hdr_out[it*nshots+is].d1 = roundf(dz*1000.0)/1000.0;
- hdr_out[it*nshots+is].d2 = roundf(dx*1000.0)/1000.0;
- hdr_out[it*nshots+is].gx = gx[it*nshots+is];
- hdr_out[it*nshots+is].gy = gy[it*nshots+is];
- hdr_out[it*nshots+is].sdepth = gz[it*nshots+is];
- }
- }
+ // fp_out = fopen("block.su", "w+");
+
+ // if (nshots>1) dz = ((float)(gz[1] - gz[0]))/1000.0;
+ // else dz = 1.0;
+
+ // for (is = 0; is < nshots; is++) {
+ // for (it = 0; it < nray; it++) {
+ // hdr_out[it*nshots+is].fldr = it+1;
+ // hdr_out[it*nshots+is].tracl = it+1;
+ // hdr_out[it*nshots+is].tracf = it+1;
+ // hdr_out[it*nshots+is].scalco = -1000;
+ // hdr_out[it*nshots+is].scalel = -1000;
+ // hdr_out[it*nshots+is].trid = 1;
+ // hdr_out[it*nshots+is].ns = nt;
+ // hdr_out[it*nshots+is].trwf = nray;
+ // hdr_out[it*nshots+is].ntr = nray;
+ // hdr_out[it*nshots+is].f1 = (((float)gz[0])/1000.0);
+ // hdr_out[it*nshots+is].f2 = (((float)gx[0])/1000.0);
+ // hdr_out[it*nshots+is].dt = ((int)(dt*1E6));
+ // hdr_out[it*nshots+is].d1 = dt;
+ // hdr_out[it*nshots+is].d2 = roundf(dx*1000.0)/1000.0;
+ // hdr_out[it*nshots+is].gx = gx[it*nshots+is];
+ // hdr_out[it*nshots+is].gy = gy[it*nshots+is];
+ // hdr_out[it*nshots+is].sdepth = gz[it*nshots+is];
+ // }
+ // }
- ret = writeData3D(fp_out, &block[0], hdr_out, nt, nray*nshots);
- if (ret < 0 ) verr("error on writing output file.");
+ // ret = writeData3D(fp_out, &block[0], hdr_out, nt, nray*nshots);
+ // if (ret < 0 ) verr("error on writing output file.");
- fclose(fp_out);
+ // fclose(fp_out);
free(image); free(hdr_out); free(block);
@@ -363,7 +396,7 @@ int main (int argc, char **argv)
return 0;
}
-void timeShift(float *data, long nsam, long nrec, float dt, float *time, float *amp, float *delay, float fmin, float fmax)
+void timeShift(float *data, long nsam, long nrec, float dt, float *time, float *amp, float delay, float fmin, float fmax)
{
long optn, iom, iomin, iomax, nfreq, ix, sign;
float omin, omax, deltom, om, tom, df, *rdata, scl;
@@ -406,7 +439,7 @@ void timeShift(float *data, long nsam, long nrec, float dt, float *time, float *
}
for (iom = iomin ; iom < iomax ; iom++) {
om = deltom*iom;
- tom = om*-1.0*(time[ix]-delay[ix]);
+ tom = om*-1.0*(time[ix]+delay);
cdatascl[ix*nfreq+iom].r = (cdata[ix*nfreq+iom].r*cos(-tom) - cdata[ix*nfreq+iom].i*sin(-tom))/(amp[ix]*amp[ix]);
cdatascl[ix*nfreq+iom].i = (cdata[ix*nfreq+iom].i*cos(-tom) + cdata[ix*nfreq+iom].r*sin(-tom))/(amp[ix]*amp[ix]);
}
@@ -443,4 +476,4 @@ void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long
for (it = 0 ; it < nsamout ; it++)
datout[ix*nsamout+it] = scl*data[ix*nsam+it];
}
-}
\ No newline at end of file
+}