From 2c48834b1eaa6a5a9c3474e0fe2256e7dc9c1d7d Mon Sep 17 00:00:00 2001
From: JanThorbecke <janth@xs4all.nl>
Date: Thu, 19 Nov 2020 09:57:38 +0100
Subject: [PATCH] plane-wave ringing effect in Gmin reduced
---
marchenko/applyMute_tshift.c | 35 +++++++++++------------------------
marchenko/marchenko.c | 6 +++++-
2 files changed, 16 insertions(+), 25 deletions(-)
diff --git a/marchenko/applyMute_tshift.c b/marchenko/applyMute_tshift.c
index a2f956e..ac69f90 100644
--- a/marchenko/applyMute_tshift.c
+++ b/marchenko/applyMute_tshift.c
@@ -70,10 +70,10 @@ void applyMute_tshift( float *data, int *mute, int smooth, int above, int Nfoc,
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++) {
+ for (j = MIN(nt,nt-tmute+2*ts+shift-smooth),l=0; j < MIN(nt,nt-tmute+shift); j++,l++) {
Nig[j] *= costaper[smooth-l-1];
}
}
@@ -93,18 +93,18 @@ void applyMute_tshift( float *data, int *mute, int smooth, int above, int Nfoc,
imute = ixpos[i];
tmute = mute[isyn*nxs+imute];
ts = tsynW[isyn*nxs+imute];
+ for (j = 0; j < MAX(0,-2*ts+tmute-shift-smooth); j++) {
+ Nig[j] = 0.0;
+ }
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 < MAX(0,-2*ts+tmute-shift-smooth-1); j++) {
- Nig[j] = 0.0;
+ for (j = MIN(nt,nt-tmute+shift),l=0; j < MIN(nt,nt-tmute+shift+smooth); j++,l++) {
+ Nig[j] *= costaper[l];
}
- for (j = nt+1-tmute+shift+smooth; j < nt; j++) {
+ for (j = MIN(nt,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++) {
- Nig[j] *= costaper[l];
- }
}
/* To Do above==2 is not yet adapated for plane-waves */
else if (above==2){//Separates the direct part of the wavefield from the coda
@@ -147,8 +147,8 @@ void applyMute_tshift( float *data, int *mute, int smooth, int above, int Nfoc,
void timeShift(float *data, int nsam, int nrec, float dt, float shift, float fmin, float fmax)
{
- int optn, iom, iomin, iomax, nfreq, ix, it;
- float omin, omax, deltom, om, tom, df, *trace, scl;
+ int optn, iom, nfreq, ix, it;
+ float deltom, om, tom, df, *trace, scl;
complex *ctrace, ctmp;
optn = optncr(nsam);
@@ -162,10 +162,6 @@ void timeShift(float *data, int nsam, int nrec, float dt, float shift, float fmi
if (trace == NULL) verr("memory allocation error for rdata");
deltom = 2.*M_PI*df;
- omin = 2.*M_PI*fmin;
- omax = 2.*M_PI*fmax;
- iomin = (int)MIN((omin/deltom), (nfreq));
- iomax = MIN((int)(omax/deltom), (nfreq));
scl = 1.0/(float)optn;
for (ix = 0; ix < nrec; ix++) {
@@ -173,16 +169,7 @@ void timeShift(float *data, int nsam, int nrec, float dt, float shift, float fmi
for (it=nsam;it<optn;it++) trace[it]=0.0;
/* Forward time-frequency FFT */
rc1fft(&trace[0], &ctrace[0], optn, -1);
-
- for (iom = 0; iom < iomin; iom++) {
- ctrace[iom].r = 0.0;
- ctrace[iom].i = 0.0;
- }
- for (iom = iomax; iom < nfreq; iom++) {
- ctrace[iom].r = 0.0;
- ctrace[iom].i = 0.0;
- }
- for (iom = iomin ; iom < iomax ; iom++) {
+ for (iom = 0 ; iom < nfreq ; iom++) {
om = deltom*iom;
tom = om*shift;
ctmp = ctrace[iom];
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index 0d0f1ae..6189de2 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -581,9 +581,13 @@ int main (int argc, char **argv)
}
/* Apply mute with window for Gmin */
if ( plane_wave==1 ) {
- applyMute_tshift(Gmin, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, 0, tsynW);
/* for plane wave with angle tshift downward */
+ if (verbose>1) vmess("Gmin planewave tshift=%f\n", tshift);
timeShift(Gmin, nts, npos, dt, tshift, fmin, fmax);
+ itmin = NINT(0.5*tshift/dt);
+ for (i=0; i<nxs; i++) tsynW[i] -= itmin;
+ applyMute_tshift(Gmin, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, 0, tsynW);
+ for (i=0; i<nxs; i++) tsynW[i] += itmin;
}
else {
applyMute(Gmin, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
--
GitLab