From 0af1642dc82a9e2e909d67ef7da7cdfe4d9a2e21 Mon Sep 17 00:00:00 2001
From: JBrackenhoff <J.A.Brackenhoff@tudelft.nl>
Date: Wed, 10 Jun 2020 19:20:30 +0200
Subject: [PATCH] HomG
---
FFTlib/Makefile | 10 +-
FFTlib/ccmfft.c | 2 +-
FFTlib/genfft.h | 3 +
Make_include_template | 9 +-
fdelmodc3D/applySource3D.c | 9 +
fdelmodc3D/boundaries3D.c | 3 +-
fdelmodc3D/elastic4dc_3D.c | 73 +-
fdelmodc3D/fdelmodc3D.c | 38 +
fdelmodc3D/fdelmodc3D.h | 6 +
fdelmodc3D/getParameters3D.c | 3 +
fdelmodc3D/readModel3D.c | 13 +-
include/genfft.h | 3 +
marchenko/applyMute_tshift.c | 4 +-
marchenko3D/applyMute3D.c | 188 ++++
marchenko3D/makeWindow3D.c | 33 +-
marchenko3D/marchenko3D.c | 140 ++-
raytime3d/getParameters3d.c | 2 +-
utils/HomG.c | 496 +++++++---
utils/HomG_backup13may2020.c | 1610 +++++++++++++++++++++++++++++++
utils/HomG_backup25mar2020.c | 867 +++++++++++++++++
utils/Makefile | 40 +-
utils/padmodel.c | 238 +++++
utils/pwshift.c | 107 +-
utils/pwshift_backup23apr2020.c | 375 +++++++
utils/testfft3d.c | 359 +++++++
utils/truncate.c | 120 +++
26 files changed, 4423 insertions(+), 328 deletions(-)
create mode 100644 utils/HomG_backup13may2020.c
create mode 100644 utils/HomG_backup25mar2020.c
create mode 100644 utils/padmodel.c
create mode 100644 utils/pwshift_backup23apr2020.c
create mode 100644 utils/testfft3d.c
create mode 100644 utils/truncate.c
diff --git a/FFTlib/Makefile b/FFTlib/Makefile
index bf53a5a..fbb7b2c 100644
--- a/FFTlib/Makefile
+++ b/FFTlib/Makefile
@@ -17,12 +17,14 @@ CSRC = \
crmfft.c \
rc2dfft.c \
cr2dfft.c \
- xt2wx.c \
+ xt2wx.c \
xt2wkx.c \
+ yxt2wkykx.c \
wkx2xt.c \
- wx2xt.c \
- pfafft.c \
- kiss_fft.c \
+ wkykx2yxt.c \
+ wx2xt.c \
+ pfafft.c \
+ kiss_fft.c \
fft_mayer.c \
wallclock_time.c \
optnumber.c \
diff --git a/FFTlib/ccmfft.c b/FFTlib/ccmfft.c
index cbb0d29..6c998d1 100644
--- a/FFTlib/ccmfft.c
+++ b/FFTlib/ccmfft.c
@@ -7,7 +7,7 @@ void dfti_status_print(MKL_LONG status);
/**
* NAME: ccmfft
*
-* DESCRIPTION: Multipple vector complex to complex FFT
+* DESCRIPTION: Multiple vector complex to complex FFT
*
* USAGE:
* void ccmfft(complex *data, int n1, int n2, int ld1, int sign)
diff --git a/FFTlib/genfft.h b/FFTlib/genfft.h
index 4c7ccc6..bc1ae21 100644
--- a/FFTlib/genfft.h
+++ b/FFTlib/genfft.h
@@ -107,6 +107,9 @@ void xt2wkx(REAL *rdata, complex *cdata, int nt, int nx, int ldr, int ldc, int x
void wx2xt(complex *cdata, REAL *rdata, int nt, int nx, int ldc, int ldr);
void wkx2xt(complex *cdata, REAL *rdata, int nt, int nx, int ldc, int ldr, int xorig);
+void yxt2wkykx(REAL *rdata, complex *cdata, long nt, long nx, long ny, long ldt, long ldx, long ldy, long xorig, long yorig);
+void wkykx2yxt(complex *cdata, REAL *rdata, long nt, long nx, long ny, long ldt, long ldx, long ldy, long xorig, long yorig);
+
/*
void dc1_fft(double *rdata, dcomplex *cdata, int n, int sign);
void cd1_fft(dcomplex *cdata, double *rdata, int n, int sign);
diff --git a/Make_include_template b/Make_include_template
index 50d0665..ab90618 100644
--- a/Make_include_template
+++ b/Make_include_template
@@ -107,8 +107,15 @@ BLAS = -L${MKLROOT}/lib/ -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm
#OPTC += -DACML440 -I$(AMDROOT)/include
#BLAS = -L$(AMDROOT)/lib -lacml -lfortran -lffio -lrt -lm
+# FOR ZFP LIBRARIES
+ZFPROOT = ${ROOT}/zfp
+ZFP = -L${ZFPROOT}/lib -lzfp
+ZFPINCL = ${ZFPROOT}/include
+OPTC += -I${ZFPINCL}
+
#LIBARIES
-LIBS= -L$L -lgenfft $(BLAS)
+LIBS= -L$L -lgenfft $(FFT) $(BLAS) $(ZFP)
+
########################################################################
diff --git a/fdelmodc3D/applySource3D.c b/fdelmodc3D/applySource3D.c
index 90f545d..4221f2b 100644
--- a/fdelmodc3D/applySource3D.c
+++ b/fdelmodc3D/applySource3D.c
@@ -217,34 +217,43 @@ long applySource3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l
if (src.type == 1) {
if (src.orient==1) { /* monopole */
txx[iy*n1*n2+ix*n1+iz] += src_ampl;
+ tyy[iy*n1*n2+ix*n1+iz] += src_ampl;
tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
}
else if (src.orient==2) { /* dipole +/- */
txx[iy*n1*n2+ix*n1+iz] += src_ampl;
+ tyy[iy*n1*n2+ix*n1+iz] += src_ampl;
tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
txx[iy*n1*n2+ix*n1+iz+1] -= src_ampl;
+ tyy[iy*n1*n2+ix*n1+iz+1] -= src_ampl;
tzz[iy*n1*n2+ix*n1+iz+1] -= src_ampl;
}
else if (src.orient==3) { /* dipole - + */
txx[iy*n1*n2+ix*n1+iz] += src_ampl;
+ tyy[iy*n1*n2+ix*n1+iz] += src_ampl;
tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
txx[iy*n1*n2+(ix-1)*n1+iz] -= src_ampl;
+ tyy[iy*n1*n2+(ix-1)*n1+iz] -= src_ampl;
tzz[iy*n1*n2+(ix-1)*n1+iz] -= src_ampl;
}
else if (src.orient==4) { /* dipole +/0/- */
if (iz > ibndz) {
txx[iy*n1*n2+ix*n1+iz-1]+= 0.5*src_ampl;
+ tyy[iy*n1*n2+ix*n1+iz-1]+= 0.5*src_ampl;
tzz[iy*n1*n2+ix*n1+iz-1]+= 0.5*src_ampl;
}
if (iz < mod.nz+ibndz-1) {
txx[iy*n1*n2+ix*n1+iz+1] -= 0.5*src_ampl;
+ tyy[iy*n1*n2+ix*n1+iz+1] -= 0.5*src_ampl;
tzz[iy*n1*n2+ix*n1+iz+1] -= 0.5*src_ampl;
}
}
else if (src.orient==5) { /* dipole + - */
txx[iy*n1*n2+ix*n1+iz] += src_ampl;
+ tyy[iy*n1*n2+ix*n1+iz] += src_ampl;
tzz[iy*n1*n2+ix*n1+iz] += src_ampl;
txx[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl;
+ tyy[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl;
tzz[iy*n1*n2+(ix+1)*n1+iz] -= src_ampl;
}
}
diff --git a/fdelmodc3D/boundaries3D.c b/fdelmodc3D/boundaries3D.c
index 12905a6..9ede67a 100644
--- a/fdelmodc3D/boundaries3D.c
+++ b/fdelmodc3D/boundaries3D.c
@@ -3084,6 +3084,7 @@ MID left mid mid
ize = mod.ieYz;
ib = (bnd.ntap+ixo-1);
+
#pragma omp for private(ix,iy,iz)
for (ix=ixo; ix<ixe; ix++) {
#pragma ivdep
@@ -3166,7 +3167,7 @@ MID left mid mid
#pragma ivdep
for (iy=iyo; iy<iye; iy++) {
for (iz=izo; iz<ize; iz++) {
- vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
+ vz[iy*n1*n2+ix*n1+iz] -= roz[iy][ix][iz]*(
c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] +
tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) +
diff --git a/fdelmodc3D/elastic4dc_3D.c b/fdelmodc3D/elastic4dc_3D.c
index f2639ad..8aa04a4 100644
--- a/fdelmodc3D/elastic4dc_3D.c
+++ b/fdelmodc3D/elastic4dc_3D.c
@@ -86,68 +86,71 @@ long elastic4dc_3D(modPar mod, srcPar src, wavPar wav, bndPar bnd, long itime, l
/* calculate vx for all grid points except on the virtual boundary*/
#pragma omp for private (ix, iy, iz) nowait schedule(guided,1)
- for (iy=mod.ioXy; iy<mod.ieXy; iy++) {
- for (ix=mod.ioXx; ix<mod.ieXx; ix++) {
-#pragma ivdep
+ for (iy=mod.ioXy; iy<mod.ieXy; iy++) {
+ for (ix=mod.ioXx; ix<mod.ieXx; ix++) {
+ #pragma ivdep
for (iz=mod.ioXz; iz<mod.ieXz; iz++) {
vx[iy*n2*n1+ix*n1+iz] -= rox[iy][ix][iz]*(
- c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] +
- txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
- txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) +
- c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
- txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
- txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
+ c1*(txx[iy*n2*n1+ix*n1+iz] - txx[iy*n2*n1+(ix-1)*n1+iz] +
+ txy[(iy+1)*n2*n1+ix*n1+iz] - txy[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+1] - txz[iy*n2*n1+ix*n1+iz]) +
+ c2*(txx[iy*n2*n1+(ix+1)*n1+iz] - txx[iy*n2*n1+(ix-2)*n1+iz] +
+ txy[(iy+2)*n2*n1+ix*n1+iz] - txy[(iy-1)*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+ix*n1+iz+2] - txz[iy*n2*n1+ix*n1+iz-1]) );
}
}
- }
+ }
/* calculate vy for all grid points except on the virtual boundary */
-#pragma omp for private (ix, iy, iz) schedule(guided,1)
+#pragma omp for private (ix, iy, iz) nowait schedule(guided,1)
for (iy=mod.ioYy; iy<mod.ieYy; iy++) {
- for (ix=mod.ioYx; ix<mod.ieYx; ix++) {
-#pragma ivdep
+ for (ix=mod.ioYx; ix<mod.ieYx; ix++) {
+ #pragma ivdep
for (iz=mod.ioYz; iz<mod.ieYz; iz++) {
vy[iy*n2*n1+ix*n1+iz] -= roy[iy][ix][iz]*(
- c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] +
- tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] +
- txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) +
- c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
- tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
- txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
+ c1*(tyy[iy*n2*n1+ix*n1+iz] - tyy[(iy-1)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+1] - tyz[iy*n2*n1+ix*n1+iz] +
+ txy[iy*n2*n1+(ix+1)*n1+iz] - txy[iy*n2*n1+ix*n1+iz]) +
+ c2*(tyy[(iy+1)*n2*n1+ix*n1+iz] - tyy[(iy-2)*n2*n1+ix*n1+iz] +
+ tyz[iy*n2*n1+ix*n1+iz+2] - tyz[iy*n2*n1+ix*n1+iz-1] +
+ txy[iy*n2*n1+(ix+2)*n1+iz] - txy[iy*n2*n1+(ix-1)*n1+iz]) );
}
}
- }
+ }
/* calculate vz for all grid points except on the virtual boundary */
-#pragma omp for private (ix, iy, iz) schedule(guided,1)
+#pragma omp for private (ix, iy, iz) schedule(guided,1)
for (iy=mod.ioZy; iy<mod.ieZy; iy++) {
- for (ix=mod.ioZx; ix<mod.ieZx; ix++) {
-#pragma ivdep
+ for (ix=mod.ioZx; ix<mod.ieZx; ix++) {
+ #pragma ivdep
for (iz=mod.ioZz; iz<mod.ieZz; iz++) {
vz[iy*n2*n1+ix*n1+iz] -= roz[iy][ix][iz]*(
- c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] +
- tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
- txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) +
- c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
- tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
- txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
+ c1*(tzz[iy*n2*n1+ix*n1+iz] - tzz[iy*n2*n1+ix*n1+iz-1] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+1)*n1+iz] - txz[iy*n2*n1+ix*n1+iz]) +
+ c2*(tzz[iy*n2*n1+ix*n1+iz+1] - tzz[iy*n2*n1+ix*n1+iz-2] +
+ tyz[(iy+1)*n2*n1+ix*n1+iz] - tyz[iy*n2*n1+ix*n1+iz] +
+ txz[iy*n2*n1+(ix+2)*n1+iz] - txz[iy*n2*n1+(ix-1)*n1+iz]) );
}
- }
- }
+ }
+ }
/* Add force source */
if (src.type > 5) {
applySource3D(mod, src, wav, bnd, itime, ixsrc, iysrc, izsrc, vx, vy, vz, tzz, tyy, txx, txz, txy, tyz, rox, roy, roz, l2m, src_nwav, verbose);
}
-
+
/* boundary condition clears velocities on boundaries */
+ // mod.ischeme=1;
boundariesP3D(mod, bnd, vx, vy, vz, tzz, tyy, txx, txz, txy, tyz, rox, roy, roz, l2m, lam, mul, itime, verbose);
+ // mod.ischeme=5;
/* calculate Txx/tzz for all grid points except on the virtual boundary */
-#pragma omp for private (ix, iy, iz, dvx, dvy, dvz) nowait schedule(guided,1)
- for (iy=mod.ioPy; iy<mod.iePy; iy++) {
- for (ix=mod.ioPx; ix<mod.iePx; ix++) {
+#pragma omp for private (ix, iy, iz) schedule(guided,1)
#pragma ivdep
+ for (ix=mod.ioPx; ix<mod.iePx; ix++) {
+ for (iy=mod.ioPy; iy<mod.iePy; iy++) {
+ #pragma ivdep
for (iz=mod.ioPz; iz<mod.iePz; iz++) {
dvx = c1*(vx[iy*n2*n1+(ix+1)*n1+iz] - vx[iy*n2*n1+ix*n1+iz]) +
c2*(vx[iy*n2*n1+(ix+2)*n1+iz] - vx[iy*n2*n1+(ix-1)*n1+iz]);
diff --git a/fdelmodc3D/fdelmodc3D.c b/fdelmodc3D/fdelmodc3D.c
index 57901b8..7cffe95 100644
--- a/fdelmodc3D/fdelmodc3D.c
+++ b/fdelmodc3D/fdelmodc3D.c
@@ -293,6 +293,7 @@ int main(int argc, char **argv)
srcPar src;
bndPar bnd;
shotPar shot;
+ FILE *fp;
float **src_nwav;
float ***rox, ***roy, ***roz, ***l2m, ***lam, ***mul;
float *tss, *tes, *tep, *p, *q, *r;
@@ -365,6 +366,43 @@ int main(int argc, char **argv)
readModel3D(mod, bnd, rox, roy, roz, l2m, lam, mul, tss, tes, tep);
+ // tss = (float *)calloc(sizem,sizeof(float));
+ // for (iz=0; iz<mod.naz; iz++) {
+ // for (iy=0; iy<mod.nay; iy++) {
+ // for (ix=0; ix<mod.nax; ix++) {
+ // tss[iy*n2*n1+ix*n1+iz] = roy[iy][ix][iz];
+ // }
+ // }
+ // }
+ // vmess("nax:%li nay:%li naz:%li",mod.nax,mod.nay,mod.naz);
+ // fp = fopen("roy.bin","w+");
+ // fwrite(tss,sizem,sizeof(float),fp);
+ // fclose(fp);
+ // for (iz=0; iz<mod.naz; iz++) {
+ // for (iy=0; iy<mod.nay; iy++) {
+ // for (ix=0; ix<mod.nax; ix++) {
+ // tss[iy*n2*n1+ix*n1+iz] = rox[iy][ix][iz];
+ // }
+ // }
+ // }
+ // vmess("nax:%li nay:%li naz:%li",mod.nax,mod.nay,mod.naz);
+ // fp = fopen("rox.bin","w+");
+ // fwrite(tss,sizem,sizeof(float),fp);
+ // fclose(fp);
+ // for (iz=0; iz<mod.naz; iz++) {
+ // for (iy=0; iy<mod.nay; iy++) {
+ // for (ix=0; ix<mod.nax; ix++) {
+ // tss[iy*n2*n1+ix*n1+iz] = roz[iy][ix][iz];
+ // }
+ // }
+ // }
+ // vmess("nax:%li nay:%li naz:%li",mod.nax,mod.nay,mod.naz);
+ // fp = fopen("roz.bin","w+");
+ // fwrite(tss,sizem,sizeof(float),fp);
+ // fclose(fp);
+ // free(tss);
+ // return 0;
+
/* read and/or define source wavelet(s) */
/* Using a random source, which can have a random length
diff --git a/fdelmodc3D/fdelmodc3D.h b/fdelmodc3D/fdelmodc3D.h
index 061df2f..c9cc9c0 100644
--- a/fdelmodc3D/fdelmodc3D.h
+++ b/fdelmodc3D/fdelmodc3D.h
@@ -173,6 +173,12 @@ typedef struct _sourcePar { /* Source Array Parameters */
long circle;
long array;
long random;
+ float Mxx;
+ float Mxy;
+ float Mxz;
+ float Myy;
+ float Myz;
+ float Mzz;
float *tbeg;
float *tend;
long multiwav;
diff --git a/fdelmodc3D/getParameters3D.c b/fdelmodc3D/getParameters3D.c
index 6a0cf32..eb81ea7 100644
--- a/fdelmodc3D/getParameters3D.c
+++ b/fdelmodc3D/getParameters3D.c
@@ -211,6 +211,7 @@ long getParameters3D(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar
dt = mod->dt;
if (!getparlong("src_type",&src->type)) src->type=1;
+ if (!getparfloat("src_Mxx",&src->Mxx)) src->Mxx=0.0;
if (!getparlong("src_orient",&src->orient)) {
src->orient=1;
if (getparlong("dipsrc",&src->orient)) src->orient=2; // for compatability with DELPHI's fdacmod
@@ -1035,6 +1036,8 @@ criteria we have imposed.*/
case 6 : fprintf(stderr,"Fx "); break;
case 7 : fprintf(stderr,"Fz "); break;
case 8 : fprintf(stderr,"P-potential"); break;
+ case 9 : fprintf(stderr,"Double-couple"); break;
+ case 10 : fprintf(stderr,"Moment tensor"); break;
}
fprintf(stderr,"\n");
if (wav->random) vmess("Wavelet has a random signature with fmax=%.2f", wav->fmax);
diff --git a/fdelmodc3D/readModel3D.c b/fdelmodc3D/readModel3D.c
index 61dcea1..401f75a 100644
--- a/fdelmodc3D/readModel3D.c
+++ b/fdelmodc3D/readModel3D.c
@@ -108,8 +108,15 @@ long readModel3D(modPar mod, bndPar bnd, float ***rox, float ***roy, float ***ro
assert(nread == TRCBYTES);
//cs = (float *)calloc(nz*ny*nx,sizeof(float));
+ cs = (float ***)alloc3float(mod);
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (iz = 0; iz < nz; iz++) {
+ cs[iy][ix][iz] = 0.0;
+ }
+ }
+ }
if (mod.ischeme>2 && mod.ischeme!=5) {
- cs = (float ***)alloc3float(mod);
fpcs = fopen( mod.file_cs, "r" );
assert( fpcs != NULL);
nread = fread(&hdr, 1, TRCBYTES, fpcs);
@@ -1727,9 +1734,7 @@ long readModel3D(modPar mod, bndPar bnd, float ***rox, float ***roy, float ***ro
free3float(cp);
free3float(ro);
- if (mod.ischeme>2 && mod.ischeme!=5) {
- free3float(cs);
- }
+ free3float(cs);
return 0;
}
diff --git a/include/genfft.h b/include/genfft.h
index 4c7ccc6..bc1ae21 100644
--- a/include/genfft.h
+++ b/include/genfft.h
@@ -107,6 +107,9 @@ void xt2wkx(REAL *rdata, complex *cdata, int nt, int nx, int ldr, int ldc, int x
void wx2xt(complex *cdata, REAL *rdata, int nt, int nx, int ldc, int ldr);
void wkx2xt(complex *cdata, REAL *rdata, int nt, int nx, int ldc, int ldr, int xorig);
+void yxt2wkykx(REAL *rdata, complex *cdata, long nt, long nx, long ny, long ldt, long ldx, long ldy, long xorig, long yorig);
+void wkykx2yxt(complex *cdata, REAL *rdata, long nt, long nx, long ny, long ldt, long ldx, long ldy, long xorig, long yorig);
+
/*
void dc1_fft(double *rdata, dcomplex *cdata, int n, int sign);
void cd1_fft(dcomplex *cdata, double *rdata, int n, int sign);
diff --git a/marchenko/applyMute_tshift.c b/marchenko/applyMute_tshift.c
index 9a15951..a2f956e 100644
--- a/marchenko/applyMute_tshift.c
+++ b/marchenko/applyMute_tshift.c
@@ -93,10 +93,10 @@ void applyMute_tshift( float *data, int *mute, int smooth, int above, int Nfoc,
imute = ixpos[i];
tmute = mute[isyn*nxs+imute];
ts = tsynW[isyn*nxs+imute];
- for (j = -2*ts+tmute-shift-smooth,l=0; j < -2*ts+tmute-shift; j++,l++) {
+ for (j = MAX(0,-2*ts+tmute-shift-smooth),l=0; j < MAX(0,-2*ts+tmute-shift); j++,l++) {
Nig[j] *= costaper[smooth-l-1];
}
- for (j = 0; j < -2*ts+tmute-shift-smooth-1; j++) {
+ for (j = 0; j < MAX(0,-2*ts+tmute-shift-smooth-1); j++) {
Nig[j] = 0.0;
}
for (j = nt+1-tmute+shift+smooth; j < nt; j++) {
diff --git a/marchenko3D/applyMute3D.c b/marchenko3D/applyMute3D.c
index 222117a..b1e0abe 100644
--- a/marchenko3D/applyMute3D.c
+++ b/marchenko3D/applyMute3D.c
@@ -11,6 +11,11 @@
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
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)
@@ -141,3 +146,186 @@ void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, l
return;
}
+void applyMute3D_tshift( 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 iter, long *tsynW)
+{
+ long i, j, l, isyn, nxys;
+ float *costaper, scl, *Nig;
+ long imute, tmute, ts;
+
+ nxys = nxs*nys;
+ Nig = (float *)malloc(nt*sizeof(float));
+
+ if (smooth) {
+ costaper = (float *)malloc(smooth*sizeof(float));
+ scl = M_PI/((float)smooth);
+ for (i=0; i<smooth; i++) {
+ costaper[i] = 0.5*(1.0+cos((i+1)*scl));
+ }
+ }
+
+ for (isyn = 0; isyn < Nfoc; isyn++) {
+ for (i = 0; i < npos; i++) {
+ if (iter % 2 == 0) {
+ for (j = 0; j < nt; j++) {
+ Nig[j] = data[isyn*nxys*nt+i*nt+j];
+ }
+ }
+ else { // reverse back in time
+ j=0;
+ Nig[j] = data[isyn*nxys*nt+i*nt+j];
+ for (j = 1; j < nt; j++) {
+ Nig[j] = data[isyn*nxys*nt+i*nt+nt-j];
+ }
+ }
+ if (above==1) {
+ 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++) {
+ Nig[j] = 0.0;
+ }
+ for (j = MAX(0,tmute-shift-smooth),l=0; j < MAX(0,tmute-shift); j++,l++) {
+ Nig[j] *= costaper[smooth-l-1];
+ }
+ }
+ else if (above==0){
+ imute = iypos[i]*nxs+ixpos[i];
+ tmute = mute[isyn*nxys+imute];
+ ts = tsynW[isyn*nxys+imute];
+ if (tmute >= nt/2) {
+ memset(&Nig[0],0, sizeof(float)*nt);
+ continue;
+ }
+ for (j = MAX(0,tmute-shift),l=0; j < MAX(0,tmute-shift+smooth); j++,l++) {
+ Nig[j] *= costaper[l];
+ }
+ for (j = MAX(0,tmute-shift+smooth+1); j < MIN(nt,nt+1-tmute+2*ts+shift-smooth); j++) {
+ Nig[j] = 0.0;
+ }
+ for (j = MIN(nt-1,nt-tmute+2*ts+shift-smooth),l=0; j < MIN(nt,nt-tmute+shift); j++,l++) {
+ Nig[j] *= costaper[smooth-l-1];
+ }
+ }
+ else if (above==-1) {
+ imute = iypos[i]*nxs+ixpos[i];
+ tmute = mute[isyn*nxys+imute];
+ ts = tsynW[isyn*nxys+imute];
+ for (j = ts+tmute-shift,l=0; j < ts+tmute-shift+smooth; j++,l++) {
+ Nig[j] *= costaper[l];
+ }
+ for (j = ts+tmute-shift+smooth; j < nt; j++) {
+ Nig[j] = 0.0;
+ }
+ }
+ else if (above==4) { //Psi gate which is the inverse of the Theta gate (above=0)
+ imute = iypos[i]*nxs+ixpos[i];
+ tmute = mute[isyn*nxys+imute];
+ ts = tsynW[isyn*nxys+imute];
+ for (j = MAX(0,-2*ts+tmute-shift-smooth),l=0; j < MAX(0,-2*ts+tmute-shift); j++,l++) {
+ Nig[j] *= costaper[smooth-l-1];
+ }
+ for (j = 0; j < MAX(0,-2*ts+tmute-shift-smooth-1); j++) {
+ Nig[j] = 0.0;
+ }
+ for (j = nt+1-tmute+shift+smooth; j < nt; j++) {
+ Nig[j] = 0.0;
+ }
+ for (j = nt-tmute+shift,l=0; j < nt-tmute+shift+smooth; j++,l++) {
+ Nig[j] *= costaper[l];
+ }
+ }
+/* To Do above==2 is not yet adapated for plane-waves */
+ else if (above==2){//Separates the direct part of the wavefield from the coda
+ imute = iypos[i]*nxs+ixpos[i];
+ tmute = mute[isyn*nxys+imute];
+ for (j = 0; j < tmute-shift-smooth; j++) {
+ data[isyn*nxys*nt+i*nt+j] = 0.0;
+ }
+ for (j = tmute-shift-smooth,l=0; j < tmute-shift; j++,l++) {
+ data[isyn*nxys*nt+i*nt+j] *= costaper[smooth-l-1];
+ }
+ for (j = tmute+shift,l=0; j < tmute+shift+smooth; j++,l++) {
+ data[isyn*nxys*nt+i*nt+j] *= costaper[l];
+ }
+ for (j = tmute+shift+smooth; j < nt; j++) {
+ data[isyn*nxys*nt+i*nt+j] = 0.0;
+ }
+ }
+
+ if (iter % 2 == 0) {
+ for (j = 0; j < nt; j++) {
+ data[isyn*nxys*nt+i*nt+j] = Nig[j];
+ }
+ }
+ else { // reverse back in time
+ j=0;
+ data[isyn*nxys*nt+i*nt+j] = Nig[j];
+ for (j = 1; j < nt; j++) {
+ data[isyn*nxys*nt+i*nt+j] = Nig[nt-j];
+ }
+ }
+ } /* end if ipos */
+ }
+
+ if (smooth) free(costaper);
+ free(Nig);
+
+ return;
+}
+
+void timeShift(float *data, long nsam, long nrec, float dt, float shift, float fmin, float fmax)
+{
+ long optn, iom, iomin, iomax, nfreq, ix, it;
+ float omin, omax, deltom, om, tom, df, *trace, scl;
+ complex *ctrace, ctmp;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+ df = 1.0/(optn*dt);
+
+ ctrace = (complex *)malloc(nfreq*sizeof(complex));
+ if (ctrace == NULL) verr("memory allocation error for ctrace");
+
+ trace = (float *)malloc(optn*sizeof(float));
+ if (trace == NULL) verr("memory allocation error for rdata");
+
+ deltom = 2.*M_PI*df;
+ omin = 2.*M_PI*fmin;
+ omax = 2.*M_PI*fmax;
+ iomin = (long)MIN((omin/deltom), (nfreq));
+ iomax = MIN((long)(omax/deltom), (nfreq));
+ scl = 1.0/(float)optn;
+
+ for (ix = 0; ix < nrec; ix++) {
+ for (it=0;it<nsam;it++) trace[it]=data[ix*nsam+it];
+ for (it=nsam;it<optn;it++) trace[it]=0.0;
+ /* Forward time-frequency FFT */
+ rc1fft(&trace[0], &ctrace[0], optn, -1);
+
+ for (iom = 0; iom < iomin; iom++) {
+ ctrace[iom].r = 0.0;
+ ctrace[iom].i = 0.0;
+ }
+ for (iom = iomax; iom < nfreq; iom++) {
+ ctrace[iom].r = 0.0;
+ ctrace[iom].i = 0.0;
+ }
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = deltom*iom;
+ tom = om*shift;
+ ctmp = ctrace[iom];
+ ctrace[iom].r = ctmp.r*cos(-tom) - ctmp.i*sin(-tom);
+ ctrace[iom].i = ctmp.i*cos(-tom) + ctmp.r*sin(-tom);
+ }
+ /* Inverse frequency-time FFT and scale result */
+ cr1fft(ctrace, trace, optn, 1);
+ for (it=0;it<nsam;it++) data[ix*nsam+it]=trace[it]*scl;
+ }
+
+
+ free(ctrace);
+ free(trace);
+
+ return;
+}
\ No newline at end of file
diff --git a/marchenko3D/makeWindow3D.c b/marchenko3D/makeWindow3D.c
index c2f6f64..80152fc 100644
--- a/marchenko3D/makeWindow3D.c
+++ b/marchenko3D/makeWindow3D.c
@@ -22,6 +22,7 @@ typedef struct _complexStruct { /* complex number */
void findShotInMute(float *xrcvMute, float xrcvShot, long nxs, long *imute);
long readData3D(FILE *fp, float *data, segy *hdrs, long n1);
+void timeDiff(float *data, long nsam, long nrec, float dt, float fmin, float fmax, long opt);
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 makeWindow3D(char *file_ray, char *file_amp, char *file_wav, float dt, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc,
@@ -31,9 +32,13 @@ void makeWindow3D(char *file_ray, char *file_amp, char *file_wav, float dt, floa
segy hdr, *hdrs_mute, *hdrs_amp;
size_t nread;
long ig, is, i, iy, ix, j, l, nfreq, ntwav;
- float *wavelet, *wavtmp, scl, *timeval, dw, *amp;
+ float *wavelet, *wavtmp, scl, *timeval, dw, *amp, fmin, fmax;
complex *cmute, *cwav;
+
+ if (!getparfloat("fmin", &fmin)) fmin = 0.0;
+ if (!getparfloat("fmax", &fmax)) fmax = 70.0;
+
/*Define parameters*/
nfreq = ntfft/2+1;
wavelet = (float *)calloc(ntfft,sizeof(float));
@@ -58,9 +63,10 @@ void makeWindow3D(char *file_ray, char *file_amp, char *file_wav, float dt, floa
readData3D(fp, wavtmp, &hdr, ntwav);
//Fit the wavelet into the same time-axis as the Marchenko scheme
for (i=0; i<(ntfft/2); i++) {
- wavelet[i] = wavtmp[i];
- wavelet[ntfft-1-i] = wavtmp[ntwav-1-i];
+ wavelet[i] = -1.0*wavtmp[i];
+ wavelet[ntfft-1-i] = -1.0*wavtmp[ntwav-1-i];
}
+ timeDiff(wavelet, ntfft, 1, dt, fmin, fmax, 1);
rc1fft(wavelet,cwav,ntfft,-1);
free(wavtmp);
free(wavelet);
@@ -129,21 +135,27 @@ void makeWindow3D(char *file_ray, char *file_amp, char *file_wav, float dt, floa
if (file_wav!=NULL) { /*Apply the wavelet to create a first arrival*/
if (file_amp != NULL) {
for (ig=0; ig<nfreq; ig++) {
- cmute[ig].r = (dt/sqrtf((float)ntfft))*(cwav[ig].r*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)-cwav[ig].i*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0))/(amp[j*ny*nx+l*nx+i]);
- cmute[ig].i = (dt/sqrtf((float)ntfft))*(cwav[ig].i*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)+cwav[ig].r*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0))/(amp[j*ny*nx+l*nx+i]);
+ // cmute[ig].r = (dt/sqrtf((float)ntfft))*(cwav[ig].r*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)-cwav[ig].i*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0))/(amp[j*ny*nx+l*nx+i]);
+ // cmute[ig].i = (dt/sqrtf((float)ntfft))*(cwav[ig].i*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)+cwav[ig].r*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0))/(amp[j*ny*nx+l*nx+i]);
+ cmute[ig].r = (dt/sqrtf((float)ntfft))*(cwav[ig].r*cos(ig*dw*timeval[j*ny*nx+l*nx+i])-cwav[ig].i*sin(ig*dw*timeval[j*ny*nx+l*nx+i]))/((amp[j*ny*nx+l*nx+i]*amp[j*ny*nx+l*nx+i]));
+ cmute[ig].i = (dt/sqrtf((float)ntfft))*(cwav[ig].i*cos(ig*dw*timeval[j*ny*nx+l*nx+i])+cwav[ig].r*sin(ig*dw*timeval[j*ny*nx+l*nx+i]))/((amp[j*ny*nx+l*nx+i]*amp[j*ny*nx+l*nx+i]));
}
}
else { /*Use the raytime only to determine the mutewindow*/
for (ig=0; ig<nfreq; ig++) {
- cmute[ig].r = (dt/sqrtf((float)ntfft))*(cwav[ig].r*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)-cwav[ig].i*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0));
- cmute[ig].i = (dt/sqrtf((float)ntfft))*(cwav[ig].i*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)+cwav[ig].r*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0));
+ // cmute[ig].r = (dt/sqrtf((float)ntfft))*(cwav[ig].r*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)-cwav[ig].i*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0));
+ // cmute[ig].i = (dt/sqrtf((float)ntfft))*(cwav[ig].i*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0)+cwav[ig].r*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0));
+ cmute[ig].r = (dt/sqrtf((float)ntfft))*(cwav[ig].r*cos(ig*dw*timeval[j*ny*nx+l*nx+i])-cwav[ig].i*sin(ig*dw*timeval[j*ny*nx+l*nx+i]));
+ cmute[ig].i = (dt/sqrtf((float)ntfft))*(cwav[ig].i*cos(ig*dw*timeval[j*ny*nx+l*nx+i])+cwav[ig].r*sin(ig*dw*timeval[j*ny*nx+l*nx+i]));
}
}
}
else {
for (ig=0; ig<nfreq; ig++) {
- cmute[ig].r = (1.0/sqrtf((float)ntfft))*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0);
- cmute[ig].i = (1.0/sqrtf((float)ntfft))*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0);
+ // cmute[ig].r = (1.0/sqrtf((float)ntfft))*cos(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0);
+ // cmute[ig].i = (1.0/sqrtf((float)ntfft))*sin(ig*dw*timeval[j*ny*nx+l*nx+i]-M_PI/4.0);
+ cmute[ig].r = (1.0/sqrtf((float)ntfft))*cos(ig*dw*timeval[j*ny*nx+l*nx+i]);
+ cmute[ig].i = (1.0/sqrtf((float)ntfft))*sin(ig*dw*timeval[j*ny*nx+l*nx+i]);
}
}
cr1fft(cmute,&tinv[j*ny*nx*ntfft+l*nx*ntfft+i*ntfft],ntfft,1);
@@ -157,5 +169,4 @@ void makeWindow3D(char *file_ray, char *file_amp, char *file_wav, float dt, floa
}
return;
-}
-
+}
\ No newline at end of file
diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index 0cd893d..1ed3108 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -44,12 +44,15 @@ long readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
long *xnx, long Nfoc, long nx, long ny, long ntfft, long mode, long *maxval, float *tinv, long hw, long verbose);
long unique_elements(float *arr, long len);
+void timeShift(float *data, long nsam, long nrec, float dt, float shift, float fmin, float fmax);
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 *ixpos, long *iypos, long npos, long shift, long *tsynW);
+void applyMute3D_tshift( 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 iter, 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,
@@ -181,7 +184,7 @@ 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_amp;
- long i, j, l, k, ret, nshots, nxshot, nyshot, Nfoc, nt, nx, ny, nts, nxs, nys, ngath;
+ long i, j, l, k, ret, nshots, nxshot, nyshot, Nfoc, nt, nx, ny, nts, nxs, nys, ngath, *itmin;
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, zfp;
@@ -192,8 +195,8 @@ int main (int argc, char **argv)
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, *G_d, dt, dx, dy, dxs, dys, scl, mem;
- float *f1plus, *f1min, *iRN, *Ni, *Nig, *trace, *Gmin, *Gplus, *HomG;
- float scale, *tmpdata;
+ float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus, *HomG;
+ float scale, *tmpdata, tplmax;
float *ixmask, *iymask, *ampscl, *Gd, *Image, dzim;
float grad2rad, px, py, src_anglex, src_angley, src_velox, src_veloy, *mutetest;
complex *Refl, *Fop;
@@ -334,10 +337,10 @@ int main (int argc, char **argv)
f1min = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
iRN = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
Ni = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
- Nig = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
G_d = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
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
+ itmin = (long *)malloc(Nfoc*sizeof(long));
trace = (float *)malloc(ntfft*sizeof(float));
tapersx = (float *)malloc(nxs*sizeof(float));
tapersy = (float *)malloc(nys*sizeof(float));
@@ -389,43 +392,19 @@ int main (int argc, char **argv)
else dzim = 1.0;
/* compute time shift for tilted plane waves */
- if (plane_wave != 0) {
- grad2rad = 17.453292e-3;
- 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++) {
- tsynW[i] = 0;
+ if (plane_wave==1) {
+ for (j = 0; j < Nfoc; j++) {
+ itmin[j] = nt;
+ for (i=0; i<nys*nxs; i++) itmin[j] = MIN (itmin[j], muteW[j*nxs*nys+i]);
+ for (i=0; i<nys*nxs; i++) tsynW[j*nxs*nys+i] = muteW[j*nxs*nys+i]-itmin[j];
+ }
+ }
+ else { /* just fill with zero's */
+ for (j = 0; j < Nfoc; j++) {
+ itmin[j]=0;
+ for (i=0; i<nxs*nys; i++) {
+ tsynW[j*nxs*nys+i] = 0;
+ }
}
}
@@ -436,13 +415,13 @@ int main (int argc, char **argv)
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;
+ 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 = nys-ntapy; j < nys; j++)
- tapersy[j] =(cos(PI*(j-(nys-ntapy))/ntapy)+1)/2.0;
+ tapersy[j] = (cos(PI*(j-(nys-ntapy))/ntapy)+1)/2.0;
}
else {
for (j = 0; j < nxs; j++) tapersx[j] = 1.0;
@@ -772,12 +751,6 @@ int main (int argc, char **argv)
Ni[l*nys*nxs*nts+i*nts+j] = -iRN[l*nys*nxs*nts+iy*nxs*nts+ix*nts+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 */
- Nig[l*nys*nxs*nts+i*nts+j] = -iRN[l*nys*nxs*nts+iy*nxs*nts+ix*nts+j];
- for (j = 1; j < nts; j++) {
- Nig[l*nys*nxs*nts+i*nts+j] = -iRN[l*nys*nxs*nts+iy*nxs*nts+ix*nts+j];
- }
- }
}
if (iter==0) energyN0[l] = energyNi;
if (verbose >=2) vmess(" - iSyn %li: Ni at iteration %li has energy %e; relative to N0 %e",
@@ -785,37 +758,24 @@ 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, 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 (plane_wave==1) {
+ applyMute3D_tshift(Ni, muteW, smooth, above, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, iter, tsynW);
+ }
+ else {
+ applyMute3D(Ni, muteW, smooth, above, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
+ }
if (iter % 2 == 0) { /* even iterations update: => f_1^-(t) */
- if (plane_wave==0) { /* follow the standard focal point scheme */
- for (l = 0; l < Nfoc; l++) {
- for (i = 0; i < npos; i++) {
- j = 0;
- f1min[l*nys*nxs*nts+i*nts+j] -= Ni[l*nys*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- f1min[l*nys*nxs*nts+i*nts+j] -= Ni[l*nys*nxs*nts+i*nts+nts-j];
- }
- }
- }
- 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++) {
- for (i = 0; i < npos; i++) {
- j = 0;
- f1min[l*nys*nxs*nts+i*nts+j] -= Nig[l*nys*nxs*nts+i*nts+j];
- Ni[l*nys*nxs*nts+i*nts+j] = Nig[l*nys*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- f1min[l*nys*nxs*nts+i*nts+j] -= Nig[l*nys*nxs*nts+i*nts+j];
- Ni[l*nys*nxs*nts+i*nts+j] = Nig[l*nys*nxs*nts+i*nts+nts-j];
- }
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ f1min[l*nys*nxs*nts+i*nts+j] -= Ni[l*nys*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ f1min[l*nys*nxs*nts+i*nts+j] -= Ni[l*nys*nxs*nts+i*nts+nts-j];
}
}
- if (above==-2) applyMute3D(f1min, muteW, smooth, 0, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
}
+ if (above==-2) applyMute3D(f1min, muteW, smooth, 0, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
}
else {/* odd iterations update: => f_1^+(t) */
for (l = 0; l < Nfoc; l++) {
@@ -848,7 +808,6 @@ int main (int argc, char **argv)
} /* end of iterations */
free(Ni);
- free(Nig);
/* compute full Green's function G = int R * f2(t) + f2(-t) */
/* use f2 as operator on R in frequency domain */
@@ -883,7 +842,7 @@ int main (int argc, char **argv)
}
}
}
- if (!plane_wave) applyMute3D(green, muteW, smooth, 4, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
+ if (plane_wave!=1) 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) {
@@ -910,7 +869,24 @@ int main (int argc, char **argv)
}
}
/* Apply mute with window for Gmin */
- if (!plane_wave) applyMute3D(Gmin, muteW, smooth, 4, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
+ if (plane_wave==1) {
+ applyMute3D_tshift(Gmin, muteW, smooth, 4, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, 0, tsynW);
+ /* for plane wave with angle shift itmin downward */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ memcpy(&trace[0],&Gmin[l*nys*nxs*nts+i*nts],nts*sizeof(float));
+ for (j = 0; j < itmin[l]; j++) {
+ Gmin[l*nys*nxs*nts+i*nts+j] = 0.0;
+ }
+ for (j = 0; j < nts-itmin[l]; j++) {
+ Gmin[l*nys*nxs*nts+i*nts+j+itmin[l]] = trace[j];
+ }
+ }
+ }
+ }
+ else {
+ applyMute3D(Gmin, muteW, smooth, 4, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
+ }
} /* end if Gmin */
/* compute downgoing Green's function G^+,+ */
@@ -938,7 +914,12 @@ int main (int argc, char **argv)
}
}
/* Apply mute with window for Gplus */
- if (!plane_wave) applyMute3D(Gplus, muteW, smooth, 4, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
+ if (plane_wave) {
+ applyMute3D_tshift(Gplus, muteW, smooth, 4, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, iter, tsynW);
+ }
+ else {
+ applyMute3D(Gplus, muteW, smooth, 4, Nfoc, nxs, nys, nts, ixpos, iypos, npos, shift, tsynW);
+ }
} /* end if Gplus */
/* Estimate the amplitude of the Marchenko Redatuming */
@@ -1339,7 +1320,6 @@ int main (int argc, char **argv)
/*================ free memory ================*/
free(hdrs_out);
- free(tapersy);
exit(0);
}
diff --git a/raytime3d/getParameters3d.c b/raytime3d/getParameters3d.c
index 5c5952e..2f14086 100644
--- a/raytime3d/getParameters3d.c
+++ b/raytime3d/getParameters3d.c
@@ -351,7 +351,7 @@ long getParameters3d(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPa
vmess("************* receiver info ***************");
vmess("*******************************************");
vmess("ntrcv = %li nrcv = %li ", rec->nt, rec->n);
- vmess("nxrcv = %li nyrcv = %li nzrcv = %li ", rec->nx, rec->ny, rec->nz);
+ vmess("nxrcv = %li nyrcv = %li nzrcv = %li nrcv = %li", rec->nx, rec->ny, rec->nz, rec->n);
vmess("dzrcv = %f dxrcv = %f dyrcv = %f ", dzrcv, dxrcv, dyrcv);
vmess("Receiver array at coordinates: ");
vmess("zmin = %f zmax = %f ", rec->zr[0]+sub_z0, rec->zr[rec->n-1]+sub_z0);
diff --git a/utils/HomG.c b/utils/HomG.c
index 7e2d6f0..c6da0bd 100755
--- a/utils/HomG.c
+++ b/utils/HomG.c
@@ -50,7 +50,7 @@ long readSnapData3D(char *filename, float *data, segy *hdrs, long nsnaps, long n
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);
-
+long dignum(long number);
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);
@@ -58,6 +58,10 @@ void getVirReczfp(char *filename, long *nxs, long *nys, long *nxr, long *nyr, lo
void readzfpdata(char *filename, float *data, long size);
void getxyzzfp(char *filename, long *sx, long *sy, long *sz, long iz, long nz);
+void depthDiff3D(float *data, long nt, long nx, long ny, float dt, float dx, float dy, float fmin, float fmax, float c, int opt);
+void pad3d_data(float *data, long nt, long nx, long ny, long ntout, long nxout, long nyout, float *datout);
+void scl_data3D(float *data, long nt, long nx, long ny, float scl, float *datout, long ntout, long nxout);
+
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);
void convol(float *data1, float *data2, float *con, long nrec, long nsam, float dt, long shift);
@@ -104,7 +108,7 @@ int main (int argc, char **argv)
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;
+ long ntr, nxr, nyr, nsr, i, l, j, k, nxvr, nyvr, nzvr, count;
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;
@@ -138,6 +142,7 @@ int main (int argc, char **argv)
/*----------------------------------------------------------------------------*
* Split the filename so the number can be changed
*----------------------------------------------------------------------------*/
+ count = dignum(numb);
if (dnumb == 0) dnumb = 1;
sprintf(fins,"z%li",numb);
fp_in = fopen(fin, "r");
@@ -146,9 +151,9 @@ int main (int argc, char **argv)
}
fclose(fp_in);
ptr = strstr(fin,fins);
- pos1 = ptr - fin + 1;
- sprintf(fbegin,"%*.*s", pos1-1, pos1-1, fin);
- sprintf(fend,"%s", fin+pos1+1);
+ pos1 = ptr - fin;
+ sprintf(fbegin,"%*.*s", pos1, pos1, fin);
+ sprintf(fend,"%s", fin+pos1+count+1);
/*----------------------------------------------------------------------------*
* Determine the amount of files to be read
@@ -184,7 +189,7 @@ int main (int argc, char **argv)
/*----------------------------------------------------------------------------*
* Determine the other sizes of the files
*----------------------------------------------------------------------------*/
- sprintf(fins,"z%li",0);
+ sprintf(fins,"z%li",numb);
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);
@@ -220,8 +225,9 @@ int main (int argc, char **argv)
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);
+ vmess("Sampling distance is : x=%.3f y=%.3f t=%.3f",dx,dy,dt);
+ vmess("Scaling of the transforms : %.3f",scl);
+ vmess("Transform operators : fmin=%.1f fmax=%.1f cp=%.1f",fmin,fmax,cp);
}
if (ntr!=ntvs) verr("number of t-samples between virtual source (%li) and virtual receivers (%li) is not equal",ntvs,ntr);
@@ -248,16 +254,14 @@ int main (int argc, char **argv)
}
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);
- }
+ if (nyvs>1) depthDiff3D(&shotdata[0], ntvs, nxvs, nyvs, dt, dx, dy, fmin, fmax, cp, 1);
+ else 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);
- }
+ if (nyvs>1) depthDiff3D(&shotdata[ntvs*nxvs*nyvs], ntvs, nxvs, nyvs, dt, dx, dy, fmin, fmax, cp, 1);
+ else depthDiff(&shotdata[ntvs*nxvs*nyvs], ntvs, nxvs, dt, dx, fmin, fmax, cp, 1);
zsrc = labs(hdr_shot[0].sdepth);
}
else if (scheme==3) {
@@ -275,9 +279,8 @@ int main (int argc, char **argv)
}
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);
- }
+ if (nyvs>1) depthDiff3D(&shotdata_jkz[l*nyvs*nxvs*ntvs], ntvs, nxvs, nyvs, dt, dx, dy, fmin, fmax, cp, 1);
+ else depthDiff(&shotdata_jkz[l*nyvs*nxvs*ntvs], ntvs, nxvs, dt, dx, fmin, fmax, cp, 1);
}
}
else if (scheme==6) {
@@ -291,31 +294,30 @@ int main (int argc, char **argv)
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);
- }
+ if (nyvs>1) depthDiff3D(&shotdata[0*nyvs*nxvs*ntvs], ntvs, nxvs, nyvs, dt, dx, dy, fmin, fmax, cp, 1);
+ else depthDiff(&shotdata[0*nyvs*nxvs*ntvs], ntvs, nxvs, dt, dx, fmin, fmax, cp, 1);
+ conjugate(&shotdata[0*nyvs*nxvs*ntvs], ntvs, nxvs*nyvs, dt);
+ if (nyvs>1) depthDiff3D(&shotdata[1*nyvs*nxvs*ntvs], ntvs, nxvs, nyvs, dt, dx, dy, fmin, fmax, cp, 1);
+ else depthDiff(&shotdata[1*nyvs*nxvs*ntvs], ntvs, nxvs, dt, dx, fmin, fmax, cp, 1);
+ conjugate(&shotdata[1*nyvs*nxvs*ntvs], ntvs, nxvs*nyvs, 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);
- }
+ if (nyvs>1) depthDiff3D(&shotdata[0*nyvs*nxvs*ntvs], ntvs, nxvs, nyvs, dt, dx, dy, fmin, fmax, cp, 1);
+ else depthDiff(&shotdata[0*nyvs*nxvs*ntvs], ntvs, nxvs, dt, dx, fmin, fmax, cp, 1);
+ if (nyvs>1) depthDiff3D(&shotdata[1*nyvs*nxvs*ntvs], ntvs, nxvs, nyvs, dt, dx, dy, fmin, fmax, cp, 1);
+ else depthDiff(&shotdata[1*nyvs*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);
+ if (nyvs>1) depthDiff3D(&shotdata[0], ntvs, nxvs, nyvs, dt, dx, dy, fmin, fmax, cp, 1);
+ else depthDiff(&shotdata[0], ntvs, nxvs, dt, dx, fmin, fmax, cp, 1);
+ for (l = 0; l < nyvs*nxvs*ntvs; l++) {
+ shotdata_jkz[l] = shotdata[l];
}
+ conjugate(&shotdata_jkz[0], ntvs, nxvs*nyvs, dt);
}
else {
if (verbose) vmess("Marchenko Homogeneous Green's function retrieval with G source");
@@ -327,12 +329,12 @@ int main (int argc, char **argv)
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));
+ conv = (float *)calloc(nyr*nxr*ntr,sizeof(float));
if (scheme==5) {
- tmp1 = (float *)calloc(nxr*ntr,sizeof(float));
- tmp2 = (float *)calloc(nxr*ntr,sizeof(float));
+ tmp1 = (float *)calloc(nyr*nxr*ntr,sizeof(float));
+ tmp2 = (float *)calloc(nyr*nxr*ntr,sizeof(float));
}
- if (scheme==6 || scheme==8 || scheme==9 || scheme==10) tmp1 = (float *)calloc(nxr*ntr,sizeof(float));
+ if (scheme==6 || scheme==8 || scheme==9 || scheme==10) tmp1 = (float *)calloc(nyr*nxr*ntr,sizeof(float));
sprintf(fins,"z%li",ir*dnumb+numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
@@ -358,157 +360,142 @@ int main (int argc, char **argv)
}
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 (nyr > 1) depthDiff3D(&rcvdata[l*nyr*nxr*ntr], ntr, nxr, nyr, dt, dx, dy, fmin, fmax, cp, 1);
+ else depthDiff(&rcvdata[l*nxr*ntr], ntr, nxr, dt, dx, fmin, fmax, cp, 1);
+ convol(&shotdata[0], &rcvdata[l*nyr*nxr*ntr], conv, nyr*nxr, ntr, dt, 0);
+ timeDiff(conv, ntr, nyr*nxr, dt, fmin, fmax, -3);
+ for (i=0; i<nyr*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]/rho;
+ Ghom[j*nxvr*nyvr*nzvr+l*nzvr+ir] += 2.0*scl*conv[i*ntr+(j+ntr/2)]/rho;
}
}
}
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;
- }
+ convol(&shotdata[0], &rcvdata[l*nyr*nxr*ntr], conv, nxr*nyr, ntr, dt, 0);
+ timeDiff(conv, ntr, nyr*nxr, dt, fmin, fmax, -3);
+ for (i=0; i<nyr*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]/rho;
+ Ghom[j*nxvr*nyvr*nzvr+l*nzvr+ir] += 2.0*scl*conv[i*ntr+(j+ntr/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;
- }
+ if (zsrc > zrcv) {
+ if (verbose > 1) vmess("Homogeneous Green's function at %li uses G source (zsrc=%li)",zrcv,zsrc);
+ if (nyr>1) depthDiff3D(&rcvdata[l*nyr*nxr*ntr], ntr, nxr, nyr, dt, dx, dy, fmin, fmax, cp, 1);
+ else depthDiff(&rcvdata[l*nyr*nxr*ntr], ntr, nxr, dt, dx, fmin, fmax, cp, 1);
+ convol(&shotdata[0], &rcvdata[l*nyr*nxr*ntr], conv, nyr*nxr, ntr, dt, 0);
+ }
+ else {
+ if (verbose > 1) vmess("Homogeneous Green's function at %li uses f_2 source (zsrc=%li)",zrcv,zsrc);
+ convol(&shotdata[ntvs*nxvs*nyvs], &rcvdata[l*nyr*nxr*ntr], conv, nyr*nxr, ntr, dt, 0);
+ }
+ timeDiff(conv, ntr, nyr*nxr, dt, fmin, fmax, -1);
+ for (i=0; i<nyr*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]/rho;
+ Ghom[j*nxvr*nyvr*nzvr+l*nzvr+ir] += 2.0*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;
- }
+ if (nyr>1) depthDiff3D(&rcvdata[l*nyr*nxr*ntr], ntr, nxr, nyr, dt, dx, dy, fmin, fmax, cp, 1);
+ else depthDiff(&rcvdata[l*nyr*nxr*ntr], ntr, nxr, dt, dx, fmin, fmax, cp, 1);
+ convol2(&shotdata[0], &rcvdata[l*nyr*nxr*ntr], conv, nyr*nxr, ntr, dt, fmin, fmax, 1);
+ for (i=0; i<nyr*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]/rho;
+ Ghom[j*nxvr*nyvr*nzvr+l*nzvr+ir] += 2.0*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;
- }
+ if (nyr>1) depthDiff3D(&rcvdata[l*nyr*nxr*ntr], ntr, nxr, nyr, dt, dx, dy, fmin, fmax, cp, 1);
+ else depthDiff(&rcvdata[l*nyr*nxr*ntr], ntr, nxr, dt, dx, fmin, fmax, cp, 1);
+ convol(&shotdata[0], &rcvdata[l*nyr*nxr*ntr], conv, nyr*nxr, ntr, dt, 0);
+ timeDiff(conv, ntr, nyr*nxr, dt, fmin, fmax, -1);
+ for (i=0; i<nyr*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]/rho;
+ Ghom[j*nxvr*nyvr*nzvr+l*nzvr+ir] += 2.0*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];
- }
+ convol(&rcvdata[l*nyr*nxr*ntr], &shotdata_jkz[0], tmp2, nyr*nxr, ntr, dt, 0);
+ if (nyr>1) depthDiff3D(&rcvdata[l*nyr*nxr*ntr], ntr, nxr, nyr, dt, dx, dy, fmin, fmax, cp, 1);
+ else depthDiff(&rcvdata[l*nyr*nxr*ntr], ntr, nxr, dt, dx, fmin, fmax, cp, 1);
+ convol(&rcvdata[l*nyr*nxr*ntr], &shotdata[0], tmp1, nyr*nxr, ntr, dt, 0);
+ for (i = 0; i < nyr*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;
- }
+ }
+ timeDiff(conv, ntr, nyr*nxr, dt, fmin, fmax, -1);
+ for (i=0; i<nyr*nxr; i++) {
+ for (j=0; j<ntr/2; j++) {
+ Ghom[(j+ntr/2)*nxvr*nyvr*nzvr+l*nzvr+ir] += 1.0*scl*conv[i*ntr+j]/rho;
+ Ghom[j*nxvr*nyvr*nzvr+l*nzvr+ir] += 1.0*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;
- }
+ if (nyr>1) depthDiff3D(&rcvdata[l*nyr*nxr*ntr], ntr, nxr, nyr, dt, dx, dy, fmin, fmax, cp, 1);
+ else depthDiff(&rcvdata[l*nyr*nxr*ntr], ntr, nxr, dt, dx, fmin, fmax, cp, 1);
+ for (is=0; is<nshots; is++) {
+ convol(&shotdata[is*nyr*nxr*ntr], &rcvdata[l*nyr*nxr*ntr], conv, nyr*nxr, ntr, dt, 0);
+ timeDiff(conv, ntr, nyr*nxr, dt, fmin, fmax, -3);
+ for (i=0; i<nyr*nxr; i++) {
+ for (j=0; j<ntr/2; j++) {
+ Ghom[is*ntr*nxvr*nyvr*nzvr+(j+ntr/2)*nxvr*nyvr*nzvr+l*nzvr+ir] += 2.0*scl*conv[i*ntr+j]/rho;
+ Ghom[is*ntr*nxvr*nyvr*nzvr+j*nxvr*nyvr*nzvr+l*nzvr+ir] += 2.0*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;
- }
+ if (nyr>1) depthDiff3D(&rcvdata[l*nyr*nxr*ntr], ntr, nxr, nyr, dt, dx, dy, fmin, fmax, cp, 1);
+ else depthDiff(&rcvdata[l*nyr*nxr*ntr], ntr, nxr, dt, dx, fmin, fmax, cp, 1);
+ for (is=0; is<nshots; is++) {
+ convol(&shotdata[is*nyr*nxr*ntr], &rcvdata[l*nyr*nxr*ntr], conv, nyr*nxr, ntr, dt, 0);
+ timeDiff(conv, ntr, nyr*nxr, dt, fmin, fmax, -1);
+ for (i=0; i<nyr*nxr; i++) {
+ for (j=0; j<ntr/2; j++) {
+ Ghom[is*ntr*nxvr*nyvr*nzvr+(j+ntr/2)*nxvr*nyvr*nzvr+l*nzvr+ir] += 2.0*scl*conv[i*ntr+j]/rho;
+ Ghom[is*ntr*nxvr*nyvr*nzvr+j*nxvr*nyvr*nzvr+l*nzvr+ir] += 2.0*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;
- }
+ convol2(&shotdata[0*nyr*nxr*ntr], &rcvdata[l*nyr*nxr*ntr], conv, nyr*nxr, ntr, dt, fmin, fmax, 1);
+ convol2(&shotdata[1*nyr*nxr*ntr], &rcvdata[(l+1)*nyr*nxr*ntr], tmp1, nyr*nxr, ntr, dt, fmin, fmax, 1);
+ for (i=0; i<nyr*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;
- }
+ convol2(&shotdata[0*nyr*nxr*ntr], &rcvdata[(l+1)*nyr*nxr*ntr], conv, nyr*nxr, ntr, dt, fmin, fmax, 1);
+ convol2(&shotdata[1*nyr*nxr*ntr], &rcvdata[l*nyr*nxr*ntr], tmp1, nyr*nxr, ntr, dt, fmin, fmax, 1);
+ for (i=0; i<nyr*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;
- }
+ convol2(&shotdata[0], &rcvdata[(l+1)*nyr*nxr*ntr], conv, nyr*nxr, ntr, dt, fmin, fmax, 2);
+ convol2(&shotdata_jkz[0], &rcvdata[l*nyr*nxr*ntr], tmp1, nyr*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;
}
}
}
@@ -538,8 +525,8 @@ int main (int argc, char **argv)
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].scalco = -1000;
+ hdr_out[is*nxvr+ix].scalel = -1000;
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;
@@ -1129,7 +1116,7 @@ void getVirRec(char *filename, long *nxs, long *nys, long *nxr, long *nyr, long
gy = hdr.gy;
gx0 = gx;
gy0 = gy;
- itrace = 1;
+ itrace = 0;
ishot = 1;
tsize = hdr.ns*sizeof(float);
@@ -1354,4 +1341,217 @@ void getxyzzfp(char *filename, long *sx, long *sy, long *sz, long iz, long nz)
}
return;
+}
+
+long dignum(long number)
+{
+ long count = 0;
+
+ /* Calculate total digits */
+ count = (number == 0) ? 1 : (log10(number) + 1);
+
+ return count;
+}
+
+void depthDiff3D(float *data, long nt, long nx, long ny, float dt, float dx, float dy, float fmin, float fmax, float c, int opt)
+{
+ long optn, iom, iomin, iomax, nfreq, ix, iy, ikx, iky, nkx, nky, ikxmax, ikymax;
+ float omin, omax, deltom, df, dkx, dky, *rdata, kx, ky, scl;
+ float kx2, ky2, kz2, kp2, kp;
+ complex *cdata, *cdatascl, kz, kzinv;
+
+ optn = optncr(nt);
+ nfreq = optncr(nt)/2+1;
+ df = 1.0/(optn*dt);
+ nkx = optncc(nx);
+ nky = optncc(ny);
+ dkx = 2.0*PI/(nkx*dx);
+ dky = 2.0*PI/(nky*dy);
+
+ cdata = (complex *)calloc(nfreq*nkx*nky,sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nkx*nky*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes in 2 directions to reach FFT lengths */
+ pad3d_data(data, nt, nx, ny, optn, nkx, nky, rdata);
+
+ /* double forward FFT */
+ yxt2wkykx(&rdata[0], &cdata[0], optn, nkx, nky, optn, nkx, nky, 0, 0);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+
+ iomin = (long)MIN((omin/deltom), nfreq);
+ iomin = MAX(iomin, 0);
+ iomax = MIN((long)(omax/deltom), nfreq);
+
+ cdatascl = (complex *)calloc(nfreq*nkx*nky,sizeof(complex));
+ if (cdatascl == NULL) verr("memory allocation error for cdatascl");
+
+ if (opt > 0) {
+ for (iom = iomin ; iom <= iomax ; iom++) {
+ kp = (iom*deltom)/c;
+ kp2 = kp*kp;
+ ikxmax = nkx/2;
+ ikymax = nky/2;
+
+ for (iky = 0; iky < ikymax; iky++) {
+ ky = iky*dky;
+ ky2 = ky*ky;
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ if (kz2<0.0) continue;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].r = cdata[iom*nky*nkx+iky*nkx+ikx].r*kz.r-cdata[iom*nky*nkx+iky*nkx+ikx].i*kz.i;
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].i = cdata[iom*nky*nkx+iky*nkx+ikx].i*kz.r+cdata[iom*nky*nkx+iky*nkx+ikx].r*kz.i;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ if (kz2<0.0) continue;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].r = cdata[iom*nky*nkx+iky*nkx+ikx].r*kz.r-cdata[iom*nky*nkx+iky*nkx+ikx].i*kz.i;
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].i = cdata[iom*nky*nkx+iky*nkx+ikx].i*kz.r+cdata[iom*nky*nkx+iky*nkx+ikx].r*kz.i;
+ }
+ }
+ for (iky = nky-ikymax+1; iky < nky; iky++) {
+ ky = (iky-nky)*dky;
+ ky2 = ky*ky;
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ if (kz2<0.0) continue;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].r = cdata[iom*nky*nkx+iky*nkx+ikx].r*kz.r-cdata[iom*nky*nkx+iky*nkx+ikx].i*kz.i;
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].i = cdata[iom*nky*nkx+iky*nkx+ikx].i*kz.r+cdata[iom*nky*nkx+iky*nkx+ikx].r*kz.i;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ if (kz2<0.0) continue;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].r = cdata[iom*nky*nkx+iky*nkx+ikx].r*kz.r-cdata[iom*nky*nkx+iky*nkx+ikx].i*kz.i;
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].i = cdata[iom*nky*nkx+iky*nkx+ikx].i*kz.r+cdata[iom*nky*nkx+iky*nkx+ikx].r*kz.i;
+ }
+ }
+
+ }
+ }
+ else if (opt < 0) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ kp = iom*deltom/c;
+ kp2 = kp*kp;
+ ikxmax = nkx/2;
+ ikymax = nky/2;
+
+ for (iky = 0; iky < ikymax; iky++) {
+ ky = iky*dky;
+ ky2 = ky*ky;
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ if (kz2<0.0) continue;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].r = cdata[iom*nky*nkx+iky*nkx+ikx].r*kzinv.r-cdata[iom*nky*nkx+iky*nkx+ikx].i*kzinv.i;
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].i = cdata[iom*nky*nkx+iky*nkx+ikx].i*kzinv.r+cdata[iom*nky*nkx+iky*nkx+ikx].r*kzinv.i;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ if (kz2<0.0) continue;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].r = cdata[iom*nky*nkx+iky*nkx+ikx].r*kzinv.r-cdata[iom*nky*nkx+iky*nkx+ikx].i*kzinv.i;
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].i = cdata[iom*nky*nkx+iky*nkx+ikx].i*kzinv.r+cdata[iom*nky*nkx+iky*nkx+ikx].r*kzinv.i;
+ }
+ }
+ for (iky = nky-ikymax+1; iky < nky; iky++) {
+ ky = (iky-nky)*dky;
+ ky2 = ky*ky;
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ if (kz2<0.0) continue;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].r = cdata[iom*nky*nkx+iky*nkx+ikx].r*kzinv.r-cdata[iom*nky*nkx+iky*nkx+ikx].i*kzinv.i;
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].i = cdata[iom*nky*nkx+iky*nkx+ikx].i*kzinv.r+cdata[iom*nky*nkx+iky*nkx+ikx].r*kzinv.i;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ if (kz2<0.0) continue;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].r = cdata[iom*nky*nkx+iky*nkx+ikx].r*kzinv.r-cdata[iom*nky*nkx+iky*nkx+ikx].i*kzinv.i;
+ cdatascl[iom*nky*nkx+iky*nkx+ikx].i = cdata[iom*nky*nkx+iky*nkx+ikx].i*kzinv.r+cdata[iom*nky*nkx+iky*nkx+ikx].r*kzinv.i;
+ }
+ }
+
+ }
+ }
+ free(cdata);
+
+ /* inverse double FFT */
+ wkykx2yxt(&cdatascl[0], &rdata[0], optn, nkx, nky, optn, nkx, nky, 0, 0);
+ /* select original samples and traces */
+ scl = 1.0;
+ scl_data3D(rdata, nt, nx, ny, scl, data, optn, nkx);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
+
+void pad3d_data(float *data, long nt, long nx, long ny, long ntout, long nxout, long nyout, float *datout)
+{
+ int it,ix,iy;
+ for (iy=0;iy<ny;iy++) {
+ for (ix=0;ix<nx;ix++) {
+ for (it=0;it<nt;it++)
+ datout[iy*nxout*ntout+ix*ntout+it]=data[iy*nx*nt+ix*nt+it];
+ for (it=nt;it<ntout;it++)
+ datout[iy*nxout*ntout+ix*ntout+it]=0.0;
+ }
+ for (ix=nx;ix<nxout;ix++) {
+ for (it=0;it<ntout;it++)
+ datout[iy*nxout*ntout+ix*ntout+it]=0.0;
+ }
+ }
+ for (iy=ny;iy<nyout;iy++) {
+ for (ix=0;ix<nxout;ix++) {
+ for (it=0;it<ntout;it++)
+ datout[iy*nxout*ntout+ix*ntout+it]=0.0;
+ }
+ }
+}
+
+void scl_data3D(float *data, long nt, long nx, long ny, float scl, float *datout, long ntout, long nxout)
+{
+ int it,ix,iy;
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (it = 0 ; it < nt ; it++) {
+ datout[iy*nx*nt+ix*nt+it] = scl*data[iy*nxout*ntout+ix*ntout+it];
+ }
+ }
+ }
}
\ No newline at end of file
diff --git a/utils/HomG_backup13may2020.c b/utils/HomG_backup13may2020.c
new file mode 100644
index 0000000..54db1d8
--- /dev/null
+++ b/utils/HomG_backup13may2020.c
@@ -0,0 +1,1610 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#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))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} 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);
+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);
+long dignum(long number);
+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 depthDiff3D(float *data, long nt, long nx, long ny, float dt, float dx, float dy, float fmin, float fmax, float c, int opt);
+void pad3d_data(float *data, long nt, long nx, long ny, long ntout, long nxout, long nyout, float *datout);
+void scl_data3D(float *data, long nt, long nx, long ny, float scl, float *datout, long ntout, long nxout);
+
+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);
+void convol(float *data1, float *data2, float *con, long nrec, long nsam, float dt, long shift);
+void corr(float *data1, float *data2, float *cov, long nrec, long nsam, float dt, long shift);
+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[] = {
+" ",
+" HomG - Calculate a Homogeneous Green's function ",
+" ",
+" authors : Joeri Brackenhoff (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" file_vr= ................. First file of the array of virtual receivers",
+" file_vs= ................. File containing the virtual source",
+" ",
+" Optional parameters: ",
+" ",
+" file_out= ................ Filename of the output",
+" numb= .................... integer number of first snapshot file",
+" dnumb= ................... integer number of increment in snapshot files",
+" zmax= .................... Integer number of maximum depth level",
+" inx= ..................... Number of sources per depth level",
+" zrcv= .................... z-coordinate of first receiver location",
+" xrcv= .................... x-coordinate of first receiver location",
+" 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",
+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 *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, count;
+ 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);
+
+ /*----------------------------------------------------------------------------*
+ * Get the parameters passed to the function
+ *----------------------------------------------------------------------------*/
+ 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 (!getparfloat("rho", &rho)) rho=1000.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;
+ if (!getparlong("numb", &numb)) numb=0;
+ if (!getparlong("dnumb", &dnumb)) dnumb=1;
+ if (!getparlong("scheme", &scheme)) scheme = 0;
+ if (!getparlong("ntmax", &ntmax)) ntmax = 0;
+ if (!getparlong("verbose", &verbose)) verbose = 0;
+ 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
+ *----------------------------------------------------------------------------*/
+ count = dignum(numb);
+ if (dnumb == 0) dnumb = 1;
+ sprintf(fins,"z%li",numb);
+ fp_in = fopen(fin, "r");
+ if (fp_in == NULL) {
+ verr("error on opening basefile=%s", fin);
+ }
+ fclose(fp_in);
+ ptr = strstr(fin,fins);
+ pos1 = ptr - fin;
+ sprintf(fbegin,"%*.*s", pos1, pos1, fin);
+ sprintf(fend,"%s", fin+pos1+count+1);
+
+ /*----------------------------------------------------------------------------*
+ * Determine the amount of files to be read
+ *----------------------------------------------------------------------------*/
+ file_det = 1;
+ nzvr=0;
+ while (file_det) {
+ 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 (nzvr == 0) {
+ verr("error on opening basefile=%s", fin);
+ }
+ else if (nzvr == 1) {
+ vmess("1 file detected");
+ file_det = 0;
+ break;
+ }
+ else {
+ vmess("%li files detected",nzvr);
+ file_det = 0;
+ break;
+ }
+ }
+ fclose(fp_in);
+ nzvr++;
+ }
+
+ if (zmax < 1) zmax=0;
+ if (zmax < nzvr && zmax > 0) nzvr=zmax;
+
+ /*----------------------------------------------------------------------------*
+ * Determine the other sizes of the files
+ *----------------------------------------------------------------------------*/
+ sprintf(fins,"z%li",numb);
+ 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);
+ }
+
+ /*----------------------------------------------------------------------------*
+ * Get the file info for the source position and read in the data
+ *----------------------------------------------------------------------------*/
+
+ nshots = 0;
+ 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);
+
+ scl = dx*dy*dt;
+
+ 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");
+ }
+ fclose(fp_shot);
+
+ 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);
+ }
+
+ 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);
+
+ 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);
+ }
+ }
+ }
+ 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==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");
+ }
+
+#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");
+ if (fp_in == NULL) {
+ verr("Danger Will Robinson");
+ }
+ fclose(fp_in);
+
+ 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;
+ }
+ }
+ }
+ }
+ 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;
+ }
+ }
+ }
+ }
+ 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;
+ }
+ }
+ }
+ }
+ }
+ 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 (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 < 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;
+ hdr_out[is*nxvr+ix].scalel = -1000;
+ 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*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 = 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], (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], (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)
+{
+ long optn, iom, iomin, iomax, nfreq, ix, sign;
+ float omin, omax, deltom, om, df, *rdata, scl;
+ complex *cdata, *cdatascl;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+ df = 1.0/(optn*dt);
+
+ cdata = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data,nsam,nrec,optn,rdata);
+
+ /* Forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata[0], &cdata[0], optn, nrec, optn, nfreq, sign);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+ iomin = (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 = deltom*iom;
+ cdatascl[ix*nfreq+iom].r = -om*cdata[ix*nfreq+iom].i;
+ cdatascl[ix*nfreq+iom].i = om*cdata[ix*nfreq+iom].r;
+ }
+ }
+ else if (opt == -1) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = 1.0/(deltom*iom);
+ cdatascl[ix*nfreq+iom].r = om*cdata[ix*nfreq+iom].i;
+ cdatascl[ix*nfreq+iom].i = -om*cdata[ix*nfreq+iom].r;
+ }
+ }
+ 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);
+
+ /* Inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&cdatascl[0], &rdata[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
+
+void depthDiff(float *data, 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;
+ float omin, omax, deltom, df, dkx, *rdata, kx, scl;
+ float kx2, kz2, kp2, kp;
+ complex *cdata, *cdatascl, kz, kzinv;
+
+ optn = optncr(nsam);
+ nfreq = optncr(nsam)/2+1;
+ df = 1.0/(optn*dt);
+ nkx = optncc(nrec);
+ dkx = 2.0*PI/(nkx*dx);
+ cdata = (complex *)malloc(nfreq*nkx*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nkx*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes in 2 directions to reach FFT lengths */
+ pad2d_data(data,nsam,nrec,optn,nkx,rdata);
+
+ /* double forward FFT */
+ xt2wkx(&rdata[0], &cdata[0], optn, nkx, optn, nkx, 0);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+
+ iomin = (long)MIN((omin/deltom), nfreq);
+ iomin = MAX(iomin, 0);
+ iomax = MIN((long)(omax/deltom), nfreq);
+
+ cdatascl = (complex *)malloc(nfreq*nkx*sizeof(complex));
+ if (cdatascl == NULL) verr("memory allocation error for cdatascl");
+
+ for (iom = 0; iom < iomin; iom++) {
+ for (ix = 0; ix < nkx; ix++) {
+ cdatascl[iom*nkx+ix].r = 0.0;
+ cdatascl[iom*nkx+ix].i = 0.0;
+ }
+ }
+ for (iom = iomax; iom < nfreq; iom++) {
+ for (ix = 0; ix < nkx; ix++) {
+ cdatascl[iom*nkx+ix].r = 0.0;
+ cdatascl[iom*nkx+ix].i = 0.0;
+ }
+ }
+ if (opt > 0) {
+ for (iom = iomin ; iom <= iomax ; iom++) {
+ kp = (iom*deltom)/c;
+ kp2 = kp*kp;
+
+ ikxmax = MIN((long)(kp/dkx), nkx/2);
+
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kz.r-cdata[iom*nkx+ikx].i*kz.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kz.r+cdata[iom*nkx+ikx].r*kz.i;
+
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nkx+ikx].r = 0.0;
+ cdatascl[iom*nkx+ikx].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kz.r-cdata[iom*nkx+ikx].i*kz.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kz.r+cdata[iom*nkx+ikx].r*kz.i;
+ }
+ }
+ }
+ else if (opt < 0) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ kp = iom*deltom/c;
+ kp2 = kp*kp;
+ ikxmax = MIN((long)(kp/dkx), nkx/2);
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kzinv.r-cdata[iom*nkx+ikx].i*kzinv.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kzinv.r+cdata[iom*nkx+ikx].r*kzinv.i;
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nkx+ikx].r = 0.0;
+ cdatascl[iom*nkx+ikx].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kzinv.r-cdata[iom*nkx+ikx].i*kzinv.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kzinv.r+cdata[iom*nkx+ikx].r*kzinv.i;
+ }
+ }
+ }
+ free(cdata);
+
+ /* inverse double FFT */
+ wkx2xt(&cdatascl[0], &rdata[0], optn, nkx, nkx, optn, 0);
+ /* select original samples and traces */
+ scl = 1.0;
+ scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
+
+void pad2d_data(float *data, long nsam, long nrec, long nsamout, long nrecout, float *datout)
+{
+ long it,ix;
+ for (ix=0;ix<nrec;ix++) {
+ for (it=0;it<nsam;it++)
+ datout[ix*nsam+it]=data[ix*nsam+it];
+ for (it=nsam;it<nsamout;it++)
+ datout[ix*nsam+it]=0.0;
+ }
+ for (ix=nrec;ix<nrecout;ix++) {
+ for (it=0;it<nsamout;it++)
+ datout[ix*nsam+it]=0.0;
+ }
+}
+
+void conjugate(float *data, long nsam, long nrec, float dt)
+{
+ long optn, nfreq, j, ix, it, sign, ntdiff;
+ float *rdata, scl;
+ complex *cdata;
+
+ optn = optncr(nsam);
+ ntdiff = optn-nsam;
+ nfreq = optn/2+1;
+
+ cdata = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data,nsam,nrec,optn,rdata);
+
+ /* Forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata[0], &cdata[0], optn, nrec, optn, nfreq, sign);
+
+ /* take complex conjugate */
+ for(ix = 0; ix < nrec; ix++) {
+ for(j = 0; j < nfreq; j++) cdata[ix*nfreq+j].i = -cdata[ix*nfreq+j].i;
+ }
+
+ /* Inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&cdata[0], &rdata[0], optn, nrec, nfreq, optn, sign);
+ for (ix = 0; ix < nrec; ix++) {
+ for (it = 0 ; it < nsam ; it++)
+ data[ix*nsam+it] = scl*rdata[ix*optn+it+ntdiff];
+ }
+ //scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ 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 = 0;
+ 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;
+}
+
+long dignum(long number)
+{
+ long count = 0;
+
+ /* Calculate total digits */
+ count = (number == 0) ? 1 : (log10(number) + 1);
+
+ return count;
+}
+
+void depthDiff3D(float *data, long nt, long nx, long ny, float dt, float dx, float dy, float fmin, float fmax, float c, int opt)
+{
+ long optn, iom, iomin, iomax, nfreq, ix, iy, ikx, iky, nkx, nky, ikxmax, ikymax;
+ float omin, omax, deltom, df, dkx, dky, *rdata, kx, ky, scl;
+ float kx2, ky2, kz2, kp2, kp;
+ complex *cdata, *cdatascl, kz, kzinv;
+
+ optn = optncr(nt);
+ nfreq = optncr(nt)/2+1;
+ df = 1.0/(optn*dt);
+ nkx = optncc(nx);
+ nky = optncc(ny);
+ dkx = 2.0*PI/(nkx*dx);
+ dky = 2.0*PI/(nky*dy);
+ cdata = (complex *)malloc(nfreq*nkx*nky*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nkx*nky*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes in 2 directions to reach FFT lengths */
+ pad3d_data(data, nt, nx, ny, optn, nkx, nky, rdata);
+
+ /* double forward FFT */
+ yxt2wkykx(&rdata[0], &cdata[0], optn, nkx, nky, optn, nkx, nky, 0, 0);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+
+ iomin = (int)MIN((omin/deltom), nfreq);
+ iomin = MAX(iomin, 0);
+ iomax = MIN((int)(omax/deltom), nfreq);
+
+ cdatascl = (complex *)malloc(nfreq*nkx*nky*sizeof(complex));
+ if (cdatascl == NULL) verr("memory allocation error for cdatascl");
+
+ for (iom = 0; iom < iomin; iom++) {
+ for (iy = 0; iy < nky; iy++) {
+ for (ix = 0; ix < nkx; ix++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ }
+ }
+ for (iom = iomax; iom < nfreq; iom++) {
+ for (iy = 0; iy < nky; iy++) {
+ for (ix = 0; ix < nkx; ix++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ }
+ }
+ if (opt > 0) {
+ for (iom = iomin ; iom <= iomax ; iom++) {
+ kp = (iom*deltom)/c;
+ kp2 = kp*kp;
+ ikxmax = MIN((int)(kp/dkx), nkx/2);
+ ikymax = MIN((int)(kp/dky), nky/2);
+
+ for (iky = 0; iky < ikymax; iky++) {
+ ky = iky*dky;
+ ky2 = ky*ky;
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kz.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kz.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kz.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kz.i;
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kz.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kz.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kz.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kz.i;
+ }
+ }
+ for (iky = ikymax; iky <= nky-ikymax+1; iky++) {
+ for (ikx = 0; ikx <= nkx; ikx++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ }
+ for (iky = nky-ikymax+1; iky < nky; iky++) {
+ ky = (iky-nky)*dky;
+ ky2 = ky*ky;
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kz.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kz.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kz.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kz.i;
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kz.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kz.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kz.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kz.i;
+ }
+ }
+
+ }
+ }
+ else if (opt < 0) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ kp = iom*deltom/c;
+ kp2 = kp*kp;
+ ikxmax = MIN((int)(kp/dkx), nkx/2);
+ ikymax = MIN((int)(kp/dky), nky/2);
+
+ for (iky = 0; iky < ikymax; iky++) {
+ ky = iky*dky;
+ ky2 = ky*ky;
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.i;
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.i;
+ }
+ }
+ for (iky = ikymax; iky <= nky-ikymax+1; iky++) {
+ for (ikx = 0; ikx <= nkx; ikx++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ }
+ for (iky = nky-ikymax+1; iky < nky; iky++) {
+ ky = (iky-nky)*dky;
+ ky2 = ky*ky;
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.i;
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.i;
+ }
+ }
+
+ }
+ }
+ free(cdata);
+
+ /* inverse double FFT */
+ wkykx2yxt(&cdatascl[0], &rdata[0], optn, nkx, nky, optn, nkx, nky, 0, 0);
+ /* select original samples and traces */
+ scl = 1.0;
+ scl_data3D(rdata, optn, nx, ny, scl, data, nt, nx);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
+
+void pad3d_data(float *data, long nt, long nx, long ny, long ntout, long nxout, long nyout, float *datout)
+{
+ int it,ix,iy;
+ for (iy=0;iy<ny;iy++) {
+ for (ix=0;ix<nx;ix++) {
+ for (it=0;it<nt;it++)
+ datout[iy*nx*nt+ix*nt+it]=data[iy*nx*nt+ix*nt+it];
+ for (it=nt;it<ntout;it++)
+ datout[iy*nx*nt+ix*nt+it]=0.0;
+ }
+ for (ix=nx;ix<nxout;ix++) {
+ for (it=0;it<ntout;it++)
+ datout[iy*nx*nt+ix*nt+it]=0.0;
+ }
+ }
+ for (iy=ny;iy<nyout;iy++) {
+ for (ix=0;ix<nxout;ix++) {
+ for (it=0;it<ntout;it++)
+ datout[iy*nx*nt+ix*nt+it]=0.0;
+ }
+ }
+}
+
+void scl_data3D(float *data, long nt, long nx, long ny, float scl, float *datout, long ntout, long nxout)
+{
+ int it,ix,iy;
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (it = 0 ; it < ntout ; it++) {
+ datout[iy*nxout*ntout+ix*ntout+it] = scl*data[iy*nx*nt+ix*nt+it];
+ }
+ }
+ }
+}
\ No newline at end of file
diff --git a/utils/HomG_backup25mar2020.c b/utils/HomG_backup25mar2020.c
new file mode 100644
index 0000000..41eab13
--- /dev/null
+++ b/utils/HomG_backup25mar2020.c
@@ -0,0 +1,867 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#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))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+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);
+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 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);
+void convol(float *data1, float *data2, float *con, long nrec, long nsam, float dt, long shift);
+void corr(float *data1, float *data2, float *cov, long nrec, long nsam, float dt, long shift);
+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);
+
+char *sdoc[] = {
+" ",
+" HomG - Calculate a Homogeneous Green's function ",
+" ",
+" authors : Joeri Brackenhoff (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" file_in= ................. First file of the array of receivers",
+" file_shot= ............... File containing the shot location",
+" ",
+" Optional parameters: ",
+" ",
+" file_out= ................ Filename of the output",
+" numb= .................... integer number of first snapshot file",
+" dnumb= ................... integer number of increment in snapshot files",
+" zmax= .................... Integer number of maximum depth level",
+" 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",
+" shift=0.0 ................ shift per shot",
+" scheme=0 ................. Scheme used for retrieval. 0=Marchenko,",
+" 1=Marchenko with multiple sources, 2=classical",
+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;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ /*----------------------------------------------------------------------------*
+ * 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 (!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("fmin", &fmin)) fmin=0.0;
+ if (!getparfloat("fmax", &fmax)) fmax=100.0;
+ if (!getparfloat("shift", &shift)) shift=0.0;
+ if (!getparlong("numb", &numb)) numb=0;
+ if (!getparlong("dnumb", &dnumb)) dnumb=1;
+ 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");
+
+ /*----------------------------------------------------------------------------*
+ * Split the filename so the number can be changed
+ *----------------------------------------------------------------------------*/
+ if (dnumb == 0) dnumb = 1;
+ sprintf(fins,"z%li",numb);
+ fp_in = fopen(fin, "r");
+ if (fp_in == NULL) {
+ verr("error on opening basefile=%s", fin);
+ }
+ fclose(fp_in);
+ ptr = strstr(fin,fins);
+ pos1 = ptr - fin + 1;
+ sprintf(fbegin,"%*.*s", pos1-1, pos1-1, fin);
+ sprintf(fend,"%s", fin+pos1+1);
+
+ /*----------------------------------------------------------------------------*
+ * Determine the amount of files to be read
+ *----------------------------------------------------------------------------*/
+ file_det = 1;
+ nzs=0;
+ while (file_det) {
+ sprintf(fins,"z%li",nzs*dnumb+numb);
+ sprintf(fin,"%s%s%s",fbegin,fins,fend);
+ fp_in = fopen(fin, "r");
+ if (fp_in == NULL) {
+ if (nzs == 0) {
+ verr("error on opening basefile=%s", fin);
+ }
+ else if (nzs == 1) {
+ vmess("1 file detected");
+ file_det = 0;
+ break;
+ }
+ else {
+ vmess("%li files detected",nzs);
+ file_det = 0;
+ break;
+ }
+ }
+ fclose(fp_in);
+ nzs++;
+ }
+
+ if (inx < 1) {
+ inx = 1;
+ }
+
+ if (zmax < 1) zmax=1;
+ if (zmax < nzs) nzs=zmax;
+
+ nxs = inx;
+ count=0;
+ npos = nxs*nzs;
+
+ if (verbose) vmess("nxs: %li, nzs: %li",nxs,nzs);
+
+ nshots = 0;
+ getFileInfo3D(fshot, &nt, &nx, &ny, &nshots, &dt, &dx, &dy, &t0, &x0, &y0, &sclsxgx, &ntraces);
+
+ if (dxrcv < 0) dxrcv=dx;
+ if (dzrcv < 0) dzrcv=dx;
+
+ // ngath zijn het aantal schoten
+ shotdata = (float *)malloc(nt*nx*nshots*sizeof(float));
+ hdr_shot = (segy *)calloc(nx*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);
+
+
+ hdr_out = (segy *)calloc(nxs,sizeof(segy));
+ Ghom = (float *)malloc(nt*npos*sizeof(float));
+
+ 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];
+ }
+ }
+ 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==1) {
+ vmess("Marchenko representation with multiple sources");
+ }
+ 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++) {
+ sprintf(fins,"z%li",ir*dnumb+numb);
+ sprintf(fin2,"%s%s%s",fbegin,fins,fend);
+ fp_in = fopen(fin2, "r");
+ if (fp_in == NULL) {
+ 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 (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];
+ }
+ }
+ 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;
+ }
+ }
+ }
+ }
+ }
+
+ 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);
+ }
+}
+ free(shotdata);
+
+ if (verbose) vmess("nxs: %li nxz: %li f1: %.7f",nxs,nzs,f1);
+
+ ntshift=0;
+
+ 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");
+ }
+ }
+
+ 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;
+ }
+ ret = writeData3D(fp_out, &Ghom[(ir+ntshift)*nxs*nzs], hdr_out, nzs, nxs);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+
+ fclose(fp_out);
+ 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;
+
+
+ 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;
+ }
+ }
+ }
+
+ /* 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;
+}
+
+void timeDiff(float *data, long nsam, long nrec, float dt, float fmin, float fmax, long opt)
+{
+ long optn, iom, iomin, iomax, nfreq, ix, sign;
+ float omin, omax, deltom, om, df, *rdata, scl;
+ complex *cdata, *cdatascl;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+ df = 1.0/(optn*dt);
+
+ cdata = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data,nsam,nrec,optn,rdata);
+
+ /* Forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata[0], &cdata[0], optn, nrec, optn, nfreq, sign);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+ iomin = (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 = deltom*iom;
+ cdatascl[ix*nfreq+iom].r = -om*cdata[ix*nfreq+iom].i;
+ cdatascl[ix*nfreq+iom].i = om*cdata[ix*nfreq+iom].r;
+ }
+ }
+ else if (opt == -1) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = 1.0/(deltom*iom);
+ cdatascl[ix*nfreq+iom].r = om*cdata[ix*nfreq+iom].i;
+ cdatascl[ix*nfreq+iom].i = -om*cdata[ix*nfreq+iom].r;
+ }
+ }
+ 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;
+ }
+ }
+ }
+ free(cdata);
+
+ /* Inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&cdatascl[0], &rdata[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
+
+void depthDiff(float *data, 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;
+ float omin, omax, deltom, df, dkx, *rdata, kx, scl;
+ float kx2, kz2, kp2, kp;
+ complex *cdata, *cdatascl, kz, kzinv;
+
+ optn = optncr(nsam);
+ nfreq = optncr(nsam)/2+1;
+ df = 1.0/(optn*dt);
+ nkx = optncc(nrec);
+ dkx = 2.0*PI/(nkx*dx);
+ cdata = (complex *)malloc(nfreq*nkx*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nkx*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes in 2 directions to reach FFT lengths */
+ pad2d_data(data,nsam,nrec,optn,nkx,rdata);
+
+ /* double forward FFT */
+ xt2wkx(&rdata[0], &cdata[0], optn, nkx, optn, nkx, 0);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+
+ iomin = (long)MIN((omin/deltom), nfreq);
+ iomin = MAX(iomin, 0);
+ iomax = MIN((long)(omax/deltom), nfreq);
+
+ cdatascl = (complex *)malloc(nfreq*nkx*sizeof(complex));
+ if (cdatascl == NULL) verr("memory allocation error for cdatascl");
+
+ for (iom = 0; iom < iomin; iom++) {
+ for (ix = 0; ix < nkx; ix++) {
+ cdatascl[iom*nkx+ix].r = 0.0;
+ cdatascl[iom*nkx+ix].i = 0.0;
+ }
+ }
+ for (iom = iomax; iom < nfreq; iom++) {
+ for (ix = 0; ix < nkx; ix++) {
+ cdatascl[iom*nkx+ix].r = 0.0;
+ cdatascl[iom*nkx+ix].i = 0.0;
+ }
+ }
+ if (opt > 0) {
+ for (iom = iomin ; iom <= iomax ; iom++) {
+ kp = (iom*deltom)/c;
+ kp2 = kp*kp;
+
+ ikxmax = MIN((long)(kp/dkx), nkx/2);
+
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kz.r-cdata[iom*nkx+ikx].i*kz.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kz.r+cdata[iom*nkx+ikx].r*kz.i;
+
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nkx+ikx].r = 0.0;
+ cdatascl[iom*nkx+ikx].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kz.r-cdata[iom*nkx+ikx].i*kz.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kz.r+cdata[iom*nkx+ikx].r*kz.i;
+ }
+ }
+ }
+ else if (opt < 0) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ kp = iom*deltom/c;
+ kp2 = kp*kp;
+ ikxmax = MIN((long)(kp/dkx), nkx/2);
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kzinv.r-cdata[iom*nkx+ikx].i*kzinv.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kzinv.r+cdata[iom*nkx+ikx].r*kzinv.i;
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nkx+ikx].r = 0.0;
+ cdatascl[iom*nkx+ikx].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nkx+ikx].r = cdata[iom*nkx+ikx].r*kzinv.r-cdata[iom*nkx+ikx].i*kzinv.i;
+ cdatascl[iom*nkx+ikx].i = cdata[iom*nkx+ikx].i*kzinv.r+cdata[iom*nkx+ikx].r*kzinv.i;
+ }
+ }
+ }
+ free(cdata);
+
+ /* inverse double FFT */
+ wkx2xt(&cdatascl[0], &rdata[0], optn, nkx, nkx, optn, 0);
+ /* select original samples and traces */
+ scl = 1.0;
+ scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
+
+void pad2d_data(float *data, long nsam, long nrec, long nsamout, long nrecout, float *datout)
+{
+ long it,ix;
+ for (ix=0;ix<nrec;ix++) {
+ for (it=0;it<nsam;it++)
+ datout[ix*nsam+it]=data[ix*nsam+it];
+ for (it=nsam;it<nsamout;it++)
+ datout[ix*nsam+it]=0.0;
+ }
+ for (ix=nrec;ix<nrecout;ix++) {
+ for (it=0;it<nsamout;it++)
+ datout[ix*nsam+it]=0.0;
+ }
+}
+void conjugate(float *data, long nsam, long nrec, float dt)
+{
+ long optn, nfreq, j, ix, it, sign, ntdiff;
+ float *rdata, scl;
+ complex *cdata;
+
+ optn = optncr(nsam);
+ ntdiff = optn-nsam;
+ nfreq = optn/2+1;
+
+ cdata = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data,nsam,nrec,optn,rdata);
+
+ /* Forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata[0], &cdata[0], optn, nrec, optn, nfreq, sign);
+
+ /* take complex conjugate */
+ for(ix = 0; ix < nrec; ix++) {
+ for(j = 0; j < nfreq; j++) cdata[ix*nfreq+j].i = -cdata[ix*nfreq+j].i;
+ }
+
+ /* Inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&cdata[0], &rdata[0], optn, nrec, nfreq, optn, sign);
+ for (ix = 0; ix < nrec; ix++) {
+ for (it = 0 ; it < nsam ; it++)
+ data[ix*nsam+it] = scl*rdata[ix*optn+it+ntdiff];
+ }
+ //scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdata);
+ free(rdata);
+
+ return;
+}
+
+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;
+}
\ No newline at end of file
diff --git a/utils/Makefile b/utils/Makefile
index 27cb7fa..966672b 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 pwshift
+ALL: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d mutesnap padmodel truncate combine combine_induced reshape_su HomG snap2shot makeR1D pwshift
SRCM = \
makemod.c \
@@ -143,6 +143,26 @@ SRCMS = mutesnap.c \
docpkge.c \
readSnapData3D.c
+SRCPM = padmodel.c \
+ getFileInfo3D.c \
+ writeData3D.c \
+ wallclock_time.c \
+ getpars.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ readSnapData3D.c
+
+SRCTR = truncate.c \
+ getFileInfo3D.c \
+ writeData3D.c \
+ wallclock_time.c \
+ getpars.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ readSnapData3D.c
+
SRCCO = combine.c \
getFileInfo3D.c \
writeData3D.c \
@@ -274,6 +294,16 @@ OBJMS = $(SRCMS:%.c=%.o)
mutesnap: $(OBJMS)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o mutesnap $(OBJMS) $(LIBS)
+OBJPM = $(SRCPM:%.c=%.o)
+
+padmodel: $(OBJPM)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o padmodel $(OBJPM) $(LIBS)
+
+OBJTR = $(SRCTR:%.c=%.o)
+
+truncate: $(OBJTR)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o truncate $(OBJTR) $(LIBS)
+
OBJCO = $(SRCCO:%.c=%.o)
combine: $(OBJCO)
@@ -309,7 +339,7 @@ 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
+install: makemod makewave extendModel fconv correigen green green3D basop syn2d mat2su ftr1d mutesnap padmodel truncate combine combine_induced reshape_su HomG snap2shot makeR1D pwshift
cp makemod $B
cp makewave $B
cp extendModel $B
@@ -322,6 +352,8 @@ install: makemod makewave extendModel fconv correigen green green3D basop syn2d
cp mat2su $B
cp ftr1d $B
cp mutesnap $B
+ cp padmodel $B
+ cp truncate $B
cp combine $B
cp combine_induced $B
cp reshape_su $B
@@ -331,7 +363,7 @@ install: makemod makewave extendModel fconv correigen green green3D basop syn2d
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) pwshift $(OBJPW)
+ 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) truncate $(OBJTR) padmodel $(OBJPM) 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 $B/pwshift
+ 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/padmodel $B/truncate $B/combine $B/combine_induced $B/reshape_su $B/HomG $B/snap2shot $B/makeR1D $B/pwshift
diff --git a/utils/padmodel.c b/utils/padmodel.c
new file mode 100644
index 0000000..a4723dd
--- /dev/null
+++ b/utils/padmodel.c
@@ -0,0 +1,238 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+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);
+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);
+
+char *sdoc[] = {
+" ",
+" padmodel - extend a model with values",
+" ",
+" The order of the padding is over x, y and z, in case corners need to be padded",
+" ",
+" authors : Joeri Brackenhoff : (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke : (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" file_in= ................. File containing the first data",
+" ",
+" Optional parameters: ",
+" ",
+" file_out=out.su .......... Filename of the output",
+" nx=0 ..................... Number of samples to pad on both the left and right",
+" ny=0 ..................... Number of samples to pad on both the front and back",
+" nz=0 ..................... Number of samples to pad on only the bottom",
+" pad=0 .................... Pad option, (=0) pad with the value at the edge, (=1) pad with a constant value",
+" value=0.0 ................ Padding value if pad=1",
+" verbose=1 ................ Give detailed information of process",
+NULL};
+
+int main (int argc, char **argv)
+{
+ FILE *fp_in, *fp_out;
+ char *fin, *fout;
+ float *indata, *outdata, value;
+ float dz, dy, dx, z0, y0, x0, scl, z0out, y0out, x0out;
+ long nt, nz, ny, nx, ntr, ix, iy, iz;
+ long ret, verbose, pad;
+ long nxpad, nypad, nzpad, nxout, nyout, nzout;
+ segy *hdr_in, *hdr_out;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ /*----------------------------------------------------------------------------*
+ * Get the parameters passed to the function
+ *----------------------------------------------------------------------------*/
+ if (!getparstring("file_in", &fin)) fin = NULL;
+ if (!getparstring("file_out", &fout)) fout = "out.su";
+ if (!getparlong("nxpad", &nxpad)) nxpad=0;
+ if (!getparlong("nypad", &nypad)) nypad=0;
+ if (!getparlong("nzpad", &nzpad)) nzpad=0;
+ if (!getparlong("pad", &pad)) pad=0;
+ if (!getparfloat("value", &value)) value=0;
+ if (!getparlong("verbose", &verbose)) verbose=1;
+ if (fin == NULL) verr("No input file given");
+
+
+ /*----------------------------------------------------------------------------*
+ * Get the file info of the data and determine the indez of the truncation
+ *----------------------------------------------------------------------------*/
+
+ getFileInfo3D(fin, &nz, &nx, &ny, &nt, &dz, &dx, &dy, &z0, &x0, &y0, &scl, &ntr);
+
+ nxout = nx + 2*nxpad;
+ nyout = ny + 2*nypad;
+ nzout = nz + nzpad;
+
+ x0out = x0 - ((float)nxpad)*dx;
+ y0out = y0 - ((float)nypad)*dx;
+ z0out = z0;
+
+ if (verbose) {
+ vmess("******************** INPUT DATA ********************");
+ vmess("Number of samples: %li, x: %li, y: %li, z: %li",nx*ny*nz,nx,ny,nz);
+ vmess("Sampling distance for x: %.3f, y: %.3f, z: %.3f",dx,dy,dz);
+ vmess("Starting distance for x: %.3f, y: %.3f, z: %.3f",x0,y0,z0);
+ vmess("Final distance for x: %.3f, y: %.3f, z: %.3f",x0+((float)(nx-1))*dx,y0+((float)(ny-1))*dy,z0+((float)(nz-1))*dz);
+ vmess("******************** PADDING ********************");
+ if (!pad) vmess("Model is padded using the edge values");
+ else vmess("model is padded using constant value %.3f",value);
+ vmess("Number of padding samples: x: %li, y: %li, z: %li",2*nxpad,2*nypad,nzpad);
+ vmess("******************** OUTPUT DATA ********************");
+ vmess("Number of samples: %li, x: %li, y: %li, z: %li",nxout*nyout*nzout,nxout,nyout,nzout);
+ vmess("Sampling distance for x: %.3f, y: %.3f, z: %.3f",dx,dy,dz);
+ vmess("Starting distance for x: %.3f, y: %.3f, z: %.3f",x0out,y0out,z0out);
+ vmess("Final distance for x: %.3f, y: %.3f, z: %.3f",x0out+((float)(nxout-1))*dx,y0out+((float)(nyout-1))*dy,z0out+((float)(nzout-1))*dz);
+ }
+
+ /*----------------------------------------------------------------------------*
+ * Allocate and read in the data
+ *----------------------------------------------------------------------------*/
+
+ indata = (float *)calloc(nx*ny*nz,sizeof(float));
+ hdr_in = (segy *)calloc(nx*ny,sizeof(segy));
+ readSnapData3D(fin, indata, hdr_in, 1, nx, ny, nz, 0, nx, 0, ny, 0, nz);
+ if (verbose) vmess("Read data");
+
+ outdata = (float *)calloc(nxout*nyout*nzout,sizeof(float));
+ hdr_out = (segy *)calloc(nxout*nyout,sizeof(segy));
+
+ /*----------------------------------------------------------------------------*
+ * Fit the original data into the new data
+ *----------------------------------------------------------------------------*/
+
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (iz = 0; iz < nz; iz++) {
+ outdata[(iy+nypad)*nxout*nzout+(ix+nxpad)*nzout+iz] = indata[iy*nx*nz+ix*nz+iz];
+ }
+ }
+ }
+ if (verbose) vmess("Original data fitted into new model");
+ free(indata); free(hdr_in);
+
+ /*----------------------------------------------------------------------------*
+ * Pad with a constant value
+ *----------------------------------------------------------------------------*/
+
+ if (pad) {
+ if (verbose) vmess("Padding with constant value %.3f",value);
+ for (iz = 0; iz < nz; iz++) {
+ for (iy = nypad; iy < ny+nypad; iy++) {
+ for (ix = 0; ix < nxpad; ix++) {
+ outdata[iy*nxout*nzout+ix*nzout+iz] = value;
+ }
+ for (ix = nxpad+nx; ix < nxout; ix++) {
+ outdata[iy*nxout*nzout+ix*nzout+iz] = value;
+ }
+ }
+ for (ix = 0; ix < nxout; ix++) {
+ for (iy = 0; iy < nypad; iy++) {
+ outdata[iy*nxout*nzout+ix*nzout+iz] = value;
+ }
+ for (iy = nypad+ny; iy < nyout; iy++) {
+ outdata[iy*nxout*nzout+ix*nzout+iz] = value;
+ }
+ }
+ }
+ for (iz = nz; iz < nzout; iz++) {
+ for (iy = 0; iy < nyout; iy++) {
+ for (ix = 0; ix < nxout; ix++) {
+ outdata[iy*nxout*nzout+ix*nzout+iz] = value;
+ }
+ }
+ }
+ if (verbose) vmess("Padded data");
+ }
+
+ /*----------------------------------------------------------------------------*
+ * Pad with edges
+ *----------------------------------------------------------------------------*/
+
+ if (!pad) {
+ if (verbose) vmess("Padding with edge value");
+ for (iz = 0; iz < nz; iz++) {
+ for (iy = nypad; iy < ny+nypad; iy++) {
+ for (ix = nxpad-1; ix >= 0; ix--) {
+ outdata[iy*nxout*nzout+ix*nzout+iz] = outdata[iy*nxout*nzout+nxpad*nzout+iz];
+ }
+ for (ix = nxpad+nx; ix < nxout; ix++) {
+ outdata[iy*nxout*nzout+ix*nzout+iz] = outdata[iy*nxout*nzout+(nxpad+nx-1)*nzout+iz];
+ }
+ }
+ for (ix = 0; ix < nxout; ix++) {
+ for (iy = nypad-1; iy >= 0; iy--) {
+ outdata[iy*nxout*nzout+ix*nzout+iz] = outdata[nypad*nxout*nzout+ix*nzout+iz];
+ }
+ for (iy = nypad+ny; iy < nyout; iy++) {
+ outdata[iy*nxout*nzout+ix*nzout+iz] = outdata[(nypad+ny-1)*nxout*nzout+ix*nzout+iz];
+ }
+ }
+ }
+ for (iz = nz; iz < nzout; iz++) {
+ for (iy = 0; iy < nyout; iy++) {
+ for (ix = 0; ix < nxout; ix++) {
+ outdata[iy*nxout*nzout+ix*nzout+iz] = outdata[iy*nxout*nzout+ix*nzout+nz-1];
+ }
+ }
+ }
+ if (verbose) vmess("Padded data");
+ }
+
+ /*----------------------------------------------------------------------------*
+ * Write out the data
+ *----------------------------------------------------------------------------*/
+
+ for (iy = 0; iy < nyout; iy++) {
+ for (ix = 0; ix < nxout; ix++) {
+ hdr_out[iy*nxout+ix].tracl = iy*nxout+ix+1;
+ hdr_out[iy*nxout+ix].trid = 130;
+ hdr_out[iy*nxout+ix].scalco = -1000;
+ hdr_out[iy*nxout+ix].scalel = -1000;
+ hdr_out[iy*nxout+ix].gx = (long)(1000.0*(x0out + ((float)ix)*dx));
+ hdr_out[iy*nxout+ix].gy = (long)(1000.0*(y0out + ((float)iy)*dy));
+ hdr_out[iy*nxout+ix].ns = nzout;
+ hdr_out[iy*nxout+ix].trwf = nxout*nyout;
+ hdr_out[iy*nxout+ix].d1 = dz;
+ hdr_out[iy*nxout+ix].d2 = dx;
+ hdr_out[iy*nxout+ix].f1 = z0out;
+ hdr_out[iy*nxout+ix].f2 = x0out;
+ }
+ }
+
+ fp_out = fopen(fout, "w+");
+
+ ret = writeData3D(fp_out, &outdata[0], hdr_out, nzout, nxout*nyout);
+ vmess("Wrote data");
+
+ fclose(fp_out);
+ free(outdata); free(hdr_out);
+ return 0;
+}
+
diff --git a/utils/pwshift.c b/utils/pwshift.c
index 54999a0..afbaa72 100644
--- a/utils/pwshift.c
+++ b/utils/pwshift.c
@@ -48,6 +48,14 @@ char *sdoc[] = {
" Optional parameters: ",
" ",
" file_out=out.su .......... Filename of the output",
+" fmin=0.0 ................. Minimum of the frequency range",
+" fmax=70.0 ................ Maximum of the frequency range",
+" src_velox=1500.0 ......... Velocity of the x-angle of the plane wave",
+" src_veloy=1500.0 ......... Velocity of the y-angle of the plane wave",
+" src_anglex=0.0 ........... x-angle of the plane wave",
+" src_angley=0.0 ........... y-angle of the plane wave",
+" numb=0 ................... Starting number of the level for the image",
+" dnumb=1 .................. Increment number of the level for the image",
" verbose=1 ................ Give detailed information of process",
NULL};
@@ -55,10 +63,10 @@ int main (int argc, char **argv)
{
FILE *fp_gmin, *fp_ray, *fp_amp, *fp_out;
char *file_gmin, *file_ray, *file_amp, *file_out, fbr[100], fer[100], fba[100], fea[100], fins[100], fin2[100], numb1[100], *ptr;
- float *gmin, *conv, *time, *amp, *image, fmin, fmax, *delay;
- float dt, dy, dx, dz, t0, y0, x0, scl;
+ float *gmin, *conv, *time, *amp, *image, fmin, fmax;
+ float dt, dy, dx, dz, t0, y0, x0, scl, *shift;
float dt_ray, dy_ray, dx_ray, t0_ray, y0_ray, x0_ray, scl_ray, px, py, src_velox, src_veloy, src_anglex, src_angley, grad2rad;
- long nshots, nt, ny, nx, ntr;
+ long nshots, nt, ny, nx, ntr, delay;
long nray, nt_ray, ny_ray, nx_ray, ntr_ray;
long verbose, ix, iy, it, iz, is, *gx, *gy, *gz, numb, dnumb, pos, nzs, file_det;
size_t ret;
@@ -147,6 +155,7 @@ int main (int argc, char **argv)
vmess("***********************************************************");
}
+ nray = 0;
sprintf(fins,"z%li",0);
sprintf(file_ray,"%s%s%s",fbr,fins,fer);
getFileInfo3D(file_ray, &nx_ray, &nt_ray, &ny_ray, &nray, &dx_ray, &dt_ray, &dy_ray, &x0_ray, &t0_ray, &y0_ray, &scl_ray, &ntr_ray);
@@ -170,7 +179,6 @@ int main (int argc, char **argv)
hdr_gmin = (segy *)calloc(nshots*nx*ny,sizeof(segy));
gmin = (float *)calloc(nshots*nx*ny*nt,sizeof(float));
hdr_out = (segy *)calloc(nray*nshots,sizeof(segy));
- delay = (float *)calloc(nx*ny,sizeof(float));
image = (float *)calloc(nray*nshots,sizeof(float));
gx = (long *)calloc(nray*nshots,sizeof(long));
@@ -178,43 +186,46 @@ int main (int argc, char **argv)
gz = (long *)calloc(nray*nshots,sizeof(long));
readSnapData3D(file_gmin, gmin, hdr_gmin, nshots, nx, ny, nt, 0, nx, 0, ny, 0, nt);
+ if (verbose) vmess("Read in Gmin data");
/*----------------------------------------------------------------------------*
* Add the delay in case the plane wave is at an angle
*----------------------------------------------------------------------------*/
+ shift = (float *)calloc(nx*ny,sizeof(float));
grad2rad = 17.453292e-3;
px = sin(src_anglex*grad2rad)/src_velox;
py = sin(src_angley*grad2rad)/src_veloy;
-
- if (py < 0.0) {
- for (iy=0; iy<ny; iy++) {
- if (px < 0.0) {
- for (ix=0; ix<nx; ix++) {
- delay[iy*nx+ix] = fabsf((nx-1-ix)*dx*px) + fabsf((ny-1-iy)*dy*py);
- }
- }
- else {
- for (ix=0; ix<nx; ix++) {
- delay[iy*nx+ix] = ix*dx*px + fabsf((ny-1-iy)*dy*py);
- }
- }
- }
- }
- else {
- for (iy=0; iy<ny; iy++) {
- if (px < 0.0) {
- for (ix=0; ix<nx; ix++) {
- delay[iy*nx+ix] = fabsf((nx-1-ix)*dx*px) + iy*dy*py;
- }
- }
- else {
- for (ix=0; ix<nx; ix++) {
- delay[iy*nx+ix] = ix*dx*px + iy*dy*py;
- }
- }
- }
- }
+ if (verbose) vmess("px value is %f and py value is %f",px,py);
+
+ // if (py < 0.0) {
+ // for (iy=0; iy<ny; iy++) {
+ // if (px < 0.0) {
+ // for (ix=0; ix<nx; ix++) {
+ // shift[iy*nx+ix] = fabsf((nx-1-ix)*dx*px) + fabsf((ny-1-iy)*dy*py);
+ // }
+ // }
+ // else {
+ // for (ix=0; ix<nx; ix++) {
+ // shift[iy*nx+ix] = ix*dx*px + fabsf((ny-1-iy)*dy*py);
+ // }
+ // }
+ // }
+ // }
+ // else {
+ // for (iy=0; iy<ny; iy++) {
+ // if (px < 0.0) {
+ // for (ix=0; ix<nx; ix++) {
+ // shift[iy*nx+ix] = fabsf((nx-1-ix)*dx*px) + iy*dy*py;
+ // }
+ // }
+ // else {
+ // for (ix=0; ix<nx; ix++) {
+ // shift[iy*nx+ix] = ix*dx*px + iy*dy*py;
+ // }
+ // }
+ // }
+ // }
/*----------------------------------------------------------------------------*
* Apply the imaging condition
@@ -235,16 +246,40 @@ int main (int argc, char **argv)
readSnapData3D(fin2, time, hdr_time, nray, 1, ny, nx, 0, 1, 0, ny, 0, nx);
sprintf(fin2,"%s%s%s",fba,fins,fea);
readSnapData3D(file_amp, amp, hdr_amp, nray, 1, ny, nx, 0, 1, 0, ny, 0, nx);
+
for (is = 0; is < nray; is++) {
gx[is*nshots+iy] = hdr_time[is*ny].sx;
gy[is*nshots+iy] = hdr_time[is*ny].sy;
gz[is*nshots+iy] = hdr_time[is*ny].sdepth;
+
+ x0_ray = ((float)gx[is*nshots+iy])/1000.0;
+ y0_ray = ((float)gy[is*nshots+iy])/1000.0;
+
+ if (py < 0.0) {
+ if (px < 0.0) {
+ delay = NINT(fabsf((x0-x0_ray)*px)/dt) + NINT(fabsf((y0-y0_ray)*py)/dt);
+ }
+ else {
+ delay = NINT((x0_ray-x0)*px/dt) + NINT(fabsf((y0-y0_ray)*py)/dt);
+ }
+ }
+ else {
+ if (px < 0.0) {
+ delay = NINT(fabsf((x0-x0_ray)*px)/dt) + NINT((y0_ray-y0)*py/dt);
+ }
+ else {
+ delay = NINT((x0_ray-x0)*px/dt) + NINT((y0_ray-y0)*py/dt);
+ }
+ }
+
+ if (delay > nt-1) delay = nt-1;
+
for (it = 0; it < nx*ny*nt; it++) {
conv[it] = gmin[iy*nx*ny*nt+it];
}
- timeShift(&conv[0],nt,nx*ny,dt,&time[is*ny*nx],&[is*ny*nx],delay,fmin,fmax);
+ timeShift(&conv[0],nt,nx*ny,dt,&time[is*ny*nx],&[is*ny*nx],shift,fmin,fmax);
for (ix = 0; ix < ny*nx; ix++) {
- image[is*nshots+iy] += conv[ix*nt]*dx*dy*dt;
+ image[is*nshots+iy] += conv[ix*nt+delay]*dx*dy*dt;
}
}
@@ -334,7 +369,7 @@ void timeShift(float *data, long nsam, long nrec, float dt, float *time, float *
}
for (iom = iomin ; iom < iomax ; iom++) {
om = deltom*iom;
- tom = om*-1.0*(time[ix]+delay[ix]);
+ tom = om*-1.0*(time[ix]-delay[ix]);
cdatascl[ix*nfreq+iom].r = (cdata[ix*nfreq+iom].r*cos(-tom) - cdata[ix*nfreq+iom].i*sin(-tom))/(amp[ix]*amp[ix]);
cdatascl[ix*nfreq+iom].i = (cdata[ix*nfreq+iom].i*cos(-tom) + cdata[ix*nfreq+iom].r*sin(-tom))/(amp[ix]*amp[ix]);
}
diff --git a/utils/pwshift_backup23apr2020.c b/utils/pwshift_backup23apr2020.c
new file mode 100644
index 0000000..5425b48
--- /dev/null
+++ b/utils/pwshift_backup23apr2020.c
@@ -0,0 +1,375 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+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);
+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 timeShift(float *data, long nsam, long nrec, float dt, float *time, float *amp, float *delay, float fmin, float fmax);
+void pad_data(float *data, long nsam, long nrec, long nsamout, float *datout);
+void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout);
+
+char *sdoc[] = {
+" ",
+" combine - Combine results into a single result ",
+" ",
+" authors : Joeri Brackenhoff : (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke : (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" file_gmin= ................. File containing the G- data",
+" file_ray= .................. File containing the ray time data",
+" file_amp= .................. File containing the ray amplitude data",
+" ",
+" Optional parameters: ",
+" ",
+" file_out=out.su .......... Filename of the output",
+" verbose=1 ................ Give detailed information of process",
+NULL};
+
+int main (int argc, char **argv)
+{
+ FILE *fp_gmin, *fp_ray, *fp_amp, *fp_out;
+ char *file_gmin, *file_ray, *file_amp, *file_out, fbr[100], fer[100], fba[100], fea[100], fins[100], fin2[100], numb1[100], *ptr;
+ float *gmin, *conv, *time, *amp, *image, fmin, fmax, *delay;
+ float dt, dy, dx, dz, t0, y0, x0, scl;
+ float dt_ray, dy_ray, dx_ray, t0_ray, y0_ray, x0_ray, scl_ray, px, py, src_velox, src_veloy, src_anglex, src_angley, grad2rad;
+ long nshots, nt, ny, nx, ntr;
+ long nray, nt_ray, ny_ray, nx_ray, ntr_ray;
+ long verbose, ix, iy, it, iz, is, *gx, *gy, *gz, numb, dnumb, pos, nzs, file_det;
+ size_t ret;
+ segy *hdr_gmin, *hdr_time, *hdr_amp, *hdr_out;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ /*----------------------------------------------------------------------------*
+ * Get the parameters passed to the function
+ *----------------------------------------------------------------------------*/
+ if (!getparstring("file_gmin", &file_gmin)) file_gmin = NULL;
+ if (!getparstring("file_ray", &file_ray)) file_ray = NULL;
+ if (!getparstring("file_amp", &file_amp)) file_amp = NULL;
+ if (!getparstring("file_out", &file_out)) file_out = "out.su";
+ if (!getparlong("verbose", &verbose)) verbose=1;
+ if (!getparfloat("fmin", &fmin)) fmin=0.0;
+ if (!getparfloat("fmax", &fmax)) fmax=70.0;
+ if (!getparfloat("src_velox", &src_velox)) src_velox=1500.0;
+ if (!getparfloat("src_veloy", &src_veloy)) src_veloy=1500.0;
+ if (!getparfloat("src_anglex", &src_anglex)) src_anglex=0.0;
+ if (!getparfloat("src_angley", &src_angley)) src_angley=0.0;
+ if (!getparlong("numb", &numb)) numb=0;
+ if (!getparlong("dnumb", &dnumb)) dnumb=1;
+ if (file_gmin == NULL) verr("Incorrect gmin input");
+ if (file_ray == NULL) verr("Incorrect ray time input");
+ if (file_amp == NULL) verr("Incorrect ray amplitude input");
+
+ /*----------------------------------------------------------------------------*
+ * Determine the position of the number in the string
+ * and split the file into beginning, middle and end
+ *----------------------------------------------------------------------------*/
+ if (dnumb < 1) dnumb = 1;
+ sprintf(numb1,"z%li",numb);
+
+ ptr = strstr(file_ray,numb1);
+ pos = ptr - file_ray + 1;
+ sprintf(fbr,"%*.*s", pos-1, pos-1, file_ray);
+ sprintf(fer,"%s", file_ray+pos+1);
+
+ ptr = strstr(file_amp,numb1);
+ pos = ptr - file_amp + 1;
+ sprintf(fba,"%*.*s", pos-1, pos-1, file_ray);
+ sprintf(fea,"%s", file_ray+pos+1);
+
+ /*----------------------------------------------------------------------------*
+ * Determine the amount of files that are present
+ *----------------------------------------------------------------------------*/
+ file_det = 1;
+ nzs=0;
+ while (file_det) { // Check for a file with the filename
+ sprintf(fins,"z%li",nzs*dnumb+numb);
+ sprintf(file_ray,"%s%s%s",fbr,fins,fer);
+ fp_ray = fopen(file_ray, "r");
+ if (fp_ray == NULL) { // If the filename does not exist
+ if (nzs == 0) { // The filename is wrong to begin with
+ verr("error on opening basefile=%s", file_ray);
+ }
+ else if (nzs == 1) { // There is only a single file
+ vmess("1 file detected");
+ file_det = 0;
+ break;
+ }
+ else { // Stop after the final file has been detected
+ vmess("%li files detected",nzs);
+ file_det = 0;
+ break;
+ }
+ }
+ fclose(fp_ray);
+ nzs++;
+ }
+
+ /*----------------------------------------------------------------------------*
+ * Read in the first two files and determine the header values
+ * of the output
+ *----------------------------------------------------------------------------*/
+ getFileInfo3D(file_gmin, &nt, &nx, &ny, &nshots, &dt, &dx, &dy, &t0, &x0, &y0, &scl, &ntr);
+
+ if (verbose) {
+ vmess("************************ Gmin info ************************");
+ vmess("Number of depth levels : %li",nshots);
+ vmess("Number of samples x: %li, y: %li, t: %li",nx,ny,nt);
+ vmess("Starting distance for x: %.3f, y: %.3f, t: %.3f",x0,y0,t0);
+ vmess("Sampling distance for x: %.3f, y: %.3f, t: %.3f",dx,dy,dt);
+ vmess("***********************************************************");
+ }
+
+ sprintf(fins,"z%li",0);
+ sprintf(file_ray,"%s%s%s",fbr,fins,fer);
+ getFileInfo3D(file_ray, &nx_ray, &nt_ray, &ny_ray, &nray, &dx_ray, &dt_ray, &dy_ray, &x0_ray, &t0_ray, &y0_ray, &scl_ray, &ntr_ray);
+
+ if (verbose) {
+ vmess("************************ ray info *************************");
+ vmess("Number of depth levels : %li",nzs);
+ vmess("Number of focal points : %li",nray);
+ vmess("Number of samples x: %li, y: %li, t: %li",nx_ray,ny_ray,nt_ray);
+ vmess("Starting distance for x: %.3f, y: %.3f, t: %.3f",x0_ray,y0_ray,t0_ray);
+ vmess("Sampling distance for x: %.3f, y: %.3f, t: %.3f",dx_ray,dy_ray,dt_ray);
+ vmess("***********************************************************");
+ }
+
+ if (nshots!=nzs) verr("The depth levels of the rays (%li) do not match those of the plane waves (%li)",nzs,nshots);
+
+ /*----------------------------------------------------------------------------*
+ * Read in a single file to determine if the header values match
+ * and allocate the data
+ *----------------------------------------------------------------------------*/
+ hdr_gmin = (segy *)calloc(nshots*nx*ny,sizeof(segy));
+ gmin = (float *)calloc(nshots*nx*ny*nt,sizeof(float));
+ hdr_out = (segy *)calloc(nray*nshots,sizeof(segy));
+ delay = (float *)calloc(nx*ny,sizeof(float));
+
+ image = (float *)calloc(nray*nshots,sizeof(float));
+ gx = (long *)calloc(nray*nshots,sizeof(long));
+ gy = (long *)calloc(nray*nshots,sizeof(long));
+ gz = (long *)calloc(nray*nshots,sizeof(long));
+
+ readSnapData3D(file_gmin, gmin, hdr_gmin, nshots, nx, ny, nt, 0, nx, 0, ny, 0, nt);
+ if (verbose) vmess("Read in Gmin data");
+
+ /*----------------------------------------------------------------------------*
+ * Add the delay in case the plane wave is at an angle
+ *----------------------------------------------------------------------------*/
+
+ grad2rad = 17.453292e-3;
+ px = sin(src_anglex*grad2rad)/src_velox;
+ py = sin(src_angley*grad2rad)/src_veloy;
+
+ if (py < 0.0) {
+ for (iy=0; iy<ny; iy++) {
+ if (px < 0.0) {
+ for (ix=0; ix<nx; ix++) {
+ delay[iy*nx+ix] = fabsf((nx-1-ix)*dx*px) + fabsf((ny-1-iy)*dy*py);
+ }
+ }
+ else {
+ for (ix=0; ix<nx; ix++) {
+ delay[iy*nx+ix] = ix*dx*px + fabsf((ny-1-iy)*dy*py);
+ }
+ }
+ }
+ }
+ else {
+ for (iy=0; iy<ny; iy++) {
+ if (px < 0.0) {
+ for (ix=0; ix<nx; ix++) {
+ delay[iy*nx+ix] = fabsf((nx-1-ix)*dx*px) + iy*dy*py;
+ }
+ }
+ else {
+ for (ix=0; ix<nx; ix++) {
+ delay[iy*nx+ix] = ix*dx*px + iy*dy*py;
+ }
+ }
+ }
+ }
+
+ /*----------------------------------------------------------------------------*
+ * Apply the imaging condition
+ *----------------------------------------------------------------------------*/
+#pragma omp parallel for schedule(static,1) default(shared) \
+ private(is,it,ix,fins,fin2,time,amp,hdr_time,hdr_amp,conv)
+ for (iy = 0; iy < nshots; iy++) {
+
+ hdr_time = (segy *)calloc(ny*nray,sizeof(segy));
+ time = (float *)calloc(nx*ny*nray,sizeof(float));
+ hdr_amp = (segy *)calloc(ny*nray,sizeof(segy));
+ amp = (float *)calloc(nx*ny*nray,sizeof(float));
+ conv = (float *)calloc(nx*ny*nt,sizeof(float));
+
+ vmess("Depth level %li out of %li",iy+1,nshots);
+ sprintf(fins,"z%li",iy*dnumb+numb);
+ sprintf(fin2,"%s%s%s",fbr,fins,fer);
+ readSnapData3D(fin2, time, hdr_time, nray, 1, ny, nx, 0, 1, 0, ny, 0, nx);
+ sprintf(fin2,"%s%s%s",fba,fins,fea);
+ readSnapData3D(file_amp, amp, hdr_amp, nray, 1, ny, nx, 0, 1, 0, ny, 0, nx);
+ for (is = 0; is < nray; is++) {
+ gx[is*nshots+iy] = hdr_time[is*ny].sx;
+ gy[is*nshots+iy] = hdr_time[is*ny].sy;
+ gz[is*nshots+iy] = hdr_time[is*ny].sdepth;
+ for (it = 0; it < nx*ny*nt; it++) {
+ conv[it] = gmin[iy*nx*ny*nt+it];
+ }
+ timeShift(&conv[0],nt,nx*ny,dt,&time[is*ny*nx],&[is*ny*nx],delay,fmin,fmax);
+ for (ix = 0; ix < ny*nx; ix++) {
+ image[is*nshots+iy] += conv[ix*nt]*dx*dy*dt;
+ }
+ }
+
+ free(hdr_time); free(time); free(hdr_amp); free(amp); free(conv);
+ }
+ free(gmin); free(hdr_gmin);
+
+ /*----------------------------------------------------------------------------*
+ * Write out the data
+ *----------------------------------------------------------------------------*/
+ fp_out = fopen(file_out, "w+");
+
+ if (nshots>1) dz = ((float)(gz[1] - gz[0]))/1000.0;
+ else dz = 1.0;
+
+ for (it = 0; it < nray; it++) {
+ hdr_out[it].fldr = 1;
+ hdr_out[it].tracl = it+1;
+ hdr_out[it].tracf = it+1;
+ hdr_out[it].scalco = -1000;
+ hdr_out[it].scalel = -1000;
+ hdr_out[it].trid = 1;
+ hdr_out[it].ns = nshots;
+ hdr_out[it].trwf = nray;
+ hdr_out[it].ntr = nray;
+ hdr_out[it].f1 = (((float)gz[0])/1000.0);
+ hdr_out[it].f2 = (((float)gx[0])/1000.0);
+ hdr_out[it].dt = ((int)(dt*1E6));
+ hdr_out[it].d1 = roundf(dz*1000.0)/1000.0;
+ hdr_out[it].d2 = roundf(dx*1000.0)/1000.0;
+ hdr_out[it].gx = gx[it*nshots];
+ hdr_out[it].gy = gy[it*nshots];
+ }
+
+ ret = writeData3D(fp_out, &image[0], hdr_out, nshots, nray);
+ if (ret < 0 ) verr("error on writing output file.");
+
+ fclose(fp_out);
+
+ free(image); free(hdr_out);
+
+ vmess("Wrote data");
+
+ return 0;
+}
+
+void timeShift(float *data, long nsam, long nrec, float dt, float *time, float *amp, float *delay, float fmin, float fmax)
+{
+ long optn, iom, iomin, iomax, nfreq, ix, sign;
+ float omin, omax, deltom, om, tom, df, *rdata, scl;
+ complex *cdata, *cdatascl;
+
+ optn = optncr(nsam);
+ nfreq = optn/2+1;
+ df = 1.0/(optn*dt);
+
+ cdata = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nrec*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes until Fourier length is reached */
+ pad_data(data,nsam,nrec,optn,rdata);
+
+ /* Forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&rdata[0], &cdata[0], optn, nrec, optn, nfreq, sign);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+ iomin = (long)MIN((omin/deltom), (nfreq));
+ 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;
+ }
+ for (iom = iomin ; iom < iomax ; iom++) {
+ om = deltom*iom;
+ tom = om*-1.0*(time[ix]-delay[ix]);
+ cdatascl[ix*nfreq+iom].r = (cdata[ix*nfreq+iom].r*cos(-tom) - cdata[ix*nfreq+iom].i*sin(-tom))/(amp[ix]*amp[ix]);
+ cdatascl[ix*nfreq+iom].i = (cdata[ix*nfreq+iom].i*cos(-tom) + cdata[ix*nfreq+iom].r*sin(-tom))/(amp[ix]*amp[ix]);
+ }
+ }
+ free(cdata);
+
+ /* Inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/(float)optn;
+ crmfft(&cdatascl[0], &rdata[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata,optn,nrec,scl,data,nsam);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
+
+void 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];
+ }
+}
\ No newline at end of file
diff --git a/utils/testfft3d.c b/utils/testfft3d.c
new file mode 100644
index 0000000..a5ec6cd
--- /dev/null
+++ b/utils/testfft3d.c
@@ -0,0 +1,359 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+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);
+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 pad3d_data(float *data, long nt, long nx, long ny, long ntout, long nxout, long nyout, float *datout);
+void scl_data3D(float *data, long nt, long nx, long ny, float scl, float *datout, long ntout, long nxout);
+void depthDiff3D(float *data, long nt, long nx, long ny, float dt, float dx, float dy, float fmin, float fmax, float c, int opt);
+
+char *sdoc[] = {
+" ",
+" testfft3d ",
+NULL};
+
+int main (int argc, char **argv)
+{
+ FILE *fp_in, *fp_out;
+ char *fin, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100], numb1[100], *direction;
+ float *indata, *rdata;
+ float dz, dy, dx, z0, y0, x0, scl, cp, fmin, fmax;
+ long nt, nz, ny, nx, ntr, ix, iy, it, is, iz, pos, file_det, nzs, dt;
+ long numb, dnumb, ret, nzmax, compact, verbose, nxyz, sx, sy, sz;
+ long nxout, nyout, nzout, ixout, iyout, izout;
+ long nf, nft, nkx, nky;
+ complex *cdata;
+ segy *hdr_in, *hdr_out;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ /*----------------------------------------------------------------------------*
+ * Get the parameters passed to the function
+ *----------------------------------------------------------------------------*/
+ if (!getparstring("file_in", &fin)) fin = NULL;
+ if (!getparstring("file_out", &fout)) fout = "out.su";
+ if (!getparlong("verbose", &verbose)) verbose=1;
+ if(!getparfloat("cp", &cp)) cp = 1500.0;
+ if(!getparfloat("fmin", &fmin)) fmin = 0.0;
+ if (fin == NULL) verr("Incorrect downgoing input");
+
+ /*----------------------------------------------------------------------------*
+ * Read in the first two files and determine the header values
+ * of the output
+ *----------------------------------------------------------------------------*/
+
+ getFileInfo3D(fin, &nt, &nx, &ny, &nz, &dz, &dx, &dy, &z0, &x0, &y0, &scl, &ntr);
+
+ if (verbose) {
+ vmess("number of time samples: %li", nt);
+ vmess("Number of virtual receivers: %li, x: %li, y: %li",nx*ny,nx,ny);
+ vmess("Starting distance for x: %.3f, y: %.3f",x0,y0);
+ vmess("Sampling distance for x: %.3f, y: %.3f t:%.3f",dx,dy,dz);
+ }
+ if(!getparfloat("fmax", &fmax)) fmax = 1.0/(2.0*dz);
+
+ /*----------------------------------------------------------------------------*
+ * Read in a single file to determine if the header values match
+ * and allocate the data
+ *----------------------------------------------------------------------------*/
+ indata = (float *)calloc(nx*ny*nt,sizeof(float));
+ hdr_in = (segy *)calloc(nx*ny,sizeof(segy));
+
+ readSnapData3D(fin, indata, hdr_in, 1, nx, ny, nt, 0, nx, 0, ny, 0, nt);
+ vmess("Read data");
+
+ depthDiff3D(&indata[0], nt, nx, ny, dz, dx, dy, fmin, fmax, cp, 1);
+ vmess("Applied transform");
+
+ // nft = loptncr(nt);
+ // nf = nft/2+1;
+ // nkx = optncc(nx);
+ // nky = optncc(ny);
+
+ // rdata = (float *)calloc(nkx*nky*nft,sizeof(float));
+ // cdata = (complex *)calloc(nkx*nky*nf,sizeof(complex));
+
+ // pad3d_data(indata, nt, nx, ny, nft, nkx, nky, rdata);
+
+ // yxt2wkykx(rdata, cdata, nft, nkx, nky, nft, nkx, nky, 0, 0);
+
+ // wkykx2yxt(cdata, rdata, nft, nkx, nky, nft, nkx, nky, 0, 0);
+
+ // scl_data3D(rdata, nft, nx, ny, 1.0, indata, nt, nx);
+
+ /*----------------------------------------------------------------------------*
+ * Write out the data
+ *----------------------------------------------------------------------------*/
+ fp_out = fopen(fout, "w+");
+
+ ret = writeData3D(fp_out, &indata[0], hdr_in, nt, nx*ny);
+ if (ret < 0 ) verr("error on writing output file.");
+
+ fclose(fp_out);
+ free(indata);
+ vmess("Wrote data");
+ return 0;
+}
+
+void depthDiff3D(float *data, long nt, long nx, long ny, float dt, float dx, float dy, float fmin, float fmax, float c, int opt)
+{
+ long optn, iom, iomin, iomax, nfreq, ix, iy, ikx, iky, nkx, nky, ikxmax, ikymax;
+ float omin, omax, deltom, df, dkx, dky, *rdata, kx, ky, scl;
+ float kx2, ky2, kz2, kp2, kp;
+ complex *cdata, *cdatascl, kz, kzinv;
+
+ optn = optncr(nt);
+ nfreq = optncr(nt)/2+1;
+ df = 1.0/(optn*dt);
+ nkx = optncc(nx);
+ nky = optncc(ny);
+ dkx = 2.0*PI/(nkx*dx);
+ dky = 2.0*PI/(nky*dy);
+ cdata = (complex *)malloc(nfreq*nkx*nky*sizeof(complex));
+ if (cdata == NULL) verr("memory allocation error for cdata");
+
+ rdata = (float *)malloc(optn*nkx*nky*sizeof(float));
+ if (rdata == NULL) verr("memory allocation error for rdata");
+
+ /* pad zeroes in 2 directions to reach FFT lengths */
+ pad3d_data(data, nt, nx, ny, optn, nkx, nky, rdata);
+
+ /* double forward FFT */
+ yxt2wkykx(&rdata[0], &cdata[0], optn, nkx, nky, optn, nkx, nky, 0, 0);
+
+ deltom = 2.*PI*df;
+ omin = 2.*PI*fmin;
+ omax = 2.*PI*fmax;
+
+ iomin = (int)MIN((omin/deltom), nfreq);
+ iomin = MAX(iomin, 0);
+ iomax = MIN((int)(omax/deltom), nfreq);
+
+ cdatascl = (complex *)malloc(nfreq*nkx*nky*sizeof(complex));
+ if (cdatascl == NULL) verr("memory allocation error for cdatascl");
+
+ for (iom = 0; iom < iomin; iom++) {
+ for (iy = 0; iy < nky; iy++) {
+ for (ix = 0; ix < nkx; ix++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ }
+ }
+ for (iom = iomax; iom < nfreq; iom++) {
+ for (iy = 0; iy < nky; iy++) {
+ for (ix = 0; ix < nkx; ix++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ }
+ }
+ if (opt > 0) {
+ for (iom = iomin ; iom <= iomax ; iom++) {
+ kp = (iom*deltom)/c;
+ kp2 = kp*kp;
+ ikxmax = MIN((int)(kp/dkx), nkx/2);
+ ikymax = MIN((int)(kp/dky), nky/2);
+
+ for (iky = 0; iky < ikymax; iky++) {
+ ky = iky*dky;
+ ky2 = ky*ky;
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kz.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kz.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kz.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kz.i;
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kz.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kz.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kz.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kz.i;
+ }
+ }
+ for (iky = ikymax; iky <= nky-ikymax+1; iky++) {
+ for (ikx = 0; ikx <= nkx; ikx++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ }
+ for (iky = nky-ikymax+1; iky < nky; iky++) {
+ ky = (iky-nky)*dky;
+ ky2 = ky*ky;
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kz.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kz.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kz.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kz.i;
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kz.r = 0.0;
+ kz.i = sqrt(kz2);
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kz.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kz.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kz.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kz.i;
+ }
+ }
+
+ }
+ }
+ else if (opt < 0) {
+ for (iom = iomin ; iom < iomax ; iom++) {
+ kp = iom*deltom/c;
+ kp2 = kp*kp;
+ ikxmax = MIN((int)(kp/dkx), nkx/2);
+ ikymax = MIN((int)(kp/dky), nky/2);
+
+ for (iky = 0; iky < ikymax; iky++) {
+ ky = iky*dky;
+ ky2 = ky*ky;
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.i;
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.i;
+ }
+ }
+ for (iky = ikymax; iky <= nky-ikymax+1; iky++) {
+ for (ikx = 0; ikx <= nkx; ikx++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ }
+ for (iky = nky-ikymax+1; iky < nky; iky++) {
+ ky = (iky-nky)*dky;
+ ky2 = ky*ky;
+ for (ikx = 0; ikx < ikxmax; ikx++) {
+ kx = ikx*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.i;
+ }
+ for (ikx = ikxmax; ikx <= nkx-ikxmax+1; ikx++) {
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = 0.0;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = 0.0;
+ }
+ for (ikx = nkx-ikxmax+1; ikx < nkx; ikx++) {
+ kx = (ikx-nkx)*dkx;
+ kx2 = kx*kx;
+ kz2 = kp2 - kx2 - ky2;
+ kzinv.r = 0.0;
+ kzinv.i = -sqrt(kz2)/kz2;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].r = cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.r-cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.i;
+ cdatascl[iom*nky*nkx+iy*nkx+ix].i = cdata[iom*nky*nkx+iy*nkx+ix].i*kzinv.r+cdata[iom*nky*nkx+iy*nkx+ix].r*kzinv.i;
+ }
+ }
+
+ }
+ }
+ free(cdata);
+
+ /* inverse double FFT */
+ wkykx2yxt(&cdatascl[0], &rdata[0], optn, nkx, nky, optn, nkx, nky, 0, 0);
+ /* select original samples and traces */
+ scl = 1.0;
+ scl_data3D(rdata, optn, nx, ny, scl, data, nt, nx);
+
+ free(cdatascl);
+ free(rdata);
+
+ return;
+}
+
+void pad3d_data(float *data, long nt, long nx, long ny, long ntout, long nxout, long nyout, float *datout)
+{
+ int it,ix,iy;
+ for (iy=0;iy<ny;iy++) {
+ for (ix=0;ix<nx;ix++) {
+ for (it=0;it<nt;it++)
+ datout[iy*nx*nt+ix*nt+it]=data[iy*nx*nt+ix*nt+it];
+ for (it=nt;it<ntout;it++)
+ datout[iy*nx*nt+ix*nt+it]=0.0;
+ }
+ for (ix=nx;ix<nxout;ix++) {
+ for (it=0;it<ntout;it++)
+ datout[iy*nx*nt+ix*nt+it]=0.0;
+ }
+ }
+ for (iy=ny;iy<nyout;iy++) {
+ for (ix=0;ix<nxout;ix++) {
+ for (it=0;it<ntout;it++)
+ datout[iy*nx*nt+ix*nt+it]=0.0;
+ }
+ }
+}
+
+void scl_data3D(float *data, long nt, long nx, long ny, float scl, float *datout, long ntout, long nxout)
+{
+ int it,ix,iy;
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (it = 0 ; it < ntout ; it++) {
+ datout[iy*nxout*ntout+ix*ntout+it] = scl*data[iy*nx*nt+ix*nt+it];
+ }
+ }
+ }
+}
\ No newline at end of file
diff --git a/utils/truncate.c b/utils/truncate.c
new file mode 100644
index 0000000..7eb757a
--- /dev/null
+++ b/utils/truncate.c
@@ -0,0 +1,120 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+#ifndef MAX
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#endif
+#ifndef MIN
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#endif
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+#ifndef COMPLEX
+typedef struct _complexStruct { /* complex number */
+ float r,i;
+} complex;
+#endif/* complex */
+
+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);
+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);
+
+char *sdoc[] = {
+" ",
+" truncate - Truncate data below a certain depth",
+" ",
+" authors : Joeri Brackenhoff : (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke : (janth@xs4all.nl)",
+" ",
+" Required parameters: ",
+"",
+" file_in= ................. File containing the first data",
+" ",
+" Optional parameters: ",
+" ",
+" file_out=out.su .......... Filename of the output",
+" ztrunc=0.0 ............... Truncation depth",
+" verbose=1 ................ Give detailed information of process",
+NULL};
+
+int main (int argc, char **argv)
+{
+ FILE *fp_in, *fp_out;
+ char *fin, *fout;
+ float *indata, ztrunc, tmp;
+ float dz, dy, dx, z0, y0, x0, scl;
+ long izt, nt, nz, ny, nx, ntr, ix, iz;
+ long ret, verbose, nxyz;
+ segy *hdr_in;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ /*----------------------------------------------------------------------------*
+ * Get the parameters passed to the function
+ *----------------------------------------------------------------------------*/
+ if (!getparstring("file_in", &fin)) fin = NULL;
+ if (!getparstring("file_out", &fout)) fout = "out.su";
+ if (!getparfloat("ztrunc", &ztrunc)) ztrunc=0.0;
+ if (!getparlong("verbose", &verbose)) verbose=1;
+ if (fin == NULL) verr("No input file given");
+
+
+ /*----------------------------------------------------------------------------*
+ * Get the file info of the data and determine the indez of the truncation
+ *----------------------------------------------------------------------------*/
+
+ getFileInfo3D(fin, &nz, &nx, &ny, &nt, &dz, &dx, &dy, &z0, &x0, &y0, &scl, &ntr);
+
+ nxyz = nx*ny*nz;
+ izt = NINT((ztrunc-z0)/dz);
+
+ if (verbose) {
+ vmess("Number of samples: %li, x: %li, y: %li, z: %li",nxyz,nx,ny,nz);
+ vmess("Sampling distance for x: %.3f, y: %.3f, z: %.3f",dx,dy,dz);
+ vmess("Starting distance for x: %.3f, y: %.3f, z: %.3f",x0,y0,z0);
+ vmess("Final distance for x: %.3f, y: %.3f, z: %.3f",x0+((float)(nx-1))*dx,y0+((float)(ny-1))*dy,z0+((float)(nz-1))*dz);
+ vmess("Truncation depth is %f (iz=%li)",ztrunc,izt);
+ }
+ if (izt>nz) verr("The truncation depth (%f) is too large for the model (zmax=%f)",ztrunc,z0+((float)(nz-1))*dz);
+ if (izt<1) verr("The truncation depth (%f) is too small for the model (zmin=%f)",ztrunc,z0);
+
+ /*----------------------------------------------------------------------------*
+ * Allocate, read in and truncate the data
+ *----------------------------------------------------------------------------*/
+
+ indata = (float *)calloc(nx*ny*nz,sizeof(float));
+ hdr_in = (segy *)calloc(nx*ny,sizeof(segy));
+ readSnapData3D(fin, indata, hdr_in, 1, nx, ny, nz, 0, nx, 0, ny, 0, nz);
+ if (verbose) vmess("Read data");
+
+ for (ix = 0; ix < ny*nx; ix++) {
+ tmp = indata[ix*nz+izt];
+ for (iz = izt+1; iz < nz; iz++) {
+ indata[ix*nz+iz] = tmp;
+ }
+ }
+ if (verbose) vmess("Truncated data");
+
+ /*----------------------------------------------------------------------------*
+ * Write out the data
+ *----------------------------------------------------------------------------*/
+ fp_out = fopen(fout, "w+");
+
+ ret = writeData3D(fp_out, &indata[0], hdr_in, nz, nx*ny);
+ vmess("Wrote data");
+
+ fclose(fp_out);
+ free(indata); free(hdr_in);
+ return 0;
+}
+
--
GitLab