diff --git a/marchenko/demo/twoD/marchenko_ray.scr b/marchenko/demo/twoD/marchenko_ray.scr
index 5df653cb6e3016cd3e359aec9ab190fe7ec8eda3..31666a832fd784119efcbe100567edce0ec5638a 100755
--- a/marchenko/demo/twoD/marchenko_ray.scr
+++ b/marchenko/demo/twoD/marchenko_ray.scr
@@ -14,13 +14,13 @@ marchenko3D file_shot=shots/refl_rp.su verbose=1 \
file_green=greensrc.su
#Create the homogeneous Green's function
-marchenko3D file_shot=shots/refl_rp.su verbose=1 file_inp=greensrc.su \
+marchenko3D file_shot=shots/refl_rp.su verbose=1 file_inp=greensrc2.su \
tap=0 niter=6 hw=8 shift=7 smooth=3 geomspread=1 ampest=1 \
file_ray=rayz_time.su file_amp=rayz_amp.su file_wav=wavefp.su \
file_homg=homg_rayz.su scheme=0 file_green=green_rayz.su
#Create the homogeneous Green's function
-marchenko3D file_shot=shots/refl_rp.su verbose=1 file_inp=greensrc.su \
+marchenko3D file_shot=shots/refl_rp.su verbose=1 file_inp=greensrc2.su \
tap=0 niter=6 hw=8 shift=7 smooth=3 geomspread=1 ampest=1 \
file_ray=rayz_time.su file_amp=rayz_amp.su file_wav=wavefp.su \
file_homg=homg_class_rayz.su scheme=1 file_green=green_class_rayz.su
diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index 47499c1093f2841be9e74670db0e7110cace6da6..0695d33765687efed7a1d0332828eb4dc76e22a0 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -143,6 +143,8 @@ char *sdoc[] = {
" file_homg= ............... output file with homogeneous Green's function ",
" file_ampscl= ............. output file with estimated amplitudes ",
" file_iter= ............... output file with -Ni(-t) for each iteration",
+" compact=0 ................ Write out homg and imag in compact format",
+" .......................... WARNING! This write-out cannot be displayed with SU"
" verbose=0 ................ silent option; >0 displays info",
" ",
" ",
@@ -160,7 +162,7 @@ int main (int argc, char **argv)
long size, n1, n2, n3, ntap, tap, dxi, dyi, ntraces, pad, *sx, *sy, *sz;
long nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
long reci, countmin, mode, n2out, n3out, verbose, ntfft;
- long iter, niter, tracf, *muteW, *tsynW, ampest, plane_wave;
+ long iter, niter, tracf, *muteW, *tsynW, ampest, plane_wave, compact;
long hw, smooth, above, shift, *ixpos, *iypos, npos, ix, iy, nzim, nxim, nyim;
long nshots_r, *isxcount, *reci_xsrc, *reci_xrcv;
float fmin, fmax, *tapersh, *tapersy, fxf, fyf, dxf, dyf, *xsrc, *ysrc, *xrcv, *yrcv, *zsyn, *zsrc, *xrcvsyn, *yrcvsyn;
@@ -169,7 +171,7 @@ int main (int argc, char **argv)
float *green, *f2p, *pmin, *G_d, dt, dx, dy, dxs, dys, scl, mem;
float *f1plus, *f1min, *iRN, *Ni, *Nig, *trace, *Gmin, *Gplus, *HomG;
float xmin, xmax, ymin, ymax, scale, tsq, Q, f0, *tmpdata;
- float *ixmask, *iymask, *ampscl, *Gd, *Image, dzim, dyim, dxim;
+ float *ixmask, *iymask, *ampscl, *Gd, *Image, dzim;
float grad2rad, p, src_angle, src_velo;
complex *Refl, *Fop;
char *file_tinv, *file_shot, *file_green, *file_iter, *file_imag, *file_homg, *file_ampscl;
@@ -225,6 +227,7 @@ int main (int argc, char **argv)
if(!getparlong("smooth", &smooth)) smooth = 5;
if(!getparlong("above", &above)) above = 0;
if(!getparlong("shift", &shift)) shift=12;
+ if(!getparlong("compact", &compact)) compact=0;
if (!getparlong("plane_wave", &plane_wave)) plane_wave = 0;
if (!getparfloat("src_angle",&src_angle)) src_angle=0.;
@@ -551,6 +554,12 @@ int main (int argc, char **argv)
vmess("number of output traces = x:%li y:%li total:%li", n2out, n3out, n2out*n3out);
vmess("number of output samples = %li", ntfft);
vmess("Size of output data/file = %.1f MB", mem);
+ if (compact>0) {
+ vmess("Save format for homg and imag = compact");
+ }
+ else {
+ vmess("Save format for homg and imag = normal");
+ }
}
@@ -887,37 +896,65 @@ int main (int argc, char **argv)
imaging3D(Image, Gmin, f1plus, nxs, nys, ntfft, dxs, dys, dt, Nfoc, verbose);
if (file_gmin==NULL) free(Gmin);
- // Set headers
- hdrs_Nfoc = (segy *)calloc(nxim*nyim,sizeof(segy));
- for (l=0; l<nyim; l++){
- for (j=0; j<nxim; j++){
- hdrs_Nfoc[l*nxim+j].ns = nzim;
- hdrs_Nfoc[l*nxim+j].fldr = 1;
- hdrs_Nfoc[l*nxim+j].tracl = 1;
- hdrs_Nfoc[l*nxim+j].tracf = l*nxim+j+1;
- hdrs_Nfoc[l*nxim+j].trid = 2;
- hdrs_Nfoc[l*nxim+j].scalco = -1000;
- hdrs_Nfoc[l*nxim+j].scalel = -1000;
- hdrs_Nfoc[l*nxim+j].sx = xsyn[l*nxim+j]*(1e3);
- hdrs_Nfoc[l*nxim+j].sy = ysyn[l*nxim+j]*(1e3);
- hdrs_Nfoc[l*nxim+j].gx = xsyn[l*nxim+j]*(1e3);
- hdrs_Nfoc[l*nxim+j].gy = ysyn[l*nxim+j]*(1e3);
- hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l*nxim+j]*(1e3);
- hdrs_Nfoc[l*nxim+j].f1 = zsyn[0];
- hdrs_Nfoc[l*nxim+j].f2 = xsyn[0];
- hdrs_Nfoc[l*nxim+j].d1 = dzim;
- hdrs_Nfoc[l*nxim+j].d2 = dxs;
- hdrs_Nfoc[l*nxim+j].dt = (int)(hdrs_Nfoc[l*nxim+j].d1*(1E6));
- hdrs_Nfoc[l*nxim+j].trwf = nxim*nyim;
- hdrs_Nfoc[l*nxim+j].ntr = nxim*nyim;
+ // Set headers and write out image
+ fp_imag = fopen(file_imag, "w+");
+
+ if (compact > 0) {
+ hdrs_Nfoc = (segy *)calloc(1,sizeof(segy));
+
+ hdrs_Nfoc[0].ns = nzim*nyim*nxim;
+ hdrs_Nfoc[0].fldr = 1;
+ hdrs_Nfoc[0].tracr = nzim;
+ hdrs_Nfoc[0].tracl = nyim;
+ hdrs_Nfoc[0].tracf = nxim;
+ hdrs_Nfoc[0].trid = 2;
+ hdrs_Nfoc[0].scalco = -1000;
+ hdrs_Nfoc[0].scalel = -1000;
+ hdrs_Nfoc[0].sx = xsyn[0]*(1e3);
+ hdrs_Nfoc[0].sy = ysyn[0]*(1e3);
+ hdrs_Nfoc[0].sdepth = zsyn[0]*(1e3);
+ hdrs_Nfoc[0].f1 = zsyn[0];
+ hdrs_Nfoc[0].f2 = xsyn[0];
+ hdrs_Nfoc[0].ungpow = xsyn[0];
+ hdrs_Nfoc[0].d1 = dzim;
+ hdrs_Nfoc[0].d2 = dxs;
+ hdrs_Nfoc[0].unscale = dys;
+ hdrs_Nfoc[0].dt = (int)(dt*(1E6));
+
+ if (fp_imag==NULL) verr("error on creating output file %s", file_imag);
+ ret = writeData3D(fp_imag, (float *)&Image[0], hdrs_Nfoc, nzim*nyim*nxim, 1);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ else {
+ hdrs_Nfoc = (segy *)calloc(nxim*nyim,sizeof(segy));
+ for (l=0; l<nyim; l++){
+ for (j=0; j<nxim; j++){
+ hdrs_Nfoc[l*nxim+j].ns = nzim;
+ hdrs_Nfoc[l*nxim+j].fldr = 1;
+ hdrs_Nfoc[l*nxim+j].tracl = 1;
+ hdrs_Nfoc[l*nxim+j].tracf = l*nxim+j+1;
+ hdrs_Nfoc[l*nxim+j].trid = 2;
+ hdrs_Nfoc[l*nxim+j].scalco = -1000;
+ hdrs_Nfoc[l*nxim+j].scalel = -1000;
+ hdrs_Nfoc[l*nxim+j].sx = xsyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].sy = ysyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].gx = xsyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].gy = ysyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].sdepth = zsyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].f1 = zsyn[0];
+ hdrs_Nfoc[l*nxim+j].f2 = xsyn[0];
+ hdrs_Nfoc[l*nxim+j].d1 = dzim;
+ hdrs_Nfoc[l*nxim+j].d2 = dxs;
+ hdrs_Nfoc[l*nxim+j].dt = (int)(dt*(1E6));
+ hdrs_Nfoc[l*nxim+j].trwf = nxim*nyim;
+ hdrs_Nfoc[l*nxim+j].ntr = nxim*nyim;
+ }
}
+ if (fp_imag==NULL) verr("error on creating output file %s", file_imag);
+ ret = writeData3D(fp_imag, (float *)&Image[0], hdrs_Nfoc, nzim, nxim*nyim);
+ if (ret < 0 ) verr("error on writing output file.");
}
- // Write out image
- fp_imag = fopen(file_imag, "w+");
- if (fp_imag==NULL) verr("error on creating output file %s", file_imag);
- ret = writeData3D(fp_imag, (float *)&Image[0], hdrs_Nfoc, nzim, nxim*nyim);
- if (ret < 0 ) verr("error on writing output file.");
fclose(fp_imag);
free(hdrs_Nfoc);
free(Image);
@@ -931,45 +968,70 @@ int main (int argc, char **argv)
sy = (long *)calloc(nxs*nys,sizeof(long));
sz = (long *)calloc(nxs*nys,sizeof(long));
- // Determine Image
+ // Determine Homogeneous Green's function
HomG = (float *)calloc(Nfoc*ntfft,sizeof(float));
homogeneousg3D(HomG, green, f2p, zsyn, nxs, nys, ntfft, dxs, dys, dt, Nfoc, sx, sy, sz, verbose);
+ // Set headers and write out the data
fp_homg = fopen(file_homg, "w+");
if (fp_homg==NULL) verr("error on creating output file %s", file_homg);
- hdrs_Nfoc = (segy *)calloc(nxim*nyim,sizeof(segy));
- for (i=0; i<ntfft; i++) {
- // Set headers
- for (l=0; l<nyim; l++){
- for (j=0; j<nxim; j++){
- hdrs_Nfoc[l*nxim+j].ns = nzim;
- hdrs_Nfoc[l*nxim+j].fldr = i+1;
- hdrs_Nfoc[l*nxim+j].tracl = 1;
- hdrs_Nfoc[l*nxim+j].tracf = l*nxim+j+1;
- hdrs_Nfoc[l*nxim+j].trid = 2;
- hdrs_Nfoc[l*nxim+j].scalco = -1000;
- hdrs_Nfoc[l*nxim+j].scalel = -1000;
- hdrs_Nfoc[l*nxim+j].sx = sx[l*nxim+j];
- hdrs_Nfoc[l*nxim+j].sy = sy[l*nxim+j];
- hdrs_Nfoc[l*nxim+j].gx = xsyn[l*nxim+j]*(1e3);
- hdrs_Nfoc[l*nxim+j].gy = ysyn[l*nxim+j]*(1e3);
- hdrs_Nfoc[l*nxim+j].sdepth = sz[l*nxim+j];
- hdrs_Nfoc[l*nxim+j].selev = -sz[l*nxim+j];
- hdrs_Nfoc[l*nxim+j].f1 = zsyn[0];
- hdrs_Nfoc[l*nxim+j].f2 = xsyn[0];
- hdrs_Nfoc[l*nxim+j].d1 = dzim;
- hdrs_Nfoc[l*nxim+j].d2 = dxs;
- hdrs_Nfoc[l*nxim+j].dt = (int)(hdrs_Nfoc[l*nxim+j].d1*(1E6));
- hdrs_Nfoc[l*nxim+j].trwf = nxim*nyim;
- hdrs_Nfoc[l*nxim+j].ntr = nxim*nyim;
- }
- }
+ if (compact > 0) {
+ hdrs_Nfoc = (segy *)calloc(nxim*nyim,sizeof(segy));
- // Write out homogeneous Green's function
- ret = writeData3D(fp_homg, (float *)&HomG[i*Nfoc], hdrs_Nfoc, nzim, nxim*nyim);
+ hdrs_Nfoc[0].ns = nzim*nyim*nxim;
+ hdrs_Nfoc[0].fldr = ntfft;
+ hdrs_Nfoc[0].tracr = nzim;
+ hdrs_Nfoc[0].tracl = nyim;
+ hdrs_Nfoc[0].tracf = nxim;
+ hdrs_Nfoc[0].trid = 2;
+ hdrs_Nfoc[0].scalco = -1000;
+ hdrs_Nfoc[0].scalel = -1000;
+ hdrs_Nfoc[0].sx = xsyn[0]*(1e3);
+ hdrs_Nfoc[0].sy = ysyn[0]*(1e3);
+ hdrs_Nfoc[0].sdepth = zsyn[0]*(1e3);
+ hdrs_Nfoc[0].f1 = zsyn[0];
+ hdrs_Nfoc[0].f2 = xsyn[0];
+ hdrs_Nfoc[0].ungpow = xsyn[0];
+ hdrs_Nfoc[0].d1 = dzim;
+ hdrs_Nfoc[0].d2 = dxs;
+ hdrs_Nfoc[0].unscale = dys;
+ hdrs_Nfoc[0].dt = (int)(dt*(1E6));
+
+ ret = writeData3D(fp_homg, (float *)&HomG[0], hdrs_Nfoc, nzim*nyim*nxim*ntfft, 1);
if (ret < 0 ) verr("error on writing output file.");
}
+ else {
+ hdrs_Nfoc = (segy *)calloc(nxim*nyim,sizeof(segy));
+ for (i=0; i<ntfft; i++) {
+ for (l=0; l<nyim; l++){
+ for (j=0; j<nxim; j++){
+ hdrs_Nfoc[l*nxim+j].ns = nzim;
+ hdrs_Nfoc[l*nxim+j].fldr = i+1;
+ hdrs_Nfoc[l*nxim+j].tracl = 1;
+ hdrs_Nfoc[l*nxim+j].tracf = l*nxim+j+1;
+ hdrs_Nfoc[l*nxim+j].trid = 2;
+ hdrs_Nfoc[l*nxim+j].scalco = -1000;
+ hdrs_Nfoc[l*nxim+j].scalel = -1000;
+ hdrs_Nfoc[l*nxim+j].sx = sx[l*nxim+j];
+ hdrs_Nfoc[l*nxim+j].sy = sy[l*nxim+j];
+ hdrs_Nfoc[l*nxim+j].gx = xsyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].gy = ysyn[l*nxim+j]*(1e3);
+ hdrs_Nfoc[l*nxim+j].sdepth = sz[l*nxim+j];
+ hdrs_Nfoc[l*nxim+j].selev = -sz[l*nxim+j];
+ hdrs_Nfoc[l*nxim+j].f1 = zsyn[0];
+ hdrs_Nfoc[l*nxim+j].f2 = xsyn[0];
+ hdrs_Nfoc[l*nxim+j].d1 = dzim;
+ hdrs_Nfoc[l*nxim+j].d2 = dxs;
+ hdrs_Nfoc[l*nxim+j].dt = (int)(hdrs_Nfoc[l*nxim+j].d1*(1E6));
+ hdrs_Nfoc[l*nxim+j].trwf = nxim*nyim;
+ hdrs_Nfoc[l*nxim+j].ntr = nxim*nyim;
+ }
+ }
+ ret = writeData3D(fp_homg, (float *)&HomG[i*Nfoc], hdrs_Nfoc, nzim, nxim*nyim);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ }
fclose(fp_homg);
free(hdrs_Nfoc);
diff --git a/utils/MuteSnap.c b/utils/MuteSnap.c
deleted file mode 100644
index b90681ec6965d7a2ff110c80c1f771ee1aee29fb..0000000000000000000000000000000000000000
--- a/utils/MuteSnap.c
+++ /dev/null
@@ -1,227 +0,0 @@
-#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) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
-
-#ifndef COMPLEX
-typedef struct _complexStruct { /* complex number */
- float r,i;
-} complex;
-#endif/* complex */
-
-int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
-int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
-double wallclock_time(void);
-int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
-int readSnapData(char *filename, float *data, segy *hdr, int ngath, int nx, int ntfft, int sx, int ex, int sz, int ez);
-int topdet(float *data, int nt);
-int farrdet(float *data, int nt, float tol);
-
-char *sdoc[] = {
-" ",
-" HomG - Calculate a Homogeneous Green's function ",
-" ",
-" authors : Joeri Brackenhoff (J.A.Brackenhoff@tudelft.nl)",
-" : Jan Thorbecke : (janth@xs4all.nl)",
-" ",
-" Required parameters: ",
-"",
-" fhom= .................... File containing the snapshot data that will be muted",
-" fsnap= ................... File containing the snapshot data that will determine the mute window",
-" ",
-" Optional parameters: ",
-" ",
-" fout=out.su .............. Filename of the output",
-" shift=5 .................. Shift from the maximum",
-" smooth=5 ................. Length of smoothing taper",
-" mode=0 ................... Determine first arrival by maximum (mode=0), first event above tol (mode=1) or by raytime (mode=2)",
-" tol=1 .................... Tolerance for the determination of first arrival if mode=1",
-" fray ..................... File containing the raytimes of the first arrivals",
-NULL};
-
-int main (int argc, char **argv)
-{
- FILE *fp_snap, *fp_hom, *fp_out;
- char *fhom, *fsnap, *fout, *fray;
- float *homdata, *snapdata, *outdata, *rtrace, *costaper, scl, tol, *timeval, dt;
- float dx, dz, z0, x0, xmin, xmax, sclsxgx, f1, f2, dxrcv, dzrcv, dxpos, offset;
- int nt, nts, nx, nxs, nxh, ntraces, ret, ix, it, is, ir, nzs, nzh, nz, ht, indrcv, shift;
- int rmt, smooth, mode, nzh1, nzs1, nxh1, nxs1, nts1, nt1;
- segy *hdr_hom, *hdr_snap, *hdrs_mute;
-
- initargs(argc, argv);
- requestdoc(1);
-
- if (!getparstring("fhom", &fhom)) fhom = NULL;
- if (!getparstring("fsnap", &fsnap)) fsnap = NULL;
- if (!getparstring("fout", &fout)) fout = "out.su";
- if (!getparstring("fray", &fray)) fray = NULL;
- if (!getparint("shift", &shift)) shift = 5;
- if (!getparint("smooth", &smooth)) smooth = 5;
- if (!getparint("mode", &mode)) mode = 0;
- if (!getparfloat("tol", &tol)) tol = 5;
- if (fhom == NULL) verr("Incorrect G_hom input");
- if (fsnap == NULL) verr("Incorrect snapshot input");
- if (mode == 2) {
- if (fray == NULL) verr("No filename for raytimes given");
- }
- if (!getparint("nxs1", &nxs1)) nxs1 = 0;
- if (!getparint("nxh1", &nxh1)) nxh1 = 0;
- if (!getparint("nzs1", &nzs1)) nzs1 = 0;
- if (!getparint("nzh1", &nzh1)) nzh1 = 0;
- if (!getparint("nts1", &nts1)) nts1 = 0;
- if (!getparint("nt1", &nt1)) nt1 = 0;
-
- if (smooth) {
- costaper = (float *)malloc(smooth*sizeof(float));
- scl = M_PI/((float)smooth);
- for (is=0; is<smooth; is++) {
- costaper[is] = 0.5*(1.0+cos((is+1)*scl));
- }
- }
-
- getFileInfo(fsnap, &nzs, &nxs, &nts, &dz, &dx, &z0, &x0, &xmin, &xmax, &sclsxgx, &ntraces);
- getFileInfo(fhom, &nzh, &nxh, &nt, &dz, &dx, &z0, &x0, &xmin, &xmax, &sclsxgx, &ntraces);
-
- if (nxh1 != 0) nxh = nxh1;
- if (nxs1 != 0) nxs = nxs1;
- if (nzh1 != 0) nzh = nzh1;
- if (nzs1 != 0) nzs = nzs1;
- if (nt1 != 0) nt = nt1;
- if (nts1 != 0) nts = nts1;
-
- if (nzs != nzh || nxs != nxh) {
- verr("sampling in x or z direction is incorrect, nzs=%d nzh=%d, nxs=%d nxh=%d",nzs,nzh,nxs,nxh);
- }
-
- vmess("nzs=%d nzh=%d, nxs=%d nxh=%d, nts=%d nt=%d",nzs,nzh,nxs,nxh,nts,nt);
-
- nz = nzh;
- nx = nxh;
-
- snapdata = (float *)malloc(nz*nx*nts*sizeof(float));
- hdr_snap = (segy *)calloc(nx*nts,sizeof(segy));
- homdata = (float *)malloc(nz*nx*nt*sizeof(float));
- hdr_hom = (segy *)calloc(nx*nt,sizeof(segy));
- ht = (int)ceil(nt/2);
- rtrace = (float *)malloc(nts*sizeof(float));
-
- if (mode != 2) {
- readSnapData(fsnap, &snapdata[0], &hdr_snap[0], nts, nx, nz, 0, nx, 0, nz);
- vmess("Read Snapshot data");
- }
- readSnapData(fhom, &homdata[0], &hdr_hom[0], nt, nx, nz, 0, nx, 0, nz);
- vmess("Read G_hom");
-
- if (mode == 0) {
- vmess("First arrival determined through maximum");
- }
- else if (mode == 1) {
- vmess("First arrival determined through tolerance (=%.4f)",tol);
- }
- else if (mode == 2) {
- vmess("First arrival determined through raytimes");
- fp_snap = fopen(fray,"r");
- if (fp_snap == NULL) {
- verr("Could not open file");
- }
- fclose(fp_snap);
- hdrs_mute = (segy *) calloc(nz,sizeof(segy));
- timeval = (float *)calloc(nz*nx,sizeof(float));
- readSnapData(fray, timeval, hdrs_mute, nz, 1, nx, 0, 1, 0, nx);
- dt = hdr_hom[0].dt/1E6;
- }
-
- for (ir = 0; ir < nz; ir++) {
- for (is = 0; is < nx; is++) {
- for (it = 0; it < nts; it++) {
- rtrace[it] = snapdata[it*nxs*nzs+is*nzs+ir];
- }
- if (mode == 0) {
- indrcv = topdet(&rtrace[0],nts);
- }
- else if (mode == 1) {
- indrcv = farrdet(&rtrace[0],nts,tol);
- }
- else if (mode == 2) {
- indrcv = (int)roundf(timeval[ir*nx+is]/dt);
- }
- rmt = MIN(nt-indrcv,indrcv)-shift;
- for (it = ht-rmt+1; it < ht; it++) {
- if (it-(ht-rmt) < smooth) {
- homdata[it*nxs*nzs+is*nzs+ir] *= costaper[it-(ht-rmt)];
- }
- else {
- homdata[it*nxs*nzs+is*nzs+ir] = 0.0;
- }
- }
- for (it = ht; it < ht+rmt; it++) {
- if (it-(ht+rmt)+smooth > 0) {
- homdata[it*nxs*nzs+is*nzs+ir] *= costaper[smooth-(it-(ht+rmt)+smooth)];
- }
- else {
- homdata[it*nxs*nzs+is*nzs+ir] = 0.0;
- }
- }
- }
- vmess("Muting Homogeneous Green's function at depth %d from %d depths",ir+1,nzs);
- }
- free(snapdata);free(hdr_snap);
-
- fp_out = fopen(fout, "w+");
-
- for (ir = 0; ir < nt; ir++) {
- ret = writeData(fp_out, &homdata[ir*nxs*nzs], &hdr_hom[ir*nx], nz, nx);
- if (ret < 0 ) verr("error on writing output file.");
- }
-
- fclose(fp_out);
- vmess("Wrote Data");
- return 0;
-}
-
-int topdet(float *data, int nt)
-{
- int it,ind;
- float maxval;
-
- maxval = data[0];
- ind = 0;
-
- for (it = 1; it < nt; it++) {
- if (fabs(data[it]) > maxval) {
- maxval = data[it];
- ind = it;
- }
- }
-
- return ind;
-}
-
-int farrdet(float *data, int nt, float tol)
-{
- int it,ind;
-
- ind = 0;
-
- for (it = 0; it < nt; it++) {
- if (fabs(data[it]) > tol) {
- ind = it;
- break;
- }
- }
-
- return ind;
-}
diff --git a/utils/combine.c b/utils/combine.c
index c08792ec069405fce004eaf74471396895c411c6..55f5211c1ad1f4066597bd9bc246bb37173855d3 100755
--- a/utils/combine.c
+++ b/utils/combine.c
@@ -27,6 +27,8 @@ 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 readCompact(char *filename, long *nx, long *ny, long *nz, long *nt, float *dx, float *dy, float *dz, long *dt,
+ float *fx, float *fy, float *fz, long *sx, long *sy, long *sz, float *scl);
char *sdoc[] = {
" ",
@@ -41,6 +43,7 @@ char *sdoc[] = {
" ",
" Optional parameters: ",
" ",
+" compact=0 ................ Save format of input data (0=normal), (1=transpose), (2=compact)",
" file_out= ................ Filename of the output",
" numb= .................... integer number of first file",
" dnumb= ................... integer number of increment in files",
@@ -51,10 +54,10 @@ 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];
- float *indata, *outdata, dt;
+ float *indata, *outdata;
float dz, dy, dx, z0, y0, x0, scl;
- long nt, nz, ny, nx, ntr, ix, iy, it, is, iz, pos, file_det, nzs;
- long numb, dnumb, ret, nzmax, transpose, verbose, nxyz, sx, sy, sz;
+ 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;
segy *hdr_in, *hdr_out;
initargs(argc, argv);
@@ -69,7 +72,7 @@ int main (int argc, char **argv)
if (!getparlong("dnumb", &dnumb)) dnumb=0;
if (!getparlong("nzmax", &nzmax)) nzmax=0;
if (!getparlong("verbose", &verbose)) verbose=0;
- if (!getparlong("transpose", &transpose)) transpose=0;
+ if (!getparlong("compact", &compact)) compact=0;
if (fin == NULL) verr("Incorrect downgoing input");
/*----------------------------------------------------------------------------*
@@ -98,6 +101,8 @@ int main (int argc, char **argv)
}
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);
@@ -121,12 +126,18 @@ int main (int argc, char **argv)
sprintf(fins,"%li",numb);
sprintf(fin2,"%s%s%s",fbegin,fins,fend);
nt = 1;
- if (transpose==0) {
+ if (compact==0) {
+ if (verbose) vmess("Save fomat is normal");
getFileInfo3D(fin2, &nz, &nx, &ny, &nt, &dz, &dx, &dy, &z0, &x0, &y0, &scl, &ntr);
}
- else {
+ else if (compact==1) {
+ if (verbose) vmess("Save fomat is transpose");
getFileInfo3D(fin2, &nx, &nz, &ny, &nt, &dz, &dx, &dy, &z0, &x0, &y0, &scl, &ntr);
}
+ else if (compact==2) {
+ if (verbose) vmess("Save fomat is compact");
+ readCompact(fin2, &nx, &ny, &nz, &nt, &dx, &dy, &dz, &dt, &x0, &y0, &z0, &sx, &sy, &sz, &scl);
+ }
nxyz = nx*ny*nzs*nz;
@@ -143,8 +154,13 @@ int main (int argc, char **argv)
*----------------------------------------------------------------------------*/
hdr_out = (segy *)calloc(nx*ny,sizeof(segy));
outdata = (float *)calloc(nxyz*nt,sizeof(float));
- hdr_in = (segy *)calloc(nx*ny*nt,sizeof(segy));
indata = (float *)calloc(nx*ny*nz*nt,sizeof(float));
+ if (compact != 2) {
+ hdr_in = (segy *)calloc(nx*ny*nt,sizeof(segy));
+ }
+ else {
+ hdr_in = (segy *)calloc(1,sizeof(segy));
+ }
/*----------------------------------------------------------------------------*
* Combine the separate files
@@ -158,7 +174,7 @@ int main (int argc, char **argv)
verr("Error opening file");
}
fclose(fp_in);
- if (transpose==0) {
+ if (compact==0) {
readSnapData3D(fin2, indata, hdr_in, nt, nx, ny, nz, 0, nx, 0, ny, 0, nz);
for (it = 0; it < nt; it++) {
for (iy = 0; iy < ny; iy++) {
@@ -170,7 +186,7 @@ int main (int argc, char **argv)
}
}
}
- else {
+ else if (compact==1) {
readSnapData3D(fin2, indata, hdr_in, nt, nz, ny, nx, 0, nz, 0, ny, 0, nx);
for (it = 0; it < nt; it++) {
for (iy = 0; iy < ny; iy++) {
@@ -182,12 +198,26 @@ int main (int argc, char **argv)
}
}
}
+ else if (compact==2) {
+ readSnapData3D(fin2, indata, hdr_in, 1, 1, 1, nx*ny*nz*nt, 0, 1, 0, 1, 0, nx*ny*nz*nt);
+ for (it = 0; it < nt; it++) {
+ for (iy = 0; iy < ny; iy++) {
+ for (ix = 0; ix < nx; ix++) {
+ for (iz = 0; iz < nz; iz++) {
+ outdata[it*ny*nx*nz*nzs+iy*nx*nz*nzs+ix*nz*nzs+iz+is] = indata[it*ny*nx*nz+iy*nx*nz+ix*nz+iz];
+ }
+ }
+ }
+ }
+ }
}
free(indata);
- sx = hdr_in[0].sx;
- sy = hdr_in[0].sy;
- sz = hdr_in[0].sdepth;
- dt = hdr_in[0].dt;
+ if (compact != 2) {
+ sx = hdr_in[0].sx;
+ sy = hdr_in[0].sy;
+ sz = hdr_in[0].sdepth;
+ dt = hdr_in[0].dt;
+ }
free(hdr_in);
/*----------------------------------------------------------------------------*
@@ -228,4 +258,47 @@ int main (int argc, char **argv)
free(outdata); free(hdr_out);
vmess("Wrote data");
return 0;
+}
+
+void readCompact(char *filename, long *nx, long *ny, long *nz, long *nt, float *dx, float *dy, float *dz, long *dt,
+ float *fx, float *fy, float *fz, long *sx, long *sy, long *sz, float *scl)
+{
+ FILE *fp;
+ segy hdr;
+ size_t nread;
+
+ 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");
+
+ *nx = hdr.tracf;
+ *ny = hdr.tracl;
+ *nz = hdr.tracr;
+ *nt = hdr.fldr;
+
+ *dx = hdr.d2;
+ *dy = hdr.unscale;
+ *dz = hdr.d1;
+ *dt = hdr.dt;
+
+ *fx = hdr.f2;
+ *fy = hdr.ungpow;
+ *fz = hdr.f1;
+
+ *sx = hdr.sx;
+ *sy = hdr.sy;
+ *sz = hdr.sdepth;
+
+ if (hdr.scalco > 0) {
+ *scl = ((float)hdr.scalco);
+ }
+ else if (hdr.scalco < 0) {
+ *scl = (-1.0/((float)hdr.scalco));
+ }
+ else {
+ *scl = 1.0;
+ }
+
+ return;
}
\ No newline at end of file