diff --git a/marchenko_applications/Makefile b/marchenko_applications/Makefile
index b3aa063811ab7468e24b3da1a9318806f93c71fb..c810a755fc2cbf61a429ef838a0e872178af576e 100644
--- a/marchenko_applications/Makefile
+++ b/marchenko_applications/Makefile
@@ -36,7 +36,6 @@ SRCH = marchenko.c \
getpars.c \
readSnapData.c \
Cost.c \
- AmpEst.c \
freqwave.c \
getParameters.c \
getModelInfo.c \
@@ -53,7 +52,8 @@ SRCH = marchenko.c \
imaging.c \
threadAffinity.c \
makeWindow.c \
- homogeneousg.c
+ homogeneousg.c \
+ AmpEstApp.c
SRCC = combine.c \
getFileInfo.c \
diff --git a/marchenko_applications/marchenko.c b/marchenko_applications/marchenko.c
index 83a1243da812489a781374d98f040bd531473a8d..1e9fdc1e9d3df87dc519eb083910e5f211e8331a 100644
--- a/marchenko_applications/marchenko.c
+++ b/marchenko_applications/marchenko.c
@@ -43,7 +43,6 @@ int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, floa
void name_ext(char *filename, char *extension);
void Cost(float *f1p, float *f1d, float *Gm, float *Gm0, double *J, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn);
void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs, int nt, int *xrcvsyn, int npossyn, int shift, int pad, int nt0);
-void AmpEst(float *f1d, float *Gd, float *ampest, int Nsyn, int nxs, int ntfft, int *ixpossyn, int npossyn, char *file_wav);
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 *nxm);
int readData(FILE *fp, float *data, segy *hdrs, int n1);
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
@@ -54,6 +53,8 @@ void iterations (complex *Refl, int nx, int nt, int nxs, int nts, float dt, floa
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 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 AmpEst(float *amp, 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 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);
@@ -209,10 +210,10 @@ int main (int argc, char **argv)
FILE *fp_out, *fp_f1plus, *fp_f1min;
FILE *fp_gmin, *fp_gplus, *fp_f2, *fp_pmin;
int i, j, l, ret, nshots, Nsyn, nt, nx, nts, nxs, ngath;
- int size, n1, n2, ntap, tap, di, ntraces, nb, ib;
+ int size, n1, n2, ntap, tap, di, ntraces, nb, ib, ampest;
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn, *synpos;
int reci, mode, ixa, ixb, n2out, verbose, ntfft;
- int iter, niter, niterh, tracf, *muteW, pad, nt0, ampest, *hmuteW, *hxnxsyn;
+ int iter, niter, niterh, tracf, *muteW, pad, nt0, *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;
@@ -268,6 +269,7 @@ int main (int argc, char **argv)
if (!getparfloat("fmax", &fmax)) fmax = 70.0;
if (!getparint("ixa", &ixa)) ixa = 0;
if (!getparint("ixb", &ixb)) ixb = ixa;
+ if (!getparint("ampest",&est)) ampest = 0;
// if (!getparint("reci", &reci)) reci = 0;
reci=0; // source-receiver reciprocity is not yet fully build into the code
if (!getparfloat("weight", &weight)) weight = 1.0;
@@ -282,7 +284,6 @@ int main (int argc, char **argv)
if(!getparint("smooth", &smooth)) smooth = 5;
if(!getparint("above", &above)) above = 0;
if(!getparint("shift", &shift)) shift=12;
- if(!getparint("ampest", &est)) ampest=0;
if(!getparint("nb", &nb)) nb=0;
if (!getparfloat("bstart", &bstart)) bstart = 1.0;
if (!getparfloat("bend", &bend)) bend = 1.0;
@@ -336,7 +337,7 @@ int main (int argc, char **argv)
ngath = 0; /* setting ngath=0 scans all traces; n2 contains maximum traces/gather */
if (file_ray!=NULL && file_tinv==NULL) {
- ret = getFileInfo(file_ray, &n2, &n1, &ngath, &d2, &d1, &f2, &f1, &xmin, &xmax, &scl, &ntraces);
+ ret = getFileInfo(file_ray, &n2, &n1, &ngath, &d1, &d2, &f2, &f1, &xmin, &xmax, &scl, &ntraces);
n1 = 1;
ntraces = n2*ngath;
scl = 0.0010;
@@ -481,14 +482,14 @@ int main (int argc, char **argv)
}
/* check consistency of header values */
- if (xrcvsyn[0] != 0 || xrcvsyn[1] != 0 ) fxs = xrcvsyn[0];
- fxs2 = fxs + (float)(nxs-1)*dxs;
dxf = (xrcvsyn[nxs-1] - xrcvsyn[0])/(float)(nxs-1);
if (NINT(dxs*1e3) != NINT(fabs(dxf)*1e3)) {
vmess("dx in hdr.d1 (%.3f) and hdr.gx (%.3f) not equal",d2, dxf);
if (dxf != 0) dxs = fabs(dxf);
vmess("dx in operator => %f", dxs);
}
+ if (xrcvsyn[0] != 0 || xrcvsyn[1] != 0 ) fxs = xrcvsyn[0];
+ fxs2 = fxs + (float)(nxs-1)*dxs;
/*================ Reading shot records ================*/
@@ -608,7 +609,7 @@ int main (int argc, char **argv)
ngath = 1;
if (file_rays!=NULL || file_cp!=NULL) {
- vmess("check");
+ WPs.wav=1;
makeWindow(WPs, file_rays, file_amps, dt, hxrcvsyn, hxsyn, hzsyn, hxnxsyn, ngath, nxs, ntfft, mode, hmuteW, hG_d, hw, verbose);
}
else {
@@ -695,6 +696,24 @@ int main (int argc, char **argv)
}
}*/
+ /* compute downgoing Green's function G^+,+ */
+ if (ampest==1) {
+ amp = (float *)calloc(Nsyn,sizeof(float));
+ AmpEst(amp,WP,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);
+ for (l = 0; l < Nsyn; l++) {
+ for (i = 0; i < nxs*nts; i++) {
+ G_d[l*nxs*nts+i] *= amp[l];
+ f2p[l*nxs*nts+i] *= amp[l];
+ f1plus[l*nxs*nts+i] *= amp[l];
+ f1min[l*nxs*nts+i] *= amp[l];
+ pmin[l*nxs*nts+i] *= amp[l];
+ }
+ }
+ }
+
+
if (niterh==0) {
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
diff --git a/marchenko_full/marchenko.c b/marchenko_full/marchenko.c
index 685a12745d894dc11de45a35b8f37a40777171fd..811513be9b553eb624c78f550691c4d318b9be82 100644
--- a/marchenko_full/marchenko.c
+++ b/marchenko_full/marchenko.c
@@ -305,7 +305,17 @@ int main (int argc, char **argv)
ngath = 0; /* setting ngath=0 scans all traces; n2 contains maximum traces/gather */
if (file_ray!=NULL && file_tinv==NULL) {
- ret = getFileInfo(file_ray, &n2, &n1, &ngath, &d2, &d1, &f2, &f1, &xmin, &xmax, &scl, &ntraces);
+ ret = getFileInfo(file_ray, &n2, &n1, &ngath, &d1, &d2, &f2, &f1, &xmin, &xmax, &scl, &ntraces);
+ n1 = 1;
+ ntraces = n2*ngath;
+ scl = 0.0010;
+ d1 = -1.0*xmin;
+ xmin = -1.0*xmax;
+ xmax = d1;
+ WP.wav = 1;
+ shot.nz = 1;
+ shot.nx = ngath;
+ shot.n = shot.nx*shot.nz;
}
else if (file_ray==NULL && file_tinv==NULL) {
getParameters(&mod, &rec, &src, &shot, &ray, verbose);
@@ -418,14 +428,14 @@ int main (int argc, char **argv)
}
/* check consistency of header values */
- if (xrcvsyn[0] != 0 || xrcvsyn[1] != 0 ) fxs = xrcvsyn[0];
- fxs2 = fxs + (float)(nxs-1)*dxs;
dxf = (xrcvsyn[nxs-1] - xrcvsyn[0])/(float)(nxs-1);
if (NINT(dxs*1e3) != NINT(fabs(dxf)*1e3)) {
vmess("dx in hdr.d1 (%.3f) and hdr.gx (%.3f) not equal",d2, dxf);
if (dxf != 0) dxs = fabs(dxf);
vmess("dx in operator => %f", dxs);
}
+ if (xrcvsyn[0] != 0 || xrcvsyn[1] != 0 ) fxs = xrcvsyn[0];
+ fxs2 = fxs + (float)(nxs-1)*dxs;
/*================ Reading shot records ================*/
diff --git a/marchenko_full/readTinvData.c b/marchenko_full/readTinvData.c
index 5a19d891a1705ff8767017f6e3674328b8174f47..e21a38b9590ec40cea5e4a285e8082516a00e710 100644
--- a/marchenko_full/readTinvData.c
+++ b/marchenko_full/readTinvData.c
@@ -110,9 +110,11 @@ int readTinvData(char *filename, WavePar WP, char *file_ray, char *file_amp, flo
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;
+ xnx[isyn]=hdrs_mute[isyn].ns;
+ if (hdrs_mute[isyn].scalco < 0) scl=-1.0/hdrs_mute[isyn].scalco;
+ else scl=hdrs_mute[isyn].scalco;
+ xsrc[isyn] = hdrs_mute[isyn].sx*scl;
+ zsrc[isyn] = hdrs_mute[isyn].sdepth*scl;
}
}
else {
diff --git a/raytime/raytime.c b/raytime/raytime.c
index 3e8fc7bc0bd2aab6bee1cb302890a817d95654e5..84697bb5ceff7fcbcd0cbe67148228980442879d 100644
--- a/raytime/raytime.c
+++ b/raytime/raytime.c
@@ -121,7 +121,7 @@ int main(int argc, char **argv)
fcoord coordsx, coordgx, Time;
icoord grid, isrc;
float Jr, *ampl, *time, *ttime, *ttime_p, cp_average;
- float dxrcv, dzrcv;
+ float dxrcv, dzrcv, dt_tmp;
segy hdr;
char filetime[1024], fileamp[1024], *method;
size_t nwrite;
@@ -314,6 +314,8 @@ private (coordgx,irec,Time,Jr)
hdr.f1 = mod.x0+rec.x[0]*mod.dx;
hdr.d2 = (shot.x[MIN(shot.n-1,1)]-shot.x[0])*mod.dx;
hdr.f2 = mod.x0+shot.x[0]*mod.dx;
+ dt_tmp = (fabs(hdr.d1*((float)hdr.scalco)));
+ hdr.dt = (unsigned short)dt_tmp;
nwrite = fwrite( &hdr, 1, TRCBYTES, fpt);
assert(nwrite == TRCBYTES);