diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index 0695d33765687efed7a1d0332828eb4dc76e22a0..4df7250c3e1234d1a7de392234f8d30b24150082 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -68,8 +68,7 @@ void synthesis3D(complex *Refl, complex *Fop, float *Top, float *iRN, long nx, l
float dx, float dy, long ntfft, long nw, long nw_low, long nw_high, long mode, long reci, long nshots, long nxsrc, long nysrc,
long *ixpos, long *iypos, long npos, double *tfft, long *isxcount, long *reci_xsrc, long *reci_xrcv, float *ixmask, long verbose);
-
-long writeDataIter3D(char *file_iter, float *data, segy *hdrs, long n1, long n2, long n3, float d2, float f2, long n2out, long Nfoc, float *xsyn, float *ysyn, float *zsyn, long *ixpos, int *iypos, long npos, long iter);
+long writeDataIter3D(char *file_iter, float *data, segy *hdrs, long n1, long n2, long n3, long Nfoc, float *xsyn, float *ysyn, float *zsyn, long *ixpos, long *iypos, long npos, long t0shift, long iter);
void imaging3D(float *Image, float *Gmin, float *f1plus, long nx, long ny, long nt, float dx, float dy, float dt, long Nfoc, long verbose);
@@ -162,8 +161,8 @@ int main (int argc, char **argv)
long size, n1, n2, n3, ntap, tap, dxi, dyi, ntraces, pad, *sx, *sy, *sz;
long nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
long reci, countmin, mode, n2out, n3out, verbose, ntfft;
- long iter, niter, tracf, *muteW, *tsynW, ampest, plane_wave, compact;
- long hw, smooth, above, shift, *ixpos, *iypos, npos, ix, iy, nzim, nxim, nyim;
+ long iter, niter, tracf, *muteW, *tsynW, ampest, plane_wave, t0shift;
+ long hw, smooth, above, shift, *ixpos, *iypos, npos, ix, iy, nzim, nxim, nyim, compact;
long nshots_r, *isxcount, *reci_xsrc, *reci_xrcv;
float fmin, fmax, *tapersh, *tapersy, fxf, fyf, dxf, dyf, *xsrc, *ysrc, *xrcv, *yrcv, *zsyn, *zsrc, *xrcvsyn, *yrcvsyn;
double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi, energyN0;
@@ -177,7 +176,7 @@ int main (int argc, char **argv)
char *file_tinv, *file_shot, *file_green, *file_iter, *file_imag, *file_homg, *file_ampscl;
char *file_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_pmin, *file_inp;
char *file_ray, *file_amp, *file_wav;
- segy *hdrs_out, *hdrs_Nfoc;
+ segy *hdrs_out, *hdrs_Nfoc, *hdrs_iter;
initargs(argc, argv);
requestdoc(1);
@@ -609,6 +608,28 @@ int main (int argc, char **argv)
hdrs_out[k*n2+i].tracl = k*n2+i+1;
}
}
+ if (file_iter != NULL) {
+ hdrs_iter = (segy *) calloc(npos,sizeof(segy));
+ if (hdrs_iter == NULL) verr("allocation for hdrs_iter");
+ for (i = 0; i < npos; i++) {
+ ix = ixpos[i];
+ iy = iypos[i];
+ hdrs_iter[i].ns = n1;
+ hdrs_iter[i].trid = 1;
+ hdrs_iter[i].dt = dt*1000000;
+ hdrs_iter[i].f1 = f1;
+ hdrs_iter[i].f2 = f2;
+ hdrs_iter[i].d1 = d1;
+ hdrs_iter[i].d2 = d2;
+ hdrs_iter[i].trwf = npos;
+ hdrs_iter[i].scalco = -1000;
+ hdrs_iter[i].gx = NINT(1000*(f2+ix*d2));
+ hdrs_iter[i].gy = NINT(1000*(f3+iy*d3));
+ hdrs_iter[i].scalel = -1000;
+ hdrs_iter[i].tracl = i+1;
+ }
+ }
+
t1 = wallclock_time();
tread = t1-t0;
@@ -646,11 +667,11 @@ int main (int argc, char **argv)
t3 = wallclock_time();
tsyn += t3 - t2;
-/* ToDo
if (file_iter != NULL) {
- writeDataIter3D(file_iter, iRN, hdrs_out, ntfft, nxs, nys, d2, f2, n2out, n3out, Nfoc, xsyn,ysyn, zsyn, ixpos, npos, iter);
+ t0shift=1;
+ writeDataIter3D(file_iter, iRN, hdrs_iter, ntfft, nxs, nys, Nfoc, xsyn, ysyn, zsyn, ixpos, iypos, npos, t0shift, iter);
}
-*/
+
/* N_k(x,t) = -N_(k-1)(x,-t) */
/* p0^-(x,t) += iRN = (R * T_d^inv)(t) */
for (l = 0; l < Nfoc; l++) {
diff --git a/marchenko3D/synthesis3D.c b/marchenko3D/synthesis3D.c
index 416542d19aaf4a4420fd9aa0b3eaafe7166d4bd0..41a66b0dfc787a2f1e24ec1044e6b5a5851507e0 100644
--- a/marchenko3D/synthesis3D.c
+++ b/marchenko3D/synthesis3D.c
@@ -168,10 +168,10 @@ long *ixpos, long *iypos, long npos, double *tfft, long *isxcount, long *reci_xs
#ifdef _OPENMP
npe = (long)omp_get_max_threads();
- /* parallelisation is over number of virtual source positions (Nfoc) */
- if (npe > Nfoc) {
- vmess("Number of OpenMP threads set to %li (was %li)", Nfoc, npe);
- omp_set_num_threads((int)Nfoc);
+ /* parallelisation is over number of shot positions (nshots) */
+ if (npe > nshots) {
+ vmess("Number of OpenMP threads set to %li (was %li)", nshots, npe);
+ omp_set_num_threads((int)nshots);
}
#endif
@@ -226,27 +226,28 @@ long *ixpos, long *iypos, long npos, double *tfft, long *isxcount, long *reci_xs
t1 = wallclock_time();
*tfft += t1 - t0;
-/* Loop over total number of shots */
if (reci == 0 || reci == 1) {
- for (k=0; k<nshots; k++) {
- if ((xsrc[k] < fxb) || (xsrc[k] > fxe) || (ysrc[k] < fyb) || (ysrc[k] > fye)) continue;
- isrc = NINT((ysrc[k] - fysb)/dys)*nxs+NINT((xsrc[k] - fxsb)/dxs);
- inx = xnx[k]; /* number of traces per shot */
/*================ SYNTHESIS ================*/
#pragma omp parallel default(none) \
- shared(iRN, dx, dy, npe, nw, verbose) \
+ shared(iRN, dx, dy, npe, nw, nshots, verbose) \
shared(Refl, Nfoc, reci, xrcv, xsrc, yrcv, ysrc, xsyn, ysyn) \
shared(fxsb, fxse, fysb, fyse, nxs, nys, nxys, dxs, dys) \
- shared(nx, ny, nxy, dysrc, dxsrc, inx, k, nfreq, nw_low, nw_high) \
- shared(Fop, size, nts, ntfft, scl, iyrcv, ixrcv, isrc) \
- private(l, ix, iy, j, m, i, sum, rtrace)
+ shared(nx, ny, nxy, dysrc, dxsrc, nfreq, nw_low, nw_high, xnx) \
+ shared(Fop, size, nts, ntfft, scl, iyrcv, ixrcv, fxb, fxe, fyb, fye) \
+ private(l, ix, iy, j, m, i, sum, rtrace, k, isrc, inx)
{ /* start of parallel region */
- sum = (complex *)malloc(nfreq*sizeof(complex));
- rtrace = (float *)calloc(ntfft,sizeof(float));
+ sum = (complex *)malloc(nfreq*sizeof(complex));
+ rtrace = (float *)calloc(ntfft,sizeof(float));
+/* Loop over total number of shots */
#pragma omp for schedule(guided,1)
+ for (k=0; k<nshots; k++) {
+ if ((xsrc[k] < fxb) || (xsrc[k] > fxe) || (ysrc[k] < fyb) || (ysrc[k] > fye)) continue;
+ isrc = NINT((ysrc[k] - fysb)/dys)*nxs+NINT((xsrc[k] - fxsb)/dxs);
+ inx = xnx[k]; /* number of traces per shot */
+
for (l = 0; l < Nfoc; l++) {
/* compute integral over receiver positions */
/* multiply R with Fop and sum over nx */
@@ -269,19 +270,16 @@ long *ixpos, long *iypos, long npos, double *tfft, long *isxcount, long *reci_xs
for (j = 0; j < nts; j++)
iRN[l*size+isrc*nts+j] += rtrace[j]*scl*dx*dy;
- } /* end of parallel Nfoc loop */
- free(sum);
- free(rtrace);
+ } /* end of Nfoc loop */
-#ifdef _OPENMP
-#pragma omp single
- npe = (long)omp_get_num_threads();
-#endif
-} /* end of parallel region */
+ if (verbose>4) vmess("*** Shot gather %li processed ***", k);
+
+ } /* end of parallel nshots (k) loop */
+ free(sum);
+ free(rtrace);
- if (verbose>4) vmess("*** Shot gather %li processed ***", k);
+} /* end of parallel region */
- } /* end of nshots (k) loop */
} /* end of if reci */
t = wallclock_time() - t0;
@@ -290,4 +288,4 @@ long *ixpos, long *iypos, long npos, double *tfft, long *isxcount, long *reci_xs
}
return;
-}
\ No newline at end of file
+}
diff --git a/marchenko3D/writeDataIter3D.c b/marchenko3D/writeDataIter3D.c
index 0820e4301e3b7a851fe7353e2f8003341a59539a..fbaecdb2f4d71bb67a5982bae554c5aee5078955 100644
--- a/marchenko3D/writeDataIter3D.c
+++ b/marchenko3D/writeDataIter3D.c
@@ -16,55 +16,53 @@
void name_ext(char *filename, char *extension);
-long writeDataIter3D(char *file_iter, float *data, segy *hdrs, long n1, long n2, long n3, float d2, float f2, long Nfoc, float *xsyn, float *ysyn, float *zsyn, long *ixpos, int *iypos, long npos, long t0shift, long iter)
+long writeDataIter3D(char *file_iter, float *data, segy *hdrs, long n1, long n2, long n3, long Nfoc, float *xsyn, float *ysyn, float *zsyn, long *ixpos, long *iypos, long npos, long t0shift, long iter)
{
FILE *fp_iter;
size_t nwrite;
- int i, l, j, k, iy, ret, tracf, size, ix;
+ int i, l, j, ret, tracf, size, ix, iy;
char number[16], filename[1024];
float *trace;
trace = (float *)malloc(n1*sizeof(float));
strcpy(filename, file_iter);
- sprintf(number,"_%03d",(iter+1));
+ sprintf(number,"_%03d",(iter));
name_ext(filename, number);
fp_iter = fopen(filename, "w+");
if (fp_iter==NULL) verr("error on creating output file %s", filename);
- tracf=1;
size=n1*n2*n3;
+ tracf=1;
for (l = 0; l < Nfoc; l++) {
- for (k = 0; k < n3; k++) {
- for (i = 0; i < n2; i++) {
+ for (i = 0; i < npos; i++) {
ix = ixpos[i]; /* select proper position */
- hdrs[k*n2+i].fldr = l+1;
- hdrs[k*n2+i].sx = NINT(xsyn[l]*1000);
- hdrs[k*n2+i].sy = NINT(ysyn[l]*1000);
- hdrs[k*n2+i].offset = (long)NINT((f2+i*d2) - xsyn[l]);
- hdrs[k*n2+i].tracf = tracf++;
- hdrs[k*n2+i].selev = NINT(zsyn[l]*1000);
- hdrs[k*n2+i].sdepth = NINT(-zsyn[l]*1000);
-
+ iy = iypos[i];
+ hdrs[i].fldr = l+1;
+ hdrs[i].sx = NINT(xsyn[l]*1000);
+ hdrs[i].sy = NINT(ysyn[l]*1000);
+ hdrs[i].tracf = tracf++;
+ hdrs[i].selev = NINT(zsyn[l]*1000);
+ hdrs[i].sdepth = NINT(-zsyn[l]*1000);
+
if (t0shift) {
- /* rotate to get t=0 in the middle */
- hdrs[k*n2+i].f1 = -n1*0.5*hdrs[k*n2+i].d1;
-
- for (j = 0; j < n1/2; j++) {
- trace[n1/2+j] = data[l*size+iy*n2*n2+ix*n1+j];
- }
- for (j = n1/2; j < n1; j++) {
- trace[j-n1/2] = data[l*size+iy*n2*n2+ix*n1+j];
+ /* rotate to get t=0 in the middle */
+ hdrs[i].f1 = -n1*0.5*hdrs[i].d1;
+
+ for (j = 0; j < n1/2; j++) {
+ trace[n1/2+j] = data[l*size+iy*n2*n1+ix*n1+j];
+ }
+ for (j = n1/2; j < n1; j++) {
+ trace[j-n1/2] = data[l*size+iy*n2*n1+ix*n1+j];
+ }
}
- }
else {
hdrs[i].f1 = 0.0;
- memcpy(&trace[0],&data[l*size+iy*n2*n2+ix*n1],n1*sizeof(float));
+ memcpy(&trace[0],&data[l*size+iy*n2*n1+ix*n1],n1*sizeof(float));
}
- nwrite = fwrite(&hdrs[k*n2+i], 1, TRCBYTES, fp_iter);
- assert(nwrite == TRCBYTES);
- nwrite = fwrite(trace, sizeof(float), n1, fp_iter);
- assert (nwrite == n1);
- }
- }
+ nwrite = fwrite(&hdrs[i], 1, TRCBYTES, fp_iter);
+ assert(nwrite == TRCBYTES);
+ nwrite = fwrite(trace, sizeof(float), n1, fp_iter);
+ assert (nwrite == n1);
+ }
}
ret = fclose(fp_iter);
if (ret < 0 ) verr("error on writing output file.");