Merge branch 'master' of https://github.com/JanThorbecke/OpenSource

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++) {

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++) {

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.

marchenko3D/demo/marchenko3D/oneD/execute.sh

 #!/bin/bash
 
-./model.sh
+#./model.sh
 
-./shotmodel.sh
+#./shotmodel.sh
 
-./makeR.sh
+#./makeR.sh
+
+./zfp.sh
 
 ./farrmod.sh
 
 ./fmute3D.sh
 
-./zfp.sh
-
 ./marchenko3d.sh

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 \

marchenko3D/demo/marchenko3D/oneD/fmute3D.sh

 #!/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

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
 

marchenko3D/demo/marchenko3D/oneD/marchenko3d.sh

 #!/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

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

marchenko3D/demo/marchenko3D/oneD/plane_horizontal/execute_plane.sh

+#!/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

marchenko3D/demo/marchenko3D/oneD/plane_horizontal/farrmod_plane.sh

+#!/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 \

marchenko3D/demo/marchenko3D/oneD/plane_horizontal/fmute3D_plane.sh

+#!/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

marchenko3D/demo/marchenko3D/oneD/plane_horizontal/marchenko3d_plane.sh

+#!/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

marchenko3D/demo/marchenko3D/oneD/plane_tilted/execute_tilted.sh

+#!/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

marchenko3D/demo/marchenko3D/oneD/plane_tilted/farrmod_tilted_neg.sh

+#!/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 \

marchenko3D/demo/marchenko3D/oneD/plane_tilted/farrmod_tilted_pos.sh

+#!/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 \

marchenko3D/demo/marchenko3D/oneD/plane_tilted/fmute3D_tilted_neg.sh

+#!/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

marchenko3D/demo/marchenko3D/oneD/plane_tilted/fmute3D_tilted_pos.sh

+#!/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

marchenko3D/demo/marchenko3D/oneD/plane_tilted/marchenko3d_tilted_neg.sh

+#!/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

marchenko3D/demo/marchenko3D/oneD/plane_tilted/marchenko3d_tilted_pos.sh

+#!/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