diff --git a/marchenko/applyMute.c b/marchenko/applyMute.c
index df98e8bc39f8b132c10b237f7e3e01f3167cdf7d..acbbeea27dece52a24b5f0ad366f2b2a1351f06d 100644
--- a/marchenko/applyMute.c
+++ b/marchenko/applyMute.c
@@ -12,11 +12,11 @@
#endif
#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
-void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *ixpos, int npos, int shift)
+void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *ixpos, int npos, int shift, int *tsynW)
{
int i, j, l, isyn;
float *costaper, scl;
- int imute, tmute;
+ int imute, tmute, ts;
if (smooth) {
costaper = (float *)malloc(smooth*sizeof(float));
@@ -31,10 +31,11 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs
for (i = 0; i < npos; i++) {
imute = ixpos[i];
tmute = mute[isyn*nxs+imute];
- for (j = 0; j < MAX(0,tmute-shift-smooth); j++) {
+ ts = tsynW[isyn*nxs+imute];
+ for (j = 0; j < MAX(0,-2*ts+tmute-shift-smooth); j++) {
data[isyn*nxs*nt+i*nt+j] = 0.0;
}
- for (j = MAX(0,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); j++,l++) {
+ for (j = MAX(0,-2*ts+tmute-shift-smooth),l=0; j < MAX(0,-2*ts+tmute-shift); j++,l++) {
data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
}
}
@@ -43,14 +44,15 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs
for (i = 0; i < npos; i++) {
imute = ixpos[i];
tmute = mute[isyn*nxs+imute];
+ ts = tsynW[isyn*nxs+imute];
if (tmute >= nt/2) {
memset(&data[isyn*nxs*nt+i*nt],0, sizeof(float)*nt);
continue;
}
- for (j = MAX(0,tmute-shift),l=0; j < MAX(0,tmute-shift+smooth); j++,l++) {
+ for (j = MAX(0,-2*ts+tmute-shift),l=0; j < MAX(0,-2*ts+tmute-shift+smooth); j++,l++) {
data[isyn*nxs*nt+i*nt+j] *= costaper[l];
}
- for (j = MAX(0,tmute-shift+smooth)+1; j < MIN(nt,nt+1-tmute+shift-smooth); j++) {
+ for (j = MAX(0,-2*ts+tmute-shift+smooth)+1; j < MIN(nt,nt+1-tmute+shift-smooth); j++) {
data[isyn*nxs*nt+i*nt+j] = 0.0;
}
for (j = MIN(nt,nt-tmute+shift-smooth),l=0; j < MIN(nt,nt-tmute+shift); j++,l++) {
@@ -62,10 +64,11 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs
for (i = 0; i < npos; i++) {
imute = ixpos[i];
tmute = mute[isyn*nxs+imute];
- for (j = MAX(0,tmute-shift),l=0; j < MAX(0,tmute-shift+smooth); j++,l++) {
+ ts = tsynW[isyn*nxs+imute];
+ for (j = MAX(0,ts+tmute-shift),l=0; j < MAX(0,ts+tmute-shift+smooth); j++,l++) {
data[isyn*nxs*nt+i*nt+j] *= costaper[l];
}
- for (j = MAX(0,tmute-shift+smooth); j < nt; j++) {
+ for (j = MAX(0,ts+tmute-shift+smooth); j < nt; j++) {
data[isyn*nxs*nt+i*nt+j] = 0.0;
}
}
diff --git a/marchenko/fmute.c b/marchenko/fmute.c
index ba4f39acb407d3dacf414096dafc0b3ab67a2c8d..4d63faf38fa8bb4cb489ee554bca07e52ce9a439 100644
--- a/marchenko/fmute.c
+++ b/marchenko/fmute.c
@@ -25,7 +25,7 @@ int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *
int readData(FILE *fp, float *data, segy *hdrs, int n1);
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
-void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *ixpos, int npos, int shift);
+void applyMute(float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *xrcvsyn, int npos, int shift, int *muteW);
double wallclock_time(void);
/*********************** self documentation **********************/
@@ -43,7 +43,7 @@ char *sdoc[] = {
" Optional parameters: ",
" ",
" file_out= ................ output file",
-" above=0 .................. mute after(0), before(1) or around(2) the maximum times of file_mute",
+" above=0 .................. mute above(1), around(0) or below(-1) the first travel times of file_tinv",
" .......................... options 4 is the inverse of 0 and -1 the inverse of 1",
" shift=0 .................. number of points above(positive) / below(negative) maximum time for mute",
" check=0 .................. plots muting window on top of file_mute: output file check.su",
@@ -64,7 +64,7 @@ int main (int argc, char **argv)
FILE *fp_in1, *fp_in2, *fp_out, *fp_chk, *fp_psline1, *fp_psline2;
int verbose, shift, k, nx1, nt1, nx2, nt2;
int ntmax, nxmax, ret, i, j, jmax, imax, above, check;
- int size, ntraces, ngath, *maxval, hw, smooth;
+ int size, ntraces, ngath, *maxval, *tsynW, hw, smooth;
int tstart, tend, scale, *xrcv;
float dt, d2, f1, f2, t0, t1, f1b, f2b, d1, d1b, d2b;
float w1, w2, dxrcv;
@@ -169,6 +169,7 @@ int main (int argc, char **argv)
/*================ initializations ================*/
maxval = (int *)calloc(nx1,sizeof(int));
+ tsynW = (int *)calloc(nx1,sizeof(int));
xrcv = (int *)calloc(nx1,sizeof(int));
if (file_out==NULL) fp_out = stdout;
@@ -289,7 +290,7 @@ int main (int argc, char **argv)
/*================ apply mute window ================*/
- applyMute(tmpdata2, maxval, smooth, above, 1, nx2, nt2, xrcv, nx2, shift);
+ applyMute(tmpdata2, maxval, smooth, above, 1, nx2, nt2, xrcv, nx2, shift, tsynW);
/*================ write result to output file ================*/
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index 9812335d5eeca1883de2be2810350e053687f58b..2d16e436f83b026b592fdf28a5bc416497bbb7a4 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -43,7 +43,7 @@ float *zsyn, int *ixpos, int npos, int iter);
void name_ext(char *filename, char *extension);
-void applyMute(float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *xrcvsyn, int npos, int shift);
+void applyMute(float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *xrcvsyn, int npos, int shift, int *muteW);
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 *ntraces);
int readData(FILE *fp, float *data, segy *hdrs, int n1);
@@ -87,6 +87,9 @@ char *sdoc[] = {
" shift=12 ................. number of points above(positive) / below(negative) travel time for mute",
" hw=8 ..................... window in time samples to look for maximum in next trace",
" smooth=5 ................. number of points to smooth mute with cosine window",
+" plane_wave=0 ............. enable plane-wave illumination function"
+" src_angle=0 .............. angle of plane source array",
+" src_velo=1500 ............ velocity to use in src_angle definition",
" REFLECTION RESPONSE CORRECTION ",
" tsq=0.0 .................. scale factor n for t^n for true amplitude recovery",
" Q=0.0 .......,............ Q correction factor",
@@ -121,16 +124,17 @@ int main (int argc, char **argv)
int size, n1, n2, ntap, tap, di, ntraces, pad;
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
int reci, countmin, mode, n2out, verbose, ntfft;
- int iter, niter, tracf, *muteW;
- int hw, smooth, above, shift, *ixpos, npos, ix;
+ int iter, niter, tracf, *muteW, *tsynW;
+ int hw, smooth, above, shift, *ixpos, npos, ix, plane_wave;
int nshots_r, *isxcount, *reci_xsrc, *reci_xrcv;
float fmin, fmax, *tapersh, *tapersy, fxf, dxf, *xsrc, *xrcv, *zsyn, *zsrc, *xrcvsyn;
double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi, energyN0;
float d1, d2, f1, f2, fxsb, fxse, ft, fx, *xsyn, dxsrc;
float *green, *f2p, *pmin, *G_d, dt, dx, dxs, scl, mem;
- float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus;
+ float *f1plus, *f1min, *iRN, *Ni, *Nig, *trace, *Gmin, *Gplus;
float xmin, xmax, scale, tsq, Q, f0;
float *ixmask;
+ float grad2rad, p, src_angle, src_velo;
complex *Refl, *Fop;
char *file_tinv, *file_shot, *file_green, *file_iter;
char *file_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_pmin;
@@ -176,6 +180,10 @@ int main (int argc, char **argv)
if(!getparint("above", &above)) above = 0;
if(!getparint("shift", &shift)) shift=12;
+ if (!getparint("plane_wave", &plane_wave)) plane_wave = 0;
+ if (!getparfloat("src_angle",&src_angle)) src_angle=0.;
+ if (!getparfloat("src_velo",&src_velo)) src_velo=1500.;
+
if (reci && ntap) vwarn("tapering influences the reciprocal result");
/*================ Reading info about shot and initial operator sizes ================*/
@@ -214,8 +222,10 @@ int main (int argc, char **argv)
f1min = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
iRN = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
Ni = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ Nig = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
G_d = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
muteW = (int *)calloc(Nfoc*nxs,sizeof(int));
+ tsynW = (int *)malloc(Nfoc*nxs*sizeof(int)); // time-shift for Giovanni's plane-wave on non-zero times
trace = (float *)malloc(ntfft*sizeof(float));
tapersy = (float *)malloc(nxs*sizeof(float));
xrcvsyn = (float *)calloc(Nfoc*nxs,sizeof(float)); // x-rcv postions of focal points
@@ -247,6 +257,28 @@ int main (int argc, char **argv)
/* reading data added zero's to the number of time samples to be the same as ntfft */
nts = ntfft;
+ /* compute time shift for tilted plane waves */
+ if (plane_wave != 0) {
+ grad2rad = 17.453292e-3;
+ p = sin(src_angle*grad2rad)/src_velo;
+ if (p < 0.0) {
+ for (i=0; i<nxs; i++) {
+ tsynW[i] = NINT(fabsf((nxs-1-i)*dxs*p)/dt);
+ }
+ }
+ else {
+ for (i=0; i<nxs; i++) {
+ tsynW[i] = NINT(i*dxs*p/dt);
+ }
+ }
+ if (Nfoc!=1) verr("For plave-wave focusing only one function can be computed at the same time");
+ }
+ else { /* just fill with zero's */
+ for (i=0; i<nxs*Nfoc; i++) {
+ tsynW[i] = 0;
+ }
+ }
+
/* define tapers to taper edges of acquisition */
if (tap == 1 || tap == 3) {
for (j = 0; j < ntap; j++)
@@ -461,6 +493,12 @@ int main (int argc, char **argv)
pmin[l*nxs*nts+i*nts+j] += iRN[l*nxs*nts+ix*nts+j];
energyNi += iRN[l*nxs*nts+ix*nts+j]*iRN[l*nxs*nts+ix*nts+j];
}
+ if (plane_wave!=0) { /* don't reverse in time */
+ Nig[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+ix*nts+j];
+ for (j = 1; j < nts; j++) {
+ Nig[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+ix*nts+j];
+ }
+ }
}
if (iter==0) energyN0 = energyNi;
if (verbose >=2) vmess(" - iSyn %d: Ni at iteration %d has energy %e; relative to N0 %e", l, iter, sqrt(energyNi),
@@ -468,29 +506,34 @@ sqrt(energyNi/energyN0));
}
/* apply mute window based on times of direct arrival (in muteW) */
- applyMute(Ni, muteW, smooth, above, Nfoc, nxs, nts, ixpos, npos, shift);
-
- /* update f2 */
- for (l = 0; l < Nfoc; l++) {
- for (i = 0; i < npos; i++) {
- j = 0;
- f2p[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- f2p[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
- }
- }
- }
+ applyMute(Ni, muteW, smooth, above, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
+ if (plane_wave!=0) applyMute(Nig, muteW, smooth, above, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
if (iter % 2 == 0) { /* even iterations update: => f_1^-(t) */
- for (l = 0; l < Nfoc; l++) {
- for (i = 0; i < npos; i++) {
- j = 0;
- f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j];
+ if (plane_wave==0) { /* follow the standard focal point scheme */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j];
+ }
}
}
- }
+ }
+ else { /* plane wave scheme */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ f1min[l*nxs*nts+i*nts+j] -= Nig[l*nxs*nts+i*nts+j];
+ Ni[l*nxs*nts+i*nts+j] = Nig[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ f1min[l*nxs*nts+i*nts+j] -= Nig[l*nxs*nts+i*nts+j];
+ Ni[l*nxs*nts+i*nts+j] = Nig[l*nxs*nts+i*nts+nts-j];
+ }
+ }
+ }
+ }
}
else {/* odd iterations update: => f_1^+(t) */
for (l = 0; l < Nfoc; l++) {
@@ -504,6 +547,17 @@ sqrt(energyNi/energyN0));
}
}
+ /* update f2 */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ f2p[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ f2p[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
+ }
+ }
+ }
+
t2 = wallclock_time();
tcopy += t2 - t3;
@@ -512,6 +566,7 @@ sqrt(energyNi/energyN0));
} /* end of iterations */
free(Ni);
+ free(Nig);
free(G_d);
/* compute full Green's function G = int R * f2(t) + f2(-t) = Pplus + Pmin */
@@ -552,7 +607,7 @@ sqrt(energyNi/energyN0));
}
}
/* Apply mute with window for Gmin */
- applyMute(Gmin, muteW, smooth, 1, Nfoc, nxs, nts, ixpos, npos, shift);
+ applyMute(Gmin, muteW, smooth, 1, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
} /* end if Gmin */
/* compute downgoing Green's function G^+,+ */
@@ -578,6 +633,9 @@ sqrt(energyNi/energyN0));
}
} /* end if Gplus */
+ free(muteW);
+ free(tsynW);
+
t2 = wallclock_time();
if (verbose) {
vmess("Total CPU-time marchenko = %.3f", t2-t0);