Skip to content

Commit

Permalink
taper in time window in primaries app
Browse files Browse the repository at this point in the history
  • Loading branch information
@JanThorbecke
JanThorbecke committed Aug 8, 2019
1 parent 87c7fb8 commit 289151c
Show file tree
Hide file tree
Showing 6 changed files with 165 additions and 56 deletions.
24 changes: 24 additions & 0 deletions marchenko/demo/oneD/iterations.scr
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,24 @@
#!/bin/bash

#second reflector at time:
# 800/1800+600/2300
# .70531400966183574878 => sample 176
#thrids reflector at time:
# 800/1800+600/2300+800/2000
# 1.10531400966183574878 sample 276

select=451

#for istart in 176 200 276 400
for istart in 176
do
(( iend = istart + 1 ))

../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
nshots=601 verbose=1 istart=$istart iend=$iend fmax=90 \
pad=1024 niter=45 smooth=10 niterskip=600 shift=20 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su

suximage < f1min_${istart}043.su clip=1 title="${istart}" &

done

8 changes: 4 additions & 4 deletions marchenko/demo/oneD/primaries.scr
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,14 +13,14 @@ export OMP_NUM_THREADS=20
#shot record to remove internal multiples
select=451

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

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

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

# for displaying results
Expand Down
24 changes: 22 additions & 2 deletions marchenko/demo/oneD/primariesFrame.scr
View file
Original file line number Diff line number Diff line change
Expand Up @@ -23,13 +23,13 @@ fconv file_in1=shotsx0.su file_in2=wave.su file_out=shotw.su
basop file_in=shotw.su choice=5 | sugain scale=-1 > DD.su

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

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

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

# Ni(t) = -1*RN0(-t)
sustack < fconvN1.su key=fldr > RN0.su
basop file_in=RN0.su choice=5 | sugain scale=-1 > N1.su

#apply mute window for samples above nts-ii
(( itmax = nts-itime+shift ))
suzero itmax=$itmax < N1.su > N1m.su


exit
#display
file=fconvN1.su
file_base=${file%.su}
sumax < $file mode=abs outpar=nep
clipref=`cat nep | awk '{print $1/2}'`
supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
label1="time (s)" label2="distance (m)" \
n1tic=2 d2=5 x1beg=0 x1end=1.5 d1num=0.4 \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_flip.eps



exit;
Expand Down
4 changes: 2 additions & 2 deletions marchenko/marchenko.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ typedef struct _complexStruct { /* complex number */
int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int nshots,
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);
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);
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);
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);

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

Expand Down Expand Up @@ -464,7 +464,7 @@ int main (int argc, char **argv)
tsyn += t3 - t2;

if (file_iter != NULL) {
writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, iter);
writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 1, iter);
}
/* N_k(x,t) = -N_(k-1)(x,-t) */
/* p0^-(x,t) += iRN = (R * T_d^inv)(t) */
Expand Down
138 changes: 99 additions & 39 deletions marchenko/marchenko_primaries.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx

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);

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);
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);
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);
int readData(FILE *fp, float *data, segy *hdrs, int n1);
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
Expand Down Expand Up @@ -102,7 +102,7 @@ int main (int argc, char **argv)
{
FILE *fp_out, *fp_rr, *fp_w;
size_t nread;
int i, j, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath;
int i, j, k, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath;
int size, n1, n2, ntap, tap, di, ntraces;
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
int reci, countmin, mode, ixa, ixb, n2out, verbose, ntfft;
Expand All @@ -114,9 +114,9 @@ int main (int argc, char **argv)
double t0, t1, t2, t3, t4, tsyn, tread, tfft, tcopy, tii;
float d1, d2, f1, f2, fxsb, fxse, ft, fx, *xsyn, dxsrc;
float *G_d, *DD, *RR, dt, dx, dxs, scl, mem;
float *rtrace, *tmpdata, *f1min, *iRN, *Ni, *trace;
float *rtrace, *tmpdata, *f1min, *f1plus, *iRN, *Ni, *trace;
float xmin, xmax, scale, tsq;
float Q, f0, *ixmask;
float Q, f0, *ixmask, *costaper;
complex *Refl, *Fop, *ctrace, *cwave;
char *file_tinv, *file_shot, *file_rr, *file_src, *file_iter;
segy *hdrs_out, hdr;
Expand Down Expand Up @@ -157,9 +157,17 @@ int main (int argc, char **argv)
if(!getparint("shift", &shift)) shift=20;
if(!getparint("smooth", &smooth)) smooth = 5;
if(!getparint("ishot", &ishot)) ishot=300;

if (reci && ntap) vwarn("tapering influences the reciprocal result");

smooth = MIN(smooth, shift);
if (smooth) {
costaper = (float *)malloc(smooth*sizeof(float));
scl = M_PI/((float)smooth);
for (i=0; i<smooth; i++) {
costaper[i] = 0.5*(1.0+cos((i+1)*scl));
}
}

/*================ Reading info about shot and initial operator sizes ================*/

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

Fop = (complex *)calloc(nxs*nw*Nfoc,sizeof(complex));
f1min = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
f1plus = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
iRN = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
Ni = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
G_d = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
Expand Down Expand Up @@ -411,20 +420,21 @@ int main (int argc, char **argv)
if ((NINT(di*dxs) != NINT(dxf)) && verbose)
vwarn("dx in receiver (%.2f) and operator (%.2f) don't match",dx,dxs);
if (verbose) {
vmess("Number of receivers in focusop = %d", nxs);
vmess("number of shots = %d", nshots);
vmess("number of receiver/shot = %d", nx);
vmess("first model position = %.2f", fxsb);
vmess("last model position = %.2f", fxse);
vmess("first source position fxf = %.2f", fxf);
vmess("source distance dxsrc = %.2f", dxsrc);
vmess("last source position = %.2f", fxf+(nshots-1)*dxsrc);
vmess("receiver distance dxf = %.2f", dxf);
vmess("Number of receivers in focusop = %d", nxs);
vmess("number of shots = %d", nshots);
vmess("number of receiver/shot = %d", nx);
vmess("first model position = %.2f", fxsb);
vmess("last model position = %.2f", fxse);
vmess("first source position fxf = %.2f", fxf);
vmess("source distance dxsrc = %.2f", dxsrc);
vmess("last source position = %.2f", fxf+(nshots-1)*dxsrc);
vmess("receiver distance dxf = %.2f", dxf);
vmess("direction of increasing traces = %d", di);
vmess("number of time samples (nt,nts) = %d (%d,%d)", ntfft, nt, nts);
vmess("time sampling = %e ", dt);
if (file_rr != NULL) vmess("RR output file = %s ", file_rr);
if (file_iter != NULL) vmess("Iterations output file = %s ", file_iter);
vmess("number of time samples fft nt nts = %d %d %d", ntfft, nt, nts);
vmess("time sampling = %e ", dt);
vmess("smoothing taper for time-window = %d ", smooth);
if (file_rr != NULL) vmess("RR output file = %s ", file_rr);
if (file_iter != NULL) vmess("Iterations output file = %s ", file_iter);
}

/*================ initializations ================*/
Expand All @@ -433,9 +443,10 @@ int main (int argc, char **argv)
else n2out = nshots;
mem = Nfoc*n2out*ntfft*sizeof(float)/1048576.0;
if (verbose) {
vmess("number of output traces = %d", n2out);
vmess("number of output samples = %d", ntfft);
vmess("Size of output data/file = %.1f MB", mem);
vmess("Time-sample range processed = %d:%d", istart, iend);
vmess("number of output traces = %d", n2out);
vmess("number of output samples = %d", ntfft);
vmess("Size of output data/file = %.1f MB", mem);
}

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

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

Expand All @@ -503,9 +514,15 @@ int main (int argc, char **argv)
G_d[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+i*nts+j];
}
/* apply mute window for samples above nts-ii */
for (j = 0; j < nts-ii+T*shift; j++) {
for (j = 0; j < nts-ii+T*shift-smooth; j++) {
G_d[l*nxs*nts+i*nts+j] = 0.0;
}
for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) {
G_d[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
}
//for (j = 0; j < nts-ii+T*shift; j++) {
// G_d[l*nxs*nts+i*nts+j] = 0.0;
//}
}
for (i = 0; i < npos; i++) {
ix = ixpos[i];
Expand All @@ -514,6 +531,10 @@ int main (int argc, char **argv)
for (j = 1; j < nts; j++) {
f1min[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j];
}
/* apply mute window for samples above nts-ii */
//for (j = 0; j < nts-ii+T*shift; j++) {
// f1min[l*nxs*nts+i*nts+j] = 0.0;
//}
}
}
}
Expand All @@ -529,12 +550,21 @@ int main (int argc, char **argv)
for (j = 1; j < nts; j++) {
G_d[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j] - f1min[l*nxs*nts+i*nts+nts-j];
}
for (j = 0; j < nts-ii+T*shift; j++) {
/* apply mute window for samples above nts-ii */
for (j = 0; j < nts-ii+T*shift-smooth; j++) {
G_d[l*nxs*nts+i*nts+j] = 0.0;
}
for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) {
G_d[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
}
//for (j = 0; j < nts-ii+T*shift; j++) {
// G_d[l*nxs*nts+i*nts+j] = 0.0;
//}
}
}
}
/*================ initialization ================*/

memcpy(Ni, G_d, Nfoc*nxs*ntfft*sizeof(float));

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

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

t3 = wallclock_time();
Expand All @@ -568,20 +598,35 @@ int main (int argc, char **argv)
}
}

if (iter % 2 == 0) { /* even iterations: => f_1^-(t) */
if (iter % 2 == 0) { /* even iterations: => f_1^+(t) */
/* apply muting for the acausal part */
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
/* apply mute window for samples after ii */
for (j = ii-T*shift; j < nts; j++) {
Ni[l*nxs*nts+i*nts+j] = 0.0;
}
for (j = 0; j < shift; j++) {
for (j = ii-T*shift+smooth, k=0; j < ii-T*shift; j++, k++) {
Ni[l*nxs*nts+i*nts+j] *= costaper[k];
}
/* apply mute window for delta function at t=0*/
for (j = 0; j < shift-smooth; j++) {
Ni[l*nxs*nts+i*nts+j] = 0.0;
}
for (j = shift-smooth, k=1; j < shift; j++, k++) {
Ni[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
}
f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
}
}
}
if (file_iter != NULL) {
writeDataIter("f1plus.su", f1plus, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter);
}
}
else {/* odd iterations: => f_1^+(t) */
else {/* odd iterations: => f_1^-(t) */
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
ix = ixpos[i];
Expand All @@ -602,17 +647,30 @@ int main (int argc, char **argv)
f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j];
}
}
/* apply mute window */
for (j = nts-shift; j < nts; j++) {
Ni[l*nxs*nts+i*nts+j] = 0.0;
}
/* apply mute window for delta function at t=0*/
for (j = nts-shift+smooth; j < nts; j++) {
Ni[l*nxs*nts+i*nts+j] = 0.0;
}
for (j = nts-shift, k=0; j < nts-shift+smooth; j++, k++) {
Ni[l*nxs*nts+i*nts+j] *= costaper[k];
}
/* apply mute window for samples above nts-ii */
for (j = 0; j < nts-ii+T*shift; j++) {
Ni[l*nxs*nts+i*nts+j] = 0.0;
}
/* important to apply this mute after updating f1min */
//for (j = 0; j < nts-ii+T*shift; j++) {
// Ni[l*nxs*nts+i*nts+j] = 0.0;
//}
/* apply mute window for samples above nts-ii */
for (j = 0; j < nts-ii+T*shift-smooth; j++) {
Ni[l*nxs*nts+i*nts+j] = 0.0;
}
for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) {
Ni[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
}
}
}
//fprintf(stderr,"f1min at %d = %f NI=%f\n", ii, f1min[(nxs/2)*nts+ii], Ni[(nxs/2)*nts+ii]);
if (file_iter != NULL) {
writeDataIter("f1min.su", f1min, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter);
}
} /* end else (iter) branch */

t2 = wallclock_time();
Expand All @@ -627,13 +685,13 @@ int main (int argc, char **argv)
}
}

/* To Do optional write intermediate RR results to file */
/* To Do? optional write intermediate RR results to file */

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

free(Ni);
free(G_d);
free(f1min);
free(f1plus);

t2 = wallclock_time();
if (verbose) {
Expand Down
Loading

0 comments on commit 289151c

Please sign in to comment.