Skip to content

Commit 289151c

Browse files
JanThorbeckeJanThorbecke
committed
taper in time window in primaries app
1 parent 87c7fb8 commit 289151c

File tree

6 files changed

+165
-56
lines changed

6 files changed

+165
-56
lines changed

‎marchenko/demo/oneD/iterations.scr‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
#!/bin/bash
2+
3+
#second reflector at time:
4+
# 800/1800+600/2300
5+
# .70531400966183574878 => sample 176
6+
#thrids reflector at time:
7+
# 800/1800+600/2300+800/2000
8+
# 1.10531400966183574878 sample 276
9+
10+
select=451
11+
12+
#for istart in 176 200 276 400
13+
for istart in 176
14+
do
15+
(( iend = istart + 1 ))
16+
17+
../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
18+
nshots=601 verbose=1 istart=$istart iend=$iend fmax=90 \
19+
pad=1024 niter=45 smooth=10 niterskip=600 shift=20 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su
20+
21+
suximage < f1min_${istart}043.su clip=1 title="${istart}" &
22+
23+
done
24+

‎marchenko/demo/oneD/primaries.scr‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -13,14 +13,14 @@ export OMP_NUM_THREADS=20
1313
#shot record to remove internal multiples
1414
select=451
1515

16-
makewave fp=20 dt=0.004 file_out=wave.su nt=1024 t0=0.0 scale=0 scfft=1
16+
makewave fp=20 dt=0.004 file_out=wave.su nt=2048 t0=0.0 scale=0 scfft=1
1717

1818
../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
19-
nshots=601 verbose=1 istart=40 fmax=90 \
20-
niter=15 niterskip=600 shift=20 file_rr=pred_rr.su T=0
19+
nshots=601 verbose=2 istart=40 iend=500 fmax=90 pad=1024 \
20+
niter=25 smooth=10 niterskip=600 shift=20 file_rr=pred_rr.su T=0
2121

2222
#for reference original shot record from Reflection matrix
23-
suwind key=fldr min=$select max=$select < shotsdx5_rp.su > shotsx0.su
23+
suwind key=fldr min=$select max=$select < shotsdx5_rp.su itmax=2047 > shotsx0.su
2424
fconv file_in1=shotsx0.su file_in2=wave.su file_out=shotw.su
2525

2626
# for displaying results

‎marchenko/demo/oneD/primariesFrame.scr‎

Lines changed: 22 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -23,13 +23,13 @@ fconv file_in1=shotsx0.su file_in2=wave.su file_out=shotw.su
2323
basop file_in=shotw.su choice=5 | sugain scale=-1 > DD.su
2424

2525
nts=1024
26-
itime=300 # select time sample to work on
26+
itime=300 # select time sample to work on 1.2 seconds
2727
shift=20
2828
#mute time nts-ii+shift to compute G_d
2929
(( itmax = nts-itime+shift ))
3030
suzero itmax=$itmax < DD.su > G_d.su
3131
#f1min = -DD(-t) = shotw.su
32-
cp shotw.su f1min.su
32+
cp shotw.su f1min0.su
3333

3434
#first iteration
3535
cp G_d.su N0.su
@@ -39,6 +39,26 @@ fconv file_in1=shotsdx5_rp.su file_in2=N0.su file_out=fconvN1.su verbose=1 fmax=
3939

4040
suwind key=fldr min=451 max=451 < fconvN1.su | suximage
4141

42+
# Ni(t) = -1*RN0(-t)
43+
sustack < fconvN1.su key=fldr > RN0.su
44+
basop file_in=RN0.su choice=5 | sugain scale=-1 > N1.su
45+
46+
#apply mute window for samples above nts-ii
47+
(( itmax = nts-itime+shift ))
48+
suzero itmax=$itmax < N1.su > N1m.su
49+
50+
51+
exit
52+
#display
53+
file=fconvN1.su
54+
file_base=${file%.su}
55+
sumax < $file mode=abs outpar=nep
56+
clipref=`cat nep | awk '{print $1/2}'`
57+
supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
58+
label1="time (s)" label2="distance (m)" \
59+
n1tic=2 d2=5 x1beg=0 x1end=1.5 d1num=0.4 \
60+
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_flip.eps
61+
4262

4363

4464
exit;

‎marchenko/marchenko.c‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -35,7 +35,7 @@ typedef struct _complexStruct { /* complex number */
3535
int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int nshots,
3636
int nx, int nxs, float fxsb, float dxs, int ntfft, int mode, float scale, float tsq, float Q, float f0, int reci, int *nshots_r, int *isxcount, int *reci_xsrc, int *reci_xrcv, float *ixmask, int verbose);
3737
int readTinvData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nfoc, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
38-
int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int iter);
38+
int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int t0shift, int iter);
3939

4040
void name_ext(char *filename, char *extension);
4141

@@ -464,7 +464,7 @@ int main (int argc, char **argv)
464464
tsyn += t3 - t2;
465465

466466
if (file_iter != NULL) {
467-
writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, iter);
467+
writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 1, iter);
468468
}
469469
/* N_k(x,t) = -N_(k-1)(x,-t) */
470470
/* p0^-(x,t) += iRN = (R * T_d^inv)(t) */

‎marchenko/marchenko_primaries.c‎

Lines changed: 99 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx
2929

3030
int readTinvData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nfoc, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
3131

32-
int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int iter);
32+
int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int t0shift, int iter);
3333
int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
3434
int readData(FILE *fp, float *data, segy *hdrs, int n1);
3535
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
@@ -102,7 +102,7 @@ int main (int argc, char **argv)
102102
{
103103
FILE *fp_out, *fp_rr, *fp_w;
104104
size_t nread;
105-
int i, j, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath;
105+
int i, j, k, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath;
106106
int size, n1, n2, ntap, tap, di, ntraces;
107107
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
108108
int reci, countmin, mode, ixa, ixb, n2out, verbose, ntfft;
@@ -114,9 +114,9 @@ int main (int argc, char **argv)
114114
double t0, t1, t2, t3, t4, tsyn, tread, tfft, tcopy, tii;
115115
float d1, d2, f1, f2, fxsb, fxse, ft, fx, *xsyn, dxsrc;
116116
float *G_d, *DD, *RR, dt, dx, dxs, scl, mem;
117-
float *rtrace, *tmpdata, *f1min, *iRN, *Ni, *trace;
117+
float *rtrace, *tmpdata, *f1min, *f1plus, *iRN, *Ni, *trace;
118118
float xmin, xmax, scale, tsq;
119-
float Q, f0, *ixmask;
119+
float Q, f0, *ixmask, *costaper;
120120
complex *Refl, *Fop, *ctrace, *cwave;
121121
char *file_tinv, *file_shot, *file_rr, *file_src, *file_iter;
122122
segy *hdrs_out, hdr;
@@ -157,9 +157,17 @@ int main (int argc, char **argv)
157157
if(!getparint("shift", &shift)) shift=20;
158158
if(!getparint("smooth", &smooth)) smooth = 5;
159159
if(!getparint("ishot", &ishot)) ishot=300;
160-
161160
if (reci && ntap) vwarn("tapering influences the reciprocal result");
162161

162+
smooth = MIN(smooth, shift);
163+
if (smooth) {
164+
costaper = (float *)malloc(smooth*sizeof(float));
165+
scl = M_PI/((float)smooth);
166+
for (i=0; i<smooth; i++) {
167+
costaper[i] = 0.5*(1.0+cos((i+1)*scl));
168+
}
169+
}
170+
163171
/*================ Reading info about shot and initial operator sizes ================*/
164172

165173
if (file_tinv != NULL) { /* G_d is read from file_tinv */
@@ -204,6 +212,7 @@ int main (int argc, char **argv)
204212

205213
Fop = (complex *)calloc(nxs*nw*Nfoc,sizeof(complex));
206214
f1min = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
215+
f1plus = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
207216
iRN = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
208217
Ni = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
209218
G_d = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
@@ -411,20 +420,21 @@ int main (int argc, char **argv)
411420
if ((NINT(di*dxs) != NINT(dxf)) && verbose)
412421
vwarn("dx in receiver (%.2f) and operator (%.2f) don't match",dx,dxs);
413422
if (verbose) {
414-
vmess("Number of receivers in focusop = %d", nxs);
415-
vmess("number of shots = %d", nshots);
416-
vmess("number of receiver/shot = %d", nx);
417-
vmess("first model position = %.2f", fxsb);
418-
vmess("last model position = %.2f", fxse);
419-
vmess("first source position fxf = %.2f", fxf);
420-
vmess("source distance dxsrc = %.2f", dxsrc);
421-
vmess("last source position = %.2f", fxf+(nshots-1)*dxsrc);
422-
vmess("receiver distance dxf = %.2f", dxf);
423+
vmess("Number of receivers in focusop = %d", nxs);
424+
vmess("number of shots = %d", nshots);
425+
vmess("number of receiver/shot = %d", nx);
426+
vmess("first model position = %.2f", fxsb);
427+
vmess("last model position = %.2f", fxse);
428+
vmess("first source position fxf = %.2f", fxf);
429+
vmess("source distance dxsrc = %.2f", dxsrc);
430+
vmess("last source position = %.2f", fxf+(nshots-1)*dxsrc);
431+
vmess("receiver distance dxf = %.2f", dxf);
423432
vmess("direction of increasing traces = %d", di);
424-
vmess("number of time samples (nt,nts) = %d (%d,%d)", ntfft, nt, nts);
425-
vmess("time sampling = %e ", dt);
426-
if (file_rr != NULL) vmess("RR output file = %s ", file_rr);
427-
if (file_iter != NULL) vmess("Iterations output file = %s ", file_iter);
433+
vmess("number of time samples fft nt nts = %d %d %d", ntfft, nt, nts);
434+
vmess("time sampling = %e ", dt);
435+
vmess("smoothing taper for time-window = %d ", smooth);
436+
if (file_rr != NULL) vmess("RR output file = %s ", file_rr);
437+
if (file_iter != NULL) vmess("Iterations output file = %s ", file_iter);
428438
}
429439

430440
/*================ initializations ================*/
@@ -433,9 +443,10 @@ int main (int argc, char **argv)
433443
else n2out = nshots;
434444
mem = Nfoc*n2out*ntfft*sizeof(float)/1048576.0;
435445
if (verbose) {
436-
vmess("number of output traces = %d", n2out);
437-
vmess("number of output samples = %d", ntfft);
438-
vmess("Size of output data/file = %.1f MB", mem);
446+
vmess("Time-sample range processed = %d:%d", istart, iend);
447+
vmess("number of output traces = %d", n2out);
448+
vmess("number of output samples = %d", ntfft);
449+
vmess("Size of output data/file = %.1f MB", mem);
439450
}
440451

441452
ixa=0; ixb=0;
@@ -484,7 +495,7 @@ int main (int argc, char **argv)
484495
vmess("*******************************************");
485496
fprintf(stderr," %s: Progress: %3d%%",xargv[0],0);
486497
}
487-
perc=(iend-istart)/10;if(!perc)perc=1;
498+
perc=(iend-istart)/100;if(!perc)perc=1;
488499

489500
/*================ start loop over number of time-samples ================*/
490501

@@ -503,9 +514,15 @@ int main (int argc, char **argv)
503514
G_d[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+i*nts+j];
504515
}
505516
/* apply mute window for samples above nts-ii */
506-
for (j = 0; j < nts-ii+T*shift; j++) {
517+
for (j = 0; j < nts-ii+T*shift-smooth; j++) {
507518
G_d[l*nxs*nts+i*nts+j] = 0.0;
508519
}
520+
for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) {
521+
G_d[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
522+
}
523+
//for (j = 0; j < nts-ii+T*shift; j++) {
524+
// G_d[l*nxs*nts+i*nts+j] = 0.0;
525+
//}
509526
}
510527
for (i = 0; i < npos; i++) {
511528
ix = ixpos[i];
@@ -514,6 +531,10 @@ int main (int argc, char **argv)
514531
for (j = 1; j < nts; j++) {
515532
f1min[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j];
516533
}
534+
/* apply mute window for samples above nts-ii */
535+
//for (j = 0; j < nts-ii+T*shift; j++) {
536+
// f1min[l*nxs*nts+i*nts+j] = 0.0;
537+
//}
517538
}
518539
}
519540
}
@@ -529,12 +550,21 @@ int main (int argc, char **argv)
529550
for (j = 1; j < nts; j++) {
530551
G_d[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j] - f1min[l*nxs*nts+i*nts+nts-j];
531552
}
532-
for (j = 0; j < nts-ii+T*shift; j++) {
553+
/* apply mute window for samples above nts-ii */
554+
for (j = 0; j < nts-ii+T*shift-smooth; j++) {
533555
G_d[l*nxs*nts+i*nts+j] = 0.0;
534556
}
557+
for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) {
558+
G_d[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
559+
}
560+
//for (j = 0; j < nts-ii+T*shift; j++) {
561+
// G_d[l*nxs*nts+i*nts+j] = 0.0;
562+
//}
535563
}
536564
}
537565
}
566+
/*================ initialization ================*/
567+
538568
memcpy(Ni, G_d, Nfoc*nxs*ntfft*sizeof(float));
539569

540570
/*================ number of Marchenko iterations ================*/
@@ -550,7 +580,7 @@ int main (int argc, char **argv)
550580
reci, nshots, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv, ixmask, verbose);
551581

552582
if (file_iter != NULL) {
553-
writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 1000*ii+iter);
583+
writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 1, 1000*ii+iter);
554584
}
555585

556586
t3 = wallclock_time();
@@ -568,20 +598,35 @@ int main (int argc, char **argv)
568598
}
569599
}
570600

571-
if (iter % 2 == 0) { /* even iterations: => f_1^-(t) */
601+
if (iter % 2 == 0) { /* even iterations: => f_1^+(t) */
572602
/* apply muting for the acausal part */
573603
for (l = 0; l < Nfoc; l++) {
574604
for (i = 0; i < npos; i++) {
605+
/* apply mute window for samples after ii */
575606
for (j = ii-T*shift; j < nts; j++) {
576607
Ni[l*nxs*nts+i*nts+j] = 0.0;
577608
}
578-
for (j = 0; j < shift; j++) {
609+
for (j = ii-T*shift+smooth, k=0; j < ii-T*shift; j++, k++) {
610+
Ni[l*nxs*nts+i*nts+j] *= costaper[k];
611+
}
612+
/* apply mute window for delta function at t=0*/
613+
for (j = 0; j < shift-smooth; j++) {
579614
Ni[l*nxs*nts+i*nts+j] = 0.0;
580615
}
616+
for (j = shift-smooth, k=1; j < shift; j++, k++) {
617+
Ni[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
618+
}
619+
f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
620+
for (j = 1; j < nts; j++) {
621+
f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
622+
}
581623
}
582624
}
625+
if (file_iter != NULL) {
626+
writeDataIter("f1plus.su", f1plus, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter);
627+
}
583628
}
584-
else {/* odd iterations: => f_1^+(t) */
629+
else {/* odd iterations: => f_1^-(t) */
585630
for (l = 0; l < Nfoc; l++) {
586631
for (i = 0; i < npos; i++) {
587632
ix = ixpos[i];
@@ -602,17 +647,30 @@ int main (int argc, char **argv)
602647
f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j];
603648
}
604649
}
605-
/* apply mute window */
606-
for (j = nts-shift; j < nts; j++) {
607-
Ni[l*nxs*nts+i*nts+j] = 0.0;
608-
}
650+
/* apply mute window for delta function at t=0*/
651+
for (j = nts-shift+smooth; j < nts; j++) {
652+
Ni[l*nxs*nts+i*nts+j] = 0.0;
653+
}
654+
for (j = nts-shift, k=0; j < nts-shift+smooth; j++, k++) {
655+
Ni[l*nxs*nts+i*nts+j] *= costaper[k];
656+
}
609657
/* apply mute window for samples above nts-ii */
610-
for (j = 0; j < nts-ii+T*shift; j++) {
611-
Ni[l*nxs*nts+i*nts+j] = 0.0;
612-
}
658+
/* important to apply this mute after updating f1min */
659+
//for (j = 0; j < nts-ii+T*shift; j++) {
660+
// Ni[l*nxs*nts+i*nts+j] = 0.0;
661+
//}
662+
/* apply mute window for samples above nts-ii */
663+
for (j = 0; j < nts-ii+T*shift-smooth; j++) {
664+
Ni[l*nxs*nts+i*nts+j] = 0.0;
665+
}
666+
for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) {
667+
Ni[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
668+
}
613669
}
614670
}
615-
//fprintf(stderr,"f1min at %d = %f NI=%f\n", ii, f1min[(nxs/2)*nts+ii], Ni[(nxs/2)*nts+ii]);
671+
if (file_iter != NULL) {
672+
writeDataIter("f1min.su", f1min, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter);
673+
}
616674
} /* end else (iter) branch */
617675

618676
t2 = wallclock_time();
@@ -627,13 +685,13 @@ int main (int argc, char **argv)
627685
}
628686
}
629687

630-
/* To Do optional write intermediate RR results to file */
688+
/* To Do? optional write intermediate RR results to file */
631689

632690
if (verbose) {
633691
if(!((iend-ii-istart)%perc)) fprintf(stderr,"\b\b\b\b%3d%%",(ii-istart)*100/(iend-istart));
634692
if((ii-istart)==10)t4=wallclock_time();
635-
if((ii-istart)==50){
636-
t4=(wallclock_time()-t4)*((iend-istart)/40.0);
693+
if((ii-istart)==20){
694+
t4=(wallclock_time()-t4)*((iend-istart)/10.0);
637695
fprintf(stderr,"\r %s: Estimated total compute time = %.2fs.\n %s: Progress: %3d%%",xargv[0],t4,xargv[0],(ii-istart)/((iend-istart)/100));
638696
}
639697
//t4=wallclock_time();
@@ -645,6 +703,8 @@ int main (int argc, char **argv)
645703

646704
free(Ni);
647705
free(G_d);
706+
free(f1min);
707+
free(f1plus);
648708

649709
t2 = wallclock_time();
650710
if (verbose) {

0 commit comments

Comments
 (0)