Skip to content
Snippets Groups Projects
Commit 5643c815 authored by JanThorbecke's avatar JanThorbecke
Browse files

preparing for GPY MME paper publication

parent 6bc2ae4d
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
marchenko/demo/mme/clean 0 → 100755
+ 4
− 0
View file @ 5643c815
#!/bin/bash
rm *.su *.bin *.eps nep line* *.asci
marchenko/demo/mme/epsModel.scr
+ 4
− 4
View file @ 5643c815
...@@ -21,13 +21,13 @@ cat << EOF3 > line3 ...@@ -21,13 +21,13 @@ cat << EOF3 > line3
EOF3 EOF3
#model #model
supsimage hbox=4 wbox=6 labelsize=12 < model10_cp.su \ supsimage hbox=4 wbox=6 labelsize=16 < model10_cp.su \
x1beg=0 x1end=1400.0 d1num=200 lstyle=vertright legend=1 threecolor=0 \ x1beg=0 x1end=1400.0 d1num=200 lstyle=vertright legend=1 threecolor=0 \
label1="depth (m)" label2="lateral distance (m)" \ label1="depth (m)" label2="lateral distance (m)" \
curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black \ curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black \
n1tic=5 x2beg=-2250 f2num=-2000 d2num=1000 x2end=2250 > model_cp_line.eps n1tic=5 x2beg=-2250 f2num=-2000 d2num=1000 x2end=2250 > model_cp_line.eps
supsimage hbox=4 wbox=6 labelsize=12 < model10_ro.su \ supsimage hbox=4 wbox=6 labelsize=16 < model10_ro.su \
x1beg=0 x1end=1400.0 d1num=200 lstyle=vertright legend=1 threecolor=0 \ x1beg=0 x1end=1400.0 d1num=200 lstyle=vertright legend=1 threecolor=0 \
label1="depth (m)" label2="lateral distance (m)" \ label1="depth (m)" label2="lateral distance (m)" \
curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black \ curve=line1,line2,line3 npair=2,2,2 curvecolor=black,black,black \
...@@ -36,12 +36,12 @@ supsimage hbox=4 wbox=6 labelsize=12 < model10_ro.su \ ...@@ -36,12 +36,12 @@ supsimage hbox=4 wbox=6 labelsize=12 < model10_ro.su \
#wavelet #wavelet
dt=0.0005 dt=0.0005
supsgraph < wavefw.su \ supsgraph < wavefw.su \
labelsize=12 d1=$dt style=normal \ labelsize=16 d1=$dt style=normal \
label1="time (s)" label2="amplitude" \ label1="time (s)" label2="amplitude" \
d1num=0.15 wbox=6 hbox=3 x1end=0.9 > wavefw.eps d1num=0.15 wbox=6 hbox=3 x1end=0.9 > wavefw.eps
sufft < wavefw.su | suamp | sugain scale=$dt | supsgraph \ sufft < wavefw.su | suamp | sugain scale=$dt | supsgraph \
labelsize=12 style=normal \ labelsize=16 style=normal \
label1="frequency (1/s)" label2="amplitude" \ label1="frequency (1/s)" label2="amplitude" \
d1num=10 wbox=6 hbox=3 x1end=125 x2end=1.1 > wavefw_freq.eps d1num=10 wbox=6 hbox=3 x1end=125 x2end=1.1 > wavefw_freq.eps
......
marchenko/demo/mme/epsPrimaries.scr
+ 12
− 8
View file @ 5643c815
...@@ -102,6 +102,10 @@ fi ...@@ -102,6 +102,10 @@ fi
# second iteration # second iteration
if [[ "$1" == "Figure4" ]]; if [[ "$1" == "Figure4" ]];
then then
if [[ ! -f "shotsx0.su" ]]; then
echo "ERR: File shotsx0.su is not yet created, please run ./epsPrimaries.scr Figure3 first."
exit
fi
suwind itmax=1023 < Mi_276001.su > N0.su suwind itmax=1023 < Mi_276001.su > N0.su
#compute R*N0 #compute R*N0
fconv file_in1=shotsx0.su file_in2=N0.su file_out=fconvN1.su verbose=1 fmax=90 fconv file_in1=shotsx0.su file_in2=N0.su file_out=fconvN1.su verbose=1 fmax=90
...@@ -182,7 +186,7 @@ echo $piter ...@@ -182,7 +186,7 @@ echo $piter
file=Mi_${istart}${piter}.su file=Mi_${istart}${piter}.su
file_base=${file%.su} file_base=${file%.su}
#possibly add suflip flip=3 to flip the time-axis #possibly add suflip flip=3 to flip the time-axis
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \ supsimage < $file hbox=6 wbox=4 labelsize=16 linewidth=0.0 \
n1tic=2 d2=5 x1beg=$nstart x1end=$ns d1=1 d1num=100 \ n1tic=2 d2=5 x1beg=$nstart x1end=$ns d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \ label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
...@@ -197,7 +201,7 @@ piter=$(printf %03d $iter) ...@@ -197,7 +201,7 @@ piter=$(printf %03d $iter)
echo $piter echo $piter
file=Mi_${istart}${piter}.su file=Mi_${istart}${piter}.su
file_base=${file%.su} file_base=${file%.su}
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \ supsimage < $file hbox=6 wbox=4 labelsize=16 linewidth=0.0 \
n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \ n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \ label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > $file_base.eps
...@@ -218,14 +222,14 @@ file=Mi_${istart}${piter}.su ...@@ -218,14 +222,14 @@ file=Mi_${istart}${piter}.su
#shift=$(echo "scale=4; ($dtrcv*($ns/2.0-1))" | bc -l) #shift=$(echo "scale=4; ($dtrcv*($ns/2.0-1))" | bc -l)
#basop choice=shift shift=$shift file_in=$file | \ #basop choice=shift shift=$shift file_in=$file | \
file_base=${file%.su} file_base=${file%.su}
suflip flip=3 < $file | supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 \ suflip flip=3 < $file | supsimage hbox=6 wbox=4 labelsize=16 linewidth=0.0 \
n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \ n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \ label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > ${file_base}flip.eps f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > ${file_base}flip.eps
file=k1min_${istart}${piter}.su file=k1min_${istart}${piter}.su
file_base=${file%.su} file_base=${file%.su}
supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 < $file\ supsimage hbox=6 wbox=4 labelsize=16 linewidth=0.0 < $file\
n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \ n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \ label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > $file_base.eps f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > $file_base.eps
...@@ -242,7 +246,7 @@ for (( istart=246; istart<=316; istart+=10 )) ...@@ -242,7 +246,7 @@ for (( istart=246; istart<=316; istart+=10 ))
do do
file=k1min_${istart}${piter}.su file=k1min_${istart}${piter}.su
file_base=${file%.su} file_base=${file%.su}
supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 < $file\ supsimage hbox=6 wbox=4 labelsize=16 linewidth=0.0 < $file\
n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \ n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \ label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > $file_base.eps f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > $file_base.eps
...@@ -259,7 +263,7 @@ echo $piter ...@@ -259,7 +263,7 @@ echo $piter
file=Mi_${istart}${piter}.su file=Mi_${istart}${piter}.su
file_base=${file%.su} file_base=${file%.su}
#possibly add suflip flip=3 to flip the time-axis #possibly add suflip flip=3 to flip the time-axis
supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \ supsimage < $file hbox=6 wbox=4 labelsize=16 linewidth=0.0 \
n1tic=2 d2=5 x1beg=$nstart x1end=$ns d1=1 d1num=100 \ n1tic=2 d2=5 x1beg=$nstart x1end=$ns d1=1 d1num=100 \
label1="time sample number" label2="lateral distance" \ label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}T.eps f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}T.eps
...@@ -296,12 +300,12 @@ file=WindowOdd276.su ...@@ -296,12 +300,12 @@ file=WindowOdd276.su
file_base=${file%.su} file_base=${file%.su}
suwind key=tracl min=1 max=1 < $file | suwind itmin=1024 | \ suwind key=tracl min=1 max=1 < $file | suwind itmin=1024 | \
supsgraph d1=1 style=normal f1=0 \ supsgraph d1=1 style=normal f1=0 \
labelsize=12 label1="time sample number" label2="amplitude" \ labelsize=14 label1="time sample number" label2="amplitude" \
d1num=100 wbox=6 hbox=3 x2end=1.2 > $file_base.eps d1num=100 wbox=6 hbox=3 x2end=1.2 > $file_base.eps
file=WindowEven276.su file=WindowEven276.su
file_base=${file%.su} file_base=${file%.su}
suwind key=tracl min=1 max=1 < $file | suwind itmax=1024 | \ suwind key=tracl min=1 max=1 < $file | suwind itmax=1024 | \
supsgraph d1=1 style=normal f1=0 \ supsgraph d1=1 style=normal f1=0 \
labelsize=12 label1="time sample number" label2="amplitude" \ labelsize=14 label1="time sample number" label2="amplitude" \
d1num=100 wbox=6 hbox=3 x2end=1.2 > $file_base.eps d1num=100 wbox=6 hbox=3 x2end=1.2 > $file_base.eps
marchenko/marchenko_primaries.c
+ 63
− 15
View file @ 5643c815
...@@ -109,6 +109,13 @@ char *sdoc[] = { ...@@ -109,6 +109,13 @@ char *sdoc[] = {
NULL}; NULL};
/**************** end self doc ***********************************/ /**************** end self doc ***********************************/
/*
* The equation and Figure numbers mentioned in the comments refer to the paper:
* "Implementation of the Marchenko Multiple Elimination algorithm",
* Thorbecke, Jan; Zhang, Lele; Wapenaar, Kees; Slob, Evert,
* Accepted for publication, todo add full reference after publication
*/
int main (int argc, char **argv) int main (int argc, char **argv)
{ {
FILE *fp_out, *fp_rr, *fp_w, *fp_ud; FILE *fp_out, *fp_rr, *fp_w, *fp_ud;
...@@ -197,7 +204,7 @@ int main (int argc, char **argv) ...@@ -197,7 +204,7 @@ int main (int argc, char **argv)
if (reci && ntap) vwarn("tapering influences the reciprocal result"); if (reci && ntap) vwarn("tapering influences the reciprocal result");
/* defines the smooth transition zone of the time-window */ /* defines the smooth transition zone of the time-window (Figure 5) */
smooth = MIN(smooth, shift); smooth = MIN(smooth, shift);
if (smooth) { if (smooth) {
...@@ -431,7 +438,7 @@ int main (int argc, char **argv) ...@@ -431,7 +438,7 @@ int main (int argc, char **argv)
} }
/*================ Defining focusing operator(s) from R ================*/ /*================ Defining focusing operator(s) from R ================*/
/* M0 = -R(ishot,-t) */ /* M0 = -R(ishot,-t) equation (3) */
/* use ishot from Refl, complex-conjugate(time reverse), scale with -1 and convolve with wavelet */ /* use ishot from Refl, complex-conjugate(time reverse), scale with -1 and convolve with wavelet */
if (file_tinv == NULL) { if (file_tinv == NULL) {
...@@ -702,7 +709,7 @@ int main (int argc, char **argv) ...@@ -702,7 +709,7 @@ int main (int argc, char **argv)
for (j = 0; j < nts; j++) { for (j = 0; j < nts; j++) {
M0[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+j]; M0[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+j];
} }
/* apply mute window for samples above nts-ii */ /* apply mute window for samples above nts-ii : the Heaviside function in equation (3) */
for (j = 0; j < MIN(nts, nts-iw+isms); j++) { for (j = 0; j < MIN(nts, nts-iw+isms); j++) {
M0[l*nxs*nts+i*nts+j] = 0.0; M0[l*nxs*nts+i*nts+j] = 0.0;
} }
...@@ -710,6 +717,7 @@ int main (int argc, char **argv) ...@@ -710,6 +717,7 @@ int main (int argc, char **argv)
M0[l*nxs*nts+i*nts+j] *= costaper[smooth-k]; M0[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
} }
} }
/* Initilisation of k1min: the first term (R) in the right hand side of equation(6) */
for (i = 0; i < npos; i++) { for (i = 0; i < npos; i++) {
ix = ixpos[i]; ix = ixpos[i];
j = 0; j = 0;
...@@ -743,6 +751,7 @@ int main (int argc, char **argv) ...@@ -743,6 +751,7 @@ int main (int argc, char **argv)
} }
} }
} }
/* If enabled this write the result of equation(3) to disk */
if (file_iter != NULL) { if (file_iter != NULL) {
writeDataIter("M0.su", M0, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii); writeDataIter("M0.su", M0, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii);
} }
...@@ -760,15 +769,17 @@ int main (int argc, char **argv) ...@@ -760,15 +769,17 @@ int main (int argc, char **argv)
t2 = wallclock_time(); t2 = wallclock_time();
/*================ construction of Mi(-t) = - \int R(x,t) Mi(t) ================*/ /*================ construction of Mi(-t) = - \int R(x,t) Mi(t) ================*/
/* synthesis process is the compute kernel in equations (4) and (5) in the Geophysics implementation paper */ /* synthesis process is the compute kernel in equations (4) and (5) (rewritten from equation (2)) */
synthesisp(Refl, Fop, Mi, RMi, nx, nt, nacq, nts, dt, xsyn, Nfoc, synthesisp(Refl, Fop, Mi, RMi, nx, nt, nacq, nts, dt, xsyn, Nfoc,
xrcv, xsrc, xnx, fxse, fxsb, dxs, dxsrc, dx, ntfft, nw, nw_low, nw_high, mode, xrcv, xsrc, xnx, fxse, fxsb, dxs, dxsrc, dx, ntfft, nw, nw_low, nw_high, mode,
reci, nshots, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv, ixmask, verbose); reci, nshots, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv, ixmask, verbose);
/* If enbled the result of equation(4) or (5) (depending on odd/even) is written to disk */
if (file_iter != NULL) { if (file_iter != NULL) {
writeDataIter(file_iter, RMi, hdrs_out, ntfft, nacq, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter+1); writeDataIter(file_iter, RMi, hdrs_out, ntfft, nacq, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter+1);
} }
/* The Energy is computed to check the convergence of the scheme see Figure 6c and Figure 8 */
if (verbose >=2) { if (verbose >=2) {
for (l = 0; l < Nfoc; l++) { for (l = 0; l < Nfoc; l++) {
if (iter % 2 == 0) { if (iter % 2 == 0) {
...@@ -790,7 +801,8 @@ int main (int argc, char **argv) ...@@ -790,7 +801,8 @@ int main (int argc, char **argv)
t3 = wallclock_time(); t3 = wallclock_time();
tsyn += t3 - t2; tsyn += t3 - t2;
/* N_k(x,t) = -N_(k-1)(x,-t) */ /* M_i(x,t) = RMi_(i-1)(x,-t) : time-reversal of the outcome of the synthesis process */
/* see equation(3) where the time reversal of R is taken as input */
for (l = 0; l < Nfoc; l++) { for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) { for (i = 0; i < npos; i++) {
j = 0; j = 0;
...@@ -808,18 +820,19 @@ int main (int argc, char **argv) ...@@ -808,18 +820,19 @@ int main (int argc, char **argv)
for (i = 0; i < npos; i++) { for (i = 0; i < npos; i++) {
ix = ixpos[i]; ix = ixpos[i];
iw = NINT((ii*dt+twplane[ix])/dt); iw = NINT((ii*dt+twplane[ix])/dt);
/* store results in kplus */ /* store results in kplus the 'plus-part' of equation (6) */
for (j = 0; j < nts; j++) { for (j = 0; j < nts; j++) {
k1plus[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+j]; k1plus[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+j];
} }
/* apply mute window for samples after ii */ /* Apply mute window for samples after ii : the Heaviside function in equation (5) */
/* This also defines v1_plus in equation (8) */
for (j = MAX(0,iw-isme); j < nts; j++) { for (j = MAX(0,iw-isme); j < nts; j++) {
Mi[l*nxs*nts+i*nts+j] = 0.0; Mi[l*nxs*nts+i*nts+j] = 0.0;
} }
for (j = MAX(0,iw-isme), k=0; j < iw-isms; j++, k++) { for (j = MAX(0,iw-isme), k=0; j < iw-isms; j++, k++) {
Mi[l*nxs*nts+i*nts+j] *= costaper[k]; Mi[l*nxs*nts+i*nts+j] *= costaper[k];
} }
/* apply mute window for delta function at t=0*/ /* Apply mute window for delta function at t=0*/
iw = NINT((twplane[ix])/dt); iw = NINT((twplane[ix])/dt);
for (j = 0; j < MAX(0,iw+shift-smooth); j++) { for (j = 0; j < MAX(0,iw+shift-smooth); j++) {
Mi[l*nxs*nts+i*nts+j] = 0.0; Mi[l*nxs*nts+i*nts+j] = 0.0;
...@@ -827,6 +840,7 @@ int main (int argc, char **argv) ...@@ -827,6 +840,7 @@ int main (int argc, char **argv)
for (j = MAX(0,iw+shift-smooth), k=1; j < MAX(0,iw+shift); j++, k++) { for (j = MAX(0,iw+shift-smooth), k=1; j < MAX(0,iw+shift); j++, k++) {
Mi[l*nxs*nts+i*nts+j] *= costaper[smooth-k]; Mi[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
} }
/* This defines v1_plus in equation (8) */
for (j = 0; j < nts; j++) { for (j = 0; j < nts; j++) {
v1plus[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+j]; v1plus[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+j];
} }
...@@ -844,11 +858,13 @@ int main (int argc, char **argv) ...@@ -844,11 +858,13 @@ int main (int argc, char **argv)
for (j = 0; j < nts; j++) { for (j = 0; j < nts; j++) {
Mi[l*nxs*nts+i*nts+j] += -DD[l*nxs*nts+ix*nts+j]; Mi[l*nxs*nts+i*nts+j] += -DD[l*nxs*nts+ix*nts+j];
} }
/* store results in kmin as defined in equation (6) */
j = 0; j = 0;
k1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+j]; k1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) { for (j = 1; j < nts; j++) {
k1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+nts-j]; k1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+nts-j];
} }
/* This collects the eliminated internal multiples and is the sum in right hand-side in equation (1) */
if (file_update != NULL) { if (file_update != NULL) {
j=0; j=0;
Mup[l*nxs*nts+i*nts+j] += DD[l*nxs*nts+ix*nts+j] - k1min[l*nxs*nts+i*nts+j]; Mup[l*nxs*nts+i*nts+j] += DD[l*nxs*nts+ix*nts+j] - k1min[l*nxs*nts+i*nts+j];
...@@ -858,11 +874,13 @@ int main (int argc, char **argv) ...@@ -858,11 +874,13 @@ int main (int argc, char **argv)
} }
} }
else { else {
/* store results in kmin as defined in equation (6) */
j = 0; j = 0;
k1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+j]; k1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) { for (j = 1; j < nts; j++) {
k1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+nts-j]; k1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+nts-j];
} }
/* Mup collects the eliminated internal multiples and is the sum in right hand-side in equation (1) */
if (file_update != NULL) { if (file_update != NULL) {
j=0; j=0;
Mup[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+j]; Mup[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+j];
...@@ -871,7 +889,8 @@ int main (int argc, char **argv) ...@@ -871,7 +889,8 @@ int main (int argc, char **argv)
} }
} }
} }
/* apply mute window for samples above nts-ii */ /* Apply mute window for samples above nts-ii : the Heaviside function in equation (4)
* This also defines u1_min, for t>=t2-epsilon in equation (7) */
iw = NINT((ii*dt+twplane[ix])/dt); iw = NINT((ii*dt+twplane[ix])/dt);
for (j = 0; j < MIN(nts,nts-iw+isms); j++) { for (j = 0; j < MIN(nts,nts-iw+isms); j++) {
Mi[l*nxs*nts+i*nts+j] = 0.0; Mi[l*nxs*nts+i*nts+j] = 0.0;
...@@ -907,7 +926,8 @@ int main (int argc, char **argv) ...@@ -907,7 +926,8 @@ int main (int argc, char **argv)
} }
} }
/* write optional u- u+ v- v+ to disk */ /* Write optional u- u+ v- v+ to disk */
/* Based on equation(6) and (7) k1min is split in vmin and umin */
if (file_vmin != NULL) { if (file_vmin != NULL) {
for (l = 0; l < Nfoc; l++) { for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) { for (i = 0; i < npos; i++) {
...@@ -952,16 +972,44 @@ int main (int argc, char **argv) ...@@ -952,16 +972,44 @@ int main (int argc, char **argv)
} }
writeDataIter(file_umin, uv, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, ii); writeDataIter(file_umin, uv, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, ii);
} }
/*
* \begin{alignat}{2}
* \label{k1p}
* k_{1,i}^+({\bf x}'_0,{\bf x}''_0,t,t_2) = & \int_{t'=0}^{+\infty} \int_{\partial \setD_0} R({\bf x}'_0,{\bf x}_0,-t')
* v_{1,i}^-({\bf x}_0,{\bf x}''_0,t-t',t_2) {\rm d}t' {\rm d} {\bf x}_0.
* \end{alignat}
*
* Equation (\ref{k1p}) can be further split in the time domain as follows;
* %
* \begin{equation} \label{uv+}
* k_{1,i}^+({\bf x}'''_0,{\bf x}''_0,t,t_2) = \begin{cases}
* v_{1,i}^+({\bf x}'''_0,{\bf x}''_0,t,t_2) & t < t_2-\varepsilon \\
* u_{1,i}^+({\bf x}'''_0,{\bf x}''_0,t,t_2) & t \geq t_2-\varepsilon
* \end{cases},
* \end{equation}
* %
*
*/
/* The equation (\ref{k1p}) above is the k1plus equivalent of k1min of equation(6)
* and the time window in equation (\ref{uv+}) defines vplus and uplus */
/* Based on these equations and k1plus is split in vplus and uplus */
if (file_vplus != NULL) { if (file_vplus != NULL) {
for (l = 0; l < Nfoc; l++) { for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) { for (i = 0; i < npos; i++) {
/* apply mute window for delta function at t=0*/ /* apply mute window for delta function at t=0*/
iw0 = NINT((twplane[ix])/dt); iw0 = NINT((twplane[ix])/dt);
for (j = 0; j < MAX(0,iw0+shift-smooth); j++) { /* to apply muting window around epsilon below t=0 uncomment this part */
uv[l*nxs*nts+i*nts+j] = 0.0; /* This (partly) removes the delta function at t=0 */
} // for (j = 0; j < MAX(0,iw0+shift-smooth); j++) {
for (j = MAX(0,iw0+shift-smooth), k=1; j < MAX(0,iw0+shift); j++, k++) { // uv[l*nxs*nts+i*nts+j] = 0.0;
uv[l*nxs*nts+i*nts+j] = k1plus[l*nxs*nts+i*nts+j]*costaper[smooth-k]; // }
// for (j = MAX(0,iw0+shift-smooth), k=1; j < MAX(0,iw0+shift); j++, k++) {
// uv[l*nxs*nts+i*nts+j] = k1plus[l*nxs*nts+i*nts+j]*costaper[smooth-k];
// }
/* this keeps the delta pulse below at at t=0 */
for (j = 0; j < iw0+shift; j++) {
uv[l*nxs*nts+i*nts+j] = k1plus[l*nxs*nts+i*nts+j];
} }
iw = NINT((ii*dt+twplane[ix])/dt); iw = NINT((ii*dt+twplane[ix])/dt);
for (j = MAX(0,iw0+shift); j < iw-isme; j++, k++) { for (j = MAX(0,iw0+shift); j < iw-isme; j++, k++) {
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment