diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index 4ed77c7e5d14e15ef267d351cdb5312deddfbbd5..3eb355ef66a09b17e5748c5a99b64fde15225bb6 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -878,17 +878,6 @@ int main (int argc, char **argv)
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];
- }
- }
- }
timeShift(Gmin, nts, npos*Nfoc, dt, tshift, fmin, fmax);
}
else {
diff --git a/raytime3d/raytime3d.c b/raytime3d/raytime3d.c
index a689661c756702dbe87911eb696e6f78aed07060..801b4dfd73450291e5ab2a22c9a2048fd486cfb5 100644
--- a/raytime3d/raytime3d.c
+++ b/raytime3d/raytime3d.c
@@ -324,7 +324,7 @@ private (is,time0,ampl,nrx,nry,nrz,nr,cp_average,i,j,k,ix,iy,iz,hdrs,tmpdata,nwr
if (ret < 0 ) verr("error on writing output file.");
if (ray.geomspread==1) {
- ret = writeData3D(fpt, &l[0], hdrs, rec.nx, rec.nz*rec.ny);
+ ret = writeData3D(fpa, &l[0], hdrs, rec.nx, rec.nz*rec.ny);
if (ret < 0 ) verr("error on writing output file.");
}
diff --git a/utils/HomG.c b/utils/HomG.c
index 52cf1864bc34b440f49e4c60e3b92129a79d9a98..87712f6e181f442a5d6f9733ee3ca3aaef55839e 100755
--- a/utils/HomG.c
+++ b/utils/HomG.c
@@ -65,6 +65,8 @@ void scl_data3D(float *data, long nt, long nx, long ny, float scl, float *datout
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 deconv(float *data1, float *data2, float *decon, long nrec, long nsam,
+ float dt, float eps, float reps, 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);
@@ -97,6 +99,9 @@ char *sdoc[] = {
" zfpr=0 ................... virtual receiver data are in SU format (=0) or zfp compressed (=1)",
" cp=1000.0 ................ Velocity at the top of the medium in m/s",
" rho=1000.0 ............... Density at the top of the medium in kg/m^3",
+" file_wav= ................ Wavelet that is deconvolved from the virtual receiver data",
+" eps=0.0 .................. Absolute stabilization factor for deconvolution",
+" reps=0.01 ................ Relative stabilization factor for deconvolution (is multiplied by the maximum of the data)",
" scheme=0 ................. Scheme for the retrieval",
" .......................... scheme=0 Marchenko homogeneous Green's function retrieval with G source",
" .......................... scheme=1 Marchenko homogeneous Green's function retrieval with f2 source",
@@ -113,17 +118,17 @@ NULL};
int main (int argc, char **argv)
{
- FILE *fp_in, *fp_shot, *fp_out;
- char *fin, *fshot, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100], *direction;
+ FILE *fp_in, *fp_shot, *fp_out, *fp_wav;
+ char *fin, *fshot, *fout, *ptr, fbegin[100], fend[100], fins[100], fin2[100], *direction, *file_wav;
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;
+ float *conv, *conv2, *decon, *tmp1, *tmp2, cp, shift, eps, reps;
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, num, isn;
- float dtr, dxr, dyr, ftr, fxr, fyr, sclr, scl;
+ float dtr, dxr, dyr, ftr, fxr, fyr, sclr, scl, *wavelet;
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;
+ segy *hdr_rcv, *hdr_out, *hdr_shot, *hdr_wav;
initargs(argc, argv);
requestdoc(1);
@@ -144,6 +149,9 @@ int main (int argc, char **argv)
if (!getparlong("dnumb", &dnumb)) dnumb=1;
if (!getparlong("scheme", &scheme)) scheme = 0;
if (!getparlong("verbose", &verbose)) verbose = 0;
+ if (!getparstring("file_wav", &file_wav)) file_wav = NULL;
+ if(!getparfloat("eps", &eps)) eps=0.00;
+ if(!getparfloat("reps", &reps)) reps=0.01;
if (!getparlong("zfps", &zfps)) zfps = 0;
if (!getparlong("zfpr", &zfpr)) zfpr = 0;
if (!getparstring("direction", &direction)) direction = "z";
@@ -281,6 +289,24 @@ int main (int argc, char **argv)
yvr = (long *)malloc(nxvr*nyvr*nzvr*sizeof(long));
zvr = (long *)malloc(nxvr*nyvr*nzvr*sizeof(long));
+ /*----------------------------------------------------------------------------*
+ * Read in the wavelet and create the gather
+ *----------------------------------------------------------------------------*/
+
+ if (file_wav) {
+ if (verbose) vmess("Reading in wavelet for deconvolution");
+ wavelet = (float *)calloc(ntvs*nxvs*nyvs,sizeof(float));
+ hdr_wav = (segy *)calloc(nxvs*nyvs,sizeof(segy));
+ readSnapData3D(file_wav, &wavelet[0], &hdr_wav[0], 1, 1, 1, ntvs, 0, 1, 0, 1, 0, ntvs);
+ for (i = 1; i < nxvs*nyvs; i++){
+ for (j = 0; j < ntvs; j++) {
+ wavelet[i*ntvs+j] = wavelet[j];
+ hdr_wav[i] = hdr_wav[0];
+ }
+ }
+ free(hdr_wav);
+ }
+
/*----------------------------------------------------------------------------*
* Get the file info for the source position
*----------------------------------------------------------------------------*/
@@ -371,6 +397,7 @@ int main (int argc, char **argv)
tmp2 = (float *)calloc(nyr*nxr*ntr,sizeof(float));
}
if (scheme==8 || scheme==9 || scheme==10) tmp1 = (float *)calloc(nyr*nxr*ntr,sizeof(float));
+ if (file_wav!=NULL) decon = (float *)calloc(nyr*nxr*ntr,sizeof(float));
sprintf(fins,"%s%li",direction,ir*dnumb+numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
@@ -439,7 +466,13 @@ int main (int argc, char **argv)
else if (scheme==3) { //Marchenko representation without time-reversal G source
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);
+ if (file_wav!=NULL){
+ deconv(&rcvdata[l*nyr*nxr*ntr], wavelet, decon, nyr*nxr, ntr, dt, eps, reps, 0);
+ convol2(&shotdata[0], decon, conv, nyr*nxr, ntr, dt, fmin, fmax, 1);
+ }
+ else {
+ 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;
@@ -542,9 +575,11 @@ int main (int argc, char **argv)
free(tmp1); free(tmp2);
}
if (scheme==8 || scheme==9 || scheme==10) free(tmp1);
+ if (file_wav!=NULL) free(decon);
}
free(shotdata);
+ if (file_wav!=NULL) free(wavelet);
if (strcmp(direction,"z") == 0) {
if (nxvr>1) dxrcv = (float)((xvr[nzvr] - xvr[0])/1000.0);
@@ -1203,7 +1238,7 @@ void getVirRec(char *filename, long *nxs, long *nys, long *nxr, long *nyr, long
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");
+ if (nread != TRCBYTES) verr("Could not read the header of the VR file");
*nxs = 1;
*nys = 1;
@@ -1655,4 +1690,105 @@ void scl_data3D(float *data, long nt, long nx, long ny, float scl, float *datout
}
}
}
+}
+
+void deconv(float *data1, float *data2, float *decon, long nrec, long nsam,
+ float dt, float eps, float reps, long shift)
+{
+ long i, j, n, optn, nfreq, sign;
+ float df, dw, om, tau, *den, scl;
+ float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1, *rdata2, maxden, leps;
+ complex *cdata1, *cdata2, *cdec, 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");
+ cdec = (complex *)malloc(nfreq*nrec*sizeof(complex));
+ if (cdec == 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");
+ den = (float *)malloc(nfreq*nrec*sizeof(float));
+ if (den == NULL) verr("memory allocation error for rdata1");
+
+ /* 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 deconvolution */
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ n = nrec*nfreq;
+ maxden=0.0;
+ for (j = 0; j < n; j++) {
+ den[j] = *p2r**p2r + *p2i**p2i;
+ maxden = MAX(den[j], maxden);
+ p2r += 2;
+ p2i += 2;
+ }
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ qr = (float *) &cdec[0].r;
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ qi = qr + 1;
+ leps = reps*maxden+eps;
+ for (j = 0; j < n; j++) {
+
+ if (fabs(*p2r)>=fabs(*p2i)) {
+ *qr = (*p2r**p1r+*p2i**p1i)/(den[j]+leps);
+ *qi = (*p2r**p1i-*p2i**p1r)/(den[j]+leps);
+ } else {
+ *qr = (*p1r**p2r+*p1i**p2i)/(den[j]+leps);
+ *qi = (*p1i**p2r-*p1r**p2i)/(den[j]+leps);
+ }
+ qr += 2;
+ qi += 2;
+ p1r += 2;
+ p1i += 2;
+ p2r += 2;
+ p2i += 2;
+ }
+ free(cdata1);
+ free(cdata2);
+ free(den);
+
+ 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 = cdec[j*nfreq+i].r*cos(om*tau) + cdec[j*nfreq+i].i*sin(om*tau);
+ tmp.i = cdec[j*nfreq+i].i*cos(om*tau) - cdec[j*nfreq+i].r*sin(om*tau);
+ cdec[j*nfreq+i] = tmp;
+ om += dw;
+ }
+ }
+ }
+
+ /* inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = dt/(float)optn;
+ crmfft(&cdec[0], &rdata1[0], optn, nrec, nfreq, optn, sign);
+ scl_data(rdata1,optn,nrec,scl,decon,nsam);
+
+ free(cdec);
+ free(rdata1);
+ free(rdata2);
+ return;
}
\ No newline at end of file
diff --git a/utils/combine_induced.c b/utils/combine_induced.c
index b1c63a8679c6c1bbb2a847fd255df9cd1ab9c43a..3a6229258614f9b9bcc43993520f6c878012843c 100755
--- a/utils/combine_induced.c
+++ b/utils/combine_induced.c
@@ -48,6 +48,7 @@ char *sdoc[] = {
" nshift=0 ................. Sample number shift",
" shift=0.0 ................ Base shift in seconds for the data",
" dtshift=0.0 .............. Shift per gather in seconds",
+" opt=0 .................... Option for muting acausal events (=0) or no muting (=1)",
" verbose=1 ................ Give information about the process (=1) or not (=0)",
NULL};
@@ -58,7 +59,7 @@ int main (int argc, char **argv)
float *indata, *outdata, *loopdata, shift, dtshift;
float dt, dz, dy, dx, t0, x0, y0, z0, scl, dxrcv, dyrcv, dzrcv;
long nt, nz, ny, nx, nxyz, ntr, ix, iy, it, is, iz, pos, file_det, nxs, nys, nzs;
- long numb, dnumb, ret, nzmax, verbose, nt_out, ishift, nshift, *sx, *sy, *sz;
+ long numb, dnumb, ret, nzmax, verbose, nt_out, ishift, nshift, *sx, *sy, *sz, opt;
segy *hdr_in, *hdr_loop, *hdr_out;
initargs(argc, argv);
@@ -74,6 +75,7 @@ int main (int argc, char **argv)
if (!getparlong("nzmax", &nzmax)) nzmax=0;
if (!getparlong("verbose", &verbose)) verbose=0;
if (!getparlong("nshift", &nshift)) nshift=0;
+ if (!getparlong("opt", &opt)) opt=0;
if (!getparfloat("shift", &shift)) shift=0.0;
if (!getparfloat("dtshift", &dtshift)) dtshift=0.0;
if (fin == NULL) verr("Incorrect downgoing input");
@@ -137,6 +139,8 @@ int main (int argc, char **argv)
vmess("Starting distance for x: %.3f, y: %.3f, z: %.3f",x0,y0,z0);
vmess("Sampling distance for x: %.3f, y: %.3f, z: %.3f",dx,dy,dz);
vmess("Number of virtual sources: %li",nzs);
+ if (opt==0) vmess("Muting is applied");
+ if (opt==1) vmess("Muting is not applied");
}
/*----------------------------------------------------------------------------*
@@ -165,6 +169,8 @@ int main (int argc, char **argv)
hdr_out = (segy *)calloc(nx*ny,sizeof(segy));
outdata = (float *)calloc(nt_out*nxyz,sizeof(float));
+ if (dtshift==0.0) dtshift=((float) nshift)*dt;
+
/*----------------------------------------------------------------------------*
* Parallel loop for reading in and combining the various shots
*----------------------------------------------------------------------------*/
@@ -189,12 +195,25 @@ int main (int argc, char **argv)
sz[is] = hdr_loop[0].sdepth;
ishift = nshift*is;
- if (verbose) vmess("Shifting %li timesteps for a total of %.3f seconds",ishift,shift+(dtshift*((float)iz)));
- for (it = ishift; it < nt_out; it++) {
- for (iy = 0; iy < ny; iy++) {
- for (ix = 0; ix < nx; ix++) {
- for (iz = 0; iz < nz; iz++) {
- outdata[it*nxyz+iy*nx*nz+ix*nz+iz] += loopdata[(it-ishift+(nt/2))*nxyz+iy*nx*nz+ix*nz+iz];
+ if (verbose) vmess("Shifting %li timesteps for a total of %.3f seconds",ishift,shift+(dtshift*((float)is)));
+ if (opt==0) {
+ for (it = ishift; it < nt_out; it++) {
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (iz = 0; iz < nz; iz++) {
+ outdata[it*nxyz+iy*nx*nz+ix*nz+iz] += loopdata[(it-ishift+(nt/2))*nxyz+iy*nx*nz+ix*nz+iz];
+ }
+ }
+ }
+ }
+ }
+ else {
+ for (it = 0; it < nt_out; it++) {
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (iz = 0; iz < nz; iz++) {
+ outdata[it*nxyz+iy*nx*nz+ix*nz+iz] += loopdata[(it-ishift+(nt/2))*nxyz+iy*nx*nz+ix*nz+iz];
+ }
}
}
}
diff --git a/utils/convhomg.c b/utils/convhomg.c
index 2bd57624d4b95dd4f68053bbaca5d862dfeb2cc7..e328f0ec666108871913735e2854ce70b4aef74f 100644
--- a/utils/convhomg.c
+++ b/utils/convhomg.c
@@ -31,6 +31,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 scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout);
void convol(float *data1, float *data2, float *con, long ntfft);
+void deconv(float *data1, float *data2, float *con, long ntfft, float eps, float reps);
char *sdoc[] = {
" ",
@@ -41,22 +42,25 @@ char *sdoc[] = {
" ",
" Required parameters: ",
"",
-" file_hom= ................. First file of the array of virtual receivers",
-" file_wav= ................. File containing the virtual source",
+" file_hom= ................. File containing the homogeneous Green's function",
+" file_wavcon= .............. File containing the wavelet that will be used for convolution",
" ",
" Optional parameters: ",
" ",
+" file_wavdecon= ........... File containing the wavelet that will be used for deconvolution",
+" eps=0.0 .................. Absolute stabilization factor for deconvolution",
+" reps=0.01 ................ Relative stabilization factor for deconvolution (is multiplied by the maximum of the data)",
" file_out= ................ Filename of the output",
NULL};
int main (int argc, char **argv)
{
FILE *fp_hom, *fp_wav, *fp_out;
- char *file_hom, *file_wav, *file_out;
+ char *file_hom, *file_wavcon, *file_wavdecon, *file_out;
long nt, nx, ny, nz, ntwav, ntfft, ntr, verbose;
long nt_wav, it, ipos;
float dt, dx, dy, dz, x0, y0, z0, scl, *Ghom, *wav, *wavelet, dt_wav;
- float *trace, *conv;
+ float *trace, *conv, *decon, *tmp, eps, reps;
segy *hdr_hom, hdr_wav;
size_t nread;
@@ -67,12 +71,15 @@ int main (int argc, char **argv)
* Get the parameters passed to the function
*----------------------------------------------------------------------------*/
if (!getparstring("file_hom", &file_hom)) file_hom = NULL;
- if (!getparstring("file_wav", &file_wav)) file_wav = NULL;
+ if (!getparstring("file_wavcon", &file_wavcon)) file_wavcon = NULL;
+ if (!getparstring("file_wavdecon", &file_wavdecon)) file_wavdecon = NULL;
if (!getparstring("file_out", &file_out)) file_out = "out.su";
+ if (!getparfloat("eps", &eps)) eps = 0.00;
+ if (!getparfloat("reps", &reps)) reps = 0.01;
if (!getparlong("verbose", &verbose)) verbose = 1;
if (file_hom==NULL) verr("Error file_hom is not given");
- if (file_wav==NULL) verr("Error file_wav is not given");
+ if (file_wavcon==NULL) verr("Error file_wav is not given");
/*----------------------------------------------------------------------------*
* Get the file info of the data and determine the indez of the truncation
@@ -91,8 +98,8 @@ int main (int argc, char **argv)
* Allocate and read in the data
*----------------------------------------------------------------------------*/
- fp_wav = fopen( file_wav, "r" );
- if (fp_wav == NULL) verr("File %s does not exist or cannot be opened", file_wav);
+ fp_wav = fopen( file_wavcon, "r" );
+ if (fp_wav == NULL) verr("File %s does not exist or cannot be opened", file_wavcon);
nread = fread( &hdr_wav, 1, TRCBYTES, fp_wav );
assert(nread == TRCBYTES);
nt_wav = hdr_wav.ns;
@@ -121,18 +128,51 @@ int main (int argc, char **argv)
wavelet = (float *)calloc(ntfft,sizeof(float));
trace = (float *)calloc(ntfft,sizeof(float));
+ tmp = (float *)calloc(ntfft,sizeof(float));
conv = (float *)calloc(ntfft,sizeof(float));
for (it = 0; it < nt_wav/2; it++) {
wavelet[it] = wav[it];
wavelet[ntfft-1-it] = wav[nt_wav-1-it];
}
- free(wav);
+
+ if (file_wavdecon!=NULL) {
+ if (verbose) vmess("Reading in deconvolution wavelet");
+ fp_wav = fopen( file_wavdecon, "r" );
+ if (fp_wav == NULL) verr("File %s does not exist or cannot be opened", file_wavdecon);
+ nread = fread( &hdr_wav, 1, TRCBYTES, fp_wav );
+ assert(nread == TRCBYTES);
+ nt_wav = hdr_wav.ns;
+ dt_wav = (float)(hdr_wav.dt/1E6);
+ wav = (float *)calloc(nt_wav,sizeof(float));
+ nread = fread(&wav[0], sizeof(float), nt_wav, fp_wav);
+ assert(nread==nt_wav);
+ fclose(fp_wav);
+
+ decon = (float *)calloc(ntfft,sizeof(float));
+
+ for (it = 0; it < nt_wav/2; it++) {
+ decon[it] = wav[it];
+ decon[ntfft-1-it] = wav[nt_wav-1-it];
+ }
+ vmess("Time sampling wavelet decon : %.3f",dt_wav);
+ if (dt_wav!=dt) vmess("WARNING! dt of HomG (%f) and deconvolution wavelet (%f) do not match.",dt,dt_wav);
+ }
+ free(wav);
+
for (ipos = 0; ipos<nx*ny*nz; ipos++) {
- for (it = 0; it < nt; it++) {
- trace[it] = Ghom[it*nx*ny*nz+ipos];
- }
- convol(trace,wavelet,conv,ntfft);
+ if (file_wavdecon!=NULL) {
+ for (it = 0; it < nt; it++) {
+ tmp[it] = Ghom[it*nx*ny*nz+ipos];
+ }
+ deconv(tmp,decon,trace,ntfft,eps,reps);
+ }
+ else{
+ for (it = 0; it < nt; it++) {
+ trace[it] = Ghom[it*nx*ny*nz+ipos];
+ }
+ }
+ convol(trace,wavelet,conv,ntfft);
for (it = 0; it < nt; it++) {
Ghom[it*nx*ny*nz+ipos] = conv[it];
}
@@ -143,6 +183,7 @@ int main (int argc, char **argv)
}
free(wavelet); free(trace); free(conv);
+ if (file_wavdecon!=NULL) free(decon);
fp_out = fopen(file_out, "w+");
@@ -158,7 +199,7 @@ int main (int argc, char **argv)
void convol(float *data1, float *data2, float *con, long ntfft)
{
- long i, j, n, nfreq, sign;
+ long i, j, nfreq, sign;
float scl;
float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1;
complex *cdata1, *cdata2, *ccon;
@@ -187,8 +228,7 @@ void convol(float *data1, float *data2, float *con, long ntfft)
p1i = p1r + 1;
p2i = p2r + 1;
qi = qr + 1;
- n = nfreq;
- for (j = 0; j < n; j++) {
+ for (j = 0; j < nfreq; j++) {
*qr = (*p2r**p1r-*p2i**p1i);
*qi = (*p2r**p1i+*p2i**p1r);
qr += 2;
@@ -212,6 +252,81 @@ void convol(float *data1, float *data2, float *con, long ntfft)
return;
}
+void deconv(float *data1, float *data2, float *con, long ntfft, float eps, float reps)
+{
+ long i, j, nfreq, sign;
+ float scl, maxden, *den, leps;
+ float *qr, *qi, *p1r, *p1i, *p2r, *p2i, *rdata1;
+ complex *cdata1, *cdata2, *ccon;
+
+ nfreq = ntfft/2+1;
+
+ den = (float *)malloc(nfreq*sizeof(float));
+
+ cdata1 = (complex *)malloc(nfreq*sizeof(complex));
+ if (cdata1 == NULL) verr("memory allocation error for cdata1");
+ cdata2 = (complex *)malloc(nfreq*sizeof(complex));
+ if (cdata2 == NULL) verr("memory allocation error for cdata2");
+ ccon = (complex *)malloc(nfreq*sizeof(complex));
+ if (ccon == NULL) verr("memory allocation error for ccov");
+
+ rdata1 = (float *)malloc(ntfft*sizeof(float));
+ if (rdata1 == NULL) verr("memory allocation error for rdata1");
+
+ /* forward time-frequency FFT */
+ sign = -1;
+ rcmfft(&data1[0], &cdata1[0], (int)ntfft, 1, (int)ntfft, (int)nfreq, (int)sign);
+ rcmfft(&data2[0], &cdata2[0], (int)ntfft, 1, (int)ntfft, (int)nfreq, (int)sign);
+
+ /* apply deconvolution */
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ maxden=0.0;
+ for (j = 0; j < nfreq; j++) {
+ den[j] = *p2r**p2r + *p2i**p2i;
+ maxden = MAX(den[j], maxden);
+ p2r += 2;
+ p2i += 2;
+ }
+ p1r = (float *) &cdata1[0];
+ p2r = (float *) &cdata2[0];
+ qr = (float *) &ccon[0].r;
+ p1i = p1r + 1;
+ p2i = p2r + 1;
+ qi = qr + 1;
+ leps = reps*maxden+eps;
+ for (j = 0; j < nfreq; j++) {
+
+ if (fabs(*p2r)>=fabs(*p2i)) {
+ *qr = (*p2r**p1r+*p2i**p1i)/(den[j]+leps);
+ *qi = (*p2r**p1i-*p2i**p1r)/(den[j]+leps);
+ } else {
+ *qr = (*p1r**p2r+*p1i**p2i)/(den[j]+leps);
+ *qi = (*p1i**p2r-*p1r**p2i)/(den[j]+leps);
+ }
+ qr += 2;
+ qi += 2;
+ p1r += 2;
+ p1i += 2;
+ p2r += 2;
+ p2i += 2;
+ }
+ free(cdata1);
+ free(cdata2);
+
+ /* inverse frequency-time FFT and scale result */
+ sign = 1;
+ scl = 1.0/((float)(ntfft));
+ crmfft(&ccon[0], &rdata1[0], (int)ntfft, 1, (int)nfreq, (int)ntfft, (int)sign);
+ scl_data(rdata1,ntfft,1,scl,con,ntfft);
+
+ free(ccon);
+ free(rdata1);
+ return;
+}
+
void scl_data(float *data, long nsam, long nrec, float scl, float *datout, long nsamout)
{
long it,ix;
diff --git a/utils/mutesnap.c b/utils/mutesnap.c
index 7afa393e19ddef9bd7f8c2354e8b6d742fff77b7..6102c1751490f45b362f9f7e615df4891f298242 100644
--- a/utils/mutesnap.c
+++ b/utils/mutesnap.c
@@ -61,7 +61,7 @@ int main (int argc, char **argv)
float *homdata, *snapdata, *outdata, *rtrace, *costaper, scl, tol, *timeval, dt;
float dxs, dys, dzs, scls, fzs, fxs, fys;
float dxh, dyh, dzh, sclh, fzh, fxh, fyh;
- float dxrcv, dyrcv, dzrcv, dxpos, offset;
+ float dxrcv, dyrcv, dzrcv, dxpos, offset, maxz;
long nts, nxs, nys, nzs, ntrs, nth, nxh, nyh, nzh, ntrh, opt;
long nxyz, nxy, ret, ix, iy, iz, it, ht, indrcv, shift, rmt, mode, smooth, verbose;
segy *hdr_hom, *hdr_snap, *hdrs_mute;
@@ -82,6 +82,7 @@ int main (int argc, char **argv)
if (!getparlong("verbose", &verbose)) verbose = 1;
if (!getparlong("opt", &opt)) opt = 0;
if (!getparfloat("tol", &tol)) tol = 5;
+ if (!getparfloat("maxz", &maxz)) maxz = 0;
if (fhom == NULL) verr("Incorrect G_hom input");
if (mode != 2) {
if (fsnap == NULL) verr("Incorrect snapshot input");
@@ -104,7 +105,7 @@ int main (int argc, char **argv)
/*----------------------------------------------------------------------------*
* Determine the parameters of the files and determine if they match
*----------------------------------------------------------------------------*/
- getFileInfo3D(fhom, &nzh, &nxh, &nyh, &nth, &dzh, &dxh, &dyh, &fzh, &fxh, &fyh, &sclh, &ntrh);
+ getFileInfo3D(fhom, &nzh, &nyh, &nxh, &nth, &dzh, &dxh, &dyh, &fzh, &fxh, &fyh, &sclh, &ntrh);
if (mode != 2) {
getFileInfo3D(fsnap, &nzs, &nxs, &nys, &nts, &dzs, &dxs, &dys, &fzs, &fxs, &fys, &scls, &ntrs);
@@ -122,6 +123,7 @@ int main (int argc, char **argv)
if (verbose) {
vmess("Number of virtual receivers: %li, nz: %li, nx: %li, ny: %li",nxyz,nzh,nxh,nyh);
vmess("Sampling distance in direction of z: %.3f, x: %.3f, y: %.3f",dzh,dxh,dyh);
+ vmess("Starting distance in direction of z: %.3f, x: %.3f, y: %.3f",fzh,fxh,fyh);
vmess("Number of time samples: %li",nth);
}
@@ -159,7 +161,7 @@ int main (int argc, char **argv)
fclose(fp_snap);
hdrs_mute = (segy *) calloc(nzh*nyh,sizeof(segy));
timeval = (float *)calloc(nxyz,sizeof(float));
- readSnapData3D(fray, timeval, hdrs_mute, nzh, 1, nyh, nxh, 0, 1, 0, nyh, 0, nxh);
+ readSnapData3D(fray, timeval, hdrs_mute, 1, nzh, nyh, nxh, 0, nzh, 0, nyh, 0, nxh);
dt = hdr_hom[0].dt/1E6;
}