From 60d2a52d0d9ce8785eb2316d64d8af24c00be3b5 Mon Sep 17 00:00:00 2001
From: JanThorbecke <janth@xs4all.nl>
Date: Thu, 7 Jan 2021 10:32:45 +0100
Subject: [PATCH] negatine dipping plane-waves
---
fdelmodc/getParameters.c | 2 ++
marchenko/applyMute_tshift.c | 24 +++++++++++++++++++-----
marchenko/marchenko.c | 34 ++++++++++++++++++++++------------
marchenko/readShotData.c | 2 +-
marchenko/readTinvData.c | 4 ++--
marchenko/synthesis.c | 4 ++--
6 files changed, 48 insertions(+), 22 deletions(-)
diff --git a/fdelmodc/getParameters.c b/fdelmodc/getParameters.c
index 61983b0..0894e7c 100644
--- a/fdelmodc/getParameters.c
+++ b/fdelmodc/getParameters.c
@@ -891,6 +891,8 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
for (is=0; is<nsrc; is++) {
src->tbeg[is] = fabsf((nsrc-is-1)*dx*p);
}
+ //rec->delay+=NINT(fabsf((nsrc-1)*dx*p)/mod->dt);
+ //mod->nt += NINT(fabsf((nsrc-1)*dx*p)/mod->dt);
}
else {
for (is=0; is<nsrc; is++) {
diff --git a/marchenko/applyMute_tshift.c b/marchenko/applyMute_tshift.c
index 1863396..158f270 100644
--- a/marchenko/applyMute_tshift.c
+++ b/marchenko/applyMute_tshift.c
@@ -51,7 +51,6 @@ void applyMute_tshift( float *data, int *mute, int smooth, int above, int Nfoc,
if (above==4) above=-4;
tmute = tmute-2*ts;
}
- //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);
@@ -67,6 +66,17 @@ void applyMute_tshift( float *data, int *mute, int smooth, int above, int Nfoc,
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+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];
+ }
+ if (nt-(tmute+2*ts)+shift+smooth > nt) {
+ for (j = 0; j < MAX(0,-(tmute+2*ts)+shift); j++) {
+ Nig[j] = 0.0;
+ }
+ for (j = MAX(0,-(tmute+2*ts)+shift),l=0; j < MAX(0,-(tmute+2*ts)+shift+smooth); j++,l++) {
+ Nig[j] *= costaper[smooth-l-1];
+ }
+ }
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];
@@ -75,10 +85,6 @@ void applyMute_tshift( float *data, int *mute, int smooth, int above, int Nfoc,
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){
if (tmute >= nt/2) {
@@ -102,6 +108,14 @@ void applyMute_tshift( float *data, int *mute, int smooth, int above, int Nfoc,
Nig[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[l];
+ }
+ for (j = MIN(nt,nt+1+tmute-shift); j < nt; j++) {
+ Nig[j] = 0.0;
+ }
+ }
}
else if (above==4) { //Psi gate which is the inverse of the Theta gate (above=0)
if (tmute >= nt/2) continue;
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index 6de7bbc..1ed7503 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -134,7 +134,7 @@ int main (int argc, char **argv)
float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus;
float xmin, xmax, scale, tsq, Q, f0;
float *ixmask;
- float grad2rad, p, src_angle, src_velo, tshift;
+ float grad2rad, p, src_angle, src_velo, tshift, tneg;
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;
@@ -262,11 +262,19 @@ 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]);
+ //itmin = nt;
+ //for (i=0; i<nxs; i++) itmin = MIN (itmin, muteW[i]);
+
+ /* for negative angles tshift is negative and */
+ if (src_angle < 0.0) {
+ itmin = NINT(tshift/dt);
+ //for (i=0; i<nxs; i++) muteW[i] = MAX(0, muteW[i]-itmin);
+ for (i=0; i<nxs; i++) muteW[i] = muteW[i]-itmin;
+ timeShift(G_d, nts, nxs, dt, tshift, fmin, fmax);
+ }
/*
- 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");
@@ -274,7 +282,7 @@ int main (int argc, char **argv)
//for (i=0; i<nxs; i++) fprintf(stderr,"i=%d window=%f\n", i, tsynW[i]*dt);
}
else { /* just fill with zero's */
- itmin=0;
+ //itmin=0;
for (i=0; i<nxs*Nfoc; i++) {
tsynW[i] = 0;
}
@@ -491,7 +499,7 @@ 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);
+ //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 */
@@ -500,8 +508,7 @@ int main (int argc, char **argv)
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);
-
+ //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++) {
@@ -510,6 +517,7 @@ int main (int argc, char **argv)
}
}
applyMute_tshift(Ni, muteW, smooth, 0, Nfoc, nxs, nts, ixpos, npos, shift, iter, tsynW);
+ //applyMute(Ni, muteW, smooth, -above, Nfoc, nxs, nts, ixpos, npos, shift, 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++) {
@@ -517,6 +525,7 @@ int main (int argc, char **argv)
}
}
applyMute_tshift(Ni, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, iter, tsynW);
+ //applyMute(Ni, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
writeDataIter("mute4.su", Ni, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, iter+1);
*/
@@ -593,6 +602,8 @@ int main (int argc, char **argv)
xrcv, xsrc, xnx, fxse, fxsb, dxs, dxsrc, dx, ntfft, nw, nw_low, nw_high, mode,
reci, nshots, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv, ixmask, verbose);
+ writeDataIter("iRN.su", iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 1, iter+1);
+
/* compute upgoing Green's G^-,+ */
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
@@ -609,10 +620,9 @@ int main (int argc, char **argv)
/* for plane wave with angle tshift downward */
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;
- //for (i=0; i<nxs; i++) tsynW[i] += itmin;
+ if (src_angle > 0.0) {
+ timeShift(Gmin, nts, npos, dt, tshift, fmin, fmax);
+ }
}
else {
applyMute(Gmin, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
diff --git a/marchenko/readShotData.c b/marchenko/readShotData.c
index 80e43fc..ed1c0bc 100644
--- a/marchenko/readShotData.c
+++ b/marchenko/readShotData.c
@@ -131,7 +131,7 @@ int nx, int nxs, float fxsb, float dxs, int ntfft, int mode, float scale, float
}
if ((sx_shot != hdr.sx) || (fldr_shot != hdr.fldr) ) break;
}
- if (verbose>2) {
+ if (verbose>3) {
vmess("finished reading shot %d (%d) with %d traces",sx_shot,igath,itrace);
}
diff --git a/marchenko/readTinvData.c b/marchenko/readTinvData.c
index 11ef1a4..c0724b1 100644
--- a/marchenko/readTinvData.c
+++ b/marchenko/readTinvData.c
@@ -99,7 +99,7 @@ int readTinvData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx
}
if ((sx_shot != hdr.sx) || (fldr_shot != hdr.fldr) ) break;
}
- if (verbose>2) {
+ if (verbose>3) {
fprintf(stderr,"finished reading shot %d (%d) with %d traces\n",sx_shot,isyn,itrace);
//disp_fileinfo(filename, nt, xnx[isyn], hdr.f1, xrcv[isyn*nxm], d1, d2, &hdr);
}
@@ -144,7 +144,7 @@ int readTinvData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx
}
}
maxval[isyn*nx+imax] = jmax;
- if (verbose >= 3) vmess("Mute max at src-trace %d is sample %d", imax, maxval[imax]);
+ if (verbose >= 2) vmess("Mute max at src-trace %d is sample %d", imax, maxval[imax]);
/* search forward in trace direction from maximum in file */
for (i = imax+1; i < nx1; i++) {
diff --git a/marchenko/synthesis.c b/marchenko/synthesis.c
index 39a3150..8d7a163 100644
--- a/marchenko/synthesis.c
+++ b/marchenko/synthesis.c
@@ -180,7 +180,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
} /* end of Nfoc loop */
- if (verbose>4) vmess("*** Shot gather %d processed ***", k);
+ if (verbose>3) vmess("*** Shot gather %d processed ***", k);
} /* end of nparallel shots (k) loop */
free(sum);
@@ -375,7 +375,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
} /* end of Nfoc loop */
- if (verbose>4) vmess("*** Shot gather %d processed ***", k);
+ if (verbose>3) vmess("*** Shot gather %d processed ***", k);
} /* end of nparallel shots (k) loop */
free(sum);
--
GitLab