From b55bc12a966160b73562a575c1e63680206d09e3 Mon Sep 17 00:00:00 2001
From: JBrackenhoff <J.A.Brackenhoff@tudelft.nl>
Date: Tue, 21 Apr 2020 17:23:46 +0200
Subject: [PATCH] plane_wave
---
fdelmodc3D/acoustic4_3D.c | 1 -
fdelmodc3D/applySource3D.c | 38 +-
fdelmodc3D/fdelmodc3D.h | 4 +
fdelmodc3D/getParameters3D.c | 88 ++-
fdelmodc3D/sourceOnSurface3D.c | 13 +-
marchenko/marchenko.c | 2 +-
marchenko3D/applyMute3D.c | 27 +-
marchenko3D/fmute3D.c | 10 +-
marchenko3D/homogeneousg3D.c | 5 +-
marchenko3D/marchenko3D.c | 433 +++++++---
marchenko3D/readShotData3Dzfp.c | 1 +
utils/HomG.c | 1300 ++++++++++++++++++++++---------
utils/Makefile | 24 +-
utils/getFileInfo3D.c | 2 +-
14 files changed, 1395 insertions(+), 553 deletions(-)
diff --git a/fdelmodc3D/acoustic4_3D.c b/fdelmodc3D/acoustic4_3D.c
index 6029d53..fb624ee 100644
--- a/fdelmodc3D/acoustic4_3D.c
+++ b/fdelmodc3D/acoustic4_3D.c
@@ -165,7 +165,6 @@ long acoustic4_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, lo
if (bnd.lef==2) mod.ioPx -= bnd.npml;
if (bnd.rig==2) mod.iePx += bnd.npml;
-
/* Add stress source */
if (src.type < 6) {
applySource3D(mod, src, wav, bnd, itime, ixsrc, iysrc, izsrc, vx, vy, vz, p, NULL, NULL, NULL, NULL, NULL, rox, roy, roz, l2m, src_nwav, verbose);
diff --git a/fdelmodc3D/applySource3D.c b/fdelmodc3D/applySource3D.c
index 74e88f3..90f545d 100644
--- a/fdelmodc3D/applySource3D.c
+++ b/fdelmodc3D/applySource3D.c
@@ -25,7 +25,7 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l
float ***rox, float ***roy, float ***roz, float ***l2m, float **src_nwav, long verbose)
{
long is0, ibndz, ibndy, ibndx;
- long isrc, ix, iy, iz, n1, n2;
+ long isrc, ix, iy, iz, n1, n2, ix0, iy0, ixe, iye;
long id1, id2, id3;
float src_ampl, time, scl, dt, sdx;
float Mxx, Myy, Mzz, Mxz, Myz, Mxy;
@@ -98,9 +98,19 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l
iy = src.y[isrc] + ibndy;
iz = src.z[isrc] + ibndz;
}
- else { /* plane wave and point sources */
- ix = ixsrc + ibndx + is0 + isrc;
- iy = iysrc + ibndy + is0 + isrc;
+ else if (src.plane) {/* plane wave sources */
+ ix = ixsrc + ibndx + src.x[isrc];
+ iy = iysrc + ibndy + src.y[isrc];
+ iz = izsrc + ibndz + src.z[isrc];
+ ix0 = ixsrc + ibndx + src.x[0];
+ ixe = ixsrc + ibndx + src.x[src.n-1];
+ iy0 = iysrc + ibndy + src.y[0];
+ iye = iysrc + ibndy + src.y[src.n-1];
+ // vmess("ix0 %li ixe %li iy0 %li iye %li",ix0,ixe,iy0,iye);
+ }
+ else { /* point sources */
+ ix = ixsrc + ibndx + isrc;
+ iy = iysrc + ibndy + isrc;
iz = izsrc + ibndz;
}
time = itime*dt - src.tbeg[isrc];
@@ -111,7 +121,7 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l
/* delay not reached or no samples left in source wavelet? */
if ( (time < 0.0) || ( (itime*dt) >= src.tend[isrc]) ) continue;
-// fprintf(stderr,"isrc=%li ix=%li iz=%li src.x=%li src.z=%li\n", isrc, ix, iz, src.x[isrc], src.z[isrc]);
+ // fprintf(stderr,"isrc=%li ix=%li iy=%li iz=%li src.x=%li src.y=%li src.z=%li\n", isrc, ix, iy, iz, src.x[isrc], src.y[isrc], src.z[isrc]);
if (!src.multiwav) { /* only one wavelet for all sources */
src_ampl = src_nwav[0][id1]*(id2-time/dt) + src_nwav[0][id2]*(time/dt-id1);
@@ -128,15 +138,19 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l
}
/* cosine squared windowing to reduce edge effects on shot arrays */
- if ( (src.n>1) && src.window) {
- scl = 1.0;
- if (isrc < src.window) {
- scl = cos(0.5*M_PI*(src.window - isrc)/src.window);
+ if (src.plane) {
+ if (src.nxwindow > 0){
+ scl = 1.0;
+ if (ix-ix0 < src.nxwindow) scl = cos(0.5*M_PI*(src.nxwindow - (ix-ix0))/src.nxwindow);
+ else if (ixe-ix < src.nxwindow) scl = cos(0.5*M_PI*(src.nxwindow - (ixe-ix))/src.nxwindow);
+ src_ampl *= scl*scl;
}
- else if (isrc > src.n-src.window+1) {
- scl = cos(0.5*M_PI*(src.window - (src.n-isrc+1))/src.window);
+ if (src.nywindow > 0){
+ scl = 1.0;
+ if (iy-iy0 < src.nywindow) scl = cos(0.5*M_PI*(src.nywindow - (iy-iy0))/src.nywindow);
+ else if (iye-iy < src.nywindow) scl = cos(0.5*M_PI*(src.nywindow - (iye-iy))/src.nywindow);
+ src_ampl *= scl*scl;
}
- src_ampl *= scl*scl;
}
/* source scaling factor to compensate for discretisation */
diff --git a/fdelmodc3D/fdelmodc3D.h b/fdelmodc3D/fdelmodc3D.h
index ed4ce7b..061df2f 100644
--- a/fdelmodc3D/fdelmodc3D.h
+++ b/fdelmodc3D/fdelmodc3D.h
@@ -161,6 +161,8 @@ typedef struct _waveletPar { /* Wavelet Parameters */
typedef struct _sourcePar { /* Source Array Parameters */
long n;
+ long nx;
+ long ny;
long type;
long orient;
long *z;
@@ -181,6 +183,8 @@ typedef struct _sourcePar { /* Source Array Parameters */
float strike;
float rake;
long distribution;
+ long nxwindow;
+ long nywindow;
long window;
long injectionrate;
long sinkdepth;
diff --git a/fdelmodc3D/getParameters3D.c b/fdelmodc3D/getParameters3D.c
index 8683bc5..6a0cf32 100644
--- a/fdelmodc3D/getParameters3D.c
+++ b/fdelmodc3D/getParameters3D.c
@@ -42,6 +42,7 @@ long getParameters3D(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar
{
long isnapmax1, isnapmax2, isnapmax, sna_nrsna;
long n1, n2, n3, nx, ny, nz, nsrc, ix, axis, ioPz, is0, optn;
+ long npxsrc, npysrc, isx0, isy0, isx, isy;
long idzshot, idxshot, idyshot, nsrctext;
long src_ix0, src_iy0, src_iz0, src_ix1, src_iy1, src_iz1;
long disable_check;
@@ -54,7 +55,7 @@ long getParameters3D(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar
float tsnap1, tsnap2, dtsnap, dxsnap, dysnap, dzsnap, dtrcv;
float xsnap1, xsnap2, ysnap1, ysnap2, zsnap1, zsnap2, xmax, ymax, zmax;
float xsrc1, xsrc2, ysrc1, ysrc2, zsrc1, zsrc2, tsrc1, tsrc2, tlength, tactive;
- float src_angle, src_velo, p, grad2rad, rdelay, scaledt;
+ float src_anglex, src_angley, src_velox, src_veloy, px, py, grad2rad, rdelay, scaledt;
float *xsrca, *ysrca, *zsrca, rrcv;
float rsrc, oxsrc, oysrc, ozsrc, dphisrc, ncsrc;
size_t nsamp;
@@ -684,10 +685,13 @@ criteria we have imposed.*/
/* number of sources per shot modeling */
-
- if (!getparlong("src_window",&src->window)) src->window=0;
- if (!getparfloat("src_angle",&src_angle)) src_angle=0.;
- if (!getparfloat("src_velo",&src_velo)) src_velo=1500.;
+ if (!getparlong("src_nxwindow",&src->nxwindow)) src->nxwindow=0;
+ if (!getparlong("src_nywindow",&src->nywindow)) src->nywindow=0;
+ src->window = 0;
+ if (!getparfloat("src_anglex",&src_anglex)) src_anglex=0.;
+ if (!getparfloat("src_angley",&src_angley)) src_angley=0.;
+ if (!getparfloat("src_velox",&src_velox)) src_velox=1500.;
+ if (!getparfloat("src_veloy",&src_veloy)) src_veloy=1500.;
if (!getparlong("distribution",&src->distribution)) src->distribution=0;
if (!getparlong("src_multiwav",&src->multiwav)) src->multiwav=0;
if (!getparfloat("amplitude", &src->amplitude)) src->amplitude=0.0;
@@ -898,14 +902,26 @@ criteria we have imposed.*/
free(selev);
}
else {
- if (src->plane) { if (!getparlong("nsrc",&nsrc)) nsrc=1;}
- else nsrc=1;
+ if (src->plane) {
+ if (!getparlong("npxsrc",&npxsrc)) npxsrc=1;
+ if (!getparlong("npysrc",&npysrc)) npysrc=1;
+ }
+ else {
+ npxsrc=1;
+ npysrc=1;
+ }
- if (nsrc > nx) {
- vwarn("Number of sources used in plane wave is larger than ");
- vwarn("number of gridpoints in X. Plane wave will be clipped to the edges of the model");
- nsrc = mod->nx;
+ if (npxsrc > nx) {
+ vwarn("Number of sources used in plane wave (%li) is larger than ",npxsrc);
+ vwarn("number of gridpoints in X (%li). Plane wave will be clipped to the edges of the model",nx);
+ npxsrc = mod->nx;
}
+ if (npysrc > ny) {
+ vwarn("Number of sources used in plane wave (%li) is larger than ",npysrc);
+ vwarn("number of gridpoints in Y (%li). Plane wave will be clipped to the edges of the model",ny);
+ npysrc = mod->ny;
+ }
+ nsrc = npxsrc*npysrc;
/* for a source defined on mutliple gridpoint calculate p delay factor */
@@ -915,25 +931,47 @@ criteria we have imposed.*/
src->tbeg = (float *)malloc(nsrc*sizeof(float));
src->tend = (float *)malloc(nsrc*sizeof(float));
grad2rad = 17.453292e-3;
- p = sin(src_angle*grad2rad)/src_velo;
- if (p < 0.0) {
- for (is=0; is<nsrc; is++) {
- src->tbeg[is] = fabsf((nsrc-is-1)*dx*p);
+ px = sin(src_anglex*grad2rad)/src_velox;
+ py = sin(src_angley*grad2rad)/src_veloy;
+ if (py < 0.0) {
+ for (isy=0; isy<npysrc; isy++) {
+ if (px < 0.0) {
+ for (isx=0; isx<npxsrc; isx++) {
+ src->tbeg[isy*npxsrc+isx] = fabsf((npysrc-isy-1)*dy*py) + fabsf((npxsrc-isx-1)*dx*px);
+ }
+ }
+ else {
+ for (isx=0; isx<npxsrc; isx++) {
+ src->tbeg[isy*npxsrc+isx] = fabsf((npysrc-isy-1)*dy*py) + isx*dx*px;
+ }
+ }
}
}
else {
- for (is=0; is<nsrc; is++) {
- src->tbeg[is] = is*dx*p;
+ for (isy=0; isy<npysrc; isy++) {
+ if (px < 0.0) {
+ for (isx=0; isx<npxsrc; isx++) {
+ src->tbeg[isy*npxsrc+isx] = isy*dy*py + fabsf((npxsrc-isx-1)*dx*px);
+ }
+ }
+ else {
+ for (isx=0; isx<npxsrc; isx++) {
+ src->tbeg[isy*npxsrc+isx] = isy*dy*py + isx*dx*px;
+ }
+ }
}
}
for (is=0; is<nsrc; is++) {
src->tend[is] = src->tbeg[is] + (wav->nt-1)*wav->dt;
}
- is0 = -1*floor((nsrc-1)/2);
- for (is=0; is<nsrc; is++) {
- src->x[is] = is0 + is;
- src->y[is] = 0;
- src->z[is] = 0;
+ isx0 = -1*floor((npxsrc-1)/2);
+ isy0 = -1*floor((npysrc-1)/2);
+ for (isy=0; isy<npysrc; isy++) {
+ for (isx=0; isx<npxsrc; isx++) {
+ src->x[isy*npxsrc+isx] = isx0 + isx;
+ src->y[isy*npxsrc+isx] = isy0 + isy;
+ src->z[isy*npxsrc+isx] = 0;
+ }
}
if (wav->random) {
@@ -962,6 +1000,8 @@ criteria we have imposed.*/
wav->nsamp[nsrc] = nsamp; /* put total number of samples in last position */
wav->nst = nsamp; /* put total number of samples in nst part */
}
+ src->nx = npxsrc;
+ src->ny = npysrc;
}
if (src->multiwav) {
@@ -1002,9 +1042,9 @@ criteria we have imposed.*/
vmess("*******************************************");
vmess("*********** source array info *************");
vmess("*******************************************");
- vmess("Areal source array is defined with %li sources.",nsrc);
+ vmess("Areal source array is defined with %li sources (x=%li, y=%li).",nsrc,npxsrc,npysrc);
vmess("Memory requirement for sources = %.2f MB.",sizeof(float)*(nsamp/(1024.0*1024.0)));
- if (src->plane) vmess("Computed p-value = %f.",p);
+ if (src->plane) vmess("Computed px-value = %f. and py-value = %f",px,py);
}
if (src->random) {
vmess("Sources are placed at random locations in domain: ");
diff --git a/fdelmodc3D/sourceOnSurface3D.c b/fdelmodc3D/sourceOnSurface3D.c
index 4d8e1c3..b577199 100644
--- a/fdelmodc3D/sourceOnSurface3D.c
+++ b/fdelmodc3D/sourceOnSurface3D.c
@@ -90,10 +90,15 @@ long storeSourceOnSurface3D(modPar mod, srcPar src, bndPar bnd,
iys = src.y[isrc] + ibndy;
izs = src.z[isrc] + ibndz;
}
- else { /* plane wave and point sources */
- ixs = ixsrc + ibndx + is0 + isrc;
- iys = iysrc + ibndy + is0 + isrc;
- izs = izsrc + ibndz;
+ else if (src.plane) {/* plane wave sources */
+ ixs = ixsrc + ibndx + src.x[isrc];
+ iys = iysrc + ibndy + src.y[isrc];
+ izs = izsrc + ibndz + src.z[isrc];
+ }
+ else { /* point sources */
+ ixs = ixsrc + ibndx + isrc;
+ iys = iysrc + ibndy + isrc;
+ izs = izsrc + ibndz;
}
/* check if there are source positions placed on the boundaries.
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index e947263..673fe02 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -268,7 +268,7 @@ int main (int argc, char **argv)
tsynW[i] = NINT(i*dxs*p/dt);
}
}
- if (Nfoc!=1) verr("For plave-wave focusing only one function can be computed at the same time");
+ if (Nfoc!=1) verr("For plane-wave focusing only one function can be computed at the same time");
}
else { /* just fill with zero's */
for (i=0; i<nxs*Nfoc; i++) {
diff --git a/marchenko3D/applyMute3D.c b/marchenko3D/applyMute3D.c
index 03097a2..222117a 100644
--- a/marchenko3D/applyMute3D.c
+++ b/marchenko3D/applyMute3D.c
@@ -12,11 +12,12 @@
#endif
#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
-void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, long nxs, long nys, long nt, long *ixpos, long *iypos, long npos, long shift)
+void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, long nxs, long nys, long nt,
+ long *ixpos, long *iypos, long npos, long shift, long *tsynW)
{
long ix, iy, i, j, l, isyn, nxys;
float *costaper, scl;
- long imute, tmute;
+ long imute, tmute, ts;
nxys = nxs*nys;
@@ -33,10 +34,11 @@ void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, l
for (i = 0; i < npos; i++) {
imute = iypos[i]*nxs+ixpos[i];
tmute = mute[isyn*nxys+imute];
- for (j = 0; j < MAX(0,tmute-shift-smooth); j++) {
+ ts = tsynW[isyn*nxys+imute];
+ for (j = 0; j < MAX(0,-2*ts+tmute-shift-smooth); j++) {
data[isyn*nxys*nt+i*nt+j] = 0.0;
}
- for (j = MAX(0,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); j++,l++) {
+ for (j = MAX(0,-2*ts+tmute-shift-smooth),l=0; j < MAX(0,-2*ts+tmute-shift); j++,l++) {
data[isyn*nxys*nt+i*nt+j] *= costaper[smooth-l-1];
}
}
@@ -45,14 +47,15 @@ void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, l
for (i = 0; i < npos; i++) {
imute = iypos[i]*nxs+ixpos[i];
tmute = mute[isyn*nxys+imute];
+ ts = tsynW[isyn*nxys+imute];
if (tmute >= nt/2) {
memset(&data[isyn*nxys*nt+i*nt],0, sizeof(float)*nt);
continue;
}
- for (j = MAX(0,tmute-shift),l=0; j < MAX(0,tmute-shift+smooth); j++,l++) {
+ for (j = MAX(0,-2*ts+tmute-shift),l=0; j < MAX(0,-2*ts+tmute-shift+smooth); j++,l++) {
data[isyn*nxys*nt+i*nt+j] *= costaper[l];
}
- for (j = MAX(0,tmute-shift+smooth)+1; j < MIN(nt,nt+1-tmute+shift-smooth); j++) {
+ for (j = MAX(0,-2*ts+tmute-shift+smooth)+1; j < MIN(nt,nt+1-tmute+shift-smooth); j++) {
data[isyn*nxys*nt+i*nt+j] = 0.0;
}
for (j = MIN(nt,nt-tmute+shift-smooth),l=0; j < MIN(nt,nt-tmute+shift); j++,l++) {
@@ -64,14 +67,15 @@ void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, l
for (i = 0; i < npos; i++) {
imute = iypos[i]*nxs+ixpos[i];
tmute = mute[isyn*nxys+imute];
+ ts = tsynW[isyn*nxys+imute];
if (tmute >= nt/2) {
memset(&data[isyn*nxys*nt+i*nt],0, sizeof(float)*nt);
continue;
}
- for (j = MAX(0,tmute+shift),l=0; j < MAX(0,tmute+shift+smooth); j++,l++) {
+ for (j = MAX(0,-2*ts+tmute+shift),l=0; j < MAX(0,-2*ts+tmute+shift+smooth); j++,l++) {
data[isyn*nxys*nt+i*nt+j] *= costaper[l];
}
- for (j = MAX(0,tmute+shift+smooth)+1; j < MIN(nt,nt+1-tmute-shift-smooth); j++) {
+ for (j = MAX(0,-2*ts+tmute+shift+smooth)+1; j < MIN(nt,nt+1-tmute-shift-smooth); j++) {
data[isyn*nxys*nt+i*nt+j] = 0.0;
}
for (j = MIN(nt,nt-tmute-shift-smooth),l=0; j < MIN(nt,nt-tmute-shift); j++,l++) {
@@ -83,10 +87,11 @@ void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, l
for (i = 0; i < npos; i++) {
imute = iypos[i]*nxs+ixpos[i];
tmute = mute[isyn*nxys+imute];
- for (j = MAX(0,tmute-shift),l=0; j < MAX(0,tmute-shift+smooth); j++,l++) {
+ ts = tsynW[isyn*nxys+imute];
+ for (j = MAX(0,ts+tmute-shift),l=0; j < MAX(0,ts+tmute-shift+smooth); j++,l++) {
data[isyn*nxys*nt+i*nt+j] *= costaper[l];
}
- for (j = MAX(0,tmute-shift+smooth); j < nt; j++) {
+ for (j = MAX(0,ts+tmute-shift+smooth); j < nt; j++) {
data[isyn*nxys*nt+i*nt+j] = 0.0;
}
}
@@ -95,6 +100,7 @@ void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, l
for (i = 0; i < npos; i++) {
imute = iypos[i]*nxs+ixpos[i];
tmute = mute[isyn*nxys+imute];
+ ts = tsynW[isyn*nxys+imute];
for (j = MAX(0,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); j++,l++) {
data[isyn*nxys*nt+i*nt+j] *= costaper[smooth-l-1];
}
@@ -113,6 +119,7 @@ void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, l
for (i = 0; i < npos; i++) {
imute = iypos[i]*nxs+ixpos[i];
tmute = mute[isyn*nxys+imute];
+ ts = tsynW[isyn*nxys+imute];
for (j = 0; j < MAX(0,tmute-shift-smooth); j++) {
data[isyn*nxys*nt+i*nt+j] = 0.0;
}
diff --git a/marchenko3D/fmute3D.c b/marchenko3D/fmute3D.c
index b8fea43..5715a5d 100644
--- a/marchenko3D/fmute3D.c
+++ b/marchenko3D/fmute3D.c
@@ -26,7 +26,8 @@ long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, fl
long readData3D(FILE *fp, float *data, segy *hdrs, long n1);
long writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2);
long disp_fileinfo3D(char *file, long n1, long n2, long n3, float f1, float f2, float f3, float d1, float d2, float d3, segy *hdrs);
-void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, long nxs, long nys, long nt, long *ixpos, long *iypos, long npos, long shift);
+void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, long nxs, long nys, long nt,
+ long *ixpos, long *iypos, long npos, long shift, long *tsynW);
double wallclock_time(void);
/*********************** self documentation **********************/
@@ -66,7 +67,7 @@ int main (int argc, char **argv)
FILE *fp_in1, *fp_in2, *fp_out, *fp_chk, *fp_psline1, *fp_psline2;
long verbose, shift, k, nx1, ny1, nt1, nx2, ny2, nt2, nxy;
long ntmax, nxmax, nymax, ret, i, j, l, jmax, ixmax, iymax, above, check;
- long size, ntraces, ngath, *maxval, hw, smooth;
+ long size, ntraces, ngath, *maxval, *tsynW, hw, smooth;
long tstart, tend, scale, *xrcv, *yrcv;
float dt, dt1, dx1, dy1, ft1, fx1, fy1, t0, t1, dt2, dx2, dy2, ft2, fx2, fy2;
float w1, w2, dxrcv, dyrcv;
@@ -176,6 +177,7 @@ int main (int argc, char **argv)
nxy = nx1*ny1;
maxval = (long *)calloc(nxy,sizeof(long));
+ tsynW = (long *)calloc(nxy,sizeof(long));
xrcv = (long *)calloc(nxy,sizeof(long));
yrcv = (long *)calloc(nxy,sizeof(long));
@@ -236,7 +238,7 @@ int main (int argc, char **argv)
dxrcv = (hdrs_in1[nxy-1].gx - hdrs_in1[0].gx)*sclsxgx/(float)(nx1-1);
dyrcv = (hdrs_in1[nxy-1].gy - hdrs_in1[0].gy)*sclsxgx/(float)(ny1-1);
ixmax = NINT(((hdrs_in1[0].sx-hdrs_in1[0].gx)*sclsxgx)/dxrcv);
- iymax = NINT(((hdrs_in1[0].sy-hdrs_in1[0].gy)*sclsxgx)/dxrcv);
+ iymax = NINT(((hdrs_in1[0].sy-hdrs_in1[0].gy)*sclsxgx)/dyrcv);
if (iymax > ny1-1) {
vmess("source of y is past array, snapping to nearest y");
iymax = ny1-1;
@@ -416,7 +418,7 @@ int main (int argc, char **argv)
/*================ apply mute window ================*/
- applyMute3D(tmpdata2, maxval, smooth, above, 1, nx2, ny2, nt2, xrcv, yrcv, nx2*ny2, shift);
+ applyMute3D(tmpdata2, maxval, smooth, above, 1, nx2, ny2, nt2, xrcv, yrcv, nx2*ny2, shift, tsynW);
/*================ write result to output file ================*/
diff --git a/marchenko3D/homogeneousg3D.c b/marchenko3D/homogeneousg3D.c
index 6f3d0f7..b0ab228 100644
--- a/marchenko3D/homogeneousg3D.c
+++ b/marchenko3D/homogeneousg3D.c
@@ -197,7 +197,7 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *f1p, float *f1
count+=1;
zsrc = hdr_inp[0].selev;
- if (verbose>2) vmess("Creating Homogeneous G at location %li out of %li",l,Nfoc);
+ if (verbose>2) vmess("Creating Homogeneous G at location %li out of %li",l+1,Nfoc);
if (scheme==0) { //Marchenko representation with G source
for (k = 0; k < ny; k++) {
@@ -365,8 +365,7 @@ void homogeneousg3D(float *HomG, float *green, float *f2p, float *f1p, float *f1
free(tmp1);
free(tmp2);
}
- if (scheme==6) free(tmp1);
- if (scheme==7) free(tmp1);
+ if (scheme==6 || scheme==8 || scheme==9 || scheme==10) free(tmp1);
}
if (scheme==5) {
free(input);
diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index 6581add..665b680 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -13,6 +13,8 @@
#include <math.h>
#include <assert.h>
#include <genfft.h>
+#include "zfpmar.h"
+#include <zfp.h>
int omp_get_max_threads(void);
int omp_get_num_threads(void);
@@ -46,7 +48,8 @@ void name_ext(char *filename, char *extension);
void convol(float *data1, float *data2, float *con, long nrec, long nsam, float dt, long shift);
-void applyMute3D(float *data, long *mute, long smooth, long above, long Nfoc, long nxs, long nys, long nt, long *xrcvsyn, long *yrcvsyn, long npos, long shift);
+void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, long nxs, long nys, long nt,
+ long *ixpos, long *iypos, long npos, long shift, long *tsynW);
long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath,
float *d1, float *d2, float *d3, float *f1, float *f2, float *f3,
@@ -62,6 +65,8 @@ long writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2);
long disp_fileinfo3D(char *file, long n1, long n2, long n3, float f1, float f2, float f3, float d1, float d2, float d3, segy *hdrs);
double wallclock_time(void);
+long zfpcompress(float* data, long nx, long ny, long nz, double tolerance, zfpmar zfpm, FILE *file);
+
void AmpEst3D(float *f1d, float *Gd, float *ampest, long Nfoc, long nxs, long nys, long ntfft, long *ixpos, long *iypos, long npos,
char *file_wav, float dx, float dy, float dt);
@@ -120,7 +125,12 @@ 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"
+" MUTE-WINDOW ",
+" plane_wave=0 ............. enable plane-wave illumination function",
+" src_anglex=0 ............. angle of the plane wave in the x-direction",
+" src_angley=0 ............. angle of the plane wave in the y-direction",
+" src_velox=0 .............. velocity of the plane wave in the x-direction",
+" src_veloy=0 .............. velocity of the plane wave in the y-direction",
" REFLECTION RESPONSE CORRECTION ",
" scale=2 .................. scale factor of R for summation of Ni with G_d (only for time shot data)",
" pad=0 .................... amount of samples to pad the reflection series",
@@ -151,18 +161,18 @@ char *sdoc[] = {
" file_gmin= ............... output file with G- ",
" file_f1plus= ............. output file with f1+ ",
" file_f1min= .............. output file with f1- ",
-" file_f2= ................. output file with f2 (=p+) ",
-" file_pplus= .............. output file with p+ ",
-" file_pmin= ............... output file with p- ",
+" file_f2= ................. output file with f2 ",
" file_ampscl= ............. output file with estimated amplitudes ",
" file_iter= ............... output file with -Ni(-t) for each iteration",
" compact=0 ................ Write out homg and imag in compact format",
" .......................... WARNING! This write-out cannot be displayed with SU",
+" zfp=0 .................... Write out the standard output in compressed zfp format",
+" tolerance=1e-7 ........... accuracy of the zfp compression,",
" verbose=0 ................ silent option; >0 displays info",
" ",
" ",
-" author : Jan Thorbecke : 2016 (j.w.thorbecke@tudelft.nl)",
-" author : Joeri Brackenhoff : 2019 (j.a.brackenhoff@tudelft.nl)",
+" author : Jan Thorbecke : 2016 (j.w.thorbecke@tudelft.nl)",
+" author 3D: Joeri Brackenhoff : 2019 (j.a.brackenhoff@tudelft.nl)",
" ",
NULL};
/**************** end self doc ***********************************/
@@ -170,27 +180,30 @@ NULL};
int main (int argc, char **argv)
{
FILE *fp_out, *fp_f1plus, *fp_f1min, *fp_imag, *fp_homg;
- FILE *fp_gmin, *fp_gplus, *fp_f2, *fp_pmin, *fp_amp;
+ FILE *fp_gmin, *fp_gplus, *fp_f2, *fp_amp;
long i, j, l, k, ret, nshots, nxshot, nyshot, Nfoc, nt, nx, ny, nts, nxs, nys, ngath;
- long size, n1, n2, n3, ntap, tap, dxi, dyi, ntraces, pad, *sx, *sy, *sz;
+ long size, n1, n2, n3, ntap, ntapx, ntapy, tap, dxi, dyi, ntraces, pad, *sx, *sy, *sz;
long nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
- long reci, mode, n2out, n3out, verbose, ntfft;
+ long reci, mode, n2out, n3out, verbose, ntfft, zfp;
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 *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;
+ float fmin, fmax, *tapershx, *tapershy, *tapersy, *tapersx, fxf, fyf, dxf, dyf, *xsrc, *ysrc, *xrcv, *yrcv, *zsyn, *zsrc, *xrcvsyn, *yrcvsyn;
+ double t0, t1, t2, t3, tsyn, tread, tfft, tcopy, energyNi, *energyN0, tolerance;
float d1, d2, d3, f1, f2, f3, fxsb, fxse, fysb, fyse, ft, fx, fy, *xsyn, *ysyn, dxsrc, dysrc;
- float *green, *f2p, *pmin, *G_d, dt, dx, dy, dxs, dys, scl, mem;
+ float *green, *f2p, *G_d, dt, dx, dy, dxs, dys, scl, mem;
float *f1plus, *f1min, *iRN, *Ni, *Nig, *trace, *Gmin, *Gplus, *HomG;
float scale, *tmpdata;
float *ixmask, *iymask, *ampscl, *Gd, *Image, dzim;
- float grad2rad, p, src_angle, src_velo, *mutetest;
+ float grad2rad, px, py, src_anglex, src_angley, src_velox, src_veloy, *mutetest;
complex *Refl, *Fop;
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_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_inp;
char *file_ray, *file_amp, *file_wav, *file_shotw, *file_shotzfp;
segy *hdrs_out, *hdrs_Nfoc, *hdrs_iter;
+ zfptop zfpt;
+ zfpmar zfpm;
+ size_t nread;
initargs(argc, argv);
requestdoc(1);
@@ -209,9 +222,7 @@ int main (int argc, char **argv)
if (!getparstring("file_f1min", &file_f1min)) file_f1min = NULL;
if (!getparstring("file_gplus", &file_gplus)) file_gplus = NULL;
if (!getparstring("file_gmin", &file_gmin)) file_gmin = NULL;
- if (!getparstring("file_pplus", &file_f2)) file_f2 = NULL;
if (!getparstring("file_f2", &file_f2)) file_f2 = NULL;
- if (!getparstring("file_pmin", &file_pmin)) file_pmin = NULL;
if (!getparstring("file_iter", &file_iter)) file_iter = NULL;
if (!getparstring("file_wav", &file_wav)) file_wav = NULL;
if (!getparstring("file_imag", &file_imag)) file_imag = NULL;
@@ -232,8 +243,12 @@ int main (int argc, char **argv)
if (!getparfloat("scale", &scale)) scale = 2.0;
if (!getparlong("tap", &tap)) tap = 0;
if (!getparlong("ntap", &ntap)) ntap = 0;
+ if (!getparlong("ntapx", &ntapx)) ntapx = 0;
+ if (!getparlong("ntapy", &ntapy)) ntapy = 0;
if (!getparlong("pad", &pad)) pad = 0;
if (!getparlong("ampest", &est)) ampest = 0;
+ if (!getparlong("zfp", &zfp)) zfp = 0;
+ if (!getpardouble("tolerance", &tolerance)) tolerance = 1e-7;
if(!getparlong("niter", &niter)) niter = 10;
if(!getparlong("hw", &hw)) hw = 8;
@@ -243,8 +258,10 @@ int main (int argc, char **argv)
if(!getparlong("compact", &compact)) compact=0;
if (!getparlong("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 (!getparfloat("src_anglex",&src_anglex)) src_anglex=0.;
+ if (!getparfloat("src_velox",&src_velox)) src_velox=1500.;
+ if (!getparfloat("src_angley",&src_angley)) src_angley=0.;
+ if (!getparfloat("src_veloy",&src_veloy)) src_veloy=1500.;
if (reci && ntap) vwarn("tapering influences the reciprocal result");
@@ -313,7 +330,6 @@ int main (int argc, char **argv)
Fop = (complex *)calloc(nys*nxs*nw*Nfoc,sizeof(complex));
green = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
f2p = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
- pmin = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
f1plus = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
f1min = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
iRN = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
@@ -323,7 +339,8 @@ int main (int argc, char **argv)
muteW = (long *)calloc(Nfoc*nys*nxs,sizeof(long));
tsynW = (long *)malloc(Nfoc*nys*nxs*sizeof(long)); // time-shift for Giovanni's plane-wave on non-zero times
trace = (float *)malloc(ntfft*sizeof(float));
- tapersy = (float *)malloc(nxs*sizeof(float));
+ tapersx = (float *)malloc(nxs*sizeof(float));
+ tapersy = (float *)malloc(nys*sizeof(float));
xrcvsyn = (float *)calloc(Nfoc*nys*nxs,sizeof(float)); // x-rcv postions of focal points
yrcvsyn = (float *)calloc(Nfoc*nys*nxs,sizeof(float)); // x-rcv postions of focal points
xsyn = (float *)malloc(Nfoc*sizeof(float)); // x-src position of focal points
@@ -335,7 +352,6 @@ int main (int argc, char **argv)
energyN0= (double *)calloc(Nfoc,sizeof(double)); // minimum energy for each focal position
Refl = (complex *)malloc(nw*ny*nx*nshots*sizeof(complex));
- tapersh = (float *)malloc(nx*sizeof(float));
xrcv = (float *)calloc(nshots*ny*nx,sizeof(float)); // x-rcv postions of shots
yrcv = (float *)calloc(nshots*ny*nx,sizeof(float)); // x-rcv postions of shots
xsrc = (float *)calloc(nshots,sizeof(float)); //x-src position of shots
@@ -375,18 +391,37 @@ int main (int argc, char **argv)
/* 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 (Nfoc!=1) verr("For plave-wave focusing only one function can be computed at the same time");
+ px = sin(src_anglex*grad2rad)/src_velox;
+ py = sin(src_angley*grad2rad)/src_veloy;
+ if (py < 0.0) {
+ for (j=0; j<nys; j++) {
+ if (px < 0.0) {
+ for (i=0; i<nxs; i++) {
+ tsynW[j*nxs+i] = NINT(fabsf((nxs-1-i)*dxs*px)/dt) + NINT(fabsf((nys-1-j)*dys*py)/dt);
+ }
+ }
+ else {
+ for (i=0; i<nxs; i++) {
+ tsynW[j*nxs+i] = NINT(i*dxs*px/dt) + NINT(fabsf((nys-1-j)*dys*py)/dt);
+ }
+ }
+ }
+ }
+ else {
+ for (j=0; j<nys; j++) {
+ if (px < 0.0) {
+ for (i=0; i<nxs; i++) {
+ tsynW[j*nxs+i] = NINT(fabsf((nxs-1-i)*dxs*px)/dt) + NINT(j*dys*py/dt);
+ }
+ }
+ else {
+ for (i=0; i<nxs; i++) {
+ tsynW[j*nxs+i] = NINT(i*dxs*px/dt) + NINT(j*dys*py/dt);
+ }
+ }
+ }
+ }
+ // if (Nfoc!=1) verr("For plane-wave focusing only one function can be computed at the same time");
}
else { /* just fill with zero's */
for (i=0; i<nys*nxs*Nfoc; i++) {
@@ -396,23 +431,30 @@ int main (int argc, char **argv)
/* define tapers to taper edges of acquisition */
if (tap == 1 || tap == 3) {
- for (j = 0; j < ntap; j++)
- tapersy[j] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
- for (j = ntap; j < nxs-ntap; j++)
+ for (j = 0; j < ntapx; j++)
+ tapersx[j] = (cos(PI*(j-ntapx)/ntapx)+1)/2.0;
+ for (j = ntapx; j < nxs-ntapx; j++)
+ tapersx[j] = 1.0;
+ for (j = nxs-ntapx; j < nxs; j++)
+ tapersx[j] =(cos(PI*(j-(nxs-ntapx))/ntapx)+1)/2.0;
+ for (j = 0; j < ntapy; j++)
+ tapersy[j] = (cos(PI*(j-ntapy)/ntapy)+1)/2.0;
+ for (j = ntapy; j < nys-ntapy; j++)
tapersy[j] = 1.0;
- for (j = nxs-ntap; j < nxs; j++)
- tapersy[j] =(cos(PI*(j-(nxs-ntap))/ntap)+1)/2.0;
+ for (j = nys-ntapy; j < nys; j++)
+ tapersy[j] =(cos(PI*(j-(nys-ntapy))/ntapy)+1)/2.0;
}
else {
- for (j = 0; j < nxs; j++) tapersy[j] = 1.0;
+ for (j = 0; j < nxs; j++) tapersx[j] = 1.0;
+ for (j = 0; j < nys; j++) tapersy[j] = 1.0;
}
if (tap == 1 || tap == 3) {
- if (verbose) vmess("Taper for operator applied ntap=%li", ntap);
+ if (verbose) vmess("Taper for operator applied nxtap=%li nytap=%li",ntapx,ntapy);
for (l = 0; l < Nfoc; l++) {
for (k = 0; k < nys; k++) {
for (i = 0; i < nxs; i++) {
for (j = 0; j < nts; j++) {
- G_d[l*nys*nxs*nts+k*nxs*nts+i*nts+j] *= tapersy[i];
+ G_d[l*nys*nxs*nts+k*nxs*nts+i*nts+j] *= tapersx[i]*tapersy[k];
}
}
}
@@ -473,32 +515,40 @@ int main (int argc, char **argv)
}
mode=1;
- tapersh = (float *)malloc(nx*sizeof(float));
+ tapershx = (float *)malloc(nx*sizeof(float));
+ tapershy = (float *)malloc(ny*sizeof(float));
if (tap == 2 || tap == 3) {
- for (j = 0; j < ntap; j++)
- tapersh[j] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
- for (j = ntap; j < nx-ntap; j++)
- tapersh[j] = 1.0;
- for (j = nx-ntap; j < nx; j++)
- tapersh[j] =(cos(PI*(j-(nx-ntap))/ntap)+1)/2.0;
+ for (j = 0; j < ntapx; j++)
+ tapershx[j] = (cos(PI*(j-ntapx)/ntapx)+1)/2.0;
+ for (j = ntapx; j < nx-ntapx; j++)
+ tapershx[j] = 1.0;
+ for (j = nx-ntapx; j < nx; j++)
+ tapershx[j] =(cos(PI*(j-(nx-ntapx))/ntapx)+1)/2.0;
+ for (j = 0; j < ntapy; j++)
+ tapershy[j] = (cos(PI*(j-ntapy)/ntapy)+1)/2.0;
+ for (j = ntapy; j < ny-ntapy; j++)
+ tapershy[j] = 1.0;
+ for (j = ny-ntapy; j < ny; j++)
+ tapershy[j] =(cos(PI*(j-(ny-ntapy))/ntapy)+1)/2.0;
}
else {
- for (j = 0; j < nx; j++) tapersh[j] = 1.0;
+ for (j = 0; j < nx; j++) tapershx[j] = 1.0;
+ for (j = 0; j < ny; j++) tapershy[j] = 1.0;
}
if (tap == 2 || tap == 3) {
- if (verbose) vmess("Taper for shots applied ntap=%li", ntap);
+ if (verbose) vmess("Taper for shots applied ntapx=%li ntapy=%li",ntapx,ntapy);
for (l = 0; l < nshots; l++) {
for (j = 1; j < nw; j++) {
for (k = 0; k < ny; k++) {
for (i = 0; i < nx; i++) {
- Refl[l*nx*ny*nw+j*nx*ny+k*nx+i].r *= tapersh[i];
- Refl[l*nx*ny*nw+j*nx*ny+k*nx+i].i *= tapersh[i];
+ Refl[l*nx*ny*nw+j*nx*ny+k*nx+i].r *= tapershx[i]*tapershy[k];
+ Refl[l*nx*ny*nw+j*nx*ny+k*nx+i].i *= tapershx[i]*tapershy[k];
}
}
}
}
}
- free(tapersh);
+ free(tapershx); free(tapershy);
/* check consistency of header values */
nxshot = unique_elements(xsrc,nshots);
@@ -567,14 +617,13 @@ int main (int argc, char **argv)
vmess("number of time samples (nt,nts) = %li (%li,%li)", ntfft, nt, nts);
vmess("frequency cutoffs = min:%.3f max:%.3f",fmin,fmax);
vmess("time sampling = %e ", dt);
- if (file_green != NULL) vmess("Green output file = %s ", file_green);
- if (file_gmin != NULL) vmess("Gmin output file = %s ", file_gmin);
- if (file_gplus != NULL) vmess("Gplus output file = %s ", file_gplus);
- if (file_pmin != NULL) vmess("Pmin output file = %s ", file_pmin);
- if (file_f2 != NULL) vmess("f2 (=pplus) output file = %s ", file_f2);
- if (file_f1min != NULL) vmess("f1min output file = %s ", file_f1min);
- if (file_f1plus != NULL)vmess("f1plus output file = %s ", file_f1plus);
- if (file_iter != NULL) vmess("Iterations output file = %s ", file_iter);
+ if (file_green != NULL) vmess("Green output file = %s ", file_green);
+ if (file_gmin != NULL) vmess("Gmin output file = %s ", file_gmin);
+ if (file_gplus != NULL) vmess("Gplus output file = %s ", file_gplus);
+ if (file_f2 != NULL) vmess("f2 output file = %s ", file_f2);
+ if (file_f1min != NULL) vmess("f1min output file = %s ", file_f1min);
+ if (file_f1plus != NULL)vmess("f1plus output file = %s ", file_f1plus);
+ if (file_iter != NULL) vmess("Iterations output file = %s ", file_iter);
}
@@ -593,12 +642,8 @@ int main (int argc, char **argv)
vmess("number of output traces = x:%li y:%li total:%li", n2out, n3out, n2out*n3out);
vmess("number of output samples = %li", ntfft);
vmess("Size of output data/file = %.1f MB", mem);
- if (compact>0) {
- vmess("Save format for homg and imag = compact");
- }
- else {
- vmess("Save format for homg and imag = normal");
- }
+ if (compact==0) vmess("Save format for homg and imag = compact");
+ else vmess("Save format for homg and imag = normal");
}
@@ -721,11 +766,9 @@ int main (int argc, char **argv)
ix = ixpos[i];
iy = iypos[i];
Ni[l*nys*nxs*nts+i*nts+j] = -iRN[l*nys*nxs*nts+iy*nxs*nts+ix*nts+j];
- pmin[l*nys*nxs*nts+i*nts+j] += iRN[l*nys*nxs*nts+iy*nxs*nts+ix*nts+j];
energyNi += iRN[l*nys*nxs*nts+ix*nts+j]*iRN[l*nys*nxs*nts+ix*nts+j];
for (j = 1; j < nts; j++) {
Ni[l*nys*nxs*nts+i*nts+j] = -iRN[l*nys*nxs*nts+iy*nxs*nts+ix*nts+nts-j];
- pmin[l*nys*nxs*nts+i*nts+j] += iRN[l*nys*nxs*nts+iy*nxs*nts+ix*nts+j];
energyNi += iRN[l*nys*nxs*nts+iy*nxs*nts+ix*nts+j]*iRN[l*nys*nxs*nts+iy*nxs*nts+ix*nts+j];
}
if (plane_wave!=0) { /* don't reverse in time */
@@ -741,8 +784,8 @@ int main (int argc, char **argv)
}
/* apply mute window based on times of direct arrival (in muteW) */
- applyMute3D(Ni, muteW, smooth, above, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift);
- if (plane_wave!=0) applyMute3D(Nig, muteW, smooth, above, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift); // ToDo add tsynW taper for tilted plane-waves
+ applyMute3D(Ni, muteW, smooth, above, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
+ if (plane_wave!=0) applyMute3D(Nig, muteW, smooth, above, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW); // ToDo add tsynW taper for tilted plane-waves
if (iter % 2 == 0) { /* even iterations update: => f_1^-(t) */
@@ -756,7 +799,7 @@ int main (int argc, char **argv)
}
}
}
- if (above==-2) applyMute3D(f1min, muteW, smooth, 0, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift);
+ if (above==-2) applyMute3D(f1min, muteW, smooth, 0, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
}
else { /* plane wave scheme */
for (l = 0; l < Nfoc; l++) {
@@ -770,7 +813,7 @@ int main (int argc, char **argv)
}
}
}
- if (above==-2) applyMute3D(f1min, muteW, smooth, 0, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift);
+ if (above==-2) applyMute3D(f1min, muteW, smooth, 0, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
}
}
else {/* odd iterations update: => f_1^+(t) */
@@ -806,22 +849,40 @@ int main (int argc, char **argv)
free(Ni);
free(Nig);
- /* compute full Green's function G = int R * f2(t) + f2(-t) = Pplus + Pmin */
+ /* compute full Green's function G = int R * f2(t) + f2(-t) */
+ /* use f2 as operator on R in frequency domain */
+ mode=1;
+ if (niter==0) {
+ synthesis3D(Refl, Fop, G_d, iRN, nx, ny, nt, nxs, nys, nts, dt, xsyn, ysyn,
+ Nfoc, xrcv, yrcv, xsrc, ysrc, xnx, fxse, fxsb, fyse, fysb, dxs, dys,
+ dxsrc, dysrc, dx, dy, ntfft, nw, nw_low, nw_high, mode, reci, nshots,
+ nxshot, nyshot, ixpos, iypos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv,
+ ixmask, verbose);
+ }
+ else {
+ synthesis3D(Refl, Fop, f2p, iRN, nx, ny, nt, nxs, nys, nts, dt, xsyn, ysyn,
+ Nfoc, xrcv, yrcv, xsrc, ysrc, xnx, fxse, fxsb, fyse, fysb, dxs, dys,
+ dxsrc, dysrc, dx, dy, ntfft, nw, nw_low, nw_high, mode, reci, nshots,
+ nxshot, nyshot, ixpos, iypos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv,
+ ixmask, verbose);
+ }
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
j = 0;
+ ix = ixpos[i];
+ iy = iypos[i];
/* set green to zero if mute-window exceeds nt/2 */
- if (muteW[l*nys*nxs+iypos[i]*nxs+ixpos[i]] >= nts/2) {
+ if (muteW[l*nys*nxs+iy*nxs+ix] >= nts/2) {
memset(&green[l*nys*nxs*nts+i*nts],0, sizeof(float)*nt);
continue;
}
- green[l*nys*nxs*nts+i*nts+j] = f2p[l*nys*nxs*nts+i*nts+j] + pmin[l*nys*nxs*nts+i*nts+j];
+ green[l*nys*nxs*nts+i*nts+j] = f2p[l*nys*nxs*nts+i*nts+j] + iRN[l*nys*nxs*nts+iy*nxs*nts+ix*nts+j];
for (j = 1; j < nts; j++) {
- green[l*nys*nxs*nts+i*nts+j] = f2p[l*nys*nxs*nts+i*nts+nts-j] + pmin[l*nys*nxs*nts+i*nts+j];
+ green[l*nys*nxs*nts+i*nts+j] = f2p[l*nys*nxs*nts+i*nts+nts-j] + iRN[l*nys*nxs*nts+iy*nxs*nts+ix*nts+j];
}
}
}
- applyMute3D(green, muteW, smooth, 4, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift);
+ if (!plane_wave) applyMute3D(green, muteW, smooth, 4, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
/* compute upgoing Green's function G^+,- */
if (file_gmin != NULL || file_imag!= NULL) {
@@ -829,20 +890,11 @@ int main (int argc, char **argv)
/* use f1+ as operator on R in frequency domain */
mode=1;
- if (niter==0) {
- synthesis3D(Refl, Fop, G_d, iRN, nx, ny, nt, nxs, nys, nts, dt, xsyn, ysyn,
- Nfoc, xrcv, yrcv, xsrc, ysrc, xnx, fxse, fxsb, fyse, fysb, dxs, dys,
- dxsrc, dysrc, dx, dy, ntfft, nw, nw_low, nw_high, mode, reci, nshots,
- nxshot, nyshot, ixpos, iypos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv,
- ixmask, verbose);
- }
- else {
- synthesis3D(Refl, Fop, f1plus, iRN, nx, ny, nt, nxs, nys, nts, dt, xsyn, ysyn,
- Nfoc, xrcv, yrcv, xsrc, ysrc, xnx, fxse, fxsb, fyse, fysb, dxs, dys,
- dxsrc, dysrc, dx, dy, ntfft, nw, nw_low, nw_high, mode, reci, nshots,
- nxshot, nyshot, ixpos, iypos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv,
- ixmask, verbose);
- }
+ synthesis3D(Refl, Fop, f1plus, iRN, nx, ny, nt, nxs, nys, nts, dt, xsyn, ysyn,
+ Nfoc, xrcv, yrcv, xsrc, ysrc, xnx, fxse, fxsb, fyse, fysb, dxs, dys,
+ dxsrc, dysrc, dx, dy, ntfft, nw, nw_low, nw_high, mode, reci, nshots,
+ nxshot, nyshot, ixpos, iypos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv,
+ ixmask, verbose);
/* compute upgoing Green's G^-,+ */
for (l = 0; l < Nfoc; l++) {
@@ -857,7 +909,7 @@ int main (int argc, char **argv)
}
}
/* Apply mute with window for Gmin */
- applyMute3D(Gmin, muteW, smooth, 4, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift);
+ if (!plane_wave) applyMute3D(Gmin, muteW, smooth, 4, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
} /* end if Gmin */
/* compute downgoing Green's function G^+,+ */
@@ -885,7 +937,7 @@ int main (int argc, char **argv)
}
}
/* Apply mute with window for Gplus */
- applyMute3D(Gplus, muteW, smooth, 4, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift);
+ if (!plane_wave) applyMute3D(Gplus, muteW, smooth, 4, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
} /* end if Gplus */
/* Estimate the amplitude of the Marchenko Redatuming */
@@ -896,7 +948,7 @@ int main (int argc, char **argv)
ampscl = (float *)calloc(Nfoc,sizeof(float));
Gd = (float *)calloc(Nfoc*nxs*nys*ntfft,sizeof(float));
memcpy(Gd,Gplus,sizeof(float)*Nfoc*nxs*nys*ntfft);
- applyMute3D(Gd, muteW, smooth, 2, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift);
+ applyMute3D(Gd, muteW, smooth, 2, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
// Determine amplitude and apply scaling
AmpEst3D(G_d, Gd, ampscl, Nfoc, nxs, nys, ntfft, ixpos, iypos, npos, file_wav, dxs, dys, dt);
@@ -904,7 +956,6 @@ int main (int argc, char **argv)
for (j=0; j<nxs*nys*nts; j++) {
green[l*nxs*nts+j] *= ampscl[l];
f2p[l*nxs*nys*nts+j] *= ampscl[l];
- pmin[l*nxs*nys*nts+j] *= ampscl[l];
f1plus[l*nxs*nys*nts+j] *= ampscl[l];
f1min[l*nxs*nys*nts+j] *= ampscl[l];
if (file_gplus != NULL) Gplus[l*nxs*nys*nts+j] *= ampscl[l];
@@ -1131,10 +1182,6 @@ int main (int argc, char **argv)
fp_f2 = fopen(file_f2, "w+");
if (fp_f2==NULL) verr("error on creating output file %s", file_f2);
}
- if (file_pmin != NULL) {
- fp_pmin = fopen(file_pmin, "w+");
- if (fp_pmin==NULL) verr("error on creating output file %s", file_pmin);
- }
if (file_f1plus != NULL) {
fp_f1plus = fopen(file_f1plus, "w+");
if (fp_f1plus==NULL) verr("error on creating output file %s", file_f1plus);
@@ -1144,6 +1191,48 @@ int main (int argc, char **argv)
if (fp_f1min==NULL) verr("error on creating output file %s", file_f1min);
}
+ if (zfp) {
+ vmess("zfp compression applied");
+ zfpt.dx = dx;
+ zfpt.dy = dy;
+ zfpt.dz = dt;
+ zfpt.ndim = 3;
+ zfpt.ns = Nfoc;
+ zfpt.scale = -1000;
+ zfpt.nt = ntfft;
+ zfpt.fmin = fmin;
+ zfpt.fmax = fmax;
+ zfpt.nz = ntfft ;
+ zfpt.tolerance = tolerance;
+ zfpt.fz = f1;
+ zfpt.fx = f2;
+ zfpt.fy = f3;
+
+ if (file_green != NULL) {
+ nread = fwrite(&zfpt, 1, TOPBYTES, fp_out);
+ assert(nread == TOPBYTES);
+ }
+ if (file_gmin != NULL) {
+ nread = fwrite(&zfpt, 1, TOPBYTES, fp_gmin);
+ assert(nread == TOPBYTES);
+ }
+ if (file_gplus != NULL) {
+ nread = fwrite(&zfpt, 1, TOPBYTES, fp_gplus);
+ assert(nread == TOPBYTES);
+ }
+ if (file_f2 != NULL) {
+ nread = fwrite(&zfpt, 1, TOPBYTES, fp_f2);
+ assert(nread == TOPBYTES);
+ }
+ if (file_f1plus != NULL) {
+ nread = fwrite(&zfpt, 1, TOPBYTES, fp_f1plus);
+ assert(nread == TOPBYTES);
+ }
+ if (file_f1min != NULL) {
+ nread = fwrite(&zfpt, 1, TOPBYTES, fp_f1min);
+ assert(nread == TOPBYTES);
+ }
+ }
tracf = 1;
for (l = 0; l < Nfoc; l++) {
@@ -1169,41 +1258,74 @@ int main (int argc, char **argv)
}
}
+ if (zfp) {
+ zfpm.gx = NINT(1000*(f2));
+ zfpm.gy = NINT(1000*(f3));
+ zfpm.sx = NINT(xsyn[l]*1000);
+ zfpm.sy = NINT(ysyn[l]*1000);
+ zfpm.sz = NINT(-zsyn[l]*1000);
+ }
+
if (file_green != NULL) {
- ret = writeData3D(fp_out, (float *)&green[l*size], hdrs_out, n1, n2*n3);
- if (ret < 0 ) verr("error on writing output file.");
+ if (zfp==0) {
+ ret = writeData3D(fp_out, (float *)&green[l*size], hdrs_out, n1, n2*n3);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ else {
+ zfpcompress((float *)&green[l*size],n2,n3,n1,tolerance,zfpm,fp_out);
+ }
}
if (file_gmin != NULL) {
- ret = writeData3D(fp_gmin, (float *)&Gmin[l*size], hdrs_out, n1, n2*n3);
- if (ret < 0 ) verr("error on writing output file.");
+ if (zfp==0) {
+ ret = writeData3D(fp_gmin, (float *)&Gmin[l*size], hdrs_out, n1, n2*n3);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ else {
+ zfpcompress((float *)&Gmin[l*size],n2,n3,n1,tolerance,zfpm,fp_gmin);
+ }
}
if (file_gplus != NULL) {
- ret = writeData3D(fp_gplus, (float *)&Gplus[l*size], hdrs_out, n1, n2*n3);
- if (ret < 0 ) verr("error on writing output file.");
+ if (zfp==0) {
+ ret = writeData3D(fp_gplus, (float *)&Gplus[l*size], hdrs_out, n1, n2*n3);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ else {
+ zfpcompress((float *)&Gplus[l*size],n2,n3,n1,tolerance,zfpm,fp_gplus);
+ }
}
if (file_f2 != NULL) {
- ret = writeData3D(fp_f2, (float *)&f2p[l*size], hdrs_out, n1, n2*n3);
- if (ret < 0 ) verr("error on writing output file.");
- }
- if (file_pmin != NULL) {
- ret = writeData3D(fp_pmin, (float *)&pmin[l*size], hdrs_out, n1, n2*n3);
- if (ret < 0 ) verr("error on writing output file.");
+ if (zfp==0) {
+ ret = writeData3D(fp_f2, (float *)&f2p[l*size], hdrs_out, n1, n2*n3);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ else {
+ zfpcompress((float *)&f2p[l*size],n2,n3,n1,tolerance,zfpm,fp_f2);
+ }
}
if (file_f1plus != NULL) {
- ret = writeData3D(fp_f1plus, (float *)&f1plus[l*size], hdrs_out, n1, n2*n3);
- if (ret < 0 ) verr("error on writing output file.");
+ if (zfp==0) {
+ ret = writeData3D(fp_f1plus, (float *)&f1plus[l*size], hdrs_out, n1, n2*n3);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ else {
+ zfpcompress((float *)&f1plus[l*size],n2,n3,n1,tolerance,zfpm,fp_f1plus);
+ }
}
if (file_f1min != NULL) {
- ret = writeData3D(fp_f1min, (float *)&f1min[l*size], hdrs_out, n1, n2*n3);
- if (ret < 0 ) verr("error on writing output file.");
+ if (zfp==0) {
+ ret = writeData3D(fp_f1min, (float *)&f1min[l*size], hdrs_out, n1, n2*n3);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ else {
+ zfpcompress((float *)&f1min[l*size],n2,n3,n1,tolerance,zfpm,fp_f1min);
+ }
}
}
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);}
- if (file_pmin != NULL) {ret += fclose(fp_pmin);}
if (file_f1plus != NULL) {ret += fclose(fp_f1plus);}
if (file_f1min != NULL) {ret += fclose(fp_f1min);}
if (ret < 0) verr("err %li on closing output file",ret);
@@ -1238,3 +1360,78 @@ long unique_elements(float *arr, long len)
}
return unique;
}
+
+long zfpcompress(float* data, long nx, long ny, long nz, double tolerance, zfpmar zfpm, FILE *file)
+{
+
+ zfp_field* field = NULL;
+ zfp_stream* zfp = NULL;
+ bitstream* stream = NULL;
+ void* fi = NULL;
+ void* fo = NULL;
+ void* buffer = NULL;
+ size_t rawsize = 0;
+ size_t zfpsize = 0;
+ size_t bufsize = 0;
+ size_t nwrite;
+ zfp_exec_policy exec = zfp_exec_serial;
+
+ zfp = zfp_stream_open(NULL);
+ field = zfp_field_alloc();
+
+ zfp_field_set_pointer(field, (void *)data);
+
+ zfp_field_set_type(field, zfp_type_float);
+ if (ny<2) zfp_field_set_size_2d(field, (uint)nz, (uint)nx);
+ else zfp_field_set_size_3d(field, (uint)nz, (uint)nx, (uint)ny);
+
+ zfp_stream_set_accuracy(zfp, tolerance);
+
+ if (!zfp_stream_set_execution(zfp, exec)) {
+ fprintf(stderr, "serial execution not available\n");
+ return EXIT_FAILURE;
+ }
+
+ bufsize = zfp_stream_maximum_size(zfp, field);
+ if (!bufsize) {
+ fprintf(stderr, "invalid compression parameters\n");
+ return EXIT_FAILURE;
+ }
+
+ buffer = malloc(bufsize);
+ if (!buffer) {
+ fprintf(stderr, "cannot allocate memory\n");
+ return EXIT_FAILURE;
+ }
+
+ stream = stream_open(buffer, bufsize);
+ if (!stream) {
+ fprintf(stderr, "cannot open compressed stream\n");
+ return EXIT_FAILURE;
+ }
+ zfp_stream_set_bit_stream(zfp, stream);
+
+ if (!zfp_stream_set_execution(zfp, exec)) {
+ fprintf(stderr, "serial execution not available\n");
+ return EXIT_FAILURE;
+ }
+
+ zfpsize = zfp_compress(zfp, field);
+ if (zfpsize == 0) {
+ fprintf(stderr, "compression failed\n");
+ return EXIT_FAILURE;
+ }
+
+ zfpm.nx = nx;
+ zfpm.ny = ny;
+ zfpm.compsize = zfpsize;
+
+ nwrite = fwrite(&zfpm, 1, MARBYTES, file);
+ assert(nwrite == MARBYTES);
+ if (fwrite(buffer, 1, zfpsize, file) != zfpsize) {
+ fprintf(stderr, "cannot write compressed file\n");
+ return EXIT_FAILURE;
+ }
+
+ return 1;
+}
\ No newline at end of file
diff --git a/marchenko3D/readShotData3Dzfp.c b/marchenko3D/readShotData3Dzfp.c
index 8de346f..bbc4ae2 100644
--- a/marchenko3D/readShotData3Dzfp.c
+++ b/marchenko3D/readShotData3Dzfp.c
@@ -97,6 +97,7 @@ long readShotData3Dzfp(char *filename, float *xrcv, float *yrcv, float *xsrc, fl
}
}
+ fclose(fp);
free(trace);
return 0;
diff --git a/utils/HomG.c b/utils/HomG.c
index 736b4f4..7e2d6f0 100755
--- a/utils/HomG.c
+++ b/utils/HomG.c
@@ -6,6 +6,8 @@
#include <math.h>
#include <assert.h>
#include <genfft.h>
+#include "zfpmar.h"
+#include <zfp.h>
#ifndef MAX
#define MAX(x,y) ((x) > (y) ? (x) : (y))
@@ -21,6 +23,25 @@ typedef struct _complexStruct { /* complex number */
} complex;
#endif/* complex */
+/*
+The schemes in this module use a variety of retrieval representations
+For more information about Green's function retrieval see:
+Brackenhoff, J., Thorbecke, J., & Wapenaar, K. (2019).
+Virtual sources and receivers in the real Earth: Considerations for practical applications.
+Journal of Geophysical Research: Solid Earth, 124, 11802– 11821.
+https://doi.org/10.1029/2019JB018485
+
+Brackenhoff, J., Thorbecke, J., and Wapenaar, K.:
+Monitoring of induced distributed double-couple sources using Marchenko-based virtual receivers.
+Solid Earth, 10, 1301–1319,
+https://doi.org/10.5194/se-10-1301-2019, 2019.
+
+Wapenaar, K., Brackenhoff, J., Thorbecke, J. et al.
+Virtual acoustics in inhomogeneous media with single-sided access.
+Sci Rep 8, 2497 (2018).
+https://doi.org/10.1038/s41598-018-20924-x
+*/
+
long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, float *d1, float *d2, float *d3,
float *f1, float *f2, float *f3, float *sclsxgxsygy, long *nxm);
double wallclock_time(void);
@@ -28,6 +49,14 @@ long writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2);
long readSnapData3D(char *filename, float *data, segy *hdrs, long nsnaps, long nx, long ny, long nz,
long sx, long ex, long sy, long ey, long sz, long ez);
void conjugate(float *data, long nsam, long nrec, float dt);
+long zfpdecompress(float* data, long nx, long ny, long nz, long comp, double tolerance, FILE *file);
+
+void getFileInfo3Dzfp(char *filename, long *n1, long *n2, long *n3, long *ngath,
+ float *d1, float *d2, float *d3, float *f1, float *f2, float *f3,
+ float *fmin, float *fmax, float *scl, long *nxm);
+void getVirReczfp(char *filename, long *nxs, long *nys, long *nxr, long *nyr, long *nt);
+void readzfpdata(char *filename, float *data, long size);
+void getxyzzfp(char *filename, long *sx, long *sy, long *sz, long iz, long nz);
void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout);
void pad_data(float *data, long nsam, long nrec, long nsamout, float *datout);
@@ -36,6 +65,8 @@ void corr(float *data1, float *data2, float *cov, long nrec, long nsam, float dt
void timeDiff(float *data, long nsam, long nrec, float dt, float fmin, float fmax, long opt);
void depthDiff(float *data, long nsam, long nrec, float dt, float dx, float fmin, float fmax, float c, long opt);
void pad2d_data(float *data, long nsam, long nrec, long nsamout, long nrecout, float *datout);
+void getVirRec(char *filename, long *nxs, long *nys, long *nxr, long *nyr, long *nt);
+void convol2(float *data1, float *data2, float *con, long nrec, long nsam, float dt, float fmin, float fmax, long opt);
char *sdoc[] = {
" ",
@@ -46,8 +77,8 @@ char *sdoc[] = {
" ",
" Required parameters: ",
"",
-" file_in= ................. First file of the array of receivers",
-" file_shot= ............... File containing the shot location",
+" file_vr= ................. First file of the array of virtual receivers",
+" file_vs= ................. File containing the virtual source",
" ",
" Optional parameters: ",
" ",
@@ -58,8 +89,8 @@ char *sdoc[] = {
" inx= ..................... Number of sources per depth level",
" zrcv= .................... z-coordinate of first receiver location",
" xrcv= .................... x-coordinate of first receiver location",
-" dzrcv= ................... z-spacing of receivers",
-" dxrcv= ................... x-spacing of receivers",
+" zfps=0 ................... virtual source data are in SU format (=0) or zfp compressed (=1)",
+" zfpr=0 ................... virtual receiver data are in SU format (=0) or zfp compressed (=1)",
" shift=0.0 ................ shift per shot",
" scheme=0 ................. Scheme used for retrieval. 0=Marchenko,",
" 1=Marchenko with multiple sources, 2=classical",
@@ -67,14 +98,17 @@ NULL};
int main (int argc, char **argv)
{
- FILE *fp_in, *fp_shot, *fp_out;
- char *fin, *fshot, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100];
- float *indata, *Ghom, *shotdata, *shotdata_jkz, rho, fmin, fmax;
- float dt, dy, dx, t0, y0, x0, xmin, xmax1, sclsxgx, f1, f2, f3, dxrcv, dyrcv, dzrcv;
- float *conv, *conv2, *tmp1, *tmp2, cp, shift;
- long nshots, nt, ny, nx, ntraces, ret, ix, iy, it, is, ir, ig, file_det, nxs, nys, nzs, verbose;
- long pos1, npos, zmax, inx, numb, dnumb, count, scheme, ntmax, ntshift, shift_num;
- segy *hdr_in, *hdr_out, *hdr_shot;
+ FILE *fp_in, *fp_shot, *fp_out;
+ char *fin, *fshot, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100];
+ float *rcvdata, *Ghom, *shotdata, *shotdata_jkz, rho, fmin, fmax;
+ float dt, dy, dx, t0, y0, x0, xmin, xmax1, sclsxgx, dxrcv, dyrcv, dzrcv;
+ float *conv, *conv2, *tmp1, *tmp2, cp, shift;
+ long nshots, ntvs, nyvs, nxvs, ntraces, ret, ix, iy, it, is, ir, ig, file_det, verbose;
+ long ntr, nxr, nyr, nsr, i, l, j, k, nxvr, nyvr, nzvr;
+ float dtr, dxr, dyr, ftr, fxr, fyr, sclr, scl;
+ long pos1, npos, zmax, numb, dnumb, scheme, ntmax, ntshift, shift_num, zfps, zfpr, size;
+ long ixr, iyr, zsrc, zrcv, *xvr, *yvr, *zvr;
+ segy *hdr_rcv, *hdr_out, *hdr_shot;
initargs(argc, argv);
requestdoc(1);
@@ -82,17 +116,12 @@ int main (int argc, char **argv)
/*----------------------------------------------------------------------------*
* Get the parameters passed to the function
*----------------------------------------------------------------------------*/
- if (!getparstring("fin", &fin)) fin = NULL;
- if (!getparstring("fshot", &fshot)) fshot = NULL;
- if (!getparstring("fout", &fout)) fout = "out.su";
+ if (!getparstring("file_vr", &fin)) fin = NULL;
+ if (!getparstring("file_vs", &fshot)) fshot = NULL;
+ if (!getparstring("file_out", &fout)) fout = "out.su";
if (!getparlong("zmax", &zmax)) zmax = 0;
- if (!getparlong("inx", &inx)) inx = 0;
- if (!getparfloat("zrcv", &f1)) f1 = 0;
- if (!getparfloat("xrcv", &f2)) f2 = 0;
- if (!getparfloat("dzrcv", &dzrcv)) dzrcv = -1;
- if (!getparfloat("dxrcv", &dxrcv)) dxrcv = -1;
if (!getparfloat("rho", &rho)) rho=1000.0;
- if (!getparfloat("cp", &cp)) cp = 1500.0;
+ if (!getparfloat("cp", &cp)) cp = 1000.0;
if (!getparfloat("fmin", &fmin)) fmin=0.0;
if (!getparfloat("fmax", &fmax)) fmax=100.0;
if (!getparfloat("shift", &shift)) shift=0.0;
@@ -101,8 +130,10 @@ int main (int argc, char **argv)
if (!getparlong("scheme", &scheme)) scheme = 0;
if (!getparlong("ntmax", &ntmax)) ntmax = 0;
if (!getparlong("verbose", &verbose)) verbose = 0;
- if (fin == NULL) verr("Incorrect f2 input");
- if (fshot == NULL) verr("Incorrect Green input");
+ if (!getparlong("zfps", &zfps)) zfps = 0;
+ if (!getparlong("zfpr", &zfpr)) zfpr = 0;
+ if (fin == NULL) verr("Incorrect vr input");
+ if (fshot == NULL) verr("Incorrect vs input");
/*----------------------------------------------------------------------------*
* Split the filename so the number can be changed
@@ -123,101 +154,186 @@ int main (int argc, char **argv)
* Determine the amount of files to be read
*----------------------------------------------------------------------------*/
file_det = 1;
- nzs=0;
+ nzvr=0;
while (file_det) {
- sprintf(fins,"z%li",nzs*dnumb+numb);
+ sprintf(fins,"z%li",nzvr*dnumb+numb);
sprintf(fin,"%s%s%s",fbegin,fins,fend);
fp_in = fopen(fin, "r");
if (fp_in == NULL) {
- if (nzs == 0) {
+ if (nzvr == 0) {
verr("error on opening basefile=%s", fin);
}
- else if (nzs == 1) {
+ else if (nzvr == 1) {
vmess("1 file detected");
file_det = 0;
break;
}
else {
- vmess("%li files detected",nzs);
+ vmess("%li files detected",nzvr);
file_det = 0;
break;
}
}
fclose(fp_in);
- nzs++;
+ nzvr++;
}
- if (inx < 1) {
- inx = 1;
- }
+ if (zmax < 1) zmax=0;
+ if (zmax < nzvr && zmax > 0) nzvr=zmax;
- if (zmax < 1) zmax=1;
- if (zmax < nzs) nzs=zmax;
-
- nxs = inx;
- count=0;
- npos = nxs*nzs;
+ /*----------------------------------------------------------------------------*
+ * Determine the other sizes of the files
+ *----------------------------------------------------------------------------*/
+ sprintf(fins,"z%li",0);
+ sprintf(fin,"%s%s%s",fbegin,fins,fend);
+ if (zfpr) getVirReczfp(fin, &nxvr, &nyvr, &nxr, &nyr, &ntr);
+ else getVirRec(fin, &nxvr, &nyvr, &nxr, &nyr, &ntr);
+
+ if (verbose) {
+ if (zfpr) vmess("Virtual receiver data are zfp compressed");
+ else vmess("Virtual receiver data are in SU format");
+ vmess("Number of virtual receivers : %li (x=%li) (y=%li) (z=%li)",nxvr*nyvr*nzvr,nxvr,nyvr,nzvr);
+ vmess("Number of samples for each receiver : x=%li y=%li t=%li",nxr,nyr,ntr);
+ }
- if (verbose) vmess("nxs: %li, nzs: %li",nxs,nzs);
+ /*----------------------------------------------------------------------------*
+ * Get the file info for the source position and read in the data
+ *----------------------------------------------------------------------------*/
nshots = 0;
- getFileInfo3D(fshot, &nt, &nx, &ny, &nshots, &dt, &dx, &dy, &t0, &x0, &y0, &sclsxgx, &ntraces);
+ if (zfps) getFileInfo3Dzfp(fshot, &ntvs, &nxvs, &nyvs, &nshots, &dt, &dx, &dy, &t0, &x0, &y0, &fmin, &fmax, &sclsxgx, &ntraces);
+ else getFileInfo3D(fshot, &ntvs, &nxvs, &nyvs, &nshots, &dt, &dx, &dy, &t0, &x0, &y0, &sclsxgx, &ntraces);
- if (dxrcv < 0) dxrcv=dx;
- if (dzrcv < 0) dzrcv=dx;
+ scl = dx*dy*dt;
- // ngath zijn het aantal schoten
- shotdata = (float *)malloc(nt*nx*nshots*sizeof(float));
- hdr_shot = (segy *)calloc(nx*nshots,sizeof(segy));
+ shotdata = (float *)malloc(ntvs*nxvs*nyvs*nshots*sizeof(float));
+ hdr_shot = (segy *)calloc(nxvs*nyvs*nshots,sizeof(segy));
fp_shot = fopen(fshot,"r");
if (fp_shot == NULL) {
verr("Could not open file");
}
- vmess("nt: %li nx: %li nshots: %li",nt,nx,nshots);
fclose(fp_shot);
- readSnapData3D(fshot, &shotdata[0], &hdr_shot[0], nshots, nx, ny, nt, 0, nx, 0, ny, 0, nt);
+ if (verbose) {
+ if (zfps) vmess("Virtual source data are zfp compressed");
+ else vmess("Virtual source data are in SU format");
+ vmess("Number of virtual sources : %li ",nshots);
+ vmess("Number of samples for each source : x=%li y=%li t=%li",nxvs,nyvs,ntvs);
+ vmess("Sampling distance is : x=%.3e y=%.3e t=%.3e",dx,dy,dt);
+ vmess("Scaling of the transforms : %.3e",scl);
+ }
- hdr_out = (segy *)calloc(nxs,sizeof(segy));
- Ghom = (float *)malloc(nt*npos*sizeof(float));
+ if (ntr!=ntvs) verr("number of t-samples between virtual source (%li) and virtual receivers (%li) is not equal",ntvs,ntr);
+ if (nxr!=nxvs) verr("number of x-samples between virtual source (%li) and virtual receivers (%li) is not equal",nxvs,nxr);
+ if (nyr!=nyvs) verr("number of y-samples between virtual source (%li) and virtual receivers (%li) is not equal",nyvs,nyr);
- if (scheme==2) {
- vmess("Classical representation");
- shotdata_jkz = (float *)malloc(nt*nx*nshots*sizeof(float));
- for (ix = 0; ix < nx; ix++) {
- for (it = 0; it < nt; it++) {
- shotdata_jkz[ix*nt+it] = shotdata[ix*nt+it];
+ size = nxr*nyr*ntr;
+
+ if (zfps) readzfpdata(fshot, &shotdata[0], size);
+ else readSnapData3D(fshot, &shotdata[0], &hdr_shot[0], nshots, nxvs, nyvs, ntvs, 0, nxvs, 0, nyvs, 0, ntvs);
+
+ hdr_out = (segy *)calloc(nxvr*nyvr,sizeof(segy));
+ Ghom = (float *)calloc(ntr*nxvr*nyvr*nzvr,sizeof(float));
+ xvr = (long *)malloc(nxvr*nyvr*nzvr*sizeof(long));
+ yvr = (long *)malloc(nxvr*nyvr*nzvr*sizeof(long));
+ zvr = (long *)malloc(nxvr*nyvr*nzvr*sizeof(long));
+
+ /*----------------------------------------------------------------------------*
+ * Get the file info for the source position
+ *----------------------------------------------------------------------------*/
+
+ if (scheme==0) {
+ if (verbose) vmess("Marchenko Homogeneous Green's function retrieval with G source");
+ }
+ else if (scheme==1) {
+ if (verbose) vmess("Marchenko Homogeneous Green's function retrieval with f2 source");
+ for (k = 0; k < nyvs; k++) {
+ depthDiff(&shotdata[k*nxvs*ntvs], ntvs, nxvs, dt, dx, fmin, fmax, cp, 1);
+ }
+ }
+ else if (scheme==2) {
+ if (verbose) vmess("Marchenko Green's function retrieval with source depending on position");
+ if (nshots<2) verr("Number of shots required is 2 (1=G, 2=f_2)");
+ for (k = 0; k < nyvs; k++) {
+ depthDiff(&shotdata[ntvs*nxvs*nyvs+k*nxvs*ntvs], ntvs, nxvs, dt, dx, fmin, fmax, cp, 1);
+ }
+ zsrc = labs(hdr_shot[0].sdepth);
+ }
+ else if (scheme==3) {
+ if (verbose) vmess("Marchenko Green's function retrieval with G source");
+ }
+ else if (scheme==4) {
+ if (verbose) vmess("Marchenko Green's function retrieval with f2 source");
+ }
+ else if (scheme==5) { // Scale the Green's functions if the classical scheme is used
+ if (verbose) vmess("Classical Homogeneous Green's function retrieval");
+ shotdata_jkz = (float *)calloc(nshots*nxvs*nyvs*ntvs,sizeof(float));
+ for (l = 0; l < nshots; l++) {
+ for (i = 0; i < nyvs*nxvs*ntvs; i++) {
+ shotdata_jkz[l*nyvs*nxvs*ntvs+i] = shotdata[l*nyvs*nxvs*ntvs+i];
+ }
+ conjugate(&shotdata_jkz[l*nyvs*nxvs*ntvs], ntvs, nxvs*nyvs, dt);
+ conjugate(&shotdata[l*nyvs*nxvs*ntvs], ntvs, nxvs*nyvs, dt);
+ for (k = 0; k < nyvs; k++) {
+ depthDiff(&shotdata_jkz[l*nyvs*nxvs*ntvs+k*nxvs*ntvs], ntvs, nxvs, dt, dx, fmin, fmax, cp, 1);
}
}
- conjugate(shotdata_jkz, nt, nx, dt);
- conjugate(shotdata, nt, nx, dt);
- depthDiff(shotdata_jkz, nt, nx, dt, dx, fmin, fmax, cp, 1);
- if (verbose) vmess("Applied jkz to source data");
}
- else if (scheme==0) {
- vmess("Marchenko representation");
+ else if (scheme==6) {
+ if (verbose) vmess("Marchenko Homogeneous Green's function retrieval with multiple sources");
+ if (verbose) vmess("Looping over %li source positions",nshots);
}
- else if (scheme==1) {
- vmess("Marchenko representation with multiple sources");
+ else if (scheme==7) {
+ if (verbose) vmess("Back propagation with multiple sources");
+ if (verbose) vmess("Looping over %li source positions",nshots);
+ }
+ else if (scheme==8) { // 0=f1p 1=f1m
+ if (verbose) vmess("f1+ redatuming");
+ if (nshots<2) verr("Not enough input for the homogeneous Green's function");
+ for (k = 0; k < nyvs; k++) {
+ depthDiff(&shotdata[0*nyvs*nxvs*ntvs+k*nxvs*ntvs], ntvs, nxvs, dt, dx, fmin, fmax, cp, 1);
+ conjugate(&shotdata[0*nyvs*nxvs*ntvs+k*nxvs*ntvs], ntvs, nxvs, dt);
+ depthDiff(&shotdata[1*nyvs*nxvs*ntvs+k*nxvs*ntvs], ntvs, nxvs, dt, dx, fmin, fmax, cp, 1);
+ conjugate(&shotdata[1*nyvs*nxvs*ntvs+k*nxvs*ntvs], ntvs, nxvs, dt);
+ }
+ }
+ else if (scheme==9) { // 0=f1p 1=f1m
+ if (verbose) vmess("f1- redatuming");
+ if (nshots<2) verr("Not enough input for the homogeneous Green's function");
+ for (k = 0; k < nyvs; k++) {
+ depthDiff(&shotdata[0*nyvs*nxvs*ntvs+k*nxvs*ntvs], ntvs, nxvs, dt, dx, fmin, fmax, cp, 1);
+ depthDiff(&shotdata[1*nyvs*nxvs*ntvs+k*nxvs*ntvs], ntvs, nxvs, dt, dx, fmin, fmax, cp, 1);
+ }
+ }
+ else if (scheme==10) {
+ if (verbose) vmess("2i IM(f1) redatuming");
+ shotdata_jkz = (float *)calloc(nshots*nxvs*nyvs*ntvs,sizeof(float));
+ for (k = 0; k < nyvs; k++) {
+ depthDiff(&shotdata[k*nxvs*ntvs] , ntvs, nxvs, dt, dx, fmin, fmax, cp, 1);
+ for (l = 0; l < nxvs*ntvs; l++) {
+ shotdata_jkz[k*nxvs*ntvs+l] = shotdata[k*nxvs*ntvs+l];
+ }
+ conjugate(&shotdata_jkz[k*nxvs*ntvs], ntvs, nxvs, dt);
+ }
+ }
+ else {
+ if (verbose) vmess("Marchenko Homogeneous Green's function retrieval with G source");
}
- else if (scheme==3) {
- vmess("Marchenko representation with multiple shot gathers");
- }
-#pragma omp parallel default(shared) \
- private(ix,it,is,indata, hdr_in,fins,fin2,fp_in,conv,ig,conv2,tmp1,tmp2)
-{
- indata = (float *)malloc(nt*nx*nxs*sizeof(float));
- hdr_in = (segy *)calloc(nx*nxs,sizeof(segy));
- conv = (float *)calloc(nx*nt,sizeof(float));
- conv2 = (float *)calloc(nx*nt,sizeof(float));
- if (scheme==2) {
- tmp1 = (float *)calloc(nx*nt,sizeof(float));
- tmp2 = (float *)calloc(nx*nt,sizeof(float));
- }
-#pragma omp for
- for (ir = 0; ir < nzs; ir++) {
+#pragma omp parallel for schedule(static,1) default(shared) \
+ private(ix,iy,k,i,j,l,it,is,rcvdata,hdr_rcv,fins,fin2,fp_in,conv,ig,tmp1,tmp2)
+ for (ir = 0; ir < nzvr; ir++) {
+
+ rcvdata = (float *)malloc(ntr*nxvr*nyvr*nxr*nyr*sizeof(float));
+ hdr_rcv = (segy *)calloc(nxvr*nyvr*nxr*nyr,sizeof(segy));
+ conv = (float *)calloc(nxr*ntr,sizeof(float));
+ if (scheme==5) {
+ tmp1 = (float *)calloc(nxr*ntr,sizeof(float));
+ tmp2 = (float *)calloc(nxr*ntr,sizeof(float));
+ }
+ if (scheme==6 || scheme==8 || scheme==9 || scheme==10) tmp1 = (float *)calloc(nxr*ntr,sizeof(float));
+
sprintf(fins,"z%li",ir*dnumb+numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
fp_in = fopen(fin2, "r");
@@ -225,319 +341,311 @@ int main (int argc, char **argv)
verr("Danger Will Robinson");
}
fclose(fp_in);
- readSnapData3D(fin2, &indata[0], &hdr_in[0], nxs, nx, ny, nt, 0, nx, 0, ny, 0, nt);
- for (is=0;is<nxs;is++) {
- if (scheme==0) { //Marchenko representation
- depthDiff(&indata[is*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
- convol(shotdata, &indata[is*nx*nt], conv, nx, nt, dt, -2);
- timeDiff(conv, nt, nx, dt, fmin, fmax, -2);
- for (ix=0; ix<nx; ix++) {
- for (it=0; it<nt/2; it++) {
- Ghom[(it+nt/2)*nxs*nzs+is*nzs+ir] += conv[ix*nt+it]/rho;
- Ghom[it*nxs*nzs+is*nzs+ir] += conv[ix*nt+(it+nt/2)]/rho;
- }
- }
- }
- else if (scheme==1) { //Marchenko representation with multiple sources
- depthDiff(&indata[is*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
- convol(shotdata, &indata[is*nx*nt], conv, nx, nt, dt, 0);
- timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
- for (ix=0; ix<nx; ix++) {
- for (it=0; it<nt/2; it++) {
- Ghom[(it+nt/2)*nxs*nzs+is*nzs+ir] += 2*conv[ix*nt+it]/rho;
- Ghom[it*nxs*nzs+is*nzs+ir] += 2*conv[ix*nt+(it+nt/2)]/rho;
- }
- }
- }
- else if (scheme==2) { //classical representation
- convol(&indata[is*nx*nt], shotdata_jkz, tmp1, nx, nt, dt, 0);
- depthDiff(&indata[is*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
- convol(&indata[is*nx*nt], shotdata, tmp2, nx, nt, dt, 0);
- //corr(&indata[is*nx*nt], shotdata, tmp2, nx, nt, dt, 0);
- for (ix = 0; ix < nx; ix++) {
- for (it = 0; it < nt; it++) {
- conv[ix*nt+it] = tmp2[ix*nt+it]+tmp1[ix*nt+it];
- }
- }
- timeDiff(conv, nt, nx, dt, fmin, fmax, -1);
- for (ix=0; ix<nx; ix++) {
- for (it=0; it<nt/2; it++) {
- Ghom[(it+nt/2)*nxs*nzs+is*nzs+ir] += conv[ix*nt+it]/rho;
- Ghom[it*nxs*nzs+is*nzs+ir] += conv[ix*nt+(it+nt/2)]/rho;
- }
+
+ if (zfpr) readzfpdata(fin2, &rcvdata[0], size);
+ else readSnapData3D(fin2, &rcvdata[0], &hdr_rcv[0], nxvr*nyvr, nxr, nyr, ntr, 0, nxr, 0, nyr, 0, ntr);
+
+ if (zfpr) getxyzzfp(fin2, xvr, yvr, zvr, ir, nzvr);
+
+ zrcv = labs(zvr[ir]);
+
+ for (l = 0; l < nxvr*nyvr; l++) {
+
+ if (!zfpr) {
+ xvr[l*nzvr+ir] = hdr_rcv[l*nxr*nyr].sx;
+ yvr[l*nzvr+ir] = hdr_rcv[l*nxr*nyr].sy;
+ zvr[l*nzvr+ir] = labs(hdr_rcv[l*nxr*nyr].sdepth);
+ }
+
+ if (scheme==0) { //Marchenko representation with G source
+ for (k = 0; k < nyr; k++) {
+ depthDiff(&rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], ntr, nxr, dt, dx, fmin, fmax, cp, 1);
+ convol(&shotdata[k*nxr*ntr], &rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], conv, nxr, ntr, dt, 0);
+ timeDiff(conv, ntr, nxr, dt, fmin, fmax, -3);
+ for (i=0; i<nxr; i++) {
+ for (j=0; j<ntr/2; j++) {
+ Ghom[(j+ntr/2)*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+j]/rho;
+ Ghom[j*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+(j+ntr/2)]/rho;
+ }
+ }
}
}
- if (scheme==3) { //Marchenko representation with multiple shot gathers
- depthDiff(&indata[is*nx*nt], nt, nx, dt, dx, fmin, fmax, cp, 1);
- for (ig=0; ig<nshots; ig++) {
- convol(&shotdata[ig*nx*nt], &indata[is*nx*nt], conv, nx, nt, dt, -2);
- timeDiff(conv, nt, nx, dt, fmin, fmax, -2);
- shift_num = ig*((long)(shift/dt));
- for (ix = 0; ix < nx; ix++) {
- for (it = nt/2+1; it < nt; it++) {
- conv[ix*nt+it] = 0.0;
- }
- for (it = shift_num; it < nt; it++) {
- conv2[ix*nt+it] = conv[ix*nt+it-shift_num];
- }
- for (it = 0; it < shift_num; it++) {
- conv2[ix*nt+it] = conv[ix*nt+nt-shift_num+it];
+ else if (scheme==1) { //Marchenko representation with f2 source
+ for (k = 0; k < nyr; k++) {
+ convol(&shotdata[k*nxr*ntr], &rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], conv, nxr, ntr, dt, 0);
+ timeDiff(conv, ntr, nxr, dt, fmin, fmax, -3);
+ for (i=0; i<nxr; i++) {
+ for (j=0; j<ntr/2; j++) {
+ Ghom[(j+ntr/2)*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+j]/rho;
+ Ghom[j*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+(j+ntr/2)]/rho;
}
- }
- for (ix=0; ix<nx; ix++) {
- Ghom[(-1+nt/2)*nxs*nzs+is*nzs+ir] += conv2[ix*nt+nt-1]/rho;
- for (it=0; it<nt/2; it++) {
- Ghom[(it+nt/2)*nxs*nzs+is*nzs+ir] += conv2[ix*nt+it]/rho;
- //Ghom[it*nxs*nzs+is*nzs+ir] += conv2[ix*nt+(it+nt/2)]/rho;
- }
- }
+ }
+ }
+ }
+ else if (scheme==2) { //Marchenko representation without time-reversal using varying sources
+ for (k = 0; k < nyr; k++) {
+ if (zsrc > zrcv) {
+ if (verbose > 1) vmess("Homogeneous Green's function at %li uses G source (zsrc=%li)",zrcv);
+ depthDiff(&rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], ntr, nxr, dt, dx, fmin, fmax, cp, 1);
+ convol(&shotdata[k*nxr*ntr], &rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], conv, nxr, ntr, dt, 0);
+ }
+ else {
+ if (verbose > 1) vmess("Homogeneous Green's function at %li uses f_2 source (zsrc=%li)",zrcv);
+ convol(&shotdata[ntvs*nxvs*nyvs+k*nxr*ntr], &rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], conv, nxr, ntr, dt, 0);
+ }
+ timeDiff(conv, ntr, nxr, dt, fmin, fmax, -1);
+ for (i=0; i<nxr; i++) {
+ for (j=0; j<ntr/2; j++) {
+ Ghom[(j+ntr/2)*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+j]/rho;
+ Ghom[j*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+(j+ntr/2)]/rho;
+ }
+ }
+ }
+ }
+ else if (scheme==3) { //Marchenko representation without time-reversal G source
+ for (k = 0; k < nyr; k++) {
+ depthDiff(&rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], ntr, nxr, dt, dx, fmin, fmax, cp, 1);
+ convol2(&shotdata[k*nxr*ntr], &rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], conv, nxr, ntr, dt, fmin, fmax, 1);
+ for (i=0; i<nxr; i++) {
+ for (j=0; j<ntr/2; j++) {
+ Ghom[(j+ntr/2)*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+j]/rho;
+ Ghom[j*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+(j+ntr/2)]/rho;
+ }
+ }
+ }
+ }
+ else if (scheme==4) { //Marchenko representation without time-reversal f2 source
+ for (k = 0; k < nyr; k++) {
+ depthDiff(&rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], ntr, nxr, dt, dx, fmin, fmax, cp, 1);
+ convol(&shotdata[k*nxr*ntr], &rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], conv, nxr, ntr, dt, 0);
+ timeDiff(conv, ntr, nxr, dt, fmin, fmax, -1);
+ for (i=0; i<nxr; i++) {
+ for (j=0; j<ntr/2; j++) {
+ Ghom[(j+ntr/2)*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+j]/rho;
+ Ghom[j*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+(j+ntr/2)]/rho;
+ }
+ }
+ }
+ }
+ else if (scheme==5) { //classical representation
+ for (k = 0; k < nyr; k++) {
+ convol(&rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], &shotdata_jkz[k*nxr*ntr], tmp2, nxr, ntr, dt, 0);
+ depthDiff(&rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], ntr, nxr, dt, dx, fmin, fmax, cp, 1);
+ convol(&rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], &shotdata[k*nxr*ntr], tmp1, nxr, ntr, dt, 0);
+ for (i = 0; i < nxr; i++) {
+ for (j = 0; j < ntr; j++) {
+ conv[i*ntr+j] = tmp1[i*ntr+j]+tmp2[i*ntr+j];
+ }
+ }
+ timeDiff(conv, ntr, nxr, dt, fmin, fmax, -1);
+ for (i=0; i<nxr; i++) {
+ for (j=0; j<ntr/2; j++) {
+ Ghom[(j+ntr/2)*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+j]/rho;
+ Ghom[j*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+(j+ntr/2)]/rho;
+ }
+ }
+ }
+ }
+ else if (scheme==6) { //Marchenko representation with multiple shot gathers
+ for (k = 0; k < nyr; k++) {
+ depthDiff(&rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], ntr, nxr, dt, dx, fmin, fmax, cp, 1);
+ for (is=0; is<nshots; is++) {
+ convol(&shotdata[is*nyr*nxr*ntr+k*nxr*ntr], &rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], conv, nxr, ntr, dt, 0);
+ timeDiff(conv, ntr, nxr, dt, fmin, fmax, -3);
+ for (i=0; i<nxr; i++) {
+ for (j=0; j<ntr/2; j++) {
+ Ghom[is*ntr*nxvr*nyvr*nzvr+(j+ntr/2)*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+j]/rho;
+ Ghom[is*ntr*nxvr*nyvr*nzvr+j*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+(j+ntr/2)]/rho;
+ }
+ }
+ }
+ }
+ }
+ else if (scheme==7) { //Marchenko representation with multiple shot gathers without time-reversal
+ for (k = 0; k < nyr; k++) {
+ depthDiff(&rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], ntr, nxr, dt, dx, fmin, fmax, cp, 1);
+ for (is=0; is<nshots; is++) {
+ convol(&shotdata[is*nyr*nxr*ntr+k*nxr*ntr], &rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], conv, nxr, ntr, dt, 0);
+ timeDiff(conv, ntr, nxr, dt, fmin, fmax, -1);
+ for (i=0; i<nxr; i++) {
+ for (j=0; j<ntr/2; j++) {
+ Ghom[is*ntr*nxvr*nyvr*nzvr+(j+ntr/2)*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+j]/rho;
+ Ghom[is*ntr*nxvr*nyvr*nzvr+j*nxvr*nyvr*nzvr+l*nzvr+ir] += scl*conv[i*ntr+(j+ntr/2)]/rho;
+ }
+ }
+ }
+ }
+ }
+ else if (scheme==8) { // f1+ redatuming 0=f1p 1=f1m
+ for (k = 0; k < nyr; k++) {
+ convol2(&shotdata[0*nyr*nxr*ntr+k*nxr*ntr], &rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], conv, nxr, ntr, dt, fmin, fmax, 1);
+ convol2(&shotdata[1*nyr*nxr*ntr+k*nxr*ntr], &rcvdata[(l+1)*nyr*nxr*ntr+k*nxr*ntr], tmp1, nxr, ntr, dt, fmin, fmax, 1);
+ for (i=0; i<nxr; i++) {
+ for (j=0; j<ntr/2; j++) {
+ Ghom[(j+ntr/2)*nxvr*nyvr*nzvr+l*nzvr+ir] -= 2.0*scl*(conv[i*ntr+j] + tmp1[i*ntr+j])/rho;
+ Ghom[j*nxvr*nyvr*nzvr+l*nzvr+ir] -= 2.0*scl*(conv[i*ntr+(j+ntr/2)] + tmp1[i*ntr+(j+ntr/2)])/rho;
+ }
+ }
+ }
+ }
+ else if (scheme==9) { // f1- redatuming 0=f1p 1=f1m
+ for (k = 0; k < nyr; k++) {
+ convol2(&shotdata[0*nyr*nxr*ntr+k*nxr*ntr], &rcvdata[(l+1)*nyr*nxr*ntr+k*nxr*ntr], conv, nxr, ntr, dt, fmin, fmax, 1);
+ convol2(&shotdata[1*nyr*nxr*ntr+k*nxr*ntr], &rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], tmp1, nxr, ntr, dt, fmin, fmax, 1);
+ for (i=0; i<nxr; i++) {
+ for (j=0; j<ntr/2; j++) {
+ Ghom[(j+ntr/2)*nxvr*nyvr*nzvr+l*nzvr+ir] -= 2.0*scl*(conv[i*ntr+j] + tmp1[i*ntr+j])/rho;
+ Ghom[j*nxvr*nyvr*nzvr+l*nzvr+ir] -= 2.0*scl*(conv[i*ntr+(j+ntr/2)] + tmp1[i*ntr+(j+ntr/2)])/rho;
+ }
+ }
+ }
+ }
+ else if (scheme==10) { // 2i IM(f1) redatuming
+ for (k = 0; k < nyr; k++) {
+ convol2(&shotdata[k*nxr*ntr], &rcvdata[(l+1)*nyr*nxr*ntr+k*nxr*ntr], conv, nxr, ntr, dt, fmin, fmax, 2);
+ convol2(&shotdata_jkz[k*nxr*ntr], &rcvdata[l*nyr*nxr*ntr+k*nxr*ntr], tmp1, nxr, ntr, dt, fmin, fmax, 2);
+ for (i=0; i<nxr; i++) {
+ for (j=0; j<ntr/2; j++) {
+ Ghom[(j+ntr/2)*nxvr*nyvr*nzvr+l*nzvr+ir] += 4.0*scl*(conv[i*ntr+j] - tmp1[i*ntr+j])/rho;
+ Ghom[j*nxvr*nyvr*nzvr+l*nzvr+ir] += 4.0*scl*(conv[i*ntr+(j+ntr/2)] - tmp1[i*ntr+(j+ntr/2)])/rho;
+ }
+ }
}
}
}
-
- count+=1;
- if (verbose) vmess("Creating Homogeneous Green's function at depth %li from %li depths",count,nzs);
- }
- free(conv); free(indata); free(hdr_in); free(conv2);
- if (scheme==2) {
- free(tmp1);free(tmp2);
+ if (verbose) vmess("Creating Homogeneous Green's function at depth %li from %li depths",ir+1,nzvr);
+ free(conv); free(rcvdata); free(hdr_rcv);
+ if (scheme==5) {
+ free(tmp1); free(tmp2);
+ }
+ if (scheme==6 || scheme==8 || scheme==9 || scheme==10) free(tmp1);
}
-}
- free(shotdata);
-
- if (verbose) vmess("nxs: %li nxz: %li f1: %.7f",nxs,nzs,f1);
- ntshift=0;
+ free(shotdata);
- if (ntmax > 0) {
- if (ntmax < nt) {
- ntshift = (nt-ntmax)/2;
- if (verbose) vmess("Time shifted %li samples",ntshift);
- nt=ntmax;
- }
- else {
- if (verbose) vmess("Max time samples larger than original samples");
- }
- }
+ if (nxvr>1) dxrcv = (float)((xvr[nzvr] - xvr[0])/1000.0);
+ else dxrcv = 1.0;
+ if (nyvr>1) dyrcv = (float)((yvr[nxvr*nzvr] - yvr[0])/1000.0);
+ else dyrcv = 1.0;
+ if (nzvr>1) dzrcv = (float)((zvr[1] - zvr[0])/1000.0);
+ else dzrcv = 1.0;
fp_out = fopen(fout, "w+");
-
- for (ir = 0; ir < nt; ir++) {
- for (ix = 0; ix < nxs; ix++) {
- hdr_out[ix].fldr = ir+1;
- hdr_out[ix].tracl = ir*nxs+ix+1;
- hdr_out[ix].tracf = ix+1;
- hdr_out[ix].scalco = hdr_shot[0].scalco;
- hdr_out[ix].scalel = hdr_shot[0].scalel;
- hdr_out[ix].sdepth = hdr_shot[0].sdepth;
- hdr_out[ix].trid = 1;
- hdr_out[ix].ns = nzs;
- hdr_out[ix].trwf = nxs;
- hdr_out[ix].ntr = hdr_out[ix].fldr*hdr_out[ix].trwf;
- hdr_out[ix].f1 = f1;
- hdr_out[ix].f2 = f2/1000;
- hdr_out[ix].dt = dt*(1E6);
- hdr_out[ix].d1 = dzrcv;
- hdr_out[ix].d2 = dxrcv;
- hdr_out[ix].sx = hdr_shot[0].sx;
- hdr_out[ix].gx = (int)roundf(f2 + (ix*hdr_out[ix].d2)*1000.0);
- hdr_out[ix].offset = (hdr_out[ix].gx - hdr_out[ix].sx)/1000.0;
+
+ for (ir = 0; ir < ntr; ir++) {
+ for (is = 0; is < nyvr; is++) {
+ for (ix = 0; ix < nxvr; ix++) {
+ hdr_out[is*nxvr+ix].fldr = ir-ntr/2;
+ hdr_out[is*nxvr+ix].tracl = ir*nyvr*nxvr+is*nxvr+ix+1;
+ hdr_out[is*nxvr+ix].tracf = is*nxvr+ix+1;
+ hdr_out[is*nxvr+ix].scalco = 1000.0;
+ hdr_out[is*nxvr+ix].scalel = 1000.0;
+ hdr_out[is*nxvr+ix].sdepth = hdr_shot[0].sdepth;
+ hdr_out[is*nxvr+ix].selev = -hdr_shot[0].sdepth;
+ hdr_out[is*nxvr+ix].trid = 1;
+ hdr_out[is*nxvr+ix].ns = nzvr;
+ hdr_out[is*nxvr+ix].trwf = nxvr*nyvr;
+ hdr_out[is*nxvr+ix].ntr = (ir+1)*hdr_out[is*nxvr+ix].trwf;
+ hdr_out[is*nxvr+ix].f1 = ((float)(zvr[0]/1000.0));
+ hdr_out[is*nxvr+ix].f2 = ((float)(xvr[0]/1000.0));
+ hdr_out[is*nxvr+ix].dt = dt*(1E6);
+ hdr_out[is*nxvr+ix].d1 = dzrcv;
+ hdr_out[is*nxvr+ix].d2 = dxrcv;
+ hdr_out[is*nxvr+ix].sx = hdr_shot[0].sx;
+ hdr_out[is*nxvr+ix].sy = hdr_shot[0].sy;
+ hdr_out[is*nxvr+ix].gx = xvr[ix*nzvr];
+ hdr_out[is*nxvr+ix].gy = yvr[is*nxvr*nzvr];
+ hdr_out[is*nxvr+ix].offset = (hdr_out[is*nxvr+ix].gx - hdr_out[is*nxvr+ix].sx)/1000.0;
+ }
}
- ret = writeData3D(fp_out, &Ghom[(ir+ntshift)*nxs*nzs], hdr_out, nzs, nxs);
+ ret = writeData3D(fp_out, &Ghom[ir*nyvr*nxvr*nzvr], hdr_out, nzvr, nxvr*nyvr);
if (ret < 0 ) verr("error on writing output file.");
}
fclose(fp_out);
+ free(hdr_shot);
return 0;
}
void convol(float *data1, float *data2, float *con, long nrec, long nsam, float dt, long shift)
{
- long i, j, n, optn, nfreq, sign;
- float df, dw, om, tau, scl;
- float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1, *rdata2;
- complex *cdata1, *cdata2, *ccon, tmp;
-
- optn = optncr(nsam);
- nfreq = optn/2+1;
+ long i, j, n, optn, nfreq, sign;
+ float df, dw, om, tau, scl;
+ float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1, *rdata2;
+ complex *cdata1, *cdata2, *ccon, tmp;
+ optn = loptncr(nsam);
+ nfreq = optn/2+1;
- cdata1 = (complex *)malloc(nfreq*nrec*sizeof(complex));
- if (cdata1 == NULL) verr("memory allocation error for cdata1");
- cdata2 = (complex *)malloc(nfreq*nrec*sizeof(complex));
- if (cdata2 == NULL) verr("memory allocation error for cdata2");
- ccon = (complex *)malloc(nfreq*nrec*sizeof(complex));
- if (ccon == NULL) verr("memory allocation error for ccov");
-
- rdata1 = (float *)malloc(optn*nrec*sizeof(float));
- if (rdata1 == NULL) verr("memory allocation error for rdata1");
- rdata2 = (float *)malloc(optn*nrec*sizeof(float));
- if (rdata2 == NULL) verr("memory allocation error for rdata2");
-
- /* pad zeroes until Fourier length is reached */
- pad_data(data1, nsam, nrec, optn, rdata1);
- pad_data(data2, nsam, nrec, optn, rdata2);
-
- /* forward time-frequency FFT */
- sign = -1;
- rcmfft(&rdata1[0], &cdata1[0], optn, nrec, optn, nfreq, sign);
- rcmfft(&rdata2[0], &cdata2[0], optn, nrec, optn, nfreq, sign);
-
- /* apply convolution */
- p1r = (float *) &cdata1[0];
- p2r = (float *) &cdata2[0];
- qr = (float *) &ccon[0].r;
- p1i = p1r + 1;
- p2i = p2r + 1;
- qi = qr + 1;
- n = nrec*nfreq;
- for (j = 0; j < n; j++) {
- *qr = (*p2r**p1r-*p2i**p1i);
- *qi = (*p2r**p1i+*p2i**p1r);
- qr += 2;
- qi += 2;
- p1r += 2;
- p1i += 2;
- p2r += 2;
- p2i += 2;
- }
- free(cdata1);
- free(cdata2);
-
- if (shift==1) {
- df = 1.0/(dt*optn);
- dw = 2*PI*df;
- tau = dt*(nsam/2);
- for (j = 0; j < nrec; j++) {
- om = 0.0;
- for (i = 0; i < nfreq; i++) {
- tmp.r = ccon[j*nfreq+i].r*cos(om*tau) + ccon[j*nfreq+i].i*sin(om*tau);
- tmp.i = ccon[j*nfreq+i].i*cos(om*tau) - ccon[j*nfreq+i].r*sin(om*tau);
- ccon[j*nfreq+i] = tmp;
- om += dw;
- }
- }
- }
- if (shift==-2) {
- for (j = 0; j < nrec; j++) {
- for (i = 0; i < nfreq; i++) {
- ccon[j*nfreq+i].r = ccon[j*nfreq+i].i;
- ccon[j*nfreq+i].i = 0.0;
- }
- }
- }
+
+ cdata1 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata1 == NULL) verr("memory allocation error for cdata1");
+ cdata2 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata2 == NULL) verr("memory allocation error for cdata2");
+ ccon = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (ccon == NULL) verr("memory allocation error for ccov");
+
+ rdata1 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata1 == NULL) verr("memory allocation error for rdata1");
+ rdata2 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata2 == NULL) verr("memory allocation error for rdata2");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data1, nsam, nrec, optn, rdata1);
+ pad_data(data2, nsam, nrec, optn, rdata2);
+
+ /* forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata1[0], &cdata1[0], (int)optn, (int)nrec, (int)optn, (int)nfreq, (int)sign);
+ rcmfft(&rdata2[0], &cdata2[0], (int)optn, (int)nrec, (int)optn, (int)nfreq, (int)sign);
+
+ /* apply convolution */
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ qr = (float *) &ccon[0].r;
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ qi = qr + 1;
+ n = nrec*nfreq;
+ for (j = 0; j < n; j++) {
+ *qr = (*p2r**p1r-*p2i**p1i);
+ *qi = (*p2r**p1i+*p2i**p1r);
+ qr += 2;
+ qi += 2;
+ p1r += 2;
+ p1i += 2;
+ p2r += 2;
+ p2i += 2;
+ }
+ free(cdata1);
+ free(cdata2);
+
+ if (shift) {
+ df = 1.0/(dt*optn);
+ dw = 2*PI*df;
+// tau = 1.0/(2.0*df);
+ tau = dt*(nsam/2);
+ for (j = 0; j < nrec; j++) {
+ om = 0.0;
+ for (i = 0; i < nfreq; i++) {
+ tmp.r = ccon[j*nfreq+i].r*cos(om*tau) + ccon[j*nfreq+i].i*sin(om*tau);
+ tmp.i = ccon[j*nfreq+i].i*cos(om*tau) - ccon[j*nfreq+i].r*sin(om*tau);
+ ccon[j*nfreq+i] = tmp;
+ om += dw;
+ }
+ }
+ }
/* inverse frequency-time FFT and scale result */
- sign = 1;
- scl = 1.0/((float)(optn));
- crmfft(&ccon[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
- scl_data(rdata1,optn,nrec,scl,con,nsam);
-
- free(ccon);
- free(rdata1);
- free(rdata2);
- return;
-}
-
-void pad_data(float *data, long nsam, long nrec, long nsamout, float *datout)
-{
- long it,ix;
- for (ix=0;ix<nrec;ix++) {
- for (it=0;it<nsam;it++)
- datout[ix*nsamout+it]=data[ix*nsam+it];
- for (it=nsam;it<nsamout;it++)
- datout[ix*nsamout+it]=0.0;
- }
-}
-
-void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout)
-{
- long it,ix;
- for (ix = 0; ix < nrec; ix++) {
- for (it = 0 ; it < nsamout ; it++)
- datout[ix*nsamout+it] = scl*data[ix*nsam+it];
- }
-}
-
-void corr(float *data1, float *data2, float *cov, long nrec, long nsam, float dt, long shift)
-{
- long i, j, n, optn, nfreq, sign;
- float df, dw, om, tau, scl;
- float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1, *rdata2;
- complex *cdata1, *cdata2, *ccov, tmp;
-
- optn = optncr(nsam);
- nfreq = optn/2+1;
-
- cdata1 = (complex *)malloc(nfreq*nrec*sizeof(complex));
- if (cdata1 == NULL) verr("memory allocation error for cdata1");
- cdata2 = (complex *)malloc(nfreq*nrec*sizeof(complex));
- if (cdata2 == NULL) verr("memory allocation error for cdata2");
- ccov = (complex *)malloc(nfreq*nrec*sizeof(complex));
- if (ccov == NULL) verr("memory allocation error for ccov");
-
- rdata1 = (float *)malloc(optn*nrec*sizeof(float));
- if (rdata1 == NULL) verr("memory allocation error for rdata1");
- rdata2 = (float *)malloc(optn*nrec*sizeof(float));
- if (rdata2 == NULL) verr("memory allocation error for rdata2");
-
- /* pad zeroes until Fourier length is reached */
- pad_data(data1, nsam, nrec, optn, rdata1);
- pad_data(data2, nsam, nrec, optn, rdata2);
-
- /* forward time-frequency FFT */
- sign = -1;
- rcmfft(&rdata1[0], &cdata1[0], optn, nrec, optn, nfreq, sign);
- rcmfft(&rdata2[0], &cdata2[0], optn, nrec, optn, nfreq, sign);
-
- /* apply correlation */
- p1r = (float *) &cdata1[0];
- p2r = (float *) &cdata2[0];
- qr = (float *) &ccov[0].r;
- p1i = p1r + 1;
- p2i = p2r + 1;
- qi = qr + 1;
- n = nrec*nfreq;
- for (j = 0; j < n; j++) {
- *qr = (*p1r * *p2r + *p1i * *p2i);
- *qi = (*p1i * *p2r - *p1r * *p2i);
- qr += 2;
- qi += 2;
- p1r += 2;
- p1i += 2;
- p2r += 2;
- p2i += 2;
- }
- free(cdata1);
- free(cdata2);
-
- /* shift t=0 to middle of time window (nsam/2)*/
- if (shift) {
- df = 1.0/(dt*optn);
- dw = 2*PI*df;
- tau = dt*(nsam/2);
-
- for (j = 0; j < nrec; j++) {
- om = 0.0;
- for (i = 0; i < nfreq; i++) {
- tmp.r = ccov[j*nfreq+i].r*cos(om*tau) + ccov[j*nfreq+i].i*sin(om*tau);
- tmp.i = ccov[j*nfreq+i].i*cos(om*tau) - ccov[j*nfreq+i].r*sin(om*tau);
- ccov[j*nfreq+i] = tmp;
- om += dw;
- }
- }
- }
-
- /* inverse frequency-time FFT and scale result */
- sign = 1;
- scl = 1.0/(float)optn;
- crmfft(&ccov[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
- scl_data(rdata1,optn,nrec,scl,cov,nsam);
-
- free(ccov);
- free(rdata1);
- free(rdata2);
- return;
+ sign = 1;
+ scl = 1.0/((float)(optn));
+ crmfft(&ccon[0], &rdata1[0], (int)optn, (int)nrec, (int)nfreq, (int)optn, (int)sign);
+ scl_data(rdata1,optn,nrec,scl,con,nsam);
+
+ free(ccon);
+ free(rdata1);
+ free(rdata2);
+ return;
}
void timeDiff(float *data, long nsam, long nrec, float dt, float fmin, float fmax, long opt)
@@ -596,12 +704,19 @@ void timeDiff(float *data, long nsam, long nrec, float dt, float fmin, float fma
cdatascl[ix*nfreq+iom].i = -om*cdata[ix*nfreq+iom].r;
}
}
- else if (opt == -2) {
+ else if (opt == -2) {
for (iom = iomin ; iom < iomax ; iom++) {
om = 4.0/(deltom*iom);
cdatascl[ix*nfreq+iom].r = om*cdata[ix*nfreq+iom].r;
cdatascl[ix*nfreq+iom].i = om*cdata[ix*nfreq+iom].i;
}
+ }
+ else if (opt == -3) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = 1.0/(deltom*iom);
+ cdatascl[ix*nfreq+iom].r = 2*om*cdata[ix*nfreq+iom].i;
+ cdatascl[ix*nfreq+iom].i = 0.0;
+ }
}
}
free(cdata);
@@ -620,7 +735,7 @@ void timeDiff(float *data, long nsam, long nrec, float dt, float fmin, float fma
void depthDiff(float *data, long nsam, long nrec, float dt, float dx, float fmin, float fmax, float c, long opt)
{
- long optn, iom, iomin, iomax, nfreq, ix, ikx, nkx, ikxmax;
+ long 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;
@@ -754,6 +869,7 @@ void pad2d_data(float *data, long nsam, long nrec, long nsamout, long nrecout, f
datout[ix*nsam+it]=0.0;
}
}
+
void conjugate(float *data, long nsam, long nrec, float dt)
{
long optn, nfreq, j, ix, it, sign, ntdiff;
@@ -791,9 +907,451 @@ void conjugate(float *data, long nsam, long nrec, float dt)
data[ix*nsam+it] = scl*rdata[ix*optn+it+ntdiff];
}
//scl_data(rdata,optn,nrec,scl,data,nsam);
-
- free(cdata);
+
+ free(cdata);
free(rdata);
+ return;
+}
+
+void convol2(float *data1, float *data2, float *con, long nrec, long nsam, float dt, float fmin, float fmax, long opt)
+{
+ long optn, iom, iomin, iomax, nfreq, ix, sign, i, j, n;
+ float omin, omax, deltom, om, df, dw, tau, scl;
+ float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1, *rdata2;
+ complex *cdata1, *cdata2, *ccon, tmp, *cdatascl;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+ df = 1.0/(optn*dt);
+
+ cdata1 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata1 == NULL) verr("memory allocation error for cdata1");
+ cdata2 = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata2 == NULL) verr("memory allocation error for cdata2");
+ ccon = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (ccon == NULL) verr("memory allocation error for ccov");
+
+ rdata1 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata1 == NULL) verr("memory allocation error for rdata1");
+ rdata2 = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata2 == NULL) verr("memory allocation error for rdata2");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data1, nsam, nrec, optn, rdata1);
+ pad_data(data2, nsam, nrec, optn, rdata2);
+
+ /* forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata1[0], &cdata1[0], (int)optn, (int)nrec, (int)optn, (int)nfreq, (int)sign);
+ rcmfft(&rdata2[0], &cdata2[0], (int)optn, (int)nrec, (int)optn, (int)nfreq, (int)sign);
+
+ /* apply convolution */
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ qr = (float *) &ccon[0].r;
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ qi = qr + 1;
+ n = nrec*nfreq;
+ for (j = 0; j < n; j++) {
+ *qr = (*p2r**p1r-*p2i**p1i);
+ *qi = (*p2r**p1i+*p2i**p1r);
+ qr += 2;
+ qi += 2;
+ p1r += 2;
+ p1i += 2;
+ p2r += 2;
+ p2i += 2;
+ }
+ free(cdata1);
+ free(cdata2);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+ iomin = (long)MIN((omin/deltom), (nfreq));
+ iomin = MAX(iomin, 1);
+ iomax = MIN((long)(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==1) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = 1.0/(deltom*iom);
+ cdatascl[ix*nfreq+iom].r = om*ccon[ix*nfreq+iom].i;
+ cdatascl[ix*nfreq+iom].i = -om*ccon[ix*nfreq+iom].r;
+ }
+ }
+ else if (opt==2) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = 1.0/(deltom*iom);
+ cdatascl[ix*nfreq+iom].r = 0.0;
+ cdatascl[ix*nfreq+iom].i = -om*ccon[ix*nfreq+iom].r;
+ }
+ }
+ }
+ free(ccon);
+
+ /* Inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&cdatascl[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata1,optn,nrec,scl,con,nsam);
+
+ free(cdatascl);
+ free(rdata1);
+ free(rdata2);
+
+ return;
+ return;
+}
+
+void pad_data(float *data, long nsam, long nrec, long nsamout, float *datout)
+{
+ long it,ix;
+ for (ix=0;ix<nrec;ix++) {
+ for (it=0;it<nsam;it++)
+ datout[ix*nsamout+it]=data[ix*nsam+it];
+ for (it=nsam;it<nsamout;it++)
+ datout[ix*nsamout+it]=0.0;
+ }
+}
+
+void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout)
+{
+ long it,ix;
+ for (ix = 0; ix < nrec; ix++) {
+ for (it = 0 ; it < nsamout ; it++)
+ datout[ix*nsamout+it] = scl*data[ix*nsam+it];
+ }
+}
+
+long zfpdecompress(float* data, long nx, long ny, long nz, long comp, double tolerance, FILE *file)
+{
+ zfp_field* field = NULL;
+ zfp_stream* zfp = NULL;
+ bitstream* stream = NULL;
+ zfp_exec_policy exec = zfp_exec_serial;
+ size_t nread, compsize;
+ void *buffer;
+
+ zfp = zfp_stream_open(NULL);
+ field = zfp_field_alloc();
+ compsize = comp;
+
+ buffer = malloc(compsize);
+ if (!buffer) {
+ fprintf(stderr, "cannot allocate memory\n");
+ return EXIT_FAILURE;
+ }
+ nread = fread((uchar*)buffer, 1, compsize, file);
+ assert(nread==compsize);
+
+ stream = stream_open(buffer, compsize);
+ if (!stream) {
+ fprintf(stderr, "cannot open compressed stream\n");
+ return EXIT_FAILURE;
+ }
+ zfp_stream_set_bit_stream(zfp, stream);
+
+ zfp_field_set_type(field, zfp_type_float);
+ if (ny<2) zfp_field_set_size_2d(field, (uint)nz, (uint)nx);
+ else zfp_field_set_size_3d(field, (uint)nz, (uint)nx, (uint)ny);
+
+ zfp_stream_set_accuracy(zfp, tolerance);
+
+ if (!zfp_stream_set_execution(zfp, exec)) {
+ fprintf(stderr, "serial execution not available\n");
+ return EXIT_FAILURE;
+ }
+
+ zfp_stream_rewind(zfp);
+
+ if (!zfp_stream_set_execution(zfp, exec)) {
+ fprintf(stderr, "serial execution not available\n");
+ return EXIT_FAILURE;
+ }
+
+ zfp_field_set_pointer(field, (void *)data);
+
+ while (!zfp_decompress(zfp, field)) {
+ /* fall back on serial decompression if execution policy not supported */
+ if (zfp_stream_execution(zfp) != zfp_exec_serial) {
+ if (!zfp_stream_set_execution(zfp, zfp_exec_serial)) {
+ fprintf(stderr, "cannot change execution policy\n");
+ return EXIT_FAILURE;
+ }
+ }
+ else {
+ fprintf(stderr, "decompression failed\n");
+ return EXIT_FAILURE;
+ }
+ }
+
+ return 1;
+}
+
+void getVirRec(char *filename, long *nxs, long *nys, long *nxr, long *nyr, long *nt)
+{
+ FILE *fp;
+ segy hdr;
+ size_t nread;
+ long sx, sy, gx, gy, end_of_file, nshots, gx0, gy0, itrace, isx, isy, ishot, tsize;
+ long ny;
+
+ end_of_file = 0;
+
+ fp = fopen( filename, "r" );
+ if ( fp == NULL ) verr("Could not open %s",filename);
+ nread = fread(&hdr, 1, TRCBYTES, fp);
+ if (nread != TRCBYTES) verr("Could not read the header of the input file");
+
+ *nxs = 1;
+ *nys = 1;
+ *nt = hdr.ns;
+ *nyr = 1;
+ *nxr = 1;
+ ny = 1;
+ sx = hdr.sx;
+ sy = hdr.sy;
+ gx = hdr.gx;
+ gy = hdr.gy;
+ gx0 = gx;
+ gy0 = gy;
+ itrace = 1;
+ ishot = 1;
+ tsize = hdr.ns*sizeof(float);
+
+ fseek(fp, 0, SEEK_SET);
+
+ while (!end_of_file) {
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) {
+ end_of_file = 1;
+ break;
+ }
+ if (hdr.gy != gy) {
+ gy = hdr.gy;
+ ny++;
+ }
+ if ((sx != hdr.sx) || (sy != hdr.sy)) {
+ end_of_file = 1;
+ break;
+ }
+ itrace++;
+ fseek(fp, tsize, SEEK_CUR);
+ }
+
+ *nyr = ny;
+ *nxr = itrace/(ny);
+ end_of_file = 0;
+ isx = 1;
+ isy = 1;
+ ny = 1;
+
+ fseek(fp, 0, SEEK_SET);
+
+ while (!end_of_file) {
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) {
+ end_of_file = 1;
+ break;
+ }
+ if (hdr.sx != sx) {
+ sx = hdr.sx;
+ ishot++;
+ }
+ if (hdr.sy != sy) {
+ sy = hdr.sy;
+ ny++;
+ }
+ fseek(fp, tsize, SEEK_CUR);
+ }
+
+ *nys = ny;
+ *nxs = ishot/(ny);
+
+ return;
+}
+
+void getFileInfo3Dzfp(char *filename, long *n1, long *n2, long *n3, long *ngath,
+ float *d1, float *d2, float *d3, float *f1, float *f2, float *f3,
+ float *fmin, float *fmax, float *scl, long *nxm)
+{
+ FILE *fp_in;
+ size_t nread;
+ zfpmar zfpm;
+ zfptop zfpt;
+ long ishot, compsize;
+
+ fp_in = fopen(filename, "r");
+ if (fp_in==NULL) {
+ fprintf(stderr,"input file %s has an error\n", fp_in);
+ perror("error in opening file: ");
+ fflush(stderr);
+ return;
+ }
+
+ nread = fread(&zfpt, 1, TOPBYTES, fp_in);
+ assert(nread == TOPBYTES);
+ *ngath = zfpt.ns;
+ *n1 = zfpt.nt;
+ *n2 = 1;
+ *n3 = 1;
+ *fmin = zfpt.fmin;
+ *fmax = zfpt.fmax;
+ *f1 = zfpt.fz;
+ *f2 = zfpt.fx;
+ *f3 = zfpt.fy;
+ *d1 = zfpt.dz;
+ *d2 = zfpt.dx;
+ *d3 = zfpt.dy;
+ *nxm = 1;
+
+ if (zfpt.scale < 0.0) *scl = 1.0/fabs((float)zfpt.scale);
+ else if (zfpt.scale == 0.0) *scl = 1.0;
+ else *scl = zfpt.scale;
+
+ compsize = 0;
+
+ for (ishot=0; ishot<zfpt.ns; ishot++) {
+ fseek(fp_in, compsize, SEEK_CUR);
+ nread = fread(&zfpm, 1, MARBYTES, fp_in);
+ assert(nread == MARBYTES);
+ *n2 = MAX(*n2,zfpm.nx);
+ *n3 = MAX(*n3,zfpm.ny);
+ compsize = zfpm.compsize;
+ }
+
+ *nxm = (*n2)*(*n3);
+
+ return;
+}
+
+void getVirReczfp(char *filename, long *nxs, long *nys, long *nxr, long *nyr, long *nt)
+{
+ FILE *fp_in;
+ size_t nread;
+ zfpmar zfpm;
+ zfptop zfpt;
+ long ishot, compsize, sy, ny;
+
+ fp_in = fopen(filename, "r");
+ if (fp_in==NULL) {
+ fprintf(stderr,"input file %s has an error\n", fp_in);
+ perror("error in opening file: ");
+ fflush(stderr);
+ return;
+ }
+
+ nread = fread(&zfpt, 1, TOPBYTES, fp_in);
+ assert(nread == TOPBYTES);
+ *nt = zfpt.nt;
+ *nxr = 1;
+ *nyr = 1;
+ sy = 123456;
+ ny = 0;
+ compsize = 0;
+
+ for (ishot=0; ishot<zfpt.ns; ishot++) {
+ fseek(fp_in, compsize, SEEK_CUR);
+ nread = fread(&zfpm, 1, MARBYTES, fp_in);
+ assert(nread == MARBYTES);
+ *nxr = MAX(*nxr,zfpm.nx);
+ *nyr = MAX(*nyr,zfpm.ny);
+ compsize = zfpm.compsize;
+ if (zfpm.sy != sy) {
+ sy = zfpm.sy;
+ ny++;
+ }
+ }
+
+ *nxs = zfpt.ns/ny;
+ *nys = ny;
+
+ return;
+}
+
+void readzfpdata(char *filename, float *data, long size)
+{
+ FILE *fp;
+ size_t nread;
+ zfpmar zfpm;
+ zfptop zfpt;
+ float tolerance;
+ long l;
+
+ 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;
+ }
+
+ fseek(fp, 0, SEEK_SET);
+ nread = fread( &zfpt, 1, TOPBYTES, fp );
+ assert(nread == TOPBYTES);
+ tolerance = zfpt.tolerance;
+
+ for (l=0; l<zfpt.ns; l++) {
+ nread = fread( &zfpm, 1, MARBYTES, fp );
+ if (nread != MARBYTES) { /* no more data in file */
+ break;
+ }
+ zfpdecompress(&data[l*size], zfpm.nx, zfpm.ny, zfpt.nt, zfpm.compsize, tolerance, fp);
+ }
+
+ fclose(fp);
+
+ return;
+}
+
+void getxyzzfp(char *filename, long *sx, long *sy, long *sz, long iz, long nz)
+{
+ FILE *fp_in;
+ size_t nread;
+ zfpmar zfpm;
+ zfptop zfpt;
+ long ishot, compsize;
+
+ fp_in = fopen(filename, "r");
+ if (fp_in==NULL) {
+ fprintf(stderr,"input file %s has an error\n", fp_in);
+ perror("error in opening file: ");
+ fflush(stderr);
+ return;
+ }
+
+ fseek(fp_in, 0, SEEK_SET);
+ nread = fread(&zfpt, 1, TOPBYTES, fp_in);
+ assert(nread == TOPBYTES);
+
+ compsize = 0;
+
+ for (ishot=0; ishot<zfpt.ns; ishot++) {
+ fseek(fp_in, compsize, SEEK_CUR);
+ nread = fread(&zfpm, 1, MARBYTES, fp_in);
+ assert(nread == MARBYTES);
+
+ sx[ishot*nz+iz] = zfpm.sx;
+ sy[ishot*nz+iz] = zfpm.sy;
+ sz[ishot*nz+iz] = labs(zfpm.sz);
+
+ compsize = zfpm.compsize;
+ }
+
return;
}
\ No newline at end of file
diff --git a/utils/Makefile b/utils/Makefile
index 4caf5ee..27cb7fa 100644
--- a/utils/Makefile
+++ b/utils/Makefile
@@ -6,7 +6,7 @@ include ../Make_include
#OPTC += -openmp
#OPTC += -g -O0
-ALL: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d mutesnap combine combine_induced reshape_su HomG snap2shot makeR1D
+ALL: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d mutesnap combine combine_induced reshape_su HomG snap2shot makeR1D pwshift
SRCM = \
makemod.c \
@@ -204,6 +204,16 @@ SRCSS = snap2shot.c \
name_ext.c \
readSnapData3D.c
+SRCPW = pwshift.c \
+ getFileInfo3D.c \
+ writeData3D.c \
+ wallclock_time.c \
+ getpars.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ readSnapData3D.c
+
OBJM = $(SRCM:%.c=%.o)
makemod: $(OBJM)
@@ -294,7 +304,12 @@ OBJSS = $(SRCSS:%.c=%.o)
snap2shot: $(OBJSS)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o snap2shot $(OBJSS) $(LIBS)
-install: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d mutesnap combine combine_induced reshape_su HomG snap2shot makeR1D
+OBJPW = $(SRCPW:%.c=%.o)
+
+pwshift: $(OBJPW)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o pwshift $(OBJPW) $(LIBS)
+
+install: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d mutesnap combine combine_induced reshape_su HomG snap2shot makeR1D pwshift
cp makemod $B
cp makewave $B
cp extendModel $B
@@ -313,9 +328,10 @@ install: makemod makewave extendModel fconv correigen green green3D basop syn2d
cp HomG $B
cp snap2shot $B
cp makeR1D $B
+ cp pwshift $B
clean:
- rm -f core $(OBJM) makemod $(OBJW) makewave $(OBJE) extendModel $(OBJF) fconv $(OBJG) $(OBJC) correigen green $(OBJG3) green3D $(OBJB) basop $(OBJJ) syn2d $(OBJS) mat2su $(OBJA) ftr1d $(OBJT) mutesnap $(OBJMS) combine $(OBJCO) makeR1D $(OBJMR) reshape_su $(OBJRS) combine_induced $(OBJCI) HomG $(OBJHG) snap2shot $(OBJSS)
+ rm -f core $(OBJM) makemod $(OBJW) makewave $(OBJE) extendModel $(OBJF) fconv $(OBJG) $(OBJC) correigen green $(OBJG3) green3D $(OBJB) basop $(OBJJ) syn2d $(OBJS) mat2su $(OBJA) ftr1d $(OBJT) mutesnap $(OBJMS) combine $(OBJCO) makeR1D $(OBJMR) reshape_su $(OBJRS) combine_induced $(OBJCI) HomG $(OBJHG) snap2shot $(OBJSS) pwshift $(OBJPW)
realclean: clean
- rm -f $B/makemod $B/makewave $B/extendModel $B/fconv $B/correigen $B/green $B/green3D $B/basop $B/syn2d $B/mat2su $B/ftr1d $B/mutesnap $B/combine $B/combine_induced $B/reshape_su $B/HomG $B/snap2shot $B/makeR1D
+ rm -f $B/makemod $B/makewave $B/extendModel $B/fconv $B/correigen $B/green $B/green3D $B/basop $B/syn2d $B/mat2su $B/ftr1d $B/mutesnap $B/combine $B/combine_induced $B/reshape_su $B/HomG $B/snap2shot $B/makeR1D $B/pwshift
diff --git a/utils/getFileInfo3D.c b/utils/getFileInfo3D.c
index eef408c..7bb6561 100644
--- a/utils/getFileInfo3D.c
+++ b/utils/getFileInfo3D.c
@@ -241,4 +241,4 @@ long disp_fileinfo3D(char *file, long n1, long n2, long n3, float f1, float f2,
vmess("*** fldr = %li sx = %li sy = %li", hdrs[0].fldr, hdrs[0].sx, hdrs[0].sy);
return 0;
-}
+}
\ No newline at end of file
--
GitLab