diff --git a/marchenko_applications/JespersRayTracer.c b/marchenko_applications/JespersRayTracer.c
index bd4f9da25c25868c7655609b6859c943bb25f212..ba6ae35845c1e6260999a0799339a8cc09dd66b5 100644
--- a/marchenko_applications/JespersRayTracer.c
+++ b/marchenko_applications/JespersRayTracer.c
@@ -123,7 +123,7 @@ int getWaveParameter(float *slowness, icoord size, float dgrid, fcoord s, fcoord
if (ray.useT2 != 0)
T2 += getdT2(x, z, so, angle, ds, nRayTmp, s, r, rayReference3D, slowness, size);
- if (ray.geomspread != 0) {
+ /*if (ray.geomspread != 0) {
if (so <= 0) {
dQdPhi = 0;
}
@@ -131,7 +131,7 @@ int getWaveParameter(float *slowness, icoord size, float dgrid, fcoord s, fcoord
greentmp = greenIntP(lengthRefRay, so, lengthRefRay, slowness, size, nRayTmp, r, s);
dQdPhi += greentmp*secondDerivativeU1(slowness, size, x, z, angle, r, s)*ds/so;
}
- }
+ }*/
}
if (ray.useT2)
diff --git a/marchenko_applications/Makefile b/marchenko_applications/Makefile
index 4130dd06511833931d999436e2c91bc754e1b3e0..7fbfacafae0f940c08260efef7a6fe3c1b43940d 100644
--- a/marchenko_applications/Makefile
+++ b/marchenko_applications/Makefile
@@ -7,7 +7,7 @@ LIBS += -L$L -lgenfft -lm $(LIBSM)
#ALL: fmute marchenko marchenko2
-ALL: fmute marchenko
+ALL: fmute marchenko_app combine
SRCJ = fmute.c \
getFileInfo.c \
@@ -51,7 +51,19 @@ SRCH = marchenko.c \
gaussGen.c \
iterations.c \
imaging.c \
- threadAffinity.c
+ threadAffinity.c \
+ makeWindow.c \
+ homogeneousg.c
+
+SRCC = combine.c \
+ getFileInfo.c \
+ writeData.c \
+ wallclock_time.c \
+ getpars.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ readSnapData.c
OBJJ = $(SRCJ:%.c=%.o)
@@ -60,20 +72,26 @@ fmute: $(OBJJ)
OBJH = $(SRCH:%.c=%.o)
-marchenko: $(OBJH) raytime.h
- $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko $(OBJH) $(LIBS)
+marchenko_app: $(OBJH) raytime.h
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko_app $(OBJH) $(LIBS)
+
+OBJC = $(SRCC:%.c=%.o)
+
+combine: $(OBJC)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o combine $(OBJC) $(LIBS)
-install: fmute marchenko
+install: fmute marchenko_app combine
cp fmute $B
- cp marchenko $B
+ cp marchenko_app $B
+ cp combine $B
# cp marchenko2 $B
clean:
- rm -f core fmute $(OBJJ) marchenko $(OBJH)
+ rm -f core fmute $(OBJJ) marchenko_app $(OBJH) combine $(OBJC)
realclean: clean
- rm -f $B/fmute $B/marchenko
+ rm -f $B/fmute $B/marchenko_app $B/combine
diff --git a/marchenko_applications/marchenko.c b/marchenko_applications/marchenko.c
index 1b26d99fb983051a93897fd26f74329e6604605a..7a975bf02705a4b8778709fa7d22bbd616943ff1 100644
--- a/marchenko_applications/marchenko.c
+++ b/marchenko_applications/marchenko.c
@@ -38,7 +38,7 @@ typedef struct _complexStruct { // complex number
int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int ngath, int nx, int nxm, int ntfft, int mode, float weight, float tsq, float Q, float f0, int verbose);
int readSnapData(char *filename, float *data, segy *hdrs, int nsnaps, int nx, int nz, int sx, int ex, int sz, int ez);
//int readTinvData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nsyn, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
-int readTinvData(char *filename, WavePar WP, char *file_ray, char *file_amp, float dt, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nsyn, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
+int readTinvData(char *filename, float dt, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nsyn, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nsyn, float *xsyn, float *zsyn, int iter);
void name_ext(char *filename, char *extension);
int Cost(float *f1p, float *f1d, float *Gm, float *Gm0, double *J, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn);
@@ -49,10 +49,12 @@ int readData(FILE *fp, float *data, segy *hdrs, int n1);
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
double wallclock_time(void);
-void iterations (complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int verbose);
-void imaging (float *Image, float *Gmin, complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *synpos, int verbose);
+int makeWindow(WavePar WP, char *file_ray, char *file_amp, float dt, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nsyn, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
+void iterations (complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *first, int verbose);
+void imaging (float *Image, WavePar WP, complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *synpos, int *first, int verbose);
+void homogeneousg(float *HomG, float *green, complex *Refl, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, float *pmin, float *f1min, float *f1plus, float *f2p, float *G_d, int *muteW, int smooth, int shift, int above, int pad, int nt0, int *synpos, int verbose);
-void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, double *tfft, int verbose);
+void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, double *tfft, int *first, int verbose);
void synthesisPosistions(int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int reci, int nshots, int *ixpossyn, int *npossyn, int verbose);
@@ -210,18 +212,19 @@ int main (int argc, char **argv)
int size, n1, n2, ntap, tap, di, ntraces, nb, ib;
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn, *synpos;
int reci, mode, ixa, ixb, n2out, verbose, ntfft;
- int iter, niter, tracf, *muteW, pad, nt0, ampest;
- int hw, smooth, above, shift, *ixpossyn, npossyn, ix;
+ int iter, niter, tracf, *muteW, pad, nt0, ampest, *hmuteW, *hxnxsyn;
+ int hw, smooth, above, shift, *ixpossyn, npossyn, ix, first=1;
float fmin, fmax, *tapersh, *tapersy, fxf, dxf, fxs2, *xsrc, *xrcv, *zsyn, *zsrc, *xrcvsyn;
+ float *hzsyn, *hxsyn, *hxrcvsyn, *hG_d, xloc, zloc, *HomG;
double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi, *J;
float d1, d2, f1, f2, fxs, ft, fx, *xsyn, dxsrc, Q, f0, *Costdet;
float *green, *f2p, *pmin, *G_d, dt, dx, dxs, scl, mem, *Image, *Image2;
float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus, *Gm0;
float xmin, xmax, weight, tsq, *Gd, *amp, bstart, bend, db, *bdet, bp, b, bmin;
- complex *Refl, *Fop;
- char *file_tinv, *file_shot, *file_green, *file_iter, *file_wav, *file_ray, *file_amp, *file_img;
- char *file_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_pmin, *wavtype, *wavtype2;
- segy *hdrs_out, *hdrs_out2;
+ complex *Refl, *Fop, *cshot;
+ char *file_tinv, *file_shot, *file_green, *file_iter, *file_wav, *file_ray, *file_amp, *file_img, *file_cp;
+ char *file_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_pmin, *wavtype, *wavtype2, *file_homg;
+ segy *hdrs_im, *hdrs_homg;
WavePar WP,WPs;
modPar mod;
recPar rec;
@@ -235,7 +238,8 @@ int main (int argc, char **argv)
tsyn = tread = tfft = tcopy = 0.0;
t0 = wallclock_time();
- if (!getparstring("file_img", &file_img)) file_img = "out.su";
+ if (!getparstring("file_img", &file_img)) file_img = "img.su";
+ if (!getparstring("file_homg", &file_homg)) file_homg = NULL;
if (!getparstring("file_shot", &file_shot)) file_shot = NULL;
if (!getparstring("file_tinv", &file_tinv)) file_tinv = NULL;
if (!getparstring("file_f1plus", &file_f1plus)) file_f1plus = NULL;
@@ -249,6 +253,7 @@ int main (int argc, char **argv)
if (!getparstring("file_wav", &file_wav)) file_wav=NULL;
if (!getparstring("file_ray", &file_ray)) file_ray=NULL;
if (!getparstring("file_amp", &file_amp)) file_amp=NULL;
+ if (!getparstring("file_cp", &file_cp)) file_cp = NULL;
if (!getparint("verbose", &verbose)) verbose = 0;
if (file_tinv == NULL && file_shot == NULL)
verr("file_tinv and file_shot cannot be both input pipe");
@@ -342,8 +347,9 @@ int main (int argc, char **argv)
xmax = mod.x0+rec.x[rec.n-1]*mod.dx;
scl = 0.0010;
ntraces = n2*ngath;
- WPs.wav = 1;
WP.wav = 1;
+ WP.xloc = -123456.0;
+ WP.zloc = -123456.0;
synpos = (int *)calloc(ngath,sizeof(int));
for (l=0; l<shot.nz; l++) {
for (j=0; j<shot.nx; j++) {
@@ -417,9 +423,15 @@ int main (int argc, char **argv)
WP.nt = ntfft;
WP.dt = dt;
- mode=-1; /* apply complex conjugate to read in data */
- readTinvData(file_tinv, WPs, file_ray, file_amp, dt, xrcvsyn, xsyn, zsyn, xnxsyn,
- Nsyn, nxs, ntfft, mode, muteW, G_d, hw, verbose);
+ if (file_ray!=NULL || file_cp!=NULL) {
+ makeWindow(WP, file_ray, file_amp, dt, xrcvsyn, xsyn, zsyn, xnxsyn,
+ Nsyn, nxs, ntfft, mode, muteW, G_d, hw, verbose);
+ }
+ else {
+ mode=-1; /* apply complex conjugate to read in data */
+ readTinvData(file_tinv, dt, xrcvsyn, xsyn, zsyn, xnxsyn,
+ Nsyn, nxs, ntfft, mode, muteW, G_d, hw, verbose);
+ }
/* reading data added zero's to the number of time samples to be the same as ntfft */
nts = ntfft;
@@ -554,6 +566,71 @@ int main (int argc, char **argv)
}
//memcpy(Ni, G_d, Nsyn*nxs*ntfft*sizeof(float));
+
+ if (file_homg!=NULL) {
+ hG_d = (float *)calloc(nxs*ntfft,sizeof(float));
+ hmuteW = (int *)calloc(nxs,sizeof(int));
+ hxrcvsyn = (float *)calloc(nxs,sizeof(float));
+ hxsyn = (float *)calloc(1,sizeof(float));
+ hzsyn = (float *)calloc(1,sizeof(float));
+ hxnxsyn = (int *)calloc(1,sizeof(int));
+ cshot = (complex *)calloc(nxs*nfreq,sizeof(complex));
+
+ if(!getparfloat("xloc", &WPs.xloc)) WPs.xloc = 0;
+ if(!getparfloat("zloc", &WPs.zloc)) WPs.zloc = 0;
+ xloc = WPs.xloc;
+ zloc = WPs.zloc;
+ ngath = 1;
+
+ makeWindow(WPs, file_ray, file_amp, dt, hxrcvsyn, hxsyn, hzsyn, hxnxsyn, ngath, nxs, ntfft, mode, hmuteW, hG_d, hw, verbose);
+
+ WPs.xloc = -123456.0;
+ WPs.zloc = -123456.0;
+
+ if (tap == 1 || tap == 3) {
+ if (verbose) vmess("Taper for operator applied ntap=%d", ntap);
+ for (i = 0; i < nxs; i++) {
+ for (j = 0; j < nts; j++) {
+ hG_d[i*nts+j] *= tapersy[i];
+ }
+ }
+ }
+
+ ngath = omp_get_max_threads();
+
+ synthesisPosistions(nx, nt, nxs, nts, dt, hxsyn, 1, xrcv, xsrc, fxs2, fxs,
+ dxs, dxsrc, dx, ixa, ixb, reci, nshots, ixpossyn, &npossyn, verbose);
+
+ iterations(Refl,nx,nt,nxs,nts,dt,hxsyn,1,xrcv,xsrc,fxs2,fxs,dxs,dxsrc,dx,ixa,ixb,
+ ntfft,nw,nw_low,nw_high,mode,reci,nshots,ixpossyn,npossyn,pmin,f1min,f1plus,
+ f2p,hG_d,hmuteW,smooth,shift,above,pad,nt0,&first,verbose);
+
+ /* compute full Green's function G = int R * f2(t) + f2(-t) = Pplus + Pmin */
+ for (i = 0; i < npossyn; i++) {
+ j = 0;
+ /* set green to zero if mute-window exceeds nt/2 */
+ if (hmuteW[ixpossyn[i]] >= nts/2) {
+ memset(&green[i*nts],0, sizeof(float)*nt);
+ continue;
+ }
+ green[i*nts+j] = f2p[i*nts+j] + pmin[i*nts+j];
+ for (j = 1; j < nts; j++) {
+ green[i*nts+j] = f2p[i*nts+nts-j] + pmin[i*nts+j];
+ }
+ }
+
+ applyMute(green, hmuteW, smooth, 4, 1, nxs, nts, ixpossyn, npossyn, shift, pad, nt0);
+
+ omp_set_num_threads(ngath);
+
+ /* Transform the green position to the frequency domain */
+ /*for (i = 0; i < npossyn; i++) {
+ rc1fft(&green[i*nts],&cshot[i*nfreq],ntfft,-1);
+ }*/
+ //free(hG_d);free(hmuteW);free(hxrcvsyn);
+ free(hmuteW);free(hxrcvsyn);
+ free(hxsyn);free(hzsyn);free(hxnxsyn);free(cshot);
+ }
/* dry-run of synthesis to get all x-positions calcalated by the integration */
synthesisPosistions(nx, nt, nxs, nts, dt, xsyn, Nsyn, xrcv, xsrc, fxs2, fxs,
@@ -562,119 +639,117 @@ int main (int argc, char **argv)
vmess("synthesisPosistions: nshots=%d npossyn=%d", nshots, npossyn);
}
-/*================ set variables for output data ================*/
-
- n1 = nts; n2 = n2out;
- f1 = ft; f2 = fxs+dxs*ixpossyn[0];
- d1 = dt;
- if (reci == 0) d2 = dxsrc;
- else if (reci == 1) d2 = dxs;
- else if (reci == 2) d2 = dx;
- hdrs_out2 = (segy *) calloc(n2,sizeof(segy));
- hdrs_out = (segy *) calloc(shot.nx,sizeof(segy));
- if (hdrs_out == NULL) verr("allocation for hdrs_out");
- size = nxs*nts;
-
- for (i = 0; i < n2; i++) {
- hdrs_out2[i].ns = n1;
- hdrs_out2[i].trid = 1;
- hdrs_out2[i].dt = dt*1000000;
- hdrs_out2[i].f1 = f1;
- hdrs_out2[i].f2 = f2;
- hdrs_out2[i].d1 = d1;
- hdrs_out2[i].d2 = d2;
- hdrs_out2[i].trwf = n2out;
- hdrs_out2[i].scalco = -1000;
- hdrs_out2[i].gx = NINT(1000*(f2+i*d2));
- hdrs_out2[i].scalel = -1000;
- hdrs_out2[i].tracl = i+1;
- }
t1 = wallclock_time();
tread = t1-t0;
iterations(Refl,nx,nt,nxs,nts,dt,xsyn,Nsyn,xrcv,xsrc,fxs2,fxs,dxs,dxsrc,dx,ixa,ixb,
ntfft,nw,nw_low,nw_high,mode,reci,nshots,ixpossyn,npossyn,pmin,f1min,f1plus,
- f2p,G_d,muteW,smooth,shift,above,pad,nt0,verbose);
-
- Image = (float *)malloc(Nsyn*sizeof(float));
- Gmin = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
-
- imaging(Image,Gmin,Refl,nx,nt,nxs,nts,dt,xsyn,Nsyn,xrcv,xsrc,fxs2,fxs,dxs,dxsrc,dx,ixa,ixb,
- ntfft,nw,nw_low,nw_high,mode,reci,nshots,ixpossyn,npossyn,pmin,f1min,f1plus,
- f2p,G_d,muteW,smooth,shift,above,pad,nt0,synpos,verbose);
-
-/*============= write output files ================*/
-
- fp_out = fopen(file_img, "w+");
-
- for (i = 0; i < shot.nx; i++) {
- hdrs_out[i].fldr = 1;
- hdrs_out[i].tracl = 1;
- hdrs_out[i].tracf = i+1;
- hdrs_out[i].scalco = -1000;
- hdrs_out[i].scalel = -1000;
- hdrs_out[i].sdepth = 0;
- hdrs_out[i].trid = 1;
- hdrs_out[i].ns = shot.nz;
- hdrs_out[i].trwf = shot.nx;
- hdrs_out[i].ntr = hdrs_out[i].fldr*hdrs_out[i].trwf;
- hdrs_out[i].f1 = zsyn[0];
- hdrs_out[i].f2 = xsyn[0];
- hdrs_out[i].dt = dt*(1E6);
- hdrs_out[i].d1 = (float)zsyn[shot.nx]-zsyn[0];
- hdrs_out[i].d2 = (float)xsyn[1]-xsyn[0];
- hdrs_out[i].sx = (int)roundf(xsyn[0] + (i*hdrs_out[i].d2)*1000.0);
- hdrs_out[i].gx = (int)roundf(xsyn[0] + (i*hdrs_out[i].d2)*1000.0);
- hdrs_out[i].offset = (hdrs_out[i].gx - hdrs_out[i].sx)/1000.0;
- }
- ret = writeData(fp_out, &Image[0], hdrs_out, shot.nz, shot.nx);
- if (ret < 0 ) verr("error on writing output file.");
+ f2p,G_d,muteW,smooth,shift,above,pad,nt0,&first,verbose);
- fclose(fp_out);
+ if (niter==0) {
+ for (l = 0; l < Nsyn; l++) {
+ for (i = 0; i < npossyn; i++) {
+ j = 0;
+ ix = ixpossyn[i];
+ f2p[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts+j];
+ f1plus[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts+j];
+ green[i*nts+j] = hG_d[ix*nts+j];
+ for (j = 1; j < nts; j++) {
+ f2p[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts+j];
+ f1plus[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts+j];
+ green[i*nts+j] = hG_d[ix*nts+nts-j];
+ }
+ }
+ }
+ }
- if (file_gmin != NULL) {
- fp_gmin = fopen(file_gmin, "w+");
- if (fp_gmin==NULL) verr("error on creating output file %s", file_gmin);
- }
- if (file_f1plus != NULL) {
- fp_f1plus = fopen(file_f1plus, "w+");
- if (fp_f1plus==NULL) verr("error on creating output file %s", file_f1plus);
- }
+ if (file_img!=NULL) {
+
+ /*================ set variables for output data ================*/
- tracf = 1;
- for (j = 0; j < Nsyn; j++) {
- //l = synpos[j];
- l = j;
- if (ixa || ixb) f2 = xsyn[l]-ixb*d2;
- else {
- if (reci) f2 = fxs;
- else f2 = fxf;
- }
+ hdrs_im = (segy *) calloc(shot.nx,sizeof(segy));
+ if (hdrs_im == NULL) verr("allocation for hdrs_out");
+ Image = (float *)calloc(Nsyn,sizeof(float));
- for (i = 0; i < n2; i++) {
- hdrs_out2[i].fldr = l+1;
- hdrs_out2[i].sx = NINT(xsyn[l]*1000);
- hdrs_out2[i].offset = (long)NINT((f2+i*d2) - xsyn[l]);
- hdrs_out2[i].tracf = tracf++;
- hdrs_out2[i].selev = NINT(zsyn[l]*1000);
- hdrs_out2[i].sdepth = NINT(-zsyn[l]*1000);
- hdrs_out2[i].f1 = f1;
- }
+ first=0;
+ imaging(Image,WPs,Refl,nx,nt,nxs,nts,dt,xsyn,Nsyn,xrcv,xsrc,fxs2,fxs,dxs,dxsrc,dx,ixa,ixb,
+ ntfft,nw,nw_low,nw_high,mode,reci,nshots,ixpossyn,npossyn,pmin,f1min,f1plus,
+ f2p,G_d,muteW,smooth,shift,above,pad,nt0,synpos,&first,verbose);
+
+ /*============= write output files ================*/
+
+ fp_out = fopen(file_img, "w+");
+
+ for (i = 0; i < shot.nx; i++) {
+ hdrs_im[i].fldr = 1;
+ hdrs_im[i].tracl = 1;
+ hdrs_im[i].tracf = i+1;
+ hdrs_im[i].scalco = -1000;
+ hdrs_im[i].scalel = -1000;
+ hdrs_im[i].sdepth = 0;
+ hdrs_im[i].trid = 1;
+ hdrs_im[i].ns = shot.nz;
+ hdrs_im[i].trwf = shot.nx;
+ hdrs_im[i].ntr = hdrs_im[i].fldr*hdrs_im[i].trwf;
+ hdrs_im[i].f1 = zsyn[0];
+ hdrs_im[i].f2 = xsyn[0];
+ hdrs_im[i].dt = dt*(1E6);
+ hdrs_im[i].d1 = (float)zsyn[shot.nx]-zsyn[0];
+ hdrs_im[i].d2 = (float)xsyn[1]-xsyn[0];
+ hdrs_im[i].sx = (int)roundf(xsyn[0] + (i*hdrs_im[i].d2));
+ hdrs_im[i].gx = (int)roundf(xsyn[0] + (i*hdrs_im[i].d2));
+ hdrs_im[i].offset = (hdrs_im[i].gx - hdrs_im[i].sx)/1000.0;
+ }
+ ret = writeData(fp_out, &Image[0], hdrs_im, shot.nz, shot.nx);
+ if (ret < 0 ) verr("error on writing output file.");
- if (file_gmin != NULL) {
- ret = writeData(fp_gmin, (float *)&Gmin[l*size], hdrs_out2, n1, n2);
- if (ret < 0 ) verr("error on writing output file.");
- }
- if (file_f1plus != NULL) {
- ret = writeData(fp_f1plus, (float *)&f1plus[l*size], hdrs_out2, n1, n2);
- if (ret < 0 ) verr("error on writing output file.");
- }
+ fclose(fp_out);
}
- ret = 0;
- if (file_gmin != NULL) {ret += fclose(fp_gmin);}
- if (file_f1plus != NULL) {ret += fclose(fp_f1plus);}
- if (ret < 0) verr("err %d on closing output file",ret);
+
+ if (file_homg!=NULL) {
+
+ /*================ set variables for output data ================*/
+
+ hdrs_homg = (segy *) calloc(shot.nx,sizeof(segy));
+ if (hdrs_homg == NULL) verr("allocation for hdrs_out");
+ HomG = (float *)calloc(Nsyn*ntfft,sizeof(float));
+
+ homogeneousg(HomG,green,Refl,nx,nt,nxs,nts,dt,xsyn,Nsyn,xrcv,xsrc,fxs2,fxs,dxs,dxsrc,dx,ixa,ixb,
+ ntfft,nw,nw_low,nw_high,mode,reci,nshots,ixpossyn,npossyn,pmin,f1min,f1plus,
+ f2p,G_d,muteW,smooth,shift,above,pad,nt0,synpos,verbose);
+
+ /*============= write output files ================*/
+
+ fp_out = fopen(file_homg, "w+");
+
+ for (j = 0; j < ntfft; j++) {
+ for (i = 0; i < shot.nx; i++) {
+ hdrs_homg[i].fldr = j+1;
+ hdrs_homg[i].tracl = j*shot.nx+i+1;
+ hdrs_homg[i].tracf = i+1;
+ hdrs_homg[i].scalco = -1000;
+ hdrs_homg[i].scalel = -1000;
+ hdrs_homg[i].sdepth = (int)(zloc*1000.0);
+ hdrs_homg[i].trid = 1;
+ hdrs_homg[i].ns = shot.nz;
+ hdrs_homg[i].trwf = shot.nx;
+ hdrs_homg[i].ntr = hdrs_homg[i].fldr*hdrs_homg[i].trwf;
+ hdrs_homg[i].f1 = zsyn[0];
+ hdrs_homg[i].f2 = xsyn[0];
+ hdrs_homg[i].dt = dt*(1E6);
+ hdrs_homg[i].d1 = (float)zsyn[shot.nx]-zsyn[0];
+ hdrs_homg[i].d2 = (float)xsyn[1]-xsyn[0];
+ hdrs_homg[i].sx = (int)roundf(xsyn[0] + (i*hdrs_homg[i].d2));
+ hdrs_homg[i].gx = (int)roundf(xsyn[0] + (i*hdrs_homg[i].d2));
+ hdrs_homg[i].offset = (hdrs_homg[i].gx - hdrs_homg[i].sx)/1000.0;
+ }
+ ret = writeData(fp_out, &HomG[j*shot.n], hdrs_homg, shot.nz, shot.nx);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+
+ fclose(fp_out);
+ }
if (verbose) {
t1 = wallclock_time();
@@ -682,7 +757,6 @@ int main (int argc, char **argv)
}
- free(hdrs_out);
free(tapersy);
exit(0);
@@ -691,7 +765,7 @@ int main (int argc, char **argv)
/*================ Convolution and Integration ================*/
-void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, double *tfft, int verbose)
+void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, double *tfft, int *first, int verbose)
{
int nfreq, size, iox, inx;
float scl;
@@ -699,7 +773,7 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
float *rtrace, idxs;
complex *sum, *ctrace;
int npe;
- static int first=1, *ixrcv;
+ static int *ixrcv;
static double t0, t1, t;
size = nxs*nts;
diff --git a/marchenko_applications/marchenko.h b/marchenko_applications/marchenko.h
index 7be349d5ef4fd158f467f9dad2707cca8f709b3e..fd627253990f4ff8ea2980be26cc184f786e95ff 100644
--- a/marchenko_applications/marchenko.h
+++ b/marchenko_applications/marchenko.h
@@ -10,7 +10,7 @@
#define WAVEPAR
typedef struct WaveParameters {
int nt, shift, inv, scfft, cm, cn, wav;
- float dt, fp, fmin, flef, frig, fmax, t0, db, scale, eps;
+ float dt, fp, fmin, flef, frig, fmax, t0, db, scale, eps, xloc, zloc;
char w[10], *file_wav;
} WavePar;
#endif
diff --git a/marchenko_applications/raytime.c b/marchenko_applications/raytime.c
index 777df3df4450773c0dc81215de444aefedc305b8..a0fa6d0d3a577a0d11fd927fd8b6b28ae176a74a 100644
--- a/marchenko_applications/raytime.c
+++ b/marchenko_applications/raytime.c
@@ -42,7 +42,7 @@ int defineSource(wavPar wav, srcPar src, modPar mod, float **src_nwav, int rever
int writeSrcRecPos(modPar *mod, recPar *rec, srcPar *src, shotPar *shot);
-int raytime(float *time, float *ampl, int *xnx, float *xrcv, float *xsrc, float *zsrc)
+int raytime(float *time, float *ampl, int *xnx, float *xrcv, float *xsrc, float *zsrc, float xloc, float zloc)
{
modPar mod;
recPar rec;
@@ -70,6 +70,15 @@ int raytime(float *time, float *ampl, int *xnx, float *xrcv, float *xsrc, float
if (!getparint("verbose",&verbose)) verbose=0;
getParameters(&mod, &rec, &src, &shot, &ray, verbose);
+ if (xloc!=-123456.0 && zloc!=-123456.0) {
+ if (verbose > 3) vmess("Setting source ray to x = %.3f, z = %.3f",xloc,zloc);
+ shot.nx = 1;
+ shot.nz = 1;
+ shot.n = 1;
+ shot.x[0] = NINT((xloc-mod.x0)/mod.dx);
+ shot.z[0] = NINT((zloc-mod.z0)/mod.dz);
+ }
+
/* allocate arrays for model parameters: the different schemes use different arrays */
n1 = mod.nz;
diff --git a/marchenko_applications/readTinvData.c b/marchenko_applications/readTinvData.c
index bceab866f43e80bc7dc1550471dee92a8f1a83ca..19927fbd3ef859ffa76730fde22fd021c57c68b5 100644
--- a/marchenko_applications/readTinvData.c
+++ b/marchenko_applications/readTinvData.c
@@ -25,276 +25,164 @@ void findShotInMute(float *xrcvMute, float xrcvShot, int nxs, int *imute);
int readSnapData(char *filename, float *data, segy *hdrs, int nsnaps, int nx, int nz, int sx, int ex, int sz, int ez);
int raytime(float *amp, float *time, int *xnx, float *xrcv, float *xsrc, float *zsrc);
-int readTinvData(char *filename, WavePar WP, char *file_ray, char *file_amp, float dt, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nsyn, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose)
+int readTinvData(char *filename, float dt, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nsyn, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose)
{
FILE *fp;
- segy hdr, *hdrs_mute, *hdrs_amp, *hdrs_wav;
+ segy hdr, *hdrs_mute;
size_t nread;
- char *file_cp;
int fldr_shot, sx_shot, itrace, one_shot, ig, isyn, i, j;
- int end_of_file, nt, gx0, gx1, nfreq, geosp;
+ int end_of_file, nt, gx0, gx1;
int nx1, jmax, imax, tstart, tend, nwav;
- float xmax, tmax, lmax, *wavelet, *wavelet2;
- float scl, scel, *trace, dxrcv, *timeval, dw, *amp;
- complex *cmute, *cwav;
+ float xmax, tmax, lmax;
+ float scl, scel, *trace, dxrcv;
- if (!getparstring("file_cp", &file_cp)) file_cp=NULL;
- if (!getparint("geomspread",&geosp)) geosp=1;
- if (file_cp==NULL) geosp=0;
- /*Check wheter the raytime is used or not*/
- if (file_ray!=NULL || file_cp!=NULL) {
- /*Define parameters*/
- nfreq = ntfft/2+1;
- wavelet = (float *)calloc(ntfft,sizeof(float));
- cwav = (complex *)malloc(nfreq*sizeof(complex));
- cmute = (complex *)malloc(nfreq*sizeof(complex));
- dw = 2*M_PI/(ntfft*dt);
+ /* Reading first header */
- /*Create wavelet using parameters or read in wavelet*/
- if (WP.wav) {
- if (WP.file_wav == NULL) {
- if (verbose>0) vmess("Modeling wavelet");
- freqwave(wavelet, WP.nt, WP.dt, WP.fp, WP.fmin, WP.flef, WP.frig, WP.fmax,
- WP.t0, WP.db, WP.shift, WP.cm, WP.cn, WP.w, WP.scale, WP.scfft, WP.inv, WP.eps, verbose);
- }
- else {
- if (verbose>0) vmess("Reading in wavelet");
- fp = fopen( WP.file_wav, "r" );
- if ( fp == NULL ) {
- perror("Error opening file containing wavelet");
- }
- fclose(fp);
- wavelet2= (float *)calloc(ntfft,sizeof(float));
- hdrs_wav = (segy *)calloc(1, sizeof(segy));
- readSnapData(WP.file_wav, wavelet2, hdrs_wav, Nsyn, 1, ntfft, 0, 1, 0, ntfft);
- nwav = hdrs_wav[0].ns/2;
- for (i=0; i<nwav; i++) {
- wavelet[i] = wavelet2[i];
- wavelet[ntfft-1-i] = wavelet2[hdrs_wav[0].ns-1-i];
- }
- }
- rc1fft(wavelet,cwav,ntfft,-1);
- free(wavelet);
- }
-
- timeval = (float *)calloc(Nsyn*nx,sizeof(float));
- amp = (float *)calloc(Nsyn*nx,sizeof(float));
-
- if (file_ray!=NULL) {
-
- /* Defining mute window using raytimes */
- vmess("Using raytime for mutewindow");
- hdrs_mute = (segy *) calloc(Nsyn,sizeof(segy));
- fp = fopen( file_ray, "r" );
- if ( fp == NULL ) {
- perror("Error opening file containing wavelet");
- }
- fclose(fp);
- readSnapData(file_ray, timeval, hdrs_mute, Nsyn, 1, nx, 0, 1, 0, nx);
-
- /*Check whether the amplitude is also used*/
- if (file_amp != NULL) {
- hdrs_amp = (segy *) calloc(Nsyn,sizeof(segy));
- fp = fopen( file_amp, "r" );
- if ( fp == NULL ) {
- perror("Error opening file containing wavelet");
- }
- fclose(fp);
- readSnapData(file_amp, amp, hdrs_amp, Nsyn, 1, nx, 0, 1, 0, nx);
- }
-
- /*Define source and receiver locations from the raytime*/
- for (isyn=0; isyn<Nsyn; isyn++) {
- for (itrace=0; itrace<nx; itrace++) {
- xrcv[isyn*nx+itrace] = (hdrs_mute[isyn].f1 + hdrs_mute[isyn].d1*((float)itrace));
- }
- xnx[isyn]=nx;
- xsrc[isyn] = hdrs_mute[isyn].sx;
- zsrc[isyn] = hdrs_mute[isyn].sdepth;
- }
- }
- else {
- raytime(timeval,amp,xnx,xrcv,xsrc,zsrc);
- }
+ if (filename == NULL) fp = stdin;
+ else fp = fopen( filename, "r" );
+ if ( fp == NULL ) {
+ fprintf(stderr,"input file %s has an error\n", filename);
+ perror("error in opening file: ");
+ fflush(stderr);
+ return -1;
+ }
- /*Determine the mutewindow*/
- for (j=0; j<Nsyn; j++) {
- for (i=0; i<nx; i++) {
- maxval[j*nx+i] = (int)roundf(timeval[j*nx+i]/dt);
- if (maxval[j*nx+i] > ntfft-1) maxval[j*nx+i] = ntfft-1;
- if (WP.wav) { /*Apply the wavelet to create a first arrival*/
- if (file_amp != NULL || geosp==1) {
- for (ig=0; ig<nfreq; ig++) {
- cmute[ig].r = (cwav[ig].r*cos(ig*dw*timeval[j*nx+i]-M_PI/4.0)-cwav[ig].i*sin(ig*dw*timeval[j*nx+i]-M_PI/4.0))/amp[j*nx+i];
- cmute[ig].i = (cwav[ig].i*cos(ig*dw*timeval[j*nx+i]-M_PI/4.0)+cwav[ig].r*sin(ig*dw*timeval[j*nx+i]-M_PI/4.0))/amp[j*nx+i];
- }
- }
- else { /*Use the raytime only to determine the mutewindow*/
- for (ig=0; ig<nfreq; ig++) {
- cmute[ig].r = cwav[ig].r*cos(ig*dw*timeval[j*nx+i]-M_PI/4.0)-cwav[ig].i*sin(ig*dw*timeval[j*nx+i]-M_PI/4.0);
- cmute[ig].i = cwav[ig].i*cos(ig*dw*timeval[j*nx+i]-M_PI/4.0)+cwav[ig].r*sin(ig*dw*timeval[j*nx+i]-M_PI/4.0);
- }
- }
- cr1fft(cmute,&tinv[j*nx*ntfft+i*ntfft],ntfft,1);
- tinv[j*nx*ntfft+i*ntfft] /= ntfft;
- }
- }
- }
- }
+ fseek(fp, 0, SEEK_SET);
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+ if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
+ else if (hdr.scalco == 0) scl = 1.0;
+ else scl = hdr.scalco;
+ if (hdr.scalel < 0) scel = 1.0/fabs(hdr.scalel);
+ else if (hdr.scalel == 0) scel = 1.0;
+ else scel = hdr.scalel;
+ fseek(fp, 0, SEEK_SET);
- if (WP.wav == 0 && file_ray==NULL && filename!=NULL) {
+ nt = hdr.ns;
+ trace = (float *)calloc(ntfft,sizeof(float));
- /* Reading first header */
+ end_of_file = 0;
+ one_shot = 1;
+ isyn = 0;
- if (filename == NULL) fp = stdin;
- else fp = fopen( filename, "r" );
- if ( fp == NULL ) {
- fprintf(stderr,"input file %s has an error\n", filename);
- perror("error in opening file: ");
- fflush(stderr);
- return -1;
- }
+ /* Read shots in file */
- fseek(fp, 0, SEEK_SET);
+ while (!end_of_file) {
+
+ /* start reading data (shot records) */
+ itrace = 0;
nread = fread( &hdr, 1, TRCBYTES, fp );
- assert(nread == TRCBYTES);
- if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
- else if (hdr.scalco == 0) scl = 1.0;
- else scl = hdr.scalco;
- if (hdr.scalel < 0) scel = 1.0/fabs(hdr.scalel);
- else if (hdr.scalel == 0) scel = 1.0;
- else scel = hdr.scalel;
- fseek(fp, 0, SEEK_SET);
-
- nt = hdr.ns;
- trace = (float *)calloc(ntfft,sizeof(float));
-
- end_of_file = 0;
- one_shot = 1;
- isyn = 0;
-
- /* Read shots in file */
+ if (nread != TRCBYTES) { /* no more data in file */
+ break;
+ }
- while (!end_of_file) {
-
- /* start reading data (shot records) */
- itrace = 0;
+ sx_shot = hdr.sx;
+ fldr_shot = hdr.fldr;
+ gx0 = hdr.gx;
+ xsrc[isyn] = sx_shot*scl;
+ zsrc[isyn] = hdr.selev*scel;
+ xnx[isyn] = 0;
+ ig = isyn*nx*ntfft;
+ while (one_shot) {
+ xrcv[isyn*nx+itrace] = hdr.gx*scl;
+ nread = fread( trace, sizeof(float), nt, fp );
+ assert (nread == hdr.ns);
+
+ /* copy trace to data array */
+ memcpy( &tinv[ig+itrace*ntfft], trace, nt*sizeof(float));
+
+ gx1 = hdr.gx;
+ itrace++;
+
+ /* read next hdr of next trace */
nread = fread( &hdr, 1, TRCBYTES, fp );
- if (nread != TRCBYTES) { /* no more data in file */
+ if (nread != TRCBYTES) {
+ one_shot = 0;
+ end_of_file = 1;
break;
}
+ if ((sx_shot != hdr.sx) || (fldr_shot != hdr.fldr) ) break;
+ }
+ if (verbose>2) {
+ 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);
+ }
- sx_shot = hdr.sx;
- fldr_shot = hdr.fldr;
- gx0 = hdr.gx;
- xsrc[isyn] = sx_shot*scl;
- zsrc[isyn] = hdr.selev*scel;
- xnx[isyn] = 0;
- ig = isyn*nx*ntfft;
- while (one_shot) {
- xrcv[isyn*nx+itrace] = hdr.gx*scl;
- nread = fread( trace, sizeof(float), nt, fp );
- assert (nread == hdr.ns);
-
- /* copy trace to data array */
- memcpy( &tinv[ig+itrace*ntfft], trace, nt*sizeof(float));
-
- gx1 = hdr.gx;
- itrace++;
-
- /* read next hdr of next trace */
- nread = fread( &hdr, 1, TRCBYTES, fp );
- if (nread != TRCBYTES) {
- one_shot = 0;
- end_of_file = 1;
- break;
- }
- if ((sx_shot != hdr.sx) || (fldr_shot != hdr.fldr) ) break;
- }
- if (verbose>2) {
- 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);
- }
-
-
- /* look for maximum in shot record to define mute window */
- /* find consistent (one event) maximum related to maximum value */
- nx1 = itrace;
- xnx[isyn]=nx1;
-
- if (file_ray==NULL) {/*Use the raytime to determine the mutewindow instead of searching*/
- /* alternative find maximum at source position */
- dxrcv = (gx1 - gx0)*scl/(float)(nx1-1);
- imax = NINT(((sx_shot-gx0)*scl)/dxrcv);
- tmax=0.0;
- jmax = 0;
- for (j = 0; j < nt; j++) {
- lmax = fabs(tinv[ig+imax*ntfft+j]);
- if (lmax > tmax) {
- jmax = j;
- tmax = lmax;
- if (lmax > xmax) {
- xmax=lmax;
- }
- }
- }
- maxval[isyn*nx+imax] = jmax;
- if (verbose >= 3) 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++) {
- tstart = MAX(0, (maxval[isyn*nx+i-1]-hw));
- tend = MIN(nt-1, (maxval[isyn*nx+i-1]+hw));
- jmax=tstart;
- tmax=0.0;
- for(j = tstart; j <= tend; j++) {
- lmax = fabs(tinv[ig+i*ntfft+j]);
- if (lmax > tmax) {
- jmax = j;
- tmax = lmax;
- }
- }
- maxval[isyn*nx+i] = jmax;
- }
- /* search backward in trace direction from maximum in file */
- for (i = imax-1; i >=0; i--) {
- tstart = MAX(0, (maxval[isyn*nx+i+1]-hw));
- tend = MIN(nt-1, (maxval[isyn*nx+i+1]+hw));
- jmax=tstart;
- tmax=0.0;
- for (j = tstart; j <= tend; j++) {
- lmax = fabs(tinv[ig+i*ntfft+j]);
- if (lmax > tmax) {
- jmax = j;
- tmax = lmax;
- }
- }
- maxval[isyn*nx+i] = jmax;
- }
- }
-
- if (itrace != 0) { /* end of shot record, but not end-of-file */
- fseek( fp, -TRCBYTES, SEEK_CUR );
- isyn++;
- }
- else {
- end_of_file = 1;
- }
+ /* look for maximum in shot record to define mute window */
+ /* find consistent (one event) maximum related to maximum value */
+ nx1 = itrace;
+ xnx[isyn]=nx1;
+
+ /* alternative find maximum at source position */
+ dxrcv = (gx1 - gx0)*scl/(float)(nx1-1);
+ imax = NINT(((sx_shot-gx0)*scl)/dxrcv);
+ tmax=0.0;
+ jmax = 0;
+ for (j = 0; j < nt; j++) {
+ lmax = fabs(tinv[ig+imax*ntfft+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ if (lmax > xmax) {
+ xmax=lmax;
+ }
+ }
+ }
+ maxval[isyn*nx+imax] = jmax;
+ if (verbose >= 3) 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++) {
+ tstart = MAX(0, (maxval[isyn*nx+i-1]-hw));
+ tend = MIN(nt-1, (maxval[isyn*nx+i-1]+hw));
+ jmax=tstart;
+ tmax=0.0;
+ for(j = tstart; j <= tend; j++) {
+ lmax = fabs(tinv[ig+i*ntfft+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ }
+ }
+ maxval[isyn*nx+i] = jmax;
+ }
+ /* search backward in trace direction from maximum in file */
+ for (i = imax-1; i >=0; i--) {
+ tstart = MAX(0, (maxval[isyn*nx+i+1]-hw));
+ tend = MIN(nt-1, (maxval[isyn*nx+i+1]+hw));
+ jmax=tstart;
+ tmax=0.0;
+ for (j = tstart; j <= tend; j++) {
+ lmax = fabs(tinv[ig+i*ntfft+j]);
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
+ }
+ }
+ maxval[isyn*nx+i] = jmax;
+ }
+
+ if (itrace != 0) { /* end of shot record, but not end-of-file */
+ fseek( fp, -TRCBYTES, SEEK_CUR );
+ isyn++;
+ }
+ else {
+ end_of_file = 1;
+ }
- /* copy trace to data array for mode=-1 */
- /* time reverse trace */
- if (mode==-1) {
- for (i = 0; i < nx1; i++) {
- memcpy( trace, &tinv[ig+i*ntfft], ntfft*sizeof(float));
- j=0;
- tinv[ig+i*ntfft+j] = trace[j];
- for (j=1; j<ntfft; j++) tinv[ig+i*ntfft+ntfft-j] = trace[j];
- }
+ /* copy trace to data array for mode=-1 */
+ /* time reverse trace */
+ if (mode==-1) {
+ for (i = 0; i < nx1; i++) {
+ memcpy( trace, &tinv[ig+i*ntfft], ntfft*sizeof(float));
+ j=0;
+ tinv[ig+i*ntfft+j] = trace[j];
+ for (j=1; j<ntfft; j++) tinv[ig+i*ntfft+ntfft-j] = trace[j];
}
}
- free(trace);
}
+ free(trace);
return 0;
}
diff --git a/marchenko_full/Makefile b/marchenko_full/Makefile
index e76db3ed2377f2a38c261de4e9f2a29812c39529..931555d0562f11bdb186f5b7f00e701ebd9b17f1 100644
--- a/marchenko_full/Makefile
+++ b/marchenko_full/Makefile
@@ -7,7 +7,7 @@ OPTC += -g -O0 -Wall
#ALL: fmute marchenko marchenko2
-ALL: fmute marchenko
+ALL: fmute marchenko_full
SRCJ = fmute.c \
getFileInfo.c \
@@ -58,20 +58,20 @@ fmute: $(OBJJ)
OBJH = $(SRCH:%.c=%.o)
-marchenko: $(OBJH) raytime.h
- $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko $(OBJH) $(LIBS)
+marchenko_full: $(OBJH) raytime.h
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko_full $(OBJH) $(LIBS)
-install: fmute marchenko
+install: fmute marchenko_full
cp fmute $B
- cp marchenko $B
+ cp marchenko_full $B
# cp marchenko2 $B
clean:
- rm -f core fmute $(OBJJ) marchenko $(OBJH)
+ rm -f core fmute $(OBJJ) marchenko_full $(OBJH)
realclean: clean
- rm -f $B/fmute $B/marchenko
+ rm -f $B/fmute $B/marchenko_full
diff --git a/marchenko_full/marchenko.c b/marchenko_full/marchenko.c
index 197efa3b208c37ca0e29e08dc8062270f7837426..26bf5dd27abfbab6231b94e7e02a66317ea8f916 100644
--- a/marchenko_full/marchenko.c
+++ b/marchenko_full/marchenko.c
@@ -49,6 +49,7 @@ int readData(FILE *fp, float *data, segy *hdrs, int n1);
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
double wallclock_time(void);
+int getParameters(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPar *ray, int verbose);
void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nsyn, float *xrcv, float *xsrc, float fxs2, float fxs, float dxs, float dxsrc, float dx, int ixa, int ixb, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpossyn, int npossyn, double *tfft, int verbose);
@@ -214,7 +215,7 @@ int main (int argc, char **argv)
double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi, *J;
float d1, d2, f1, f2, fxs, ft, fx, *xsyn, dxsrc, Q, f0, *Costdet;
float *green, *f2p, *pmin, *G_d, dt, dx, dxs, scl, mem;
- float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus, *Gm0;
+ float *f1plus, *f1min, *iRN, *Ni, *Gmin, *Gplus, *Gm0;
float xmin, xmax, weight, tsq, *Gd, *amp, bstart, bend, db, *bdet, bp, b, bmin;
complex *Refl, *Fop;
char *file_tinv, *file_shot, *file_green, *file_iter, *file_wav, *file_ray, *file_amp;
@@ -367,7 +368,6 @@ int main (int argc, char **argv)
Ni = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
G_d = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
muteW = (int *)calloc(Nsyn*nxs,sizeof(int));
- trace = (float *)malloc(ntfft*sizeof(float));
ixpossyn = (int *)malloc(nxs*sizeof(int));
xrcvsyn = (float *)calloc(Nsyn*nxs,sizeof(float));
xsyn = (float *)malloc(Nsyn*sizeof(float));
@@ -978,7 +978,7 @@ for (ib=0; ib<=nb; ib++) {
}
if (file_f1plus != NULL) {
/* rotate to get t=0 in the middle */
- for (i = 0; i < n2; i++) {
+ /*for (i = 0; i < n2; i++) {
hdrs_out[i].f1 = -n1*0.5*dt;
memcpy(&trace[0],&f1plus[l*size+i*nts],nts*sizeof(float));
for (j = 0; j < n1/2; j++) {
@@ -987,13 +987,13 @@ for (ib=0; ib<=nb; ib++) {
for (j = n1/2; j < n1; j++) {
f1plus[l*size+i*nts+j-n1/2] = trace[j];
}
- }
+ }*/
ret = writeData(fp_f1plus, (float *)&f1plus[l*size], hdrs_out, n1, n2);
if (ret < 0 ) verr("error on writing output file.");
}
if (file_f1min != NULL) {
/* rotate to get t=0 in the middle */
- for (i = 0; i < n2; i++) {
+ /*for (i = 0; i < n2; i++) {
hdrs_out[i].f1 = -n1*0.5*dt;
memcpy(&trace[0],&f1min[l*size+i*nts],nts*sizeof(float));
for (j = 0; j < n1/2; j++) {
@@ -1002,7 +1002,7 @@ for (ib=0; ib<=nb; ib++) {
for (j = n1/2; j < n1; j++) {
f1min[l*size+i*nts+j-n1/2] = trace[j];
}
- }
+ }*/
ret = writeData(fp_f1min, (float *)&f1min[l*size], hdrs_out, n1, n2);
if (ret < 0 ) verr("error on writing output file.");
}