diff --git a/marchenko_applications/homogeneousg.c b/marchenko_applications/homogeneousg.c
index 4f5bfa3510ea6a26bc947ab688ee478fea79919b..eaa1449def02db5aecda0779a572a0d769971f3b 100644
--- a/marchenko_applications/homogeneousg.c
+++ b/marchenko_applications/homogeneousg.c
@@ -26,17 +26,22 @@ void corr(float *data1, float *data2, float *cov, int nrec, int nsam, float dt,
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 kxwfilter(float *data, int nt, int nx, float dt, float dx, float fmin, float fmax, float angle, float cp, float perc);
+void timeDiff(float *data, int nsam, int nrec, float dt, float fmin, float fmax, int opt);
+void depthDiff(float *data, int nsam, int nrec, float dt, float dx, float fmin, float fmax, float c, int opt);
+void pad2d_data(float *data, int nsam, int nrec, int nsamout, int nrecout, float *datout);
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)
{
- int i, j, l, ret, scheme;
+ int i, j, l, ret, scheme,count=0;
int iter, niter, ix, kxwfilt;
float scl, *conv, *greenf2, fmin, fmax, alpha, cp, perc;
+ float *greenjkz, *greenf2jkz, *tmp1, *tmp2;
double t0, t2, tfft;
FILE *fp;
if (!getparint("scheme", &scheme)) scheme = 0;
if (!getparint("kxwfilt", &kxwfilt)) kxwfilt = 0;
+ if (!getparint("niter",&niter)) niter=10;
if (!getparfloat("fmin", &fmin)) fmin = 0.0;
if (!getparfloat("fmax", &fmax)) fmax = 100.0;
if (!getparfloat("alpha", &alpha)) alpha = 65.0;
@@ -50,37 +55,71 @@ void homogeneousg(float *HomG, float *green, complex *Refl, int nx, int nt, int
if (scheme==1) {
if (verbose) vmess("Classical Homogeneous Green's function retrieval");
- greenf2 = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
+ greenf2 = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
+ greenf2jkz = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
+ greenjkz = (float *)calloc(nxs*ntfft,sizeof(float));
for (l = 0; l < Nsyn; l++) {
- for (i = 0; i < npossyn; i++) {
- j = 0;
- /* set green to zero if mute-window exceeds nt/2 */
- if (muteW[l*nxs+ixpossyn[i]] >= nts/2) {
- memset(&greenf2[l*nxs*nts+i*nts],0, sizeof(float)*nt);
- continue;
- }
- greenf2[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+j] + pmin[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- greenf2[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+nts-j] + pmin[l*nxs*nts+i*nts+j];
+ if (niter<1) {
+ for (i = 0; i < npossyn; i++) {
+ j = 0;
+ ix = ixpossyn[i];
+ greenf2[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts+j];
+ greenf2jkz[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts+j];
+ for (j = 1; j < nts; j++) {
+ greenf2[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts+nts-j];
+ greenf2jkz[l*nxs*nts+i*nts+j] = G_d[l*nxs*nts+ix*nts+nts-j];
+ }
}
+ }
+ else {
+ for (i = 0; i < npossyn; i++) {
+ j = 0;
+ /* set green to zero if mute-window exceeds nt/2 */
+ if (muteW[l*nxs+ixpossyn[i]] >= nts/2) {
+ memset(&greenf2[l*nxs*nts+i*nts],0, sizeof(float)*nt);
+ memset(&greenf2jkz[l*nxs*nts+i*nts],0, sizeof(float)*nt);
+ continue;
+ }
+ greenf2[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+j] + pmin[l*nxs*nts+i*nts+j];
+ greenf2jkz[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+j] + pmin[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ greenf2[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+nts-j] + pmin[l*nxs*nts+i*nts+j];
+ greenf2jkz[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+nts-j] + pmin[l*nxs*nts+i*nts+j];
+ }
+ }
}
+ depthDiff(&greenf2jkz[l*nxs*nts], ntfft, nshots, dt, dx, fmin, fmax, cp, 1);
}
applyMute(greenf2, muteW, smooth, 4, Nsyn, nxs, nts, ixpossyn, npossyn, shift, pad, nt0);
+ applyMute(greenf2jkz, muteW, smooth, 4, Nsyn, nxs, nts, ixpossyn, npossyn, shift, pad, nt0);
+
+ for (i = 0; i < npossyn; i++) {
+ for (j = 0; j < nts; j++) {
+ greenjkz[i*nts+j] = green[i*nts+j];
+ }
+ }
+ depthDiff(greenjkz, ntfft, nshots, dt, dx, fmin, fmax, cp, 1);
}
else {
if (verbose) vmess("Marchenko Homogeneous Green's function retrieval");
}
#pragma omp parallel default(shared) \
- private(i,j,conv)
+ private(i,j,conv,tmp1,tmp2)
{
- conv = (float *)calloc(nxs*ntfft,sizeof(float));
+ conv = (float *)calloc(nxs*ntfft,sizeof(float));
+ if (scheme==1) {
+ tmp1 = (float *)calloc(nxs*ntfft,sizeof(float));
+ tmp2 = (float *)calloc(nxs*ntfft,sizeof(float));
+ }
#pragma omp for
for (l = 0; l < Nsyn; l++) {
- if (verbose > 2) vmess("Creating Homogeneous G at location %d out of %d",l+1,Nsyn);
+ count+=1;
+
+ if (verbose > 2) vmess("Creating Homogeneous G at location %d out of %d",count,Nsyn);
if (scheme==0) { //Marchenko representation
convol(green, &f2p[l*nxs*nts], conv, nxs, nts, dt, 0);
if (kxwfilt) {
@@ -96,10 +135,17 @@ void homogeneousg(float *HomG, float *green, complex *Refl, int nx, int nt, int
}
}
else if (scheme==1) { //classical representation
- corr(green, &greenf2[l*nxs*nts], conv, nxs, nts, dt, 0);
+ corr(&greenf2[l*nxs*nts], greenjkz, tmp1, nxs, nts, dt, 0);
+ corr(&greenf2jkz[l*nxs*nts], green, tmp2, nxs, nts, dt, 0);
+ for (i = 0; i < npossyn; i++) {
+ for (j = 0; j < nts; j++) {
+ conv[i*nts+j] = -(tmp2[i*nts+j]-tmp1[i*nts+j]);
+ }
+ }
+ timeDiff(conv, ntfft, nshots, dt, fmin, fmax, -1);
for (i=0; i<npossyn; i++) {
j=0;
- HomG[(j+nts/2)*Nsyn+synpos[l]] += scl*(conv[i*nts+j] + conv[i*nts+j]);
+ HomG[(j+nts/2)*Nsyn+synpos[l]] += scl*(conv[i*nts+j]);
for (j=1; j<nts/2; j++) {
HomG[(j+nts/2)*Nsyn+synpos[l]] += scl*conv[i*nts+j];
HomG[j*Nsyn+synpos[l]] += scl*conv[i*nts+(j+nts/2)];
@@ -108,8 +154,9 @@ void homogeneousg(float *HomG, float *green, complex *Refl, int nx, int nt, int
}
}
free(conv);
+ if (scheme==1) free(tmp1);free(tmp2);
}
- if (scheme==1) free(greenf2);
+ if (scheme==1) free(greenf2);free(greenf2jkz);
t2 = wallclock_time();
if (verbose) {
@@ -199,3 +246,211 @@ void corr(float *data1, float *data2, float *cov, int nrec, int nsam, float dt,
free(rdata2);
return;
}
+
+void timeDiff(float *data, int nsam, int nrec, float dt, float fmin, float fmax, int opt)
+{
+ int optn, iom, iomin, iomax, nfreq, ix, sign;
+ float omin, omax, deltom, om, df, *rdata, scl;
+ complex *cdata, *cdatascl;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+ df = 1.0/(optn*dt);
+
+ cdata = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data,nsam,nrec,optn,rdata);
+
+ /* Forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata[0], &cdata[0], optn, nrec, optn, nfreq, sign);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+ iomin = (int)MIN((omin/deltom), (nfreq));
+ iomin = MAX(iomin, 1);
+ iomax = MIN((int)(omax/deltom), (nfreq));
+
+ cdatascl = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdatascl == NULL) verr("memory allocation error for cdatascl");
+
+ for (ix = 0; ix < nrec; ix++) {
+ for (iom = 0; iom < iomin; iom++) {
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
+ for (iom = iomax; iom < nfreq; iom++) {
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
+ if (opt > 0) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = deltom*iom;
+ cdatascl[ix*nfreq+iom].r = -om*cdata[ix*nfreq+iom].i;
+ cdatascl[ix*nfreq+iom].i = om*cdata[ix*nfreq+iom].r;
+ }
+ }
+ else if (opt < 0) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = 1.0/(deltom*iom);
+ cdatascl[ix*nfreq+iom].r = om*cdata[ix*nfreq+iom].i;
+ cdatascl[ix*nfreq+iom].i = -om*cdata[ix*nfreq+iom].r;
+ }
+ }
+ }
+ free(cdata);
+
+ /* Inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&cdatascl[0], &rdata[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
+
+void depthDiff(float *data, int nsam, int nrec, float dt, float dx, float fmin, float fmax, float c, int opt)
+{
+ int optn, iom, iomin, iomax, nfreq, ix, ikx, nkx, ikxmax;
+ float omin, omax, deltom, df, dkx, *rdata, kx, scl;
+ float kx2, kz2, kp2, kp;
+ complex *cdata, *cdatascl, kz, kzinv;
+
+ optn = optncr(nsam);
+ nfreq = optncr(nsam)/2+1;
+ df = 1.0/(optn*dt);
+ nkx = optncc(nrec);
+ dkx = 2.0*PI/(nkx*dx);
+ cdata = (complex *)malloc(nfreq*nkx*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nkx*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes in 2 directions to reach FFT lengths */
+ pad2d_data(data,nsam,nrec,optn,nkx,rdata);
+
+ /* double forward FFT */
+ xt2wkx(&rdata[0], &cdata[0], optn, nkx, optn, nkx, 0);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+
+ iomin = (int)MIN((omin/deltom), nfreq);
+ iomin = MAX(iomin, 0);
+ iomax = MIN((int)(omax/deltom), nfreq);
+
+ cdatascl = (complex *)malloc(nfreq*nkx*sizeof(complex));
+ if (cdatascl == NULL) verr("memory allocation error for cdatascl");
+
+ for (iom = 0; iom < iomin; iom++) {
+ for (ix = 0; ix < nkx; ix++) {
+ cdatascl[iom*nkx+ix].r = 0.0;
+ cdatascl[iom*nkx+ix].i = 0.0;
+ }
+ }
+ for (iom = iomax; iom < nfreq; iom++) {
+ for (ix = 0; ix < nkx; ix++) {
+ cdatascl[iom*nkx+ix].r = 0.0;
+ cdatascl[iom*nkx+ix].i = 0.0;
+ }
+ }
+ if (opt > 0) {
+ for (iom = iomin ; iom <= iomax ; iom++) {
+ kp = (iom*deltom)/c;
+ kp2 = kp*kp;
+
+ ikxmax = MIN((int)(kp/dkx), nkx/2);
+
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kz.r-cdata[iom*nkx+ikx].i*kz.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kz.r+cdata[iom*nkx+ikx].r*kz.i;
+
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nkx+ikx].r = 0.0;
+ cdatascl[iom*nkx+ikx].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kz.r-cdata[iom*nkx+ikx].i*kz.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kz.r+cdata[iom*nkx+ikx].r*kz.i;
+ }
+ }
+ }
+ else if (opt < 0) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ kp = iom*deltom/c;
+ kp2 = kp*kp;
+ ikxmax = MIN((int)(kp/dkx), nkx/2);
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kzinv.r-cdata[iom*nkx+ikx].i*kzinv.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kzinv.r+cdata[iom*nkx+ikx].r*kzinv.i;
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nkx+ikx].r = 0.0;
+ cdatascl[iom*nkx+ikx].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kzinv.r-cdata[iom*nkx+ikx].i*kzinv.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kzinv.r+cdata[iom*nkx+ikx].r*kzinv.i;
+ }
+ }
+ }
+ free(cdata);
+
+ /* inverse double FFT */
+ wkx2xt(&cdatascl[0], &rdata[0], optn, nkx, nkx, optn, 0);
+ /* select original samples and traces */
+ scl = 1.0;
+ scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
+
+void pad2d_data(float *data, int nsam, int nrec, int nsamout, int nrecout, float *datout)
+{
+ int it,ix;
+ for (ix=0;ix<nrec;ix++) {
+ for (it=0;it<nsam;it++)
+ datout[ix*nsam+it]=data[ix*nsam+it];
+ for (it=nsam;it<nsamout;it++)
+ datout[ix*nsam+it]=0.0;
+ }
+ for (ix=nrec;ix<nrecout;ix++) {
+ for (it=0;it<nsamout;it++)
+ datout[ix*nsam+it]=0.0;
+ }
+}
diff --git a/marchenko_applications/marchenko.c b/marchenko_applications/marchenko.c
index 415ac0d32a49e2344029975a69f78dc1c62a0a7e..7cd6708985fc40a7f92c982e57e2847e201d051d 100644
--- a/marchenko_applications/marchenko.c
+++ b/marchenko_applications/marchenko.c
@@ -240,7 +240,7 @@ int main (int argc, char **argv)
tsyn = tread = tfft = tcopy = 0.0;
t0 = wallclock_time();
- if (!getparstring("file_img", &file_img)) file_img = "img.su";
+ if (!getparstring("file_img", &file_img)) file_img = NULL;
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;
@@ -263,6 +263,8 @@ 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 (file_img == NULL && file_homg == NULL)
+ verr("file_img and file_homg cannot both be empty");
if (!getparstring("file_green", &