Skip to content
Snippets Groups Projects
Commit 7d935a8c authored by JanThorbecke's avatar JanThorbecke
Browse files

adding plane-wave processing to basic marchenko scheme

parent 8af0254e
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
MDD/Makefile
+ 1
− 1
View file @ 7d935a8c
......@@ -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
......
MDD/deconvolve.c
+ 13
− 1
View file @ 7d935a8c
......@@ -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 */
......
MDD/mdd.c
+ 6
− 1
View file @ 7d935a8c
......@@ -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
......
marchenko/Makefile
+ 22
− 1
View file @ 7d935a8c
......@@ -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)
......
marchenko/applyMute_tshift.c 0 → 100644
+ 202
− 0
View file @ 7d935a8c
#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;
}
marchenko/demo/oneD/marchenko.scr
+ 2
− 2
View file @ 7d935a8c
......@@ -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
......
marchenko/marchenko.c
+ 65
− 70
View file @ 7d935a8c
......@@ -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);}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment