From c45d1e7f646618b192a68743b5cf7afa9adf2623 Mon Sep 17 00:00:00 2001
From: JanThorbecke <janth@xs4all.nl>
Date: Thu, 10 Dec 2020 13:02:47 +0100
Subject: [PATCH] plane-wave option: time-windows dipping planes
---
marchenko/applyMute.c | 13 ++-
marchenko/applyMute_tshift.c | 194 +++++++++++++++++---------------
marchenko/marchenko.c | 41 +++++--
marchenko/marchenko_primaries.c | 19 ++--
4 files changed, 155 insertions(+), 112 deletions(-)
diff --git a/marchenko/applyMute.c b/marchenko/applyMute.c
index 544f0b8..782df17 100644
--- a/marchenko/applyMute.c
+++ b/marchenko/applyMute.c
@@ -52,10 +52,10 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs
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[l];
}
- for (j = MAX(0,-2*ts+tmute-shift); j < MIN(nt,nt-tmute+shift); j++) {
+ for (j = MAX(0,-2*ts+tmute-shift); j < MIN(nt,nt+1-tmute+shift+2*ts); j++) {
data[isyn*nxs*nt+i*nt+j] = 0.0;
}
- for (j = MIN(nt,nt-tmute+shift),l=0; j < MIN(nt,nt-tmute+shift+smooth); j++,l++) {
+ for (j = MIN(nt,nt+1-tmute+shift+2*ts),l=0; j < MIN(nt,nt+1-tmute+shift+2*ts+smooth); j++,l++) {
data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
}
}
@@ -64,16 +64,17 @@ 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,tmute-2*ts-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];
}
- for (j = MIN(nt,nt-tmute+shift),l=0; j < MIN(nt,nt-tmute+shift+smooth); j++,l++) {
+ for (j = MIN(nt,nt+1-tmute+shift+2*ts),l=0; j < MIN(nt,nt+1-tmute+shift+2*ts+smooth); j++,l++) {
data[isyn*nxs*nt+i*nt+j] *= costaper[l];
}
- for (j = MIN(nt,nt-tmute+shift+smooth); j < nt; j++) {
+ for (j = MIN(nt,nt+1-tmute+shift+2*ts+smooth); j < nt; j++) {
data[isyn*nxs*nt+i*nt+j] = 0.0;
}
}
diff --git a/marchenko/applyMute_tshift.c b/marchenko/applyMute_tshift.c
index ac69f90..1863396 100644
--- a/marchenko/applyMute_tshift.c
+++ b/marchenko/applyMute_tshift.c
@@ -26,7 +26,7 @@ void applyMute_tshift( float *data, int *mute, int smooth, int above, int Nfoc,
if (smooth) {
costaper = (float *)malloc(smooth*sizeof(float));
- scl = M_PI/((float)smooth);
+ scl = M_PI/((float)smooth+1);
for (i=0; i<smooth; i++) {
costaper[i] = 0.5*(1.0+cos((i+1)*scl));
}
@@ -36,150 +36,166 @@ void applyMute_tshift( float *data, int *mute, int smooth, int above, int Nfoc,
for (isyn = 0; isyn < Nfoc; isyn++) {
for (i = 0; i < npos; i++) {
+ imute = ixpos[i];
+ tmute = mute[isyn*nxs+imute];
+ ts = tsynW[isyn*nxs+imute];
+ //fprintf(stderr,"i=%d tmute=%d ts=%d\n", i, tmute, ts);
+ for (j = 0; j < nt; j++) {
+ Nig[j] = data[isyn*nxs*nt+i*nt+j];
+ }
if (iter % 2 == 0) {
- for (j = 0; j < nt; j++) {
- Nig[j] = data[isyn*nxs*nt+i*nt+j];
- }
+ if (above==0) above=0;
}
- else { // reverse back in time
- j=0;
- Nig[j] = data[isyn*nxs*nt+i*nt+j];
- for (j = 1; j < nt; j++) {
- Nig[j] = data[isyn*nxs*nt+i*nt+nt-j];
- }
+ else { // switch angle
+ if (above==0) above=-1;
+ if (above==4) above=-4;
+ tmute = tmute-2*ts;
}
- if (above==1) {
- imute = ixpos[i];
- tmute = mute[isyn*nxs+imute];
- ts = tsynW[isyn*nxs+imute];
- for (j = 0; j < MAX(0,tmute-shift-smooth); j++) {
- Nig[j] = 0.0;
+ //fprintf(stderr,"above = %d\n", above);
+ if (above==-1){ /* the above=0 implementation for plane-waves at odd iterations */
+ if (tmute >= nt/2) {
+ memset(&Nig[0],0, sizeof(float)*nt);
+ continue;
}
+ /* positive time axis mute below plane-wave first arrivals */
+ /* at the zero crossing the smooth transition zone is not handled correctly
+ * this can be implemented but makes it more complicated, better solution is to
+ * rotate Nig with the time axis in the middle, something still to do */
for (j = MAX(0,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); j++,l++) {
+ Nig[j] *= costaper[l];
+ }
+ for (j = MAX(0,tmute-shift); j < MIN(nt,nt+1-tmute+shift-2*ts); j++) {
+ Nig[j] = 0.0;
+ }
+ if (tmute-shift < 0) {
+ for (j = MIN(nt,nt+1+tmute-shift-smooth),l=0; j < MIN(nt,nt+1+tmute-shift); j++,l++) {
+ Nig[j] *= costaper[l];
+ }
+ for (j = MIN(nt,nt+1+tmute-shift); j < nt; j++) {
+ Nig[j] = 0.0;
+ }
+ }
+ /* negative time axis */
+ for (j = MIN(nt,nt+1-(tmute+2*ts)+shift),l=0; j < MIN(nt,nt+1-(tmute+2*ts)+shift+smooth); j++,l++) {
Nig[j] *= costaper[smooth-l-1];
}
}
else if (above==0){
- imute = ixpos[i];
- tmute = mute[isyn*nxs+imute];
- ts = tsynW[isyn*nxs+imute];
if (tmute >= nt/2) {
memset(&Nig[0],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,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); j++,l++) {
Nig[j] *= costaper[l];
}
- for (j = MAX(0,tmute-shift+smooth); j < MIN(nt,nt-tmute+2*ts+shift-smooth); j++) {
+ for (j = MAX(0,tmute-shift); j < MIN(nt,nt+1-tmute+shift+2*ts); j++) {
Nig[j] = 0.0;
}
- for (j = MIN(nt,nt-tmute+2*ts+shift-smooth),l=0; j < MIN(nt,nt-tmute+shift); j++,l++) {
+ for (j = MIN(nt,nt+1-tmute+shift+2*ts),l=0; j < MIN(nt,nt+1-tmute+shift+2*ts+smooth); j++,l++) {
Nig[j] *= costaper[smooth-l-1];
}
+ if (nt-tmute+shift+2*ts > nt) {
+ for (j = 0; j < MAX(0,-tmute+shift+2*ts); j++) {
+ Nig[j] = 0.0;
+ }
+ for (j = MAX(0,-tmute+shift+2*ts),l=0; j < MAX(0,-tmute+shift+2*ts+smooth); j++,l++) {
+ Nig[j] *= costaper[smooth-l-1];
+ }
+ }
}
- else if (above==-1) {
- imute = ixpos[i];
- tmute = mute[isyn*nxs+imute];
- ts = tsynW[isyn*nxs+imute];
- //ts = tsynW[isyn*nxs+ixpos[npos-1]];
- for (j = ts+tmute-shift,l=0; j < ts+tmute-shift+smooth; j++,l++) {
+ else if (above==4) { //Psi gate which is the inverse of the Theta gate (above=0)
+ if (tmute >= nt/2) continue;
+ for (j = MAX(0,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); j++,l++) {
+ Nig[j] *= costaper[smooth-l-1];
+ }
+ for (j = MIN(nt,nt+1-tmute+shift+2*ts-smooth),l=0; j < MIN(nt,nt+1-tmute+shift+2*ts); j++,l++) {
Nig[j] *= costaper[l];
}
- for (j = ts+tmute-shift+smooth; j < nt; j++) {
+ for (j = MIN(nt,nt+1-tmute+shift+2*ts); j < nt; j++) {
Nig[j] = 0.0;
}
+ if (nt-tmute+shift+2*ts >= nt) {
+ for (j = MAX(0,-tmute+shift+2*ts-smooth),l=0; j < MAX(0,-tmute+shift+2*ts); j++,l++) {
+ Nig[j] *= costaper[l];
+ }
+ for (j = MAX(0,-tmute+shift+2*ts); j < MAX(0,tmute-shift-smooth); j++) {
+ Nig[j] = 0.0;
+ }
+ }
+ else {
+ for (j = 0; j < MAX(0,tmute-shift-smooth); j++) {
+ Nig[j] = 0.0;
+ }
+ }
}
- else if (above==4) { //Psi gate which is the inverse of the Theta gate (above=0)
- 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++) {
+ else if (above==-4) { //Psi gate which is the inverse of the Theta gate (above=-1)
+ if (tmute >= nt/2) continue;
+ for (j = 0; j < MAX(0,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++) {
+ for (j = MAX(0,tmute-shift-smooth), l=0; j < MAX(0,tmute-shift); j++,l++) {
Nig[j] *= costaper[smooth-l-1];
}
- for (j = MIN(nt,nt-tmute+shift),l=0; j < MIN(nt,nt-tmute+shift+smooth); j++,l++) {
- Nig[j] *= costaper[l];
- }
- for (j = MIN(nt,nt-tmute+shift+smooth); j < nt; j++) {
- Nig[j] = 0.0;
- }
- }
-/* 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
- imute = ixpos[i];
- tmute = mute[isyn*nxs+imute];
- for (j = 0; j < tmute-shift-smooth; j++) {
- data[isyn*nxs*nt+i*nt+j] = 0.0;
- }
- for (j = tmute-shift-smooth,l=0; j < tmute-shift; j++,l++) {
- data[isyn*nxs*nt+i*nt+j] *= costaper[smooth-l-1];
+ if (tmute-shift-smooth < 0) {
+ for (j = MIN(nt,nt+1+tmute-shift-smooth),l=0; j < MIN(nt,nt+1+tmute-shift); j++,l++) {
+ Nig[j] *= costaper[smooth-l-1];
+ }
}
- for (j = tmute+shift,l=0; j < tmute+shift+smooth; j++,l++) {
- data[isyn*nxs*nt+i*nt+j] *= costaper[l];
+ for (j = MIN(nt,nt+1-(tmute+2*ts)+shift-smooth),l=0; j < MIN(nt,nt+1-(tmute+2*ts)+shift); j++,l++) {
+ Nig[j] *= costaper[l];
}
- for (j = tmute+shift+smooth; j < nt; j++) {
- data[isyn*nxs*nt+i*nt+j] = 0.0;
+ for (j = MIN(nt,nt+1-(tmute+2*ts)+shift),l=0; j < MIN(nt,nt+1+tmute-shift-smooth); j++,l++) {
+ Nig[j] = 0.0;
}
- }
+ }
- if (iter % 2 == 0) {
- for (j = 0; j < nt; j++) {
- data[isyn*nxs*nt+i*nt+j] = Nig[j];
- }
- }
- else { // reverse back in time
- j=0;
+ for (j = 0; j < nt; j++) {
data[isyn*nxs*nt+i*nt+j] = Nig[j];
- for (j = 1; j < nt; j++) {
- data[isyn*nxs*nt+i*nt+j] = Nig[nt-j];
- }
}
} /* end if ipos */
}
if (smooth) free(costaper);
- free(Nig);
+ free(Nig);
return;
}
void timeShift(float *data, int nsam, int nrec, float dt, float shift, float fmin, float fmax)
{
- int optn, iom, nfreq, ix, it;
- float deltom, om, tom, df, *trace, scl;
- complex *ctrace, ctmp;
+ int optn, iom, nfreq, ix, it;
+ float deltom, om, tom, df, *trace, scl;
+ complex *ctrace, ctmp;
- optn = optncr(nsam);
- nfreq = optn/2+1;
- df = 1.0/(optn*dt);
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+ df = 1.0/(optn*dt);
- ctrace = (complex *)malloc(nfreq*sizeof(complex));
- if (ctrace == NULL) verr("memory allocation error for ctrace");
+ ctrace = (complex *)malloc(nfreq*sizeof(complex));
+ if (ctrace == NULL) verr("memory allocation error for ctrace");
- trace = (float *)malloc(optn*sizeof(float));
- if (trace == NULL) verr("memory allocation error for rdata");
+ trace = (float *)malloc(optn*sizeof(float));
+ if (trace == NULL) verr("memory allocation error for rdata");
- deltom = 2.*M_PI*df;
+ deltom = 2.*M_PI*df;
scl = 1.0/(float)optn;
- for (ix = 0; ix < nrec; ix++) {
+ for (ix = 0; ix < nrec; ix++) {
for (it=0;it<nsam;it++) trace[it]=data[ix*nsam+it];
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 < nfreq ; iom++) {
- om = deltom*iom;
- tom = om*shift;
- ctmp = ctrace[iom];
- ctrace[iom].r = ctmp.r*cos(-tom) - ctmp.i*sin(-tom);
- ctrace[iom].i = ctmp.i*cos(-tom) + ctmp.r*sin(-tom);
- }
+ /* Forward time-frequency FFT */
+ rc1fft(&trace[0], &ctrace[0], optn, -1);
+ for (iom = 0 ; iom < nfreq ; iom++) {
+ om = deltom*iom;
+ tom = om*shift;
+ ctmp = ctrace[iom];
+ ctrace[iom].r = ctmp.r*cos(-tom) - ctmp.i*sin(-tom);
+ ctrace[iom].i = ctmp.i*cos(-tom) + ctmp.r*sin(-tom);
+ }
/* Inverse frequency-time FFT and scale result */
cr1fft(ctrace, trace, optn, 1);
for (it=0;it<nsam;it++) data[ix*nsam+it]=trace[it]*scl;
- }
+ }
free(ctrace);
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index 6189de2..6de7bbc 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -255,7 +255,6 @@ int main (int argc, char **argv)
/* compute time shift for tilted plane waves */
if (plane_wave==1) {
- itmin = nt;
/* compute time shift for shifted plane waves */
grad2rad = 17.453292e-3;
p = sin(src_angle*grad2rad)/src_velo;
@@ -263,10 +262,15 @@ int main (int argc, char **argv)
/* compute mute window for plane waves */
//for (i=0; i<nxs; i++) fprintf(stderr,"i=%d window=%d\n", i, muteW[i]);
+ itmin = nt;
for (i=0; i<nxs; i++) itmin = MIN (itmin, muteW[i]);
+/*
for (i=0; i<nxs; i++) tsynW[i] = muteW[i]-itmin;
+*/
+ for (i=0; i<nxs; i++) tsynW[i] = NINT(i*dxs*p/dt);
+ //for (i=0; i<nxs; i++) tsynW[i] = 0.0;
if (Nfoc!=1) verr("For plane-wave focusing only one function can be computed at the same time");
- //fprintf(stderr,"itmin=%d\n", itmin);
+ //fprintf(stderr,"itmin=%d tshift=%f =%d \n", itmin, tshift, NINT(tshift/dt));
//for (i=0; i<nxs; i++) fprintf(stderr,"i=%d window=%f\n", i, tsynW[i]*dt);
}
else { /* just fill with zero's */
@@ -487,13 +491,34 @@ int main (int argc, char **argv)
if (verbose >=2) vmess(" - iSyn %d: Ni at iteration %d has energy %e; relative to N0 %e", l, iter, sqrt(energyNi), sqrt(energyNi/energyN0[l]));
}
+// writeDataIter("bmute.su", Ni, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, iter+1);
/* apply mute window based on times of direct arrival (in muteW) */
+
if ( plane_wave==1 ) { /* use a-symmetric shift for plane waves with non-zero angles */
- applyMute_tshift(Ni, muteW, smooth, above, Nfoc, nxs, nts, ixpos, npos, shift, iter, tsynW);
+ applyMute_tshift(Ni, muteW, smooth, 0, Nfoc, nxs, nts, ixpos, npos, shift, iter, tsynW);
}
else {
applyMute(Ni, muteW, smooth, -above, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
}
+// writeDataIter("amute.su", Ni, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, iter+1);
+
+/*
+ // for testing time-windows with dipping plane waves
+ for (i = 0; i < npos; i++) {
+ for (j = 0; j < nts; j++) {
+ Ni[i*nts+j] = 1.0;
+ }
+ }
+ applyMute_tshift(Ni, muteW, smooth, 0, Nfoc, nxs, nts, ixpos, npos, shift, iter, tsynW);
+ writeDataIter("mute0.su", Ni, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, iter+1);
+ for (i = 0; i < npos; i++) {
+ for (j = 0; j < nts; j++) {
+ Ni[i*nts+j] = 1.0;
+ }
+ }
+ applyMute_tshift(Ni, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, iter, tsynW);
+ writeDataIter("mute4.su", Ni, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, iter+1);
+*/
if (iter % 2 == 0) { /* even iterations update: => f_1^-(t) */
for (l = 0; l < Nfoc; l++) {
@@ -582,12 +607,12 @@ int main (int argc, char **argv)
/* Apply mute with window for Gmin */
if ( plane_wave==1 ) {
/* for plane wave with angle tshift downward */
- if (verbose>1) vmess("Gmin planewave tshift=%f\n", tshift);
+ if (verbose>1) vmess("Gmin planewave tshift=%f", tshift);
+ applyMute_tshift(Gmin, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, 1, tsynW);
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;
+ //for (i=0; i<nxs; i++) tsynW[i] -= itmin;
+ //for (i=0; i<nxs; i++) tsynW[i] += itmin;
}
else {
applyMute(Gmin, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
@@ -617,7 +642,7 @@ int main (int argc, char **argv)
}
/* Apply mute with window for Gplus */
if ( plane_wave==1 ) {
- applyMute_tshift(Gplus, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, 1, tsynW);
+ applyMute_tshift(Gplus, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, 0, tsynW);
}
else {
applyMute(Gplus, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
diff --git a/marchenko/marchenko_primaries.c b/marchenko/marchenko_primaries.c
index 0e4843c..9dcc6a4 100644
--- a/marchenko/marchenko_primaries.c
+++ b/marchenko/marchenko_primaries.c
@@ -1001,16 +1001,17 @@ int main (int argc, char **argv)
iw0 = NINT((twplane[ix])/dt);
/* to apply muting window around epsilon below t=0 uncomment this part */
/* This (partly) removes the delta function at t=0 */
- // for (j = 0; j < MAX(0,iw0+shift-smooth); j++) {
- // uv[l*nxs*nts+i*nts+j] = 0.0;
- // }
- // 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];
+ for (j = 0; j < MAX(0,iw0+shift-smooth); j++) {
+ uv[l*nxs*nts+i*nts+j] = 0.0;
+ }
+ 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 and around 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);
for (j = MAX(0,iw0+shift); j < iw-isme; j++, k++) {
uv[l*nxs*nts+i*nts+j] = k1plus[l*nxs*nts+i*nts+j];
--
GitLab