From 7d935a8c33790fc002cf5f2bd8e8fd29d8e9497e Mon Sep 17 00:00:00 2001
From: JanThorbecke <janth@xs4all.nl>
Date: Thu, 7 May 2020 08:03:47 +0200
Subject: [PATCH] adding plane-wave processing to basic marchenko scheme
---
MDD/Makefile | 2 +-
MDD/deconvolve.c | 14 ++-
MDD/mdd.c | 7 +-
marchenko/Makefile | 23 +++-
marchenko/applyMute_tshift.c | 202 ++++++++++++++++++++++++++++++
marchenko/demo/oneD/marchenko.scr | 4 +-
marchenko/marchenko.c | 135 ++++++++++----------
7 files changed, 311 insertions(+), 76 deletions(-)
create mode 100644 marchenko/applyMute_tshift.c
diff --git a/MDD/Makefile b/MDD/Makefile
index 561b309..523a1f8 100644
--- a/MDD/Makefile
+++ b/MDD/Makefile
@@ -12,7 +12,7 @@ ALLINC = -I.
#General BLAS library
#LIBS += $(BLAS)
-#CFLAGS += -I$(MKLROOT)/include
+CFLAGS += -I$(MKLROOT)/include
#LIBS += -lblas -llapack -L$L -lgenfft $(LIBSM) -lc -lm
all: mdd
diff --git a/MDD/deconvolve.c b/MDD/deconvolve.c
index 147fc7b..13dc809 100644
--- a/MDD/deconvolve.c
+++ b/MDD/deconvolve.c
@@ -129,11 +129,23 @@ private(iw, iwnA, iwnB, iwAB, iwBB)
/* cblas_cgemm(CblasRowMajor,CblasNoTrans, CblasConjTrans, NA, NB, nshots, &alpha.r,
&cA[iwnA].r, NA,
&cB[iwnB].r, NB, &beta.r,
- &cC[iwAB].r, NC); */
+ ocC[iwAB].r, NC); */
+/*
+ for (i=0; i<nshots; i++) {
+ for (j=0; j<nshots; j++) {
+ for (k=0; k<nstationB; k++) {
+ cC[iwAB+j*nshots+i].r += cA[iwnA+k*nshots+i].r*cB[iwnB+j*nshots+k].r - cA[iwnA+k*nshots+i].i*cB[iwnB+j*nshots+k].i;
+ cC[iwAB+j*nshots+i].i += cA[iwnA+k*nshots+i].r*cB[iwnB+j*nshots+k].i + cA[iwnA+k*nshots+i].i*cB[iwnB+j*nshots+k].r;
+ }
+ }
+ }
+*/
+
cgemm_(transa, transb, &NA, &NB, &nshots, &alpha.r,
&cA[iwnA].r, &NA,
&cB[iwnB].r, &NB, &beta.r,
&cC[iwAB].r, &NC);
+
// memcpy(&cC[iwAB].r, &cB[iwnA].r, sizeof(float)*2*nstationA*nshots);
}
else if (mdd==1) { /* Multi Dimensional deconvolution */
diff --git a/MDD/mdd.c b/MDD/mdd.c
index 5da9828..92d1976 100644
--- a/MDD/mdd.c
+++ b/MDD/mdd.c
@@ -170,7 +170,12 @@ int main (int argc, char **argv)
if (!getparfloat("cjB", &cjB)) cjB = 1.;
#ifdef _OPENMP
- npes = atoi(getenv("OMP_NUM_THREADS"));
+ npes = omp_get_max_threads();
+ /* parallelisation is over number of shot positions (nshots) */
+ if (npes == 0) {
+ vmess("Number of OpenMP threads set to %d (was %d)", 1, npes);
+ omp_set_num_threads(1);
+ }
assert(npes != 0);
if (verbose) fprintf(stderr,"Number of OpenMP thread's is %d\n", npes);
#else
diff --git a/marchenko/Makefile b/marchenko/Makefile
index cd1c541..61b0b11 100644
--- a/marchenko/Makefile
+++ b/marchenko/Makefile
@@ -7,7 +7,7 @@ include ../Make_include
#OPTC += -g -traceback #-check-pointers=rw -check-pointers-undimensioned
-ALL: fmute marchenko marchenko_primaries
+ALL: fmute marchenko marchenko_primaries marchenko_tshift
SRCJ = fmute.c \
getFileInfo.c \
@@ -27,6 +27,22 @@ SRCH = marchenko.c \
readShotData.c \
readTinvData.c \
applyMute.c \
+ applyMute_tshift.c \
+ writeData.c \
+ writeDataIter.c \
+ wallclock_time.c \
+ name_ext.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ getpars.c
+
+SRCT = marchenko_tshift.c \
+ getFileInfo.c \
+ readData.c \
+ readShotData.c \
+ readTinvData.c \
+ applyMute_tshift.c \
writeData.c \
writeDataIter.c \
wallclock_time.c \
@@ -62,6 +78,11 @@ OBJH = $(SRCH:%.c=%.o)
marchenko: $(OBJH)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko $(OBJH) $(LIBS)
+OBJT = $(SRCT:%.c=%.o)
+
+marchenko_tshift: $(OBJT)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko_tshift $(OBJT) $(LIBS)
+
OBJP = $(SRCP:%.c=%.o)
marchenko_primaries: $(OBJP)
diff --git a/marchenko/applyMute_tshift.c b/marchenko/applyMute_tshift.c
new file mode 100644
index 0000000..9a15951
--- /dev/null
+++ b/marchenko/applyMute_tshift.c
@@ -0,0 +1,202 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <assert.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+
+void applyMute_tshift( float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *ixpos, int npos, int shift, int iter, int *tsynW)
+{
+ int i, j, l, isyn;
+ float *costaper, scl, *Nig;
+ int imute, tmute, ts;
+
+ 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));
+ }
+ }
+
+ Nig = (float *)malloc(nt*sizeof(float));
+
+ for (isyn = 0; isyn < Nfoc; isyn++) {
+ for (i = 0; i < npos; i++) {
+ if (iter % 2 == 0) {
+ for (j = 0; j < nt; j++) {
+ Nig[j] = data[isyn*nxs*nt+i*nt+j];
+ }
+ }
+ 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];
+ }
+ }
+ 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;
+ }
+ for (j = MAX(0,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); 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++) {
+ Nig[j] *= costaper[l];
+ }
+ for (j = MAX(0,tmute-shift+smooth+1); j < MIN(nt,nt+1-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++) {
+ 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++) {
+ Nig[j] *= costaper[l];
+ }
+ for (j = ts+tmute-shift+smooth; j < nt; 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 = -2*ts+tmute-shift-smooth,l=0; j < -2*ts+tmute-shift; j++,l++) {
+ Nig[j] *= costaper[smooth-l-1];
+ }
+ for (j = 0; j < -2*ts+tmute-shift-smooth-1; j++) {
+ Nig[j] = 0.0;
+ }
+ for (j = nt+1-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
+ 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];
+ }
+ for (j = tmute+shift,l=0; j < tmute+shift+smooth; j++,l++) {
+ data[isyn*nxs*nt+i*nt+j] *= costaper[l];
+ }
+ for (j = tmute+shift+smooth; j < nt; j++) {
+ data[isyn*nxs*nt+i*nt+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;
+ 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);
+
+ return;
+}
+
+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;
+ complex *ctrace, ctmp;
+
+ 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");
+
+ trace = (float *)malloc(optn*sizeof(float));
+ 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++) {
+ 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 < 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++) {
+ 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);
+ free(trace);
+
+ return;
+}
diff --git a/marchenko/demo/oneD/marchenko.scr b/marchenko/demo/oneD/marchenko.scr
index 071443f..6975bb3 100755
--- a/marchenko/demo/oneD/marchenko.scr
+++ b/marchenko/demo/oneD/marchenko.scr
@@ -8,8 +8,8 @@ export MKL_NUM_THREADS=1
fmute file_shot=iniFocus_rp.su file_out=p0plus.su above=-1 shift=-8 verbose=1 check=0 hw=8
#apply the Marchenko algorithm
-marchenko file_shot=shotsdx5_rp.su file_tinv=p0plus.su nshots=901 verbose=2 \
- tap=0 niter=8 hw=8 shift=12 smooth=3 \
+../../marchenko file_shot=shotsdx5_rp.su file_tinv=p0plus.su nshots=901 verbose=2 \
+ tap=0 niter=8 hw=8 shift=12 smooth=3 file_iter=iter.su \
file_green=pgreen.su file_gplus=Gplus0.su file_gmin=Gmin0.su \
file_f1plus=f1plus0.su file_f1min=f1min0.su file_f2=f2.su
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index 673fe02..0738c46 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -40,6 +40,9 @@ int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, floa
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, int *muteW);
+void applyMute_tshift( float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *ixpos, int npos, int shift, int iter, int *tsynW);
+void timeShift(float *data, int nsam, int nrec, float dt, float shift, float fmin, float fmax);
+
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);
@@ -83,9 +86,7 @@ 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",
+" plane_wave=0 ............. enable plane-wave illumination function",
" REFLECTION RESPONSE CORRECTION ",
" tsq=0.0 .................. scale factor n for t^n for true amplitude recovery",
" Q=0.0 .......,............ Q correction factor",
@@ -121,17 +122,16 @@ int main (int argc, char **argv)
int size, n1, n2, ntap, tap, di, ntraces, pad, rotate;
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
int reci, countmin, mode, n2out, verbose, ntfft;
- int iter, niter, tracf, *muteW, *tsynW;
+ int iter, niter, tracf, *muteW, *tsynW, itmin;
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, *Nig, *trace, *Gmin, *Gplus;
+ float *f1plus, *f1min, *iRN, *Ni, *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;
@@ -156,10 +156,7 @@ int main (int argc, char **argv)
if (!getparint("verbose", &verbose)) verbose = 0;
if (file_tinv == NULL && file_shot == NULL)
verr("file_tinv and file_shot cannot be both input pipe");
- if (!getparstring("file_green", &file_green)) {
- if (verbose) vwarn("parameter file_green not found, assume pipe");
- file_green = NULL;
- }
+ if (!getparstring("file_green", &file_green)) file_green = NULL;
if (!getparfloat("fmin", &fmin)) fmin = 0.0;
if (!getparfloat("fmax", &fmax)) fmax = 70.0;
if (!getparint("reci", &reci)) reci = 0;
@@ -179,8 +176,6 @@ int main (int argc, char **argv)
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");
@@ -220,7 +215,6 @@ 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
@@ -255,22 +249,14 @@ int main (int argc, char **argv)
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 (plane_wave==1) {
+ itmin = nt;
+ for (i=0; i<nxs; i++) itmin = MIN (itmin, muteW[i]);
+ for (i=0; i<nxs; i++) tsynW[i] = muteW[i]-itmin;
if (Nfoc!=1) verr("For plane-wave focusing only one function can be computed at the same time");
}
else { /* just fill with zero's */
+ itmin=0;
for (i=0; i<nxs*Nfoc; i++) {
tsynW[i] = 0;
}
@@ -465,7 +451,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, 1, iter);
+ writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, iter+1);
}
/* N_k(x,t) = -N_(k-1)(x,-t) */
/* p0^-(x,t) += iRN = (R * T_d^inv)(t) */
@@ -482,49 +468,30 @@ 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[Nfoc] = energyNi;
- if (verbose >=2) vmess(" - iSyn %d: Ni at iteration %d has energy %e; relative to N0 %e", l, iter, sqrt(energyNi),
-sqrt(energyNi/energyN0[Nfoc]));
+ if (verbose >=2) vmess(" - iSyn %d: Ni at iteration %d has energy %e; relative to N0 %e", l, iter, sqrt(energyNi), sqrt(energyNi/energyN0[Nfoc]));
}
/* apply mute window based on times of direct arrival (in muteW) */
- 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 ( 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);
+ }
+ else {
+ applyMute(Ni, muteW, smooth, above, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
+ }
if (iter % 2 == 0) { /* even iterations update: => f_1^-(t) */
- 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];
- }
- }
- }
- if (above==-2) applyMute(f1min, muteW, smooth, 0, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
- }
- 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];
- }
+ 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 (above==-2) applyMute(f1min, muteW, smooth, 0, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
- }
+ }
+ if (above==-2) applyMute(f1min, muteW, smooth, 0, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
}
else {/* odd iterations update: => f_1^+(t) */
for (l = 0; l < Nfoc; l++) {
@@ -537,6 +504,7 @@ sqrt(energyNi/energyN0[Nfoc]));
}
}
}
+ //writeDataIter("Ni.su", Ni, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, iter+1);
/* update f2 */
for (l = 0; l < Nfoc; l++) {
@@ -557,7 +525,6 @@ sqrt(energyNi/energyN0[Nfoc]));
} /* end of iterations */
free(Ni);
- free(Nig);
free(G_d);
/* compute full Green's function G = int R * f2(t) + f2(-t) = Pplus + Pmin */
@@ -579,7 +546,7 @@ sqrt(energyNi/energyN0[Nfoc]));
/* compute upgoing Green's function G^+,- */
if (file_gmin != NULL) {
- Gmin = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ Gmin = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
/* use f1+ as operator on R in frequency domain */
mode=1;
@@ -599,7 +566,25 @@ sqrt(energyNi/energyN0[Nfoc]));
}
}
/* Apply mute with window for Gmin */
- applyMute(Gmin, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
+ if ( plane_wave==1 ) {
+ applyMute_tshift(Gmin, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, 0, tsynW);
+ /* for plane wave with angle shift itmin downward */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ memcpy(&trace[0],&Gmin[l*nxs*nts+i*nts],nts*sizeof(float));
+ for (j = 0; j < itmin; j++) {
+ Gmin[l*nxs*nts+i*nts+j] = 0.0;
+ }
+ for (j = 0; j < nts-itmin; j++) {
+ Gmin[l*nxs*nts+i*nts+j+itmin] = trace[j];
+ }
+ }
+ }
+ //timeShift(Gmin, nts, npos, dt, itmin*dt, fmin, fmax);
+ }
+ else {
+ applyMute(Gmin, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
+ }
} /* end if Gmin */
/* compute downgoing Green's function G^+,+ */
@@ -624,7 +609,12 @@ sqrt(energyNi/energyN0[Nfoc]));
}
}
/* Apply mute with window for Gplus */
- applyMute(Gplus, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
+ if ( plane_wave==1 ) {
+ applyMute_tshift(Gplus, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, 1, tsynW);
+ }
+ else {
+ applyMute(Gplus, muteW, smooth, 4, Nfoc, nxs, nts, ixpos, npos, shift, tsynW);
+ }
} /* end if Gplus */
free(muteW);
@@ -643,8 +633,10 @@ sqrt(energyNi/energyN0[Nfoc]));
/*================ write output files ================*/
- fp_out = fopen(file_green, "w+");
- if (fp_out==NULL) verr("error on creating output file %s", file_green);
+ if (file_green != NULL) {
+ fp_out = fopen(file_green, "w+");
+ if (fp_out==NULL) verr("error on creating output file %s", file_green);
+ }
if (file_gmin != NULL) {
fp_gmin = fopen(file_gmin, "w+");
if (fp_gmin==NULL) verr("error on creating output file %s", file_gmin);
@@ -686,8 +678,10 @@ sqrt(energyNi/energyN0[Nfoc]));
hdrs_out[i].f1 = f1;
}
- ret = writeData(fp_out, (float *)&green[l*size], hdrs_out, n1, n2);
- if (ret < 0 ) verr("error on writing output file.");
+ if (file_green != NULL) {
+ ret = writeData(fp_out, (float *)&green[l*size], hdrs_out, n1, n2);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
if (file_gmin != NULL) {
ret = writeData(fp_gmin, (float *)&Gmin[l*size], hdrs_out, n1, n2);
@@ -740,7 +734,8 @@ sqrt(energyNi/energyN0[Nfoc]));
if (ret < 0 ) verr("error on writing output file.");
}
}
- ret = fclose(fp_out);
+ ret=0;
+ if (file_green != NULL) {ret += fclose(fp_out);}
if (file_gplus != NULL) {ret += fclose(fp_gplus);}
if (file_gmin != NULL) {ret += fclose(fp_gmin);}
if (file_f2 != NULL) {ret += fclose(fp_f2);}
--
GitLab