diff --git a/fdacrtmc/extractMigrationSnapshots.c b/fdacrtmc/extractMigrationSnapshots.c
index 796682146a9df9ab3e819cb0c4fc0e7dbff1430c..18cf7436fd75c5f95a705dd7829f58248d18dc1c 100644
--- a/fdacrtmc/extractMigrationSnapshots.c
+++ b/fdacrtmc/extractMigrationSnapshots.c
@@ -60,7 +60,6 @@ int storeVerticalParticleVelocitySnapshot(modPar *mod, wavPar *wav, migPar *mig)
else farr=(float*)malloc(mig->sizem*sizeof(float));
// Store Vertical Particle Velocity Field - Interpolate to Tzz(P)-Grid
- mig->wav[mig->it].vz=(float*)malloc(mig->sizem*sizeof(float));
for(ix=0,ix1=mod->ioZx+mig->x1;ix<mig->nx;ix++,ix1+=mig->skipdx){
for(iz=0,iz1=mod->ioZz+mig->z1;iz<mig->nz;iz++,iz1+=mig->skipdz){
farr[ix*mig->nz+iz]=(wav->vz[ix1*mod->naz+iz1]+wav->vz[ix1*mod->naz+iz1+1])/2.0;
@@ -77,9 +76,6 @@ int storeVerticalParticleVelocitySnapshot(modPar *mod, wavPar *wav, migPar *mig)
int extractMigrationSnapshots(modPar *mod, wavPar *wav, migPar *mig, decompPar *decomp){
size_t ix, ix1, iz, iz1;
-//if(mig->it<100){mig->it++;return(0);}
-
-
switch(mig->mode){
case 1: /* Conventional Migration - Store Pressure Field (Tzz) */
// P (Tzz)
@@ -90,8 +86,10 @@ int extractMigrationSnapshots(modPar *mod, wavPar *wav, migPar *mig, decompPar *
storePressureSnapshot(mod,wav,mig);
switch(mig->orient){ //Migration Orientation
case 0: //Image Wavefields not travelling in the same direction
+ // Vertical Particle Velocity - Interpolate to Tzz(P)-Grid
+ storeVerticalParticleVelocitySnapshot(mod,wav,mig);
// Horizontal Particle Velocity - Interpolate to Tzz(P)-Grid
- storeVerticalParticleVelocitySnapshot(mod,wav,mig); //TODO: What does this case actually do?
+ storeHorizontalParticleVelocitySnapshot(mod,wav,mig); //TODO: What does this case actually do?
break;
case 1: //Up-Down Imaging
// Vertical Particle Velocity - Interpolate to Tzz(P)-Grid
diff --git a/fdacrtmc/rtmImagingCondition.c b/fdacrtmc/rtmImagingCondition.c
index 01bca46ec6b692263f2afde5468dfed571ea5668..2e07f4c529dc1e4c990729a73b7cf6498c88b42e 100644
--- a/fdacrtmc/rtmImagingCondition.c
+++ b/fdacrtmc/rtmImagingCondition.c
@@ -66,6 +66,8 @@ int rtmImagingCondition(modPar *mod, wavPar *wav, migPar *mig, decompPar *decomp
}
}
}
+ free(mig->wav[mig->it].vx);
+ free(mig->wav[mig->it].vz);
break;
case 1: // Up-Down Imaging
for(ix=0,ix1=mod->ioPx+mig->x1;ix<mig->nx;ix++,ix1+=mig->skipdx){
@@ -78,6 +80,7 @@ int rtmImagingCondition(modPar *mod, wavPar *wav, migPar *mig, decompPar *decomp
if(pzf*pzb>0) mig->image[ix*mig->nz+iz]+=mig->dt*(wav->tzz[ix1*mod->naz+iz1]*mig->wav[mig->it].tzz[ix*mig->nz+iz]);
}
}
+ free(mig->wav[mig->it].vz);
break;
case 2: // Left-Right Imaging
for(ix=0,ix1=mod->ioPx+mig->x1;ix<mig->nx;ix++,ix1+=mig->skipdx){
@@ -90,6 +93,7 @@ int rtmImagingCondition(modPar *mod, wavPar *wav, migPar *mig, decompPar *decomp
if(pxf*pxb<0) mig->image[ix*mig->nz+iz]+=mig->dt*(wav->tzz[ix1*mod->naz+iz1]*mig->wav[mig->it].tzz[ix*mig->nz+iz]);
}
}
+ free(mig->wav[mig->it].vx);
break;
case 3: //Normal Imaging
//Not yet implemented
@@ -103,8 +107,6 @@ int rtmImagingCondition(modPar *mod, wavPar *wav, migPar *mig, decompPar *decomp
break;
}
free(mig->wav[mig->it].tzz);
- free(mig->wav[mig->it].vx);
- free(mig->wav[mig->it].vz);
break;
case 3: /* Wavefield decomposition zero-lag pressure cross-correlation */
mig->it--;break;
diff --git a/marchenko/Makefile b/marchenko/Makefile
index 7fadd2cca106c94b3ae012eb6aaa78128c749097..848b6a8624741950e794ba9e25da1e43e5ab8213 100644
--- a/marchenko/Makefile
+++ b/marchenko/Makefile
@@ -4,9 +4,10 @@ include ../Make_include
#LIBS += -L$L -lgenfft -lm $(LIBSM)
#OPTC += -g -O0 -Wall
+#OPTC += -g -traceback #-check-pointers=rw -check-pointers-undimensioned
-ALL: fmute marchenko #marchenko_primaries
+ALL: fmute marchenko marchenko_primaries
SRCJ = fmute.c \
getFileInfo.c \
@@ -35,20 +36,20 @@ SRCH = marchenko.c \
docpkge.c \
getpars.c
-#SRCP = marchenko_primaries.c \
-# synthesis.c \
-# getFileInfo.c \
-# readData.c \
-# readShotData.c \
-# readTinvData.c \
-# writeData.c \
-# writeDataIter.c \
-# wallclock_time.c \
-# name_ext.c \
-# verbosepkg.c \
-# atopkge.c \
-# docpkge.c \
-# getpars.c
+SRCP = marchenko_primaries.c \
+ synthesis.c \
+ getFileInfo.c \
+ readData.c \
+ readShotData.c \
+ readTinvData.c \
+ writeData.c \
+ writeDataIter.c \
+ wallclock_time.c \
+ name_ext.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ getpars.c
OBJJ = $(SRCJ:%.c=%.o)
diff --git a/marchenko/demo/oneD/primaries.scr b/marchenko/demo/oneD/primaries.scr
index f5a82f2ebb581b3ad54f6bea1ec0e3e7f0550389..0308f62f8fb23872b18928d039a224dc6316fc55 100755
--- a/marchenko/demo/oneD/primaries.scr
+++ b/marchenko/demo/oneD/primaries.scr
@@ -8,16 +8,18 @@
#cd $SLURM_SUBMIT_DIR
export PATH=$HOME/src/OpenSource/bin:$PATH:
-export OMP_NUM_THREADS=20
+export OMP_NUM_THREADS=40
#shot record to remove internal multiples
-select=451
+select=51
makewave fp=20 dt=0.004 file_out=wave.su nt=2048 t0=0.0 scale=0 scfft=1
../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
- nshots=601 verbose=2 istart=40 iend=500 fmax=90 pad=1024 \
- niter=25 smooth=10 niterskip=600 shift=20 file_rr=pred_rr.su T=0
+ nshots=901 verbose=2 istart=40 iend=500 fmax=90 pad=1024 \
+ niter=20 smooth=10 niterskip=1 shift=20 file_rr=pred_rr.su T=0
+
+exit;
#for reference original shot record from Reflection matrix
suwind key=fldr min=$select max=$select < shotsdx5_rp.su itmax=2047 > shotsx0.su
diff --git a/marchenko/marchenko_primaries.c b/marchenko/marchenko_primaries.c
new file mode 100644
index 0000000000000000000000000000000000000000..ee5c2750d484863ca9d4b1e63727c2dd7d6f003c
--- /dev/null
+++ b/marchenko/marchenko_primaries.c
@@ -0,0 +1,759 @@
+#include "par.h"
+#include "segy.h"
+#include <time.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <assert.h>
+#include <genfft.h>
+
+int omp_get_max_threads(void);
+int omp_get_num_threads(void);
+void omp_set_num_threads(int num_threads);
+
+#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 readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int nshots, int nx, int nxs, float fxsb, float dxs, int ntfft, int mode, float scale, float tsq, float Q, float f0, int reci, int *nshots_r, int *isxcount, int *reci_xsrc, int *reci_xrcv, float *ixmask, int verbose);
+
+int readTinvData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nfoc, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
+
+int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int t0shift, int iter);
+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 readData(FILE *fp, float *data, segy *hdrs, int n1);
+int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
+int disp_fileinfo(char *file, int n1, int n2, float f1, float f2, float d1, float d2, segy *hdrs);
+double wallclock_time(void);
+
+void synthesis(complex *Refl, complex *Fop, float *Top, float *RNi, int nx, int nt, int nxs, int nts, float dt, float *xsyn, int
+Nfoc, float *xrcv, float *xsrc, int *xnx, float fxse, float fxsb, float dxs, float dxsrc, float dx, int ntfft, int
+nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int npos, double *tfft, int *isxcount, int
+*reci_xsrc, int *reci_xrcv, float *ixmask, int verbose);
+
+void synthesisPositions(int nx, int nt, int nxs, int nts, float dt, float *xsyn, int Nfoc, float *xrcv, float *xsrc, int *xnx, float fxse, float fxsb, float dxs, float dxsrc, float dx, int nshots, int *ixpos, int *npos, int *isxcount, int countmin, int reci, int verbose);
+
+
+
+/*********************** self documentation **********************/
+char *sdoc[] = {
+" ",
+" MARCHENKO_primaries - Iterative primary reflections retrieval",
+" ",
+" marchenko_primaries file_tinv= file_shot= [optional parameters]",
+" ",
+" Required parameters: ",
+" ",
+" file_tinv= ............... shot-record from R to remove internal multiples",
+" file_shot= ............... Reflection response: R",
+" ",
+" Optional parameters: ",
+" ",
+" INTEGRATION ",
+" ishot=nshots/2 ........... shot number(s) to remove internal multiples ",
+" file_src=spike ........... convolve ishot(s) with source wavelet",
+" file_tinv= ............... use file_tinv to remove internal multiples",
+" COMPUTATION",
+" tap=0 .................... lateral taper focusing(1), shot(2) or both(3)",
+" ntap=0 ................... number of taper points at boundaries",
+" fmin=0 ................... minimum frequency in the Fourier transform",
+" fmax=70 .................. maximum frequency in the Fourier transform",
+" file_src= ................ optional source wavelet to convolve selected ishot's",
+" MARCHENKO ITERATIONS ",
+" niter=22 ................. number of iterations to inialize and restart",
+" niterec=2 ................ number of iterations in recursive part of the time-samples",
+" niterskip=50 ............. restart scheme each niterskip time-samples with niter iterations",
+" istart=20 ................ start sample of iterations for primaries",
+" iend=nt .................. end sample of iterations for primaries",
+" MUTE-WINDOW ",
+" shift=20 ................. number of points to account for wavelet (epsilon in papers)",
+" smooth=shift/2 ........... number of points to smooth mute with cosine window",
+" REFLECTION RESPONSE CORRECTION ",
+" tsq=0.0 .................. scale factor n for t^n for true amplitude recovery",
+" Q=0.0 .......,............ Q correction factor",
+" f0=0.0 ................... ... for Q correction factor",
+" scale=2 .................. scale factor of R for summation of Ni with G_d",
+" pad=0 .................... amount of samples to pad the reflection series",
+" reci=0 ................... 1; add receivers as shots 2; only use receivers as shot positions",
+" countmin=0 ............... 0.3*nxrcv; minumum number of reciprocal traces for a contribution",
+" OUTPUT DEFINITION ",
+" file_rr= ................. output file with primary only shot record",
+" file_iter= ............... output file with -Ni(-t) for each iteration",
+" T=0 ...................... :1 compute transmission-losses compensated primaries ",
+" verbose=0 ................ silent option; >0 displays info",
+" ",
+" ",
+" author : Lele Zhang & Jan Thorbecke : 2019 ",
+" ",
+NULL};
+/**************** end self doc ***********************************/
+
+int main (int argc, char **argv)
+{
+ FILE *fp_out, *fp_rr, *fp_w;
+ size_t nread;
+ int i, j, k, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath;
+ int size, n1, n2, ntap, tap, di, ntraces;
+ int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
+ int reci, countmin, mode, ixa, ixb, n2out, verbose, ntfft;
+ int iter, niter, niterec, recur, niterskip, niterrun, tracf, *muteW;
+ int hw, ii, ishot, istart, iend;
+ int smooth, *ixpos, npos, ix, m, pad, T, perc;
+ int nshots_r, *isxcount, *reci_xsrc, *reci_xrcv, shift;
+ float fmin, fmax, *tapersh, *tapersy, fxf, dxf, *xsrc, *xrcv, *zsyn, *zsrc, *xrcvsyn;
+ double t0, t1, t2, t3, t4, tsyn, tread, tfft, tcopy, tii;
+ float d1, d2, f1, f2, fxsb, fxse, ft, fx, *xsyn, dxsrc;
+ float *G_d, *DD, *RR, dt, dx, dxs, scl, mem;
+ float *rtrace, *tmpdata, *f1min, *f1plus, *RNi, *Ni, *trace;
+ float xmin, xmax, scale, tsq;
+ float Q, f0, *ixmask, *costaper;
+ complex *Refl, *Fop, *ctrace, *cwave;
+ char *file_tinv, *file_shot, *file_rr, *file_src, *file_iter;
+ segy *hdrs_out, hdr;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ tsyn = tread = tfft = tcopy = tii = 0.0;
+ t0 = wallclock_time();
+
+ if (!getparstring("file_shot", &file_shot)) file_shot = NULL;
+ if (!getparstring("file_tinv", &file_tinv)) file_tinv = NULL;
+ if(!getparstring("file_src", &file_src)) file_src = NULL;
+ if (!getparstring("file_rr", &file_rr)) verr("parameter file_rr not found");
+ if (!getparstring("file_iter", &file_iter)) file_iter = NULL;
+
+ if (!getparint("verbose", &verbose)) verbose = 0;
+ if (!getparfloat("fmin", &fmin)) fmin = 0.0;
+ if (!getparfloat("fmax", &fmax)) fmax = 70.0;
+ if (!getparint("reci", &reci)) reci = 0;
+ reci=0; // source-receiver reciprocity is not yet fully build into the code
+ if (!getparfloat("scale", &scale)) scale = 2.0;
+ if (!getparfloat("Q", &Q)) Q = 0.0;
+ if (!getparfloat("tsq", &tsq)) tsq = 0.0;
+ if (!getparfloat("f0", &f0)) f0 = 0.0;
+ if (!getparint("tap", &tap)) tap = 0;
+ if (!getparint("ntap", &ntap)) ntap = 0;
+ if (!getparint("pad", &pad)) pad = 0;
+ if (!getparint("T", &T)) T = 0;
+ if (T>0) T=-1;
+ else T=1;
+
+
+ if(!getparint("niter", &niter)) niter = 22;
+ if(!getparint("niterec", &niterec)) niterec = 2;
+ if(!getparint("niterskip", &niterskip)) niterskip = 50;
+ if(!getparint("hw", &hw)) hw = 15;
+ if(!getparint("shift", &shift)) shift=20;
+ if(!getparint("smooth", &smooth)) smooth = shift/2;
+ if(!getparint("ishot", &ishot)) ishot=300;
+ if (reci && ntap) vwarn("tapering influences the reciprocal result");
+
+ smooth = MIN(smooth, shift);
+ 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));
+ }
+ }
+
+/*================ Reading info about shot and initial operator sizes ================*/
+
+ if (file_tinv != NULL) { /* G_d is read from file_tinv */
+ ngath = 0; /* setting ngath=0 scans all traces; n2 contains maximum traces/gather */
+ ret = getFileInfo(file_tinv, &n1, &n2, &ngath, &d1, &d2, &f1, &f2, &xmin, &xmax, &scl, &ntraces);
+ Nfoc = ngath;
+ nxs = n2;
+ nts = n1;
+ dxs = d2;
+ fxsb = f2;
+ }
+
+ ngath = 0; /* setting ngath=0 scans all traces; nx contains maximum traces/gather */
+ ret = getFileInfo(file_shot, &nt, &nx, &ngath, &d1, &dx, &ft, &fx, &xmin, &xmax, &scl, &ntraces);
+ nshots = ngath;
+
+ if (!getparfloat("dt", &dt)) dt = d1;
+ if(!getparint("istart", &istart)) istart=20;
+ if(!getparint("iend", &iend)) iend=nt;
+ iend = MIN(iend, nt-shift-1);
+ istart = MIN(MAX(1,istart),iend);
+
+ if (file_tinv == NULL) {/* 'G_d' is one of the shot records */
+ if(!getparint("ishot", &ishot)) ishot=1+(nshots-1)/2;
+ ishot -= 1; /* shot numbering starts at 0 */
+ Nfoc = 1;
+ nxs = nx;
+ nts = nt;
+ dxs = dx;
+ fxsb = fx;
+ }
+ assert (nxs >= nshots);
+
+ ntfft = optncr(MAX(nt+pad, nts+pad));
+ nfreq = ntfft/2+1;
+ nw_low = (int)MIN((fmin*ntfft*dt), nfreq-1);
+ nw_low = MAX(nw_low, 1);
+ nw_high = MIN((int)(fmax*ntfft*dt), nfreq-1);
+ nw = nw_high - nw_low + 1;
+ scl = 1.0/((float)ntfft);
+ if (!getparint("countmin", &countmin)) countmin = 0.3*nx;
+
+/*================ Allocating all data arrays ================*/
+
+ Fop = (complex *)calloc(nxs*nw*Nfoc,sizeof(complex));
+ f1min = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ f1plus = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ RNi = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ Ni = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ G_d = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ DD = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ RR = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ trace = (float *)malloc(ntfft*sizeof(float));
+ ixpos = (int *)malloc(nxs*sizeof(int));
+ xrcvsyn = (float *)calloc(Nfoc*nxs,sizeof(float));
+ xsyn = (float *)malloc(Nfoc*sizeof(float));
+ zsyn = (float *)malloc(Nfoc*sizeof(float));
+ xnxsyn = (int *)calloc(Nfoc,sizeof(int));
+
+ Refl = (complex *)malloc(nw*nx*nshots*sizeof(complex));
+ xsrc = (float *)calloc(nshots,sizeof(float));
+ zsrc = (float *)calloc(nshots,sizeof(float));
+ xrcv = (float *)calloc(nshots*nx,sizeof(float));
+ xnx = (int *)calloc(nshots,sizeof(int));
+
+ if (reci!=0) {
+ reci_xsrc = (int *)malloc((nxs*nxs)*sizeof(int));
+ reci_xrcv = (int *)malloc((nxs*nxs)*sizeof(int));
+ isxcount = (int *)calloc(nxs,sizeof(int));
+ ixmask = (float *)calloc(nxs,sizeof(float));
+ }
+
+/*================ Read focusing operator(s) ================*/
+
+ if (file_tinv != NULL) { /* G_d is named DD */
+ muteW = (int *)calloc(Nfoc*nxs,sizeof(int));
+ mode=-1; /* apply complex conjugate to read in data */
+ readTinvData(file_tinv, xrcvsyn, xsyn, zsyn, xnxsyn, Nfoc, nxs, ntfft,
+ mode, muteW, DD, hw, verbose);
+ /* reading data added zero's to the number of time samples to be the same as ntfft */
+ nts = ntfft;
+ free(muteW);
+
+ /* check consistency of header values */
+ if (xrcvsyn[0] != 0 || xrcvsyn[1] != 0 ) fxsb = xrcvsyn[0];
+ fxse = fxsb + (float)(nxs-1)*dxs;
+ dxf = (xrcvsyn[nxs-1] - xrcvsyn[0])/(float)(nxs-1);
+ if (NINT(dxs*1e3) != NINT(fabs(dxf)*1e3)) {
+ vmess("dx in hdr.d1 (%.3f) and hdr.gx (%.3f) not equal",d2, dxf);
+ if (dxf != 0) dxs = fabs(dxf);
+ vmess("dx in operator => %f", dxs);
+ }
+ }
+
+/* ========================= Opening optional wavelet file ====================== */
+
+ cwave = (complex *)calloc(ntfft,sizeof(complex));
+ if (file_src != NULL){
+ if (verbose) vmess("Reading wavelet from file %s.", file_src);
+
+ fp_w = fopen(file_src, "r");
+ if (fp_w == NULL) verr("error on opening input file_src=%s", file_src);
+ nread = fread(&hdr, 1, TRCBYTES, fp_w);
+ assert (nread == TRCBYTES);
+ tmpdata = (float *)malloc(hdr.ns*sizeof(float));
+ nread = fread(tmpdata, sizeof(float), hdr.ns, fp_w);
+ assert (nread == hdr.ns);
+ fclose(fp_w);
+
+ /* Check dt wavelet is the same as reflection data */
+ if ( rint(dt*1000) != rint(hdr.dt/1000.0) ) {
+ vwarn("file_src dt != dt of file_shot %e != %e", hdr.dt/1e6, dt);
+ }
+ rtrace = (float *)calloc(ntfft,sizeof(float));
+
+ /* add zero's to the number of time samples to be the same as ntfft */
+ /* Add zero-valued samples in middle */
+ if (hdr.ns <= ntfft) {
+ for (i = 0; i < hdr.ns/2; i++) {
+ rtrace[i] = tmpdata[i];
+ rtrace[ntfft-1-i] = tmpdata[hdr.ns-1-i];
+ }
+ }
+ else {
+ vwarn("file_src has more samples than reflection data: truncated to ntfft = %d samples in middle are removed ", ntfft);
+ for (i = 0; i < ntfft/2; i++) {
+ rtrace[i] = tmpdata[i];
+ rtrace[ntfft-1-i] = tmpdata[hdr.ns-1-i];
+ }
+ }
+
+ rc1fft(rtrace, cwave, ntfft, -1);
+ free(tmpdata);
+ free(rtrace);
+ }
+ else {
+ for (i = 0; i < nfreq; i++) cwave[i].r = 1.0;
+ }
+
+
+/*================ Reading shot records ================*/
+
+ mode=1;
+ readShotData(file_shot, xrcv, xsrc, zsrc, xnx, Refl, nw, nw_low, ngath, nx, nx, fxsb, dxs, ntfft,
+ mode, scale, tsq, Q, f0, reci, &nshots_r, isxcount, reci_xsrc, reci_xrcv, ixmask,verbose);
+
+ if (tap == 2 || tap == 3) {
+ tapersh = (float *)malloc(nx*sizeof(float));
+ for (j = 0; j < ntap; j++)
+ tapersh[j] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
+ for (j = ntap; j < nx-ntap; j++)
+ tapersh[j] = 1.0;
+ for (j = nx-ntap; j < nx; j++)
+ tapersh[j] =(cos(PI*(j-(nx-ntap))/ntap)+1)/2.0;
+ if (verbose) vmess("Taper for shots applied ntap=%d", ntap);
+ for (l = 0; l < nshots; l++) {
+ for (j = 1; j < nw; j++) {
+ for (i = 0; i < nx; i++) {
+ Refl[l*nx*nw+j*nx+i].r *= tapersh[i];
+ Refl[l*nx*nw+j*nx+i].i *= tapersh[i];
+ }
+ }
+ }
+ free(tapersh);
+ }
+
+ /* check consistency of header values */
+ fxf = xsrc[0];
+ if (nx > 1) dxf = (xrcv[nx-1] - xrcv[0])/(float)(nx-1);
+ else dxf = d2;
+ if (NINT(dx*1e3) != NINT(fabs(dxf)*1e3)) {
+ vmess("dx in hdr.d1 (%.3f) and hdr.gx (%.3f) not equal",dx, dxf);
+ if (dxf != 0) dx = fabs(dxf);
+ else verr("gx hdrs not set");
+ vmess("dx used => %f", dx);
+ }
+
+ dxsrc = (float)xsrc[1] - xsrc[0];
+ if (dxsrc == 0) {
+ vwarn("sx hdrs are not filled in!!");
+ dxsrc = dx;
+ }
+
+/*================ Defining focusing operator(s) from R ================*/
+/* G_d = -R(ishot,-t) */
+
+ /* use ishot from Refl, complex-conjugate(time reverse), scale with -1 and convolve with wavelet */
+ if (file_tinv == NULL) {
+ if (verbose) vmess("Selecting G_d from Refl of %s", file_shot);
+ nts = ntfft;
+
+ scl = 1.0/((float)2.0*ntfft);
+ rtrace = (float *)calloc(ntfft,sizeof(float));
+ ctrace = (complex *)calloc(nfreq+1,sizeof(complex));
+ for (i = 0; i < xnx[ishot]; i++) {
+ for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
+ ctrace[j].r = Refl[ishot*nw*nx+m*nx+i].r*cwave[j].r + Refl[ishot*nw*nx+m*nx+i].i*cwave[j].i;
+ ctrace[j].i = -Refl[ishot*nw*nx+m*nx+i].i*cwave[j].r + Refl[ishot*nw*nx+m*nx+i].r*cwave[j].i;;
+ }
+ /* transfrom result back to time domain */
+ cr1fft(ctrace, rtrace, ntfft, 1);
+ for (j = 0; j < nts; j++) {
+ DD[0*nxs*nts+i*nts+j] = -1.0*scl*rtrace[j];
+ }
+ }
+ free(ctrace);
+ free(rtrace);
+
+ xsyn[0] = xsrc[ishot];
+ zsyn[0] = zsrc[ishot];
+ xnxsyn[0] = xnx[ishot];
+ fxse = fxsb = xrcv[0];
+ /* check consistency of header values */
+ for (l=0; l<nshots; l++) {
+ for (i = 0; i < nx; i++) {
+ fxsb = MIN(xrcv[l*nx+i],fxsb);
+ fxse = MAX(xrcv[l*nx+i],fxse);
+ }
+ }
+ dxf = dx;
+ dxs = dx;
+ }
+
+ /* define tapers to taper edges of acquisition */
+ if (tap == 1 || tap == 3) {
+ tapersy = (float *)malloc(nxs*sizeof(float));
+ for (j = 0; j < ntap; j++)
+ tapersy[j] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
+ for (j = ntap; j < nxs-ntap; j++)
+ tapersy[j] = 1.0;
+ for (j = nxs-ntap; j < nxs; j++)
+ tapersy[j] =(cos(PI*(j-(nxs-ntap))/ntap)+1)/2.0;
+ if (verbose) vmess("Taper for operator applied ntap=%d", ntap);
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < nxs; i++) {
+ for (j = 0; j < nts; j++) {
+ DD[l*nxs*nts+i*nts+j] *= tapersy[i];
+ }
+ }
+ }
+ free(tapersy);
+ }
+
+/*================ Check the size of the files ================*/
+
+ if (NINT(dxsrc/dx)*dx != NINT(dxsrc)) {
+ vwarn("source (%.2f) and receiver step (%.2f) don't match",dxsrc,dx);
+ if (reci == 2) vwarn("step used from operator (%.2f) ",dxs);
+ }
+ di = NINT(dxf/dxs);
+ if ((NINT(di*dxs) != NINT(dxf)) && verbose)
+ vwarn("dx in receiver (%.2f) and operator (%.2f) don't match",dx,dxs);
+ if (verbose) {
+ vmess("Number of receivers in focusop = %d", nxs);
+ vmess("number of shots = %d", nshots);
+ vmess("number of receiver/shot = %d", nx);
+ vmess("first model position = %.2f", fxsb);
+ vmess("last model position = %.2f", fxse);
+ vmess("first source position fxf = %.2f", fxf);
+ vmess("source distance dxsrc = %.2f", dxsrc);
+ vmess("last source position = %.2f", fxf+(nshots-1)*dxsrc);
+ vmess("receiver distance dxf = %.2f", dxf);
+ vmess("direction of increasing traces = %d", di);
+ vmess("number of time samples fft nt nts = %d %d %d", ntfft, nt, nts);
+ vmess("time sampling = %e ", dt);
+ vmess("smoothing taper for time-window = %d ", smooth);
+ if (file_rr != NULL) vmess("RR output file = %s ", file_rr);
+ if (file_iter != NULL) vmess("Iterations output file = %s ", file_iter);
+ }
+
+/*================ initializations ================*/
+
+ if (reci) n2out = nxs;
+ else n2out = nshots;
+ mem = Nfoc*n2out*ntfft*sizeof(float)/1048576.0;
+ if (verbose) {
+ vmess("Time-sample range processed = %d:%d", istart, iend);
+ vmess("number of output traces = %d", n2out);
+ vmess("number of output samples = %d", ntfft);
+ vmess("Size of output data/file = %.1f MB", mem);
+ }
+
+ ixa=0; ixb=0;
+
+ synthesisPositions(nx, nt, nxs, nts, dt, xsyn, Nfoc, xrcv, xsrc, xnx, fxse, fxsb,
+ dxs, dxsrc, dx, nshots, ixpos, &npos, isxcount, countmin, reci, verbose);
+
+ if (verbose) {
+ vmess("synthesisPositions: nshots=%d npos=%d", nshots, npos);
+ }
+
+/*================ set variables for output data ================*/
+
+ n1 = nts;
+ n2 = n2out;
+ f1 = ft;
+ f2 = fxsb+dxs*ixpos[0];
+ d1 = dt;
+ if (reci == 0) d2 = dxsrc;
+ else if (reci == 1) d2 = dxs;
+ else if (reci == 2) d2 = dx;
+
+ hdrs_out = (segy *) calloc(n2,sizeof(segy));
+ if (hdrs_out == NULL) verr("allocation for hdrs_out");
+ size = nxs*nts;
+
+ for (i = 0; i < n2; i++) {
+ hdrs_out[i].ns = n1;
+ hdrs_out[i].trid = 1;
+ hdrs_out[i].dt = dt*1000000;
+ hdrs_out[i].f1 = f1;
+ hdrs_out[i].f2 = f2;
+ hdrs_out[i].d1 = d1;
+ hdrs_out[i].d2 = d2;
+ hdrs_out[i].trwf = n2out;
+ hdrs_out[i].scalco = -1000;
+ hdrs_out[i].gx = NINT(1000*(f2+i*d2));
+ hdrs_out[i].scalel = -1000;
+ hdrs_out[i].tracl = i+1;
+ }
+ t1 = wallclock_time();
+ tread = t1-t0;
+ if(verbose) {
+ vmess("*******************************************");
+ vmess("***** Computing Marchenko for all steps****");
+ vmess("*******************************************");
+ fprintf(stderr," %s: Progress: %3d%%",xargv[0],0);
+ }
+ perc=(iend-istart)/100;if(!perc)perc=1;
+
+/*================ start loop over number of time-samples ================*/
+
+ for (ii=istart; ii<iend; ii++) {
+
+/*================ initialization ================*/
+
+ /* once every 'niterskip' time-steps start from fresh G_d and do niter (~20) iterations */
+ if ( ((ii-istart)%niterskip==0) || (ii==istart) ) {
+ niterrun=niter;
+ recur=0;
+ if (verbose>1) vmess("Doing %d iterations to reset recursion at time-sample %d\n",niterrun,ii);
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < nxs; i++) {
+ for (j = 0; j < nts; j++) {
+ G_d[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+i*nts+j];
+ }
+ /* apply mute window for samples above nts-ii */
+ for (j = 0; j < MIN(nts, nts-ii+T*shift-smooth); j++) {
+ G_d[l*nxs*nts+i*nts+j] = 0.0;
+ }
+ for (j = nts-ii+T*shift-smooth, k=1; j < MIN(nts, nts-ii+T*shift); j++, k++) {
+ G_d[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
+ }
+ }
+ for (i = 0; i < npos; i++) {
+ ix = ixpos[i];
+ j = 0;
+ f1min[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+j];
+ for (j = 1; j < nts; j++) {
+ f1min[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j];
+ }
+ /* apply mute window for samples above nts-ii */
+ //for (j = 0; j < nts-ii+T*shift; j++) {
+ // f1min[l*nxs*nts+i*nts+j] = 0.0;
+ //}
+ }
+ }
+ if (file_iter != NULL) {
+ writeDataIter("G_d.su", G_d, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii);
+ }
+ }
+ else { /* use f1min from previous iteration as starting point and do niterec iterations */
+ niterrun=niterec;
+ recur=1;
+ if (verbose>1) vmess("Doing %d iterations using previous result at time-sample %d",niterrun,ii);
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j=0;
+ ix = ixpos[i];
+ G_d[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts] - f1min[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ G_d[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j] - f1min[l*nxs*nts+i*nts+nts-j];
+ }
+ /* apply mute window for samples above nts-ii */
+ for (j = 0; j < MIN(nts,nts-ii+T*shift-smooth); j++) {
+ G_d[l*nxs*nts+i*nts+j] = 0.0;
+ }
+ for (j = nts-ii+T*shift-smooth, k=1; j < MIN(nts, nts-ii+T*shift); j++, k++) {
+ G_d[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
+ }
+ }
+ }
+ }
+/*================ initialization ================*/
+
+ memcpy(Ni, G_d, Nfoc*nxs*ntfft*sizeof(float));
+
+/*================ number of Marchenko iterations ================*/
+
+ for (iter=0; iter<niterrun; iter++) {
+
+ t2 = wallclock_time();
+
+/*================ construction of Ni(-t) = - \int R(x,t) Ni(t) ================*/
+
+ synthesis(Refl, Fop, Ni, RNi, nx, nt, nxs, nts, dt, xsyn, Nfoc,
+ xrcv, xsrc, xnx, fxse, fxsb, dxs, dxsrc, dx, ntfft, nw, nw_low, nw_high, mode,
+ reci, nshots, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv, ixmask, verbose);
+
+ if (file_iter != NULL) {
+ writeDataIter(file_iter, RNi, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter);
+ }
+
+ t3 = wallclock_time();
+ tsyn += t3 - t2;
+
+ /* N_k(x,t) = -N_(k-1)(x,-t) */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ j = 0;
+ ix = ixpos[i];
+ Ni[l*nxs*nts+i*nts+j] = -RNi[l*nxs*nts+ix*nts+j];
+ for (j = 1; j < nts; j++) {
+ Ni[l*nxs*nts+i*nts+j] = -RNi[l*nxs*nts+ix*nts+nts-j];
+ }
+ }
+ }
+
+ if (iter % 2 == 0) { /* even iterations: => f_1^+(t) */
+ /* apply muting for the acausal part */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ /* apply mute window for samples after ii */
+ for (j = MAX(0,ii-T*shift+smooth); j < nts; j++) {
+ Ni[l*nxs*nts+i*nts+j] = 0.0;
+ }
+ for (j = MAX(0,ii-T*shift+smooth), k=0; j < ii; j++, k++) {
+ Ni[l*nxs*nts+i*nts+j] *= costaper[k];
+ }
+ /* apply mute window for delta function at t=0*/
+ for (j = 0; j < shift-smooth; j++) {
+ Ni[l*nxs*nts+i*nts+j] = 0.0;
+ }
+ for (j = MAX(0,shift-smooth), k=1; j < shift; j++, k++) {
+ Ni[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
+ }
+ f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
+ }
+ }
+ }
+/*
+ if (file_iter != NULL) {
+ writeDataIter("f1plus.su", f1plus, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter);
+ }
+*/
+ }
+ else {/* odd iterations: => f_1^-(t) */
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ ix = ixpos[i];
+ if (recur==1) { /* use f1min from previous iteration */
+ for (j = 0; j < nts; j++) {
+ Ni[l*nxs*nts+i*nts+j] += DD[l*nxs*nts+ix*nts+j];
+ }
+ j = 0;
+ f1min[l*nxs*nts+i*nts+j] = -Ni[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ f1min[l*nxs*nts+i*nts+j] = -Ni[l*nxs*nts+i*nts+nts-j];
+ }
+ }
+ else {
+ j = 0;
+ f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j];
+ }
+ }
+ /* apply mute window for delta function at t=0*/
+ for (j = nts-shift+smooth; j < nts; j++) {
+ Ni[l*nxs*nts+i*nts+j] = 0.0;
+ }
+ for (j = nts-shift, k=0; j < MIN(nts, nts-shift+smooth); j++, k++) {
+ Ni[l*nxs*nts+i*nts+j] *= costaper[k];
+ }
+ /* apply mute window for samples above nts-ii */
+ /* important to apply this mute after updating f1min */
+ //for (j = 0; j < nts-ii+T*shift; j++) {
+ // Ni[l*nxs*nts+i*nts+j] = 0.0;
+ //}
+ /* apply mute window for samples above nts-ii */
+ for (j = 0; j < MIN(nts,nts-ii+T*shift-smooth); j++) {
+ Ni[l*nxs*nts+i*nts+j] = 0.0;
+ }
+ for (j = nts-ii+T*shift-smooth, k=1; j < MIN(nts,nts-ii+T*shift); j++, k++) {
+ Ni[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
+ }
+ }
+ }
+ if (file_iter != NULL) {
+ writeDataIter("f1min.su", f1min, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter);
+ }
+ } /* end else (iter) branch */
+ if (file_iter != NULL) {
+ writeDataIter("Ni.su", Ni, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter);
+ }
+
+ t2 = wallclock_time();
+ tcopy += t2 - t3;
+ if (verbose>2) vmess("*** Iteration %d finished ***", iter);
+
+ } /* end of iterations */
+
+ for (l = 0; l < Nfoc; l++) {
+ for (i = 0; i < npos; i++) {
+ RR[l*nxs*nts+i*nts+ii] = f1min[l*nxs*nts+i*nts+ii];
+ }
+ }
+
+ /* To Do? optional write intermediate RR results to file */
+
+ if (verbose) {
+ if(!((iend-ii-istart)%perc)) fprintf(stderr,"\b\b\b\b%3d%%",(ii-istart)*100/(iend-istart));
+ if((ii-istart)==10)t4=wallclock_time();
+ if((ii-istart)==20){
+ t4=(wallclock_time()-t4)*((iend-istart)/10.0);
+ fprintf(stderr,"\r %s: Estimated total compute time = %.2fs.\n %s: Progress: %3d%%",xargv[0],t4,xargv[0],(ii-istart)/((iend-istart)/100.0));
+ }
+ //t4=wallclock_time();
+ tii=(t4-t1)*((float)(iend-istart)/(ii-istart+1.0))-(t4-t1);
+ //vmess("Remaining compute time at time-sample %d = %.2f s.",ii, tii);
+ }
+
+ } /* end of time iterations ii */
+
+ free(Ni);
+ free(G_d);
+ free(f1min);
+ free(f1plus);
+
+ t2 = wallclock_time();
+ if (verbose) {
+ fprintf(stderr,"\b\b\b\b%3d%%\n",100);
+ vmess("Total CPU-time marchenko = %.3f", t2-t0);
+ vmess("with CPU-time synthesis = %.3f", tsyn);
+ vmess("with CPU-time copy array = %.3f", tcopy);
+ vmess(" CPU-time fft data = %.3f", tfft);
+ vmess("and CPU-time read data = %.3f", tread);
+ }
+
+/*================ write output files ================*/
+
+ fp_rr = fopen(file_rr, "w+");
+ if (fp_rr==NULL) verr("error on creating output file %s", file_rr);
+
+ tracf = 1;
+ for (l = 0; l < Nfoc; l++) {
+ if (reci) f2 = fxsb;
+ else f2 = fxf;
+
+ for (i = 0; i < n2; i++) {
+ hdrs_out[i].fldr = l+1;
+ hdrs_out[i].sx = NINT(xsyn[l]*1000);
+ hdrs_out[i].offset = (long)NINT((f2+i*d2) - xsyn[l]);
+ hdrs_out[i].tracf = tracf++;
+ hdrs_out[i].selev = NINT(zsyn[l]*1000);
+ hdrs_out[i].sdepth = NINT(-zsyn[l]*1000);
+ hdrs_out[i].f1 = f1;
+ }
+ ret = writeData(fp_rr, (float *)&RR[l*size], hdrs_out, n1, n2);
+ if (ret < 0 ) verr("error on writing output file.");
+ }
+ ret = fclose(fp_rr);
+ if (ret < 0) verr("err %d on closing output file",ret);
+
+ if (verbose) {
+ t1 = wallclock_time();
+ vmess("and CPU-time write data = %.3f", t1-t2);
+ }
+
+/*================ free memory ================*/
+
+ free(hdrs_out);
+
+ exit(0);
+}
+
+
diff --git a/marchenko/synthesis.c b/marchenko/synthesis.c
index 828da0749d5d9158ad5f14314e999a6437e4c5e6..e975d8746ffa1b44fcf34a15f9b3fa0db03772fb 100644
--- a/marchenko/synthesis.c
+++ b/marchenko/synthesis.c
@@ -95,7 +95,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
/* transform muted Ni (Top) to frequency domain, input for next iteration */
for (l = 0; l < Nfoc; l++) {
/* set Fop to zero, so new operator can be defined within ixpos points */
- memset(&Fop[l*nxs*nw].r, 0, nxs*nw*2*sizeof(float));
+ memset(&Fop[l*nxs*nw], 0, nxs*nw*sizeof(complex));
for (i = 0; i < npos; i++) {
rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
ix = ixpos[i];
@@ -111,7 +111,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
first=0;
for (l = 0; l < Nfoc; l++) {
/* set Fop to zero, so new operator can be defined within all ix points */
- memset(&Fop[l*nxs*nw].r, 0, nxs*nw*2*sizeof(float));
+ memset(&Fop[l*nxs*nw], 0, nxs*nw*sizeof(complex));
for (i = 0; i < nxs; i++) {
rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
for (iw=0; iw<nw; iw++) {
diff --git a/marchenko3D/Makefile b/marchenko3D/Makefile
index 303417a0bd011ab60240be798e91db9ea1732c27..800660b3aeb6130fd25017c6244e4192e7b6e5bd 100644
--- a/marchenko3D/Makefile
+++ b/marchenko3D/Makefile
@@ -7,10 +7,11 @@ LIBS += -L$L -lgenfft -lm $(LIBSM)
#ALL: fmute marchenko marchenko2
-ALL: marchenko3D fmute3D
+ALL: marchenko3D fmute3D TWtransform
SRCT = marchenko3D.c \
getFileInfo3D.c \
+ getFileInfo3DW.c \
readData3D.c \
readShotData3D.c \
readTinvData3D.c \
@@ -41,6 +42,16 @@ SRCJ3 = fmute3D.c \
docpkge.c \
getpars.c
+SRCTW = TWtransform.c \
+ getFileInfo3D.c \
+ writeData3D.c \
+ wallclock_time.c \
+ name_ext.c \
+ verbosepkg.c \
+ atopkge.c \
+ docpkge.c \
+ getpars.c
+
OBJT = $(SRCT:%.c=%.o)
marchenko3D: $(OBJT)
@@ -51,14 +62,18 @@ OBJJ3 = $(SRCJ3:%.c=%.o)
fmute3D: $(OBJJ3)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o fmute3D $(OBJJ3) $(LIBS)
-OBJMS = $(SRCMS:%.c=%.o)
+OBJTW = $(SRCTW:%.c=%.o)
+
+TWtransform: $(OBJTW)
+ $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o TWtransform $(OBJTW) $(LIBS)
-install: marchenko3D fmute3D
+install: marchenko3D fmute3D TWtransform
cp marchenko3D $B
cp fmute3D $B
+ cp TWtransform $B
clean:
- rm -f core marchenko3D $(OBJT) fmute3D $(OBJJ3)
+ rm -f core marchenko3D $(OBJT) fmute3D $(OBJJ3) TWtransform $(OBJTW)
realclean: clean
- rm -f $B/marchenko2 $B/marchenko3D $B/fmute3D
+ rm -f $B/marchenko3D $B/fmute3D $B/TWtransform
diff --git a/marchenko3D/TWtransform.c b/marchenko3D/TWtransform.c
new file mode 100644
index 0000000000000000000000000000000000000000..4c8dba0b0698e29c93e6f503315c25f68b980bfe
--- /dev/null
+++ b/marchenko3D/TWtransform.c
@@ -0,0 +1,187 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "segy.h"
+#include <assert.h>
+
+typedef struct { /* complex number */
+ float r,i;
+} complex;
+
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+
+int optncr(int n);
+void cc1fft(complex *data, int n, int sign);
+void rc1fft(float *rdata, complex *cdata, int n, int sign);
+long writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2);
+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);
+
+/*********************** self documentation **********************/
+char *sdoc[] = {
+" ",
+" TWtransform - Transform data from uncompressed time domain to compressed frequency domain",
+" ",
+" TWtransform file_T= file_W= [optional parameters]",
+" ",
+" Required parameters: ",
+" ",
+" file_T= .................. File containing the uncompressed time domain data",
+" file_W= .................. Output for the compressed frequency domain data",
+" ",
+" Optional parameters: ",
+" ",
+" verbose=0 ................ silent option; >0 displays info",
+" fmin=0 ................... minimum frequency in the output",
+" fmax=70 .................. maximum frequency in the output",
+" mode=1 ................... sign of the frequency transform",
+" ",
+" ",
+" author : Joeri Brackenhoff : (j.a.brackenhoff@tudelft.nl)",
+" author : Jan Thorbecke : (j.w.thorbecke@tudelft.nl)",
+" ",
+NULL};
+/**************** end self doc ***********************************/
+
+int main (int argc, char **argv)
+{
+ FILE *fp, *fp_out;
+ segy hdr;
+ size_t nread;
+ long fldr_shot, sx_shot, sy_shot, itrace, one_shot, igath, iw;
+ long end_of_file, nt, ntfft, nxy, nx, ny, nshots, ret, ntraces;
+ long k, l, m, j, nfreq, nw_low, nw_high, nw, mode, verbose;
+ float scl, scel, *trace, dt, dx, dy, ft, fx, fy, fmin, fmax, *cdata;
+ char *file_T, *file_W;
+ complex *ctrace;
+
+ initargs(argc, argv);
+ requestdoc(1);
+
+ if (!getparstring("file_T", &file_T)) file_T = "in.su";
+ if (file_T==NULL) verr("No file_in is given");
+ if (!getparstring("file_W", &file_W)) file_W = NULL;
+ if (file_W==NULL) verr("No file_W is given");
+ if (!getparfloat("fmin", &fmin)) fmin = 0;
+ if (!getparfloat("fmax", &fmax)) fmax = 70;
+ if (!getparfloat("weight", &scl)) scl = 2.0;
+ if (!getparlong("mode", &mode)) mode = 1;
+ if (!getparlong("verbose", &verbose)) verbose = 1;
+
+ nshots = 0;
+ ret = getFileInfo3D(file_T, &nt, &nx, &ny, &nshots, &dt, &dx, &dy, &ft, &fx, &fy, &scl, &ntraces);
+
+ ntfft = loptncr(nt);
+ nfreq = ntfft/2+1;
+ nw_low = (long)MIN((fmin*ntfft*dt), nfreq-1);
+ nw_low = MAX(nw_low, 1);
+ nw_high = MIN((long)(fmax*ntfft*dt), nfreq-1);
+ nw = nw_high - nw_low + 1;
+
+ nxy = nx*ny;
+
+ /* Reading first header */
+
+ if (file_T == NULL) fp = stdin;
+ else fp = fopen( file_T, "r" );
+ if ( fp == NULL ) {
+ fprintf(stderr,"input file %s has an error\n", file_T);
+ perror("error in opening file: ");
+ fflush(stderr);
+ return -1;
+ }
+
+ if (file_W == NULL) fp_out = stdin;
+ else fp_out = fopen( file_W, "w+" );
+ if ( fp_out == NULL ) {
+ fprintf(stderr,"input file %s has an error\n", file_W);
+ perror("error in opening file: ");
+ fflush(stderr);
+ return -1;
+ }
+
+ fseek(fp, 0, SEEK_SET);
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+
+ fseek(fp, 0, SEEK_SET);
+
+ nt = hdr.ns;
+ dt = hdr.dt/(1E6);
+
+ trace = (float *)calloc(ntfft,sizeof(float));
+ ctrace = (complex *)malloc(ntfft*sizeof(complex));
+ cdata = (float *)calloc(nw*2,sizeof(float));
+
+ end_of_file = 0;
+ one_shot = 1;
+ igath = 0;
+
+ /* Read shots in file */
+
+ while (!end_of_file) {
+
+ /* start reading data (shot records) */
+ itrace = 0;
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) { /* no more data in file */
+ break;
+ }
+
+ sx_shot = hdr.sx;
+ sy_shot = hdr.sy;
+ fldr_shot = hdr.fldr;
+ while (one_shot) {
+
+ nread = fread( trace, sizeof(float), nt, fp );
+ assert (nread == hdr.ns);
+
+ /* transform to frequency domain */
+ if (ntfft > hdr.ns)
+ memset( &trace[nt-1], 0, sizeof(float)*(ntfft-nt) );
+
+ rc1fft(trace,ctrace,(int)ntfft,-1);
+ for (iw=0; iw<nw; iw++) {
+ cdata[iw*2] = scl*ctrace[nw_low+iw].r;
+ cdata[(iw*2)+1] = scl*mode*ctrace[nw_low+iw].i;
+ }
+ itrace++;
+
+ hdr.ep = ntfft;
+ hdr.ns = 2*nw;
+ hdr.unscale = fmin;
+ hdr.ungpow = fmax;
+
+ ret = writeData3D(fp_out, (float *)&cdata[0], &hdr, 2*nw, 1);
+ if (ret < 0 ) verr("error on writing output file.");
+
+ /* read next hdr of next trace */
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) {
+ one_shot = 0;
+ end_of_file = 1;
+ break;
+ }
+ if ((sx_shot != hdr.sx) || (sy_shot != hdr.sy) || (fldr_shot != hdr.fldr)) break;
+ }
+ if (verbose) {
+ vmess("finished reading shot x=%li y=%li (%li) with %li traces",sx_shot,sy_shot,igath,itrace);
+ }
+
+ if (itrace != 0) { /* end of shot record */
+ fseek( fp, -TRCBYTES, SEEK_CUR );
+ igath++;
+ }
+ else {
+ end_of_file = 1;
+ }
+ }
+
+ free(ctrace);
+ free(trace);
+
+ exit(0);
+}
\ No newline at end of file
diff --git a/marchenko3D/getFileInfo3DW.c b/marchenko3D/getFileInfo3DW.c
new file mode 100644
index 0000000000000000000000000000000000000000..f7e02203f79ac2916a8b0e60646e81f3b0dae161
--- /dev/null
+++ b/marchenko3D/getFileInfo3DW.c
@@ -0,0 +1,237 @@
+#define _FILE_OFFSET_BITS 64
+#define _LARGEFILE_SOURCE
+
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <errno.h>
+#include <math.h>
+#include "segy.h"
+
+#define MAX(x,y) ((x) > (y) ? (x) : (y))
+#define MIN(x,y) ((x) < (y) ? (x) : (y))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
+
+/**
+* gets sizes, sampling and min/max values of a SU file
+*
+* AUTHOR:
+* Jan Thorbecke (janth@xs4all.nl)
+* The Netherlands
+**/
+
+void vmess(char *fmt, ...);
+void verr(char *fmt, ...);
+int optncr(int n);
+
+long getFileInfo3DW(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 *sclsxgxsygy, long *nxm)
+{
+ FILE *fp;
+ size_t nread, data_sz;
+ off_t bytes, ret, trace_sz, ntraces;
+ long sx_shot, sy_shot, gx_start, gx_end, gy_start, gy_end, itrace, one_shot, igath, end_of_file, fldr_shot;
+ long verbose=1, igy, nsx, nsy;
+ float scl, *trace, dxsrc, dxrcv, dysrc, dyrcv;
+ segy hdr;
+
+ if (filename == NULL) { /* read from input pipe */
+ *n1=0;
+ *n2=0;
+ *n3=0;
+ return -1; /* Input pipe */
+ }
+ else fp = fopen( filename, "r" );
+ if (fp == NULL) verr("File %s does not exist or cannot be opened", filename);
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+ ret = fseeko( fp, 0, SEEK_END );
+ if (ret<0) perror("fseeko");
+ bytes = ftello( fp );
+
+ *n1 = hdr.ep;
+ if ( (hdr.trid == 1) && (hdr.dt != 0) ) {
+ *d1 = ((float) hdr.dt)*1.e-6;
+ *f1 = ((float) hdr.delrt)/1000.;
+ }
+ else {
+ *d1 = hdr.d1;
+ *f1 = hdr.f1;
+ }
+ *f2 = hdr.f2;
+ *f3 = hdr.gy;
+ *fmin = hdr.unscale;
+ *fmax = hdr.ungpow;
+
+ data_sz = sizeof(float)*(hdr.ns);
+ trace_sz = sizeof(float)*(hdr.ns)+TRCBYTES;
+ ntraces = (long) (bytes/trace_sz);
+
+ if (hdr.scalco < 0) scl = 1.0/fabs((float)hdr.scalco);
+ else if (hdr.scalco == 0) scl = 1.0;
+ else scl = hdr.scalco;
+
+ *sclsxgxsygy = scl;
+ /* check to find out number of traces in shot gather */
+
+ one_shot = 1;
+ itrace = 1;
+ igy = 1;
+ fldr_shot = hdr.fldr;
+ sx_shot = hdr.sx;
+ sy_shot = hdr.sy;
+ gx_start = hdr.gx;
+ gy_start = hdr.gy;
+ gy_end = gy_start;
+ trace = (float *)malloc(hdr.ns*sizeof(float));
+ fseeko( fp, TRCBYTES, SEEK_SET );
+
+ while (one_shot) {
+ nread = fread( trace, sizeof(float), hdr.ns, fp );
+ assert (nread == hdr.ns);
+ if (hdr.gy != gy_end) {
+ gy_end = hdr.gy;
+ igy++;
+ }
+ gx_end = hdr.gx;
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread==0) break;
+ if ((sx_shot != hdr.sx) || (sy_shot != hdr.sy) || (fldr_shot != hdr.fldr) ) break;
+ itrace++;
+ }
+
+ if (itrace>1) {
+ *n2 = itrace/igy;
+ *n3 = igy;
+ if (*n2>1) {
+ dxrcv = (float)(gx_end - gx_start)/(float)(*n2-1);
+ }
+ else {
+ dxrcv = 1.0/scl;
+ }
+ if (*n3>1) {
+ dyrcv = (float)(gy_end - gy_start)/(float)(*n3-1);
+ }
+ else {
+ dyrcv = 1.0/scl;
+ }
+ *d2 = fabs(dxrcv)*scl;
+ *d3 = fabs(dyrcv)*scl;
+ if (NINT(dxrcv*1e3) != NINT(fabs(hdr.d2)*1e3)) {
+ if (dxrcv != 0) *d2 = fabs(dxrcv)*scl;
+ else *d2 = hdr.d2;
+ }
+ }
+ else {
+ *n2 = MAX(hdr.trwf, 1);
+ *n3 = 1;
+ *d2 = hdr.d2;
+ *d3 = 1.0;
+ dxrcv = hdr.d2;
+ dyrcv = 0.0;
+ }
+
+/* check if the total number of traces (ntraces) is correct */
+
+/* expensive way to find out how many gathers there are */
+
+// fprintf(stderr, "ngath = %li dxrcv=%f d2=%f scl=%f \n", *ngath, dxrcv, *d2, scl);
+ if (*ngath == 0) {
+ *n2 = 0;
+ *n3 = 0;
+
+ end_of_file = 0;
+ one_shot = 1;
+ igath = 0;
+ fseeko( fp, 0, SEEK_SET );
+ dxrcv = *d2;
+ dyrcv = *d3;
+
+ while (!end_of_file) {
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) { break; }
+ fldr_shot = hdr.fldr;
+ sx_shot = hdr.sx;
+ gx_start = hdr.gx;
+ gx_end = hdr.gx;
+ sy_shot = hdr.sy;
+ gy_start = hdr.gy;
+ gy_end = hdr.gy;
+
+ itrace = 1;
+ igy = 1;
+ while (one_shot) {
+ fseeko( fp, data_sz, SEEK_CUR );
+ if (hdr.gx != gx_end) dxrcv = MIN(dxrcv,labs(hdr.gx-gx_end));
+ if (hdr.gy != gy_end) {
+ igy++;
+ gy_end = hdr.gy;
+ dyrcv = MIN(dyrcv,labs(hdr.gy-gy_end));
+ }
+ gx_end = hdr.gx;
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) {
+ one_shot = 0;
+ end_of_file = 1;
+ break;
+ }
+ if ((sx_shot != hdr.sx) || (sy_shot != hdr.sy) || (fldr_shot != hdr.fldr)) break;
+ itrace++;
+ }
+ if (itrace>1) {
+ *n2 = MAX(itrace/igy,*n2);
+ *n3 = igy;
+ if (*n2>1) {
+ dxrcv = (float)(gx_end - gx_start)/(float)(*n2-1);
+ }
+ else {
+ dxrcv = 1.0/scl;
+ }
+ if (*n3>1) {
+ dyrcv = (float)(gy_end - gy_start)/(float)(*n3-1);
+ }
+ else {
+ dyrcv = 1.0/scl;
+ }
+ dxsrc = (float)(hdr.sx - sx_shot)*scl;
+ dysrc = (float)(hdr.sy - sy_shot)*scl;
+ }
+ else {
+ *n2 = MAX(MAX(hdr.trwf, 1),*n2);
+ *n3 = 1;
+ *d2 = hdr.d2;
+ *d3 = 1.0;
+ dxrcv = hdr.d2/scl;
+ dyrcv = 1.0/scl;
+ }
+ if (verbose>1) {
+ fprintf(stderr," . Scanning shot %li (%li) with %li traces dxrcv=%.2f dxsrc=%.2f %li %li dyrcv=%.2f dysrc=%.2f %li %li\n",sx_shot,igath,itrace,dxrcv*scl,dxsrc,gx_end,gx_start,dyrcv*scl,dysrc,gy_end,gy_start);
+ }
+ if (itrace != 0) { /* end of shot record */
+ fseeko( fp, -TRCBYTES, SEEK_CUR );
+ igath++;
+ }
+ else {
+ end_of_file = 1;
+ }
+ }
+ *ngath = igath;
+ *d2 = dxrcv*scl;
+ *d3 = dyrcv*scl;
+ }
+ else {
+ /* read last trace header */
+
+ fseeko( fp, -trace_sz, SEEK_END );
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ *ngath = ntraces/((*n2)*(*n3));
+ }
+// *nxm = NINT((*xmax-*xmin)/dxrcv)+1;
+ *nxm = (long)ntraces;
+
+ fclose( fp );
+ free(trace);
+
+ return 0;
+}
\ No newline at end of file
diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index 622c5dd6589a4dd4bdf321b5ce0dcdb7fea610ff..02047cb1c5bb45d0209753d62ddfcafe6a347c7f 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -46,7 +46,11 @@ void convol(float *data1, float *data2, float *con, long nrec, long nsam, float
void applyMute3D(float *data, long *mute, long smooth, long above, long Nfoc, long nxs, long nys, long nt, long *xrcvsyn, long *yrcvsyn, long npos, long shift);
-long getFileInfo3D(char *filename, long *n1, long *n2, long *n3, long *ngath, float *d1, float *d2, float *d3, float *f1, float *f2, float *f3,
+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);
+long getFileInfo3DW(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 *sclsxgxsygy, long *nxm);
long readData3D(FILE *fp, float *data, segy *hdrs, long n1);
long writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2);
@@ -175,7 +179,7 @@ int main (int argc, char **argv)
complex *Refl, *Fop;
char *file_tinv, *file_shot, *file_green, *file_iter, *file_imag, *file_homg, *file_ampscl;
char *file_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_pmin, *file_inp;
- char *file_ray, *file_amp, *file_wav;
+ char *file_ray, *file_amp, *file_wav, *file_shotw;
segy *hdrs_out, *hdrs_Nfoc, *hdrs_iter;
initargs(argc, argv);
@@ -185,6 +189,8 @@ int main (int argc, char **argv)
t0 = wallclock_time();
if (!getparstring("file_shot", &file_shot)) file_shot = NULL;
+ if (!getparstring("file_shotw", &file_shotw)) file_shotw = NULL;
+ if (file_shot==NULL && file_shotw==NULL) verr("No input for the shot data given");
if (!getparstring("file_tinv", &file_tinv)) file_tinv = NULL;
if (!getparstring("file_ray", &file_ray)) file_ray = NULL;
if (!getparstring("file_amp", &file_amp)) file_amp = NULL;
@@ -203,8 +209,8 @@ int main (int argc, char **argv)
if (file_homg!=NULL && file_inp==NULL) verr("Cannot create HomG if no file_inp is given");
if (!getparstring("file_ampscl", &file_ampscl)) file_ampscl = NULL;
if (!getparlong("verbose", &verbose)) verbose = 0;
- if (file_tinv == NULL && file_shot == NULL)
- verr("file_tinv and file_shot cannot be both input pipe");
+ if (file_tinv == NULL && file_shot == NULL && file_shotw == NULL)
+ verr("file_tinv, file_shotw and file_shot cannot be both input pipe");
if (!getparstring("file_green", &file_green)) {
if (verbose) vwarn("parameter file_green not found, assume pipe");
file_green = NULL;
@@ -265,9 +271,16 @@ int main (int argc, char **argv)
fxsb = f2;
fysb = f3;
}
+ if (verbose) vmess("Retrieved file info of the first arrivals");
ngath = 0; /* setting ngath=0 scans all traces; nx contains maximum traces/gather */
- ret = getFileInfo3D(file_shot, &nt, &nx, &ny, &ngath, &d1, &dx, &dy, &ft, &fx, &fy, &scl, &ntraces);
+ if (file_shot!=NULL) {
+ ret = getFileInfo3D(file_shot, &nt, &nx, &ny, &ngath, &d1, &dx, &dy, &ft, &fx, &fy, &scl, &ntraces);
+ }
+ else if (file_shotw!=NULL) {
+ ret = getFileInfo3DW(file_shotw, &nt, &nx, &ny, &ngath, &d1, &dx, &dy, &ft, &fx, &fy, &fmin, &fmax, &scl, &ntraces);
+ }
+ if (verbose) vmess("Retrieved file info of the shot data");
nshots = ngath;
assert (nxs*nys >= nshots);
@@ -334,6 +347,7 @@ int main (int argc, char **argv)
readTinvData3D(file_tinv, xrcvsyn, yrcvsyn, xsyn, ysyn, zsyn, xnxsyn, Nfoc,
nxs, nys, ntfft, mode, muteW, G_d, hw, verbose);
}
+ if (verbose) vmess("Read in first arrivals");
/* reading data added zero's to the number of time samples to be the same as ntfft */
nts = ntfft;
@@ -430,9 +444,17 @@ int main (int argc, char **argv)
/*================ Reading shot records ================*/
- mode=1;
- readShotData3D(file_shot, xrcv, yrcv, xsrc, ysrc, zsrc, xnx, Refl, nw,
+ if (file_shot!=NULL) {
+ mode=1;
+ readShotData3D(file_shot, xrcv, yrcv, xsrc, ysrc, zsrc, xnx, Refl, nw,
nw_low, nshots, nx, ny, ntfft, mode, scale, verbose);
+ }
+ else {
+ mode=0;
+ readShotData3D(file_shotw, xrcv, yrcv, xsrc, ysrc, zsrc, xnx, Refl, nw,
+ nw_low, nshots, nx, ny, ntfft, mode, scale, verbose);
+ }
+ mode=1;
tapersh = (float *)malloc(nx*sizeof(float));
if (tap == 2 || tap == 3) {
@@ -934,12 +956,12 @@ int main (int argc, char **argv)
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 = roundf(zsyn[0]*100.0)/100.0;
- hdrs_Nfoc[0].f2 = roundf(xsyn[0]*100.0)/100.0;
- hdrs_Nfoc[0].ungpow = roundf(xsyn[0]*100.0)/100.0;
- hdrs_Nfoc[0].d1 = roundf(dzim*100.0)/100.0;
- hdrs_Nfoc[0].d2 = roundf(dxs*100.0)/100.0;
- hdrs_Nfoc[0].unscale = roundf(dys*100.0)/100.0;
+ hdrs_Nfoc[0].f1 = roundf(zsyn[0]*1000.0)/1000.0;
+ hdrs_Nfoc[0].f2 = roundf(xsyn[0]*1000.0)/1000.0;
+ hdrs_Nfoc[0].ungpow = roundf(ysyn[0]*1000.0)/1000.0;
+ hdrs_Nfoc[0].d1 = roundf(dzim*1000.0)/1000.0;
+ hdrs_Nfoc[0].d2 = roundf(dxs*1000.0)/1000.0;
+ hdrs_Nfoc[0].unscale = roundf(dys*1000.0)/1000.0;
hdrs_Nfoc[0].dt = (int)(dt*(1E6));
if (fp_imag==NULL) verr("error on creating output file %s", file_imag);
@@ -1011,12 +1033,12 @@ int main (int argc, char **argv)
hdrs_Nfoc[0].sx = sx[0];
hdrs_Nfoc[0].sy = sy[0];
hdrs_Nfoc[0].sdepth = sz[0];
- 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].f1 = roundf(zsyn[0]*100.0)/1000.0;
+ hdrs_Nfoc[0].f2 = roundf(xsyn[0]*1000.0)/1000.0;
+ hdrs_Nfoc[0].ungpow = roundf(ysyn[0]*1000.0)/1000.0;
+ hdrs_Nfoc[0].d1 = roundf(dzim*1000.0)/1000.0;
+ hdrs_Nfoc[0].d2 = roundf(dxs*1000.0)/1000.0;
+ hdrs_Nfoc[0].unscale = roundf(dys*1000.0)/1000.0;
hdrs_Nfoc[0].dt = (int)(dt*(1E6));
ret = writeData3D(fp_homg, (float *)&HomG[0], hdrs_Nfoc, nzim*nyim*nxim*ntfft, 1);
diff --git a/marchenko3D/readShotData3D.c b/marchenko3D/readShotData3D.c
index 2dc337f5c61b8a4a82fbc4ac7450d12e6bc6e615..3c477da2690f23fe5549defbcfa7b4d0b7e02002 100644
--- a/marchenko3D/readShotData3D.c
+++ b/marchenko3D/readShotData3D.c
@@ -16,7 +16,9 @@ int optncr(int n);
void cc1fft(complex *data, int n, int sign);
void rc1fft(float *rdata, complex *cdata, int n, int sign);
-long readShotData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc, long *xnx, complex *cdata, long nw, long nw_low, long nshots, long nx, long ny, long ntfft, long mode, float scale, long verbose)
+long readShotData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc,
+ long *xnx, complex *cdata, long nw, long nw_low, long nshots, long nx, long ny, long ntfft,
+ long mode, float scale, long verbose)
{
FILE *fp;
segy hdr;
@@ -56,8 +58,14 @@ long readShotData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
nt = hdr.ns;
dt = hdr.dt/(1E6);
- trace = (float *)calloc(ntfft,sizeof(float));
- ctrace = (complex *)malloc(ntfft*sizeof(complex));
+ if (mode==0){
+ trace = (float *)calloc(2*nw,sizeof(float));
+ ctrace = (complex *)malloc(nw*sizeof(complex));
+ }
+ else {
+ trace = (float *)calloc(ntfft,sizeof(float));
+ ctrace = (complex *)malloc(ntfft*sizeof(complex));
+ }
isx = (long *)malloc((nxy*nshots)*sizeof(long));
igx = (long *)malloc((nxy*nshots)*sizeof(long));
isy = (long *)malloc((nxy*nshots)*sizeof(long));
@@ -94,17 +102,30 @@ long readShotData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
xrcv[igath*nxy+itrace] = hdr.gx*scl;
yrcv[igath*nxy+itrace] = hdr.gy*scl;
- nread = fread( trace, sizeof(float), nt, fp );
- assert (nread == hdr.ns);
+ if (mode==0) {
+ nread = fread( trace, sizeof(float), 2*nw, fp );
+ assert (nread == hdr.ns);
- /* transform to frequency domain */
- if (ntfft > hdr.ns)
- memset( &trace[nt-1], 0, sizeof(float)*(ntfft-nt) );
-
- rc1fft(trace,ctrace,(int)ntfft,-1);
- for (iw=0; iw<nw; iw++) {
- cdata[igath*nxy*nw+iw*nxy+itrace].r = scale*ctrace[nw_low+iw].r;
- cdata[igath*nxy*nw+iw*nxy+itrace].i = scale*mode*ctrace[nw_low+iw].i;
+ for (iw=0; iw<nw; iw++) {
+ cdata[igath*nxy*nw+iw*nxy+itrace].r = trace[(iw*2)];
+ cdata[igath*nxy*nw+iw*nxy+itrace].i = trace[(iw*2)+1];
+ }
+ //nread = fread(&cdata[igath*nxy*nw+iw*nxy+itrace].r, sizeof(float), 2*nw, fp);
+ //assert (nread == hdr.ns);
+ }
+ else {
+ nread = fread( trace, sizeof(float), nt, fp );
+ assert (nread == hdr.ns);
+
+ /* transform to frequency domain */
+ if (ntfft > hdr.ns)
+ memset( &trace[nt-1], 0, sizeof(float)*(ntfft-nt) );
+
+ rc1fft(trace,ctrace,(int)ntfft,-1);
+ for (iw=0; iw<nw; iw++) {
+ cdata[igath*nxy*nw+iw*nxy+itrace].r = scale*ctrace[nw_low+iw].r;
+ cdata[igath*nxy*nw+iw*nxy+itrace].i = scale*mode*ctrace[nw_low+iw].i;
+ }
}
itrace++;
xnx[igath]+=1;
diff --git a/raytime3d/getParameters3d.c b/raytime3d/getParameters3d.c
index 24c70664c7b9ed01245e4cf9f2444a431c146397..5c5952ee05a3bcfaffdc946259ccd05d33950d74 100644
--- a/raytime3d/getParameters3d.c
+++ b/raytime3d/getParameters3d.c
@@ -346,12 +346,12 @@ long getParameters3d(modPar *mod, recPar *rec, srcPar *src, shotPar *shot, rayPa
if (rec->n) {
dyrcv = rec->yr[rec->nx]-rec->yr[0];
dxrcv = rec->xr[MIN(1,rec->n-1)]-rec->xr[0];
- dzrcv = rec->zr[MIN(1,rec->n-1)]-rec->zr[0];
+ dzrcv = rec->zr[rec->n-1]-rec->zr[0];
vmess("*******************************************");
vmess("************* receiver info ***************");
vmess("*******************************************");
vmess("ntrcv = %li nrcv = %li ", rec->nt, rec->n);
- vmess("nxrcv = %li nyrcv = %li ", rec->nx, rec->ny);
+ vmess("nxrcv = %li nyrcv = %li nzrcv = %li ", rec->nx, rec->ny, rec->nz);
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/raytime3d/raytime3d.h b/raytime3d/raytime3d.h
index 14867d9687f70841a2e460043d533ae31df8757f..2c9f308e489e23603f3a4d66876cff791273c841 100644
--- a/raytime3d/raytime3d.h
+++ b/raytime3d/raytime3d.h
@@ -35,6 +35,7 @@ typedef struct _receiverPar { /* Receiver Parameters */
long n;
long nx;
long ny;
+ long nz;
long nt;
long max_nrec;
long *z;
diff --git a/raytime3d/recvPar3D.c b/raytime3d/recvPar3D.c
index 0ceb2dfcfc98b8b61294cbbb1e85c4f8410b5ed4..3065ae6ce64935d9e30c4ae526d17063f9b5f92e 100644
--- a/raytime3d/recvPar3D.c
+++ b/raytime3d/recvPar3D.c
@@ -28,13 +28,13 @@ long recvPar3D(recPar *rec, float sub_x0, float sub_y0, float sub_z0,
float dx, float dy, float dz, long nx, long ny, long nz)
{
float *xrcv1, *xrcv2, *yrcv1, *yrcv2, *zrcv1, *zrcv2;
- long i, ix, iy, ir, verbose;
+ long i, ix, iy, iz, ir, verbose;
float dxrcv, dyrcv, dzrcv, *dxr, *dyr, *dzr;
float rrcv, dphi, oxrcv, oyrcv, ozrcv, arcv;
double circ, h, a, b, e, s, xr, yr, zr, dr, srun, phase;
float xrange, yrange, zrange, sub_x1, sub_y1, sub_z1;
long Nx1, Nx2, Ny1, Ny2, Nz1, Nz2, Ndx, Ndy, Ndz, iarray, nrec, nh;
- long nxrcv, nyrcv, nzrcv, ncrcv, nrcv, ntrcv, *nlxrcv, *nlyrcv;
+ long nxrcv, nyrcv, nzrcv, ncrcv, nrcv, ntrcv, *nlxrcv, *nlyrcv, *nlzrcv;
float *xrcva, *yrcva, *zrcva;
char* rcv_txt;
FILE *fp;
@@ -267,21 +267,45 @@ long recvPar3D(recPar *rec, float sub_x0, float sub_y0, float sub_z0,
nrcv = 0;
nlxrcv=(long *)malloc(Nx1*sizeof(long));
nlyrcv=(long *)malloc(Nx1*sizeof(long));
+ nlzrcv=(long *)malloc(Nx1*sizeof(long));
for (iarray=0; iarray<Nx1; iarray++) {
xrange = (xrcv2[iarray]-xrcv1[iarray]);
yrange = (yrcv2[iarray]-yrcv1[iarray]);
zrange = (zrcv2[iarray]-zrcv1[iarray]);
- if (dxr[iarray] != 0.0 && dyr[iarray] != 0.0) {
+ if (dxr[iarray] != 0.0 && dyr[iarray] != 0.0 && dzr[iarray] != 0.0) {
nlxrcv[iarray] = NINT(fabs(xrange/dxr[iarray]))+1;
nlyrcv[iarray] = NINT(fabs(yrange/dyr[iarray]))+1;
+ nlzrcv[iarray] = NINT(fabs(zrange/dzr[iarray]))+1;
+ }
+ else if (dxr[iarray] != 0.0 && dyr[iarray] != 0.0) {
+ nlxrcv[iarray] = NINT(fabs(xrange/dxr[iarray]))+1;
+ nlyrcv[iarray] = NINT(fabs(yrange/dyr[iarray]))+1;
+ nlzrcv[iarray] = 1;
+ }
+ else if (dxr[iarray] != 0.0 && dzr[iarray] != 0.0) {
+ nlxrcv[iarray] = NINT(fabs(xrange/dxr[iarray]))+1;
+ nlyrcv[iarray] = 1;
+ nlzrcv[iarray] = NINT(fabs(zrange/dzr[iarray]))+1;
+ }
+ else if (dyr[iarray] != 0.0 && dzr[iarray] != 0.0) {
+ nlxrcv[iarray] = 1;
+ nlyrcv[iarray] = NINT(fabs(yrange/dyr[iarray]))+1;
+ nlzrcv[iarray] = NINT(fabs(zrange/dzr[iarray]))+1;
}
else if (dxr[iarray] != 0.0) {
nlxrcv[iarray] = NINT(fabs(xrange/dxr[iarray]))+1;
nlyrcv[iarray] = 1;
+ nlzrcv[iarray] = 1;
}
else if (dyr[iarray] != 0.0) {
nlxrcv[iarray] = 1;
nlyrcv[iarray] = NINT(fabs(yrange/dyr[iarray]))+1;
+ nlzrcv[iarray] = 1;
+ }
+ else if (dzr[iarray] != 0.0) {
+ nlxrcv[iarray] = 1;
+ nlyrcv[iarray] = 1;
+ nlzrcv[iarray] = NINT(fabs(zrange/dzr[iarray]))+1;
}
else {
if (dzr[iarray] == 0) {
@@ -289,11 +313,13 @@ long recvPar3D(recPar *rec, float sub_x0, float sub_y0, float sub_z0,
}
nlxrcv[iarray] = NINT(fabs(zrange/dzr[iarray]))+1;
nlyrcv[iarray] = NINT(fabs(zrange/dzr[iarray]))+1;
+ nlzrcv[iarray] = NINT(fabs(zrange/dzr[iarray]))+1;
}
- nrcv+=nlyrcv[iarray]*nlxrcv[iarray];
+ nrcv+=nlzrcv[iarray]*nlyrcv[iarray]*nlxrcv[iarray];
}
rec->nx=nlxrcv[0];
rec->ny=nlyrcv[0];
+ rec->nz=nlzrcv[0];
/* Calculate Number of Receivers */
if (verbose) vmess("Total number of linear array receivers: %li",nrcv);
@@ -309,6 +335,7 @@ long recvPar3D(recPar *rec, float sub_x0, float sub_y0, float sub_z0,
free(dzr);
free(nlxrcv);
free(nlyrcv);
+ free(nlzrcv);
}
rec->max_nrec+=nrcv;
}
@@ -448,10 +475,12 @@ long recvPar3D(recPar *rec, float sub_x0, float sub_y0, float sub_z0,
yrange = (yrcv2[iarray]-yrcv1[iarray]);
zrange = (zrcv2[iarray]-zrcv1[iarray]);
if (dxr[iarray] != 0.0) {
- nrcv = nlyrcv[iarray]*nlxrcv[iarray];
+ nrcv = nlzrcv[iarray]*nlyrcv[iarray]*nlxrcv[iarray];
dxrcv = dxr[iarray];
- dyrcv = yrange/(nlyrcv[iarray]-1);
- dzrcv = zrange/(nlxrcv[iarray]-1);
+ if (nlyrcv[iarray]<2) dyrcv=0.0;
+ else dyrcv = yrange/(nlyrcv[iarray]-1);
+ if (nlzrcv[iarray]<2) dzrcv=0.0;
+ else dzrcv = zrange/(nlzrcv[iarray]-1);
if (dyrcv != dyr[iarray]) {
vwarn("For receiver array %li: calculated dyrcv=%f given=%f", iarray, dyrcv, dyr[iarray]);
vwarn("The calculated receiver distance %f is used", dyrcv);
@@ -462,10 +491,12 @@ long recvPar3D(recPar *rec, float sub_x0, float sub_y0, float sub_z0,
}
}
else if (dyr[iarray] != 0.0) {
- nrcv = nlyrcv[iarray]*nlxrcv[iarray];
- dxrcv = xrange/(nlxrcv[iarray]-1);
+ nrcv = nlzrcv[iarray]*nlyrcv[iarray]*nlxrcv[iarray];
+ if (nlxrcv[iarray]<2) dxrcv=0.0;
+ else dxrcv = xrange/(nlxrcv[iarray]-1);
dyrcv = dyr[iarray];
- dzrcv = zrange/(nlxrcv[iarray]-1);
+ if (nlzrcv[iarray]<2) dzrcv=0.0;
+ else dzrcv = zrange/(nlzrcv[iarray]-1);
if (dxrcv != dxr[iarray]) {
vwarn("For receiver array %li: calculated dxrcv=%f given=%f", iarray, dxrcv, dxr[iarray]);
vwarn("The calculated receiver distance %f is used", dxrcv);
@@ -479,10 +510,10 @@ long recvPar3D(recPar *rec, float sub_x0, float sub_y0, float sub_z0,
if (dzr[iarray] == 0) {
verr("For receiver array %li: receiver distance dzrcv is not given", iarray);
}
- nrcv = nlyrcv[iarray]*nlxrcv[iarray];
+ nrcv = nlzrcv[iarray]*nlyrcv[iarray]*nlxrcv[iarray];
dxrcv = xrange/(nrcv-1);
dyrcv = yrange/(nrcv-1);
- dzrcv = dzr[iarray];
+ dzrcv = zrange/(nrcv-1);
if (dxrcv != dxr[iarray]) {
vwarn("For receiver array %li: calculated dxrcv=%f given=%f", iarray, dxrcv, dxr[iarray]);
vwarn("The calculated receiver distance %f is used", dxrcv);
@@ -494,18 +525,20 @@ long recvPar3D(recPar *rec, float sub_x0, float sub_y0, float sub_z0,
}
// calculate coordinates
- for (iy=0; iy<nlyrcv[iarray]; iy++) {
- for (ix=0; ix<nlxrcv[iarray]; ix++) {
- rec->xr[nrec]=xrcv1[iarray]-sub_x0+ix*dxrcv;
- rec->yr[nrec]=yrcv1[iarray]-sub_y0+iy*dyrcv;
- rec->zr[nrec]=zrcv1[iarray]-sub_z0+ix*dzrcv;
-
- rec->x[nrec]=NINT((rec->xr[nrec])/dx);
- rec->y[nrec]=NINT((rec->yr[nrec])/dy);
- rec->z[nrec]=NINT((rec->zr[nrec])/dz);
- nrec++;
- }
- }
+ for (iz=0; iz<nlzrcv[iarray]; iz++) {
+ for (iy=0; iy<nlyrcv[iarray]; iy++) {
+ for (ix=0; ix<nlxrcv[iarray]; ix++) {
+ rec->xr[nrec]=xrcv1[iarray]-sub_x0+ix*dxrcv;
+ rec->yr[nrec]=yrcv1[iarray]-sub_y0+iy*dyrcv;
+ rec->zr[nrec]=zrcv1[iarray]-sub_z0+iz*dzrcv;
+
+ rec->x[nrec]=NINT((rec->xr[nrec])/dx);
+ rec->y[nrec]=NINT((rec->yr[nrec])/dy);
+ rec->z[nrec]=NINT((rec->zr[nrec])/dz);
+ nrec++;
+ }
+ }
+ }
}
free(xrcv1);
free(xrcv2);
@@ -518,8 +551,9 @@ long recvPar3D(recPar *rec, float sub_x0, float sub_y0, float sub_z0,
free(dzr);
free(nlxrcv);
free(nlyrcv);
+ free(nlzrcv);
}
rec->n=rec->max_nrec;
return 0;
-}
\ No newline at end of file
+}
diff --git a/utils/combine.c b/utils/combine.c
index 55f5211c1ad1f4066597bd9bc246bb37173855d3..4c9da58af9762e7123d964c4714d139bb52e4ad4 100755
--- a/utils/combine.c
+++ b/utils/combine.c
@@ -34,8 +34,8 @@ char *sdoc[] = {
" ",
" combine - Combine results into a single result ",
" ",
-" authors : Joeri Brackenhoff : (J.A.Brackenhoff@tudelft.nl)",
-" : Jan Thorbecke : (janth@xs4all.nl)",
+" authors : Joeri Brackenhoff : (J.A.Brackenhoff@tudelft.nl)",
+" : Jan Thorbecke : (janth@xs4all.nl)",
" ",
" Required parameters: ",
"",
@@ -44,10 +44,11 @@ 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",
-" nzmax= ................... Maximum number of files read",
+" file_out=out.su .......... Filename of the output",
+" numb=0 ................... integer number of first file",
+" dnumb=1 .................. integer number of increment in files",
+" nzmax=0 .................. Maximum number of files read",
+" verbose=1 ................ Give detailed information of process",
NULL};
int main (int argc, char **argv)
@@ -69,9 +70,9 @@ int main (int argc, char **argv)
if (!getparstring("file_in", &fin)) fin = NULL;
if (!getparstring("file_out", &fout)) fout = "out.su";
if (!getparlong("numb", &numb)) numb=0;
- if (!getparlong("dnumb", &dnumb)) dnumb=0;
+ if (!getparlong("dnumb", &dnumb)) dnumb=1;
if (!getparlong("nzmax", &nzmax)) nzmax=0;
- if (!getparlong("verbose", &verbose)) verbose=0;
+ if (!getparlong("verbose", &verbose)) verbose=1;
if (!getparlong("compact", &compact)) compact=0;
if (fin == NULL) verr("Incorrect downgoing input");
@@ -176,6 +177,9 @@ int main (int argc, char **argv)
fclose(fp_in);
if (compact==0) {
readSnapData3D(fin2, indata, hdr_in, nt, nx, ny, nz, 0, nx, 0, ny, 0, nz);
+ if (dz == 1.0) {
+ if (is==1) dz = hdr_in[0].f1-z0;
+ }
for (it = 0; it < nt; it++) {
for (iy = 0; iy < ny; iy++) {
for (ix = 0; ix < nx; ix++) {
@@ -188,6 +192,9 @@ int main (int argc, char **argv)
}
else if (compact==1) {
readSnapData3D(fin2, indata, hdr_in, nt, nz, ny, nx, 0, nz, 0, ny, 0, nx);
+ if (dz == 1.0) {
+ if (is==1) dz = hdr_in[0].f1-z0;
+ }
for (it = 0; it < nt; it++) {
for (iy = 0; iy < ny; iy++) {
for (ix = 0; ix < nx; ix++) {
@@ -200,6 +207,9 @@ 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);
+ if (dz == 1.0) {
+ if (is==1) dz = hdr_in[0].f1-z0;
+ }
for (it = 0; it < nt; it++) {
for (iy = 0; iy < ny; iy++) {
for (ix = 0; ix < nx; ix++) {
@@ -239,11 +249,11 @@ int main (int argc, char **argv)
hdr_out[iy*nx+ix].ns = nz*nzs;
hdr_out[iy*nx+ix].trwf = nx*ny;
hdr_out[iy*nx+ix].ntr = hdr_out[iy*nx+ix].fldr*hdr_out[iy*nx+ix].trwf;
- hdr_out[iy*nx+ix].f1 = z0;
- hdr_out[iy*nx+ix].f2 = x0;
+ hdr_out[iy*nx+ix].f1 = roundf(z0*1000.0)/1000.0;
+ hdr_out[iy*nx+ix].f2 = roundf(x0*1000.0)/1000.0;
hdr_out[iy*nx+ix].dt = dt;
- hdr_out[iy*nx+ix].d1 = dz;
- hdr_out[iy*nx+ix].d2 = dx;
+ hdr_out[iy*nx+ix].d1 = roundf(dz*1000.0)/1000.0;
+ hdr_out[iy*nx+ix].d2 = roundf(dx*1000.0)/1000.0;
hdr_out[iy*nx+ix].sx = sx;
hdr_out[iy*nx+ix].gx = (int)roundf(x0 + (ix*dx))*1000;
hdr_out[iy*nx+ix].sy = sy;
diff --git a/utils/mutesnap.c b/utils/mutesnap.c
index 8676a680e401b102cd9f37e47d574f711c6e3137..037aae7dbed6892b003d5e2810ed87e979b19834 100644
--- a/utils/mutesnap.c
+++ b/utils/mutesnap.c
@@ -163,8 +163,8 @@ int main (int argc, char **argv)
/*----------------------------------------------------------------------------*
* Apply the muting to the data
*----------------------------------------------------------------------------*/
- for (iy = 0; iy < nyh; iy++) {
- for (iz = 0; iz < nzh; iz++) {
+ for (iz = 0; iz < nzh; iz++) {
+ for (iy = 0; iy < nyh; iy++) {
for (ix = 0; ix < nxh; ix++) {
if (mode != 2) {
for (it = 0; it < nth; it++) {
@@ -183,17 +183,17 @@ int main (int argc, char **argv)
rmt = MAX(MIN(nth-indrcv,indrcv)-shift-smooth,0);
for (it = ht-rmt+1; it < ht+1; it++) {
if (it-(ht-rmt+1) < smooth) {
- homdata[it*nxh*nzh+ix*nzh+iz] *= costaper[it-(ht-rmt+1)];
- homdata[(nth-it)*nxh*nzh+ix*nzh+iz] *= costaper[it-(ht-rmt+1)];
+ homdata[it*nyh*nxh*nzh+iy*nxh*nzh+ix*nzh+iz] *= costaper[it-(ht-rmt+1)];
+ homdata[(nth-it)*nyh*nxh*nzh+iy*nxh*nzh+ix*nzh+iz] *= costaper[it-(ht-rmt+1)];
}
else{
- homdata[it*nxh*nzh+ix*nzh+iz] = 0.0;
- homdata[(nth-it)*nxh*nzh+ix*nzh+iz] = 0.0;
+ homdata[it*nyh*nxh*nzh+iy*nxh*nzh+ix*nzh+iz] = 0.0;
+ homdata[(nth-it)*nyh*nxh*nzh+iy*nxh*nzh+ix*nzh+iz] = 0.0;
}
}
}
- if (verbose) vmess("Muting Homogeneous Green's function at depth %li from %li depths",iz+1,nzh);
}
+ if (verbose) vmess("Muting Homogeneous Green's function at depth %li from %li depths",iz+1,nzh);
}
free(rtrace);
if (smooth) free(costaper);
diff --git a/utils/syn2d.c b/utils/syn2d.c
index 2e87e9e328784546ec256521a03703b94ab555a3..f8c01b91e6322445e4a0d7e6f8d4e5b70c6fc68f 100644
--- a/utils/syn2d.c
+++ b/utils/syn2d.c
@@ -242,7 +242,6 @@ int main (int argc, char **argv)
if (verbose>=2) vmess("dimensions used: %d x %d",ntmax,nxmax);
}
if (!getparfloat("dt", &dt)) dt = d1;
- fprintf(stderr,"dt=%e\n", dt);
size = ntmax * nxmax;
xrcv = (float *)malloc(nxmax*sizeof(float));
@@ -501,6 +500,7 @@ int main (int argc, char **argv)
for (i = 0; i < n2; i++) {
hdrs_out[i].ns = n1;
+ hdrs_out[i].dt = (unsigned short)(dt*1000000);
hdrs_out[i].trid = hdrs_in[i].trid;
hdrs_out[i].f1 = f1;
hdrs_out[i].f2 = f2;
@@ -510,7 +510,7 @@ int main (int argc, char **argv)
hdrs_out[i].trwf = n2out;
hdrs_out[i].scalco = -1000;
hdrs_out[i].sx = NINT(xsyn[l]*1000);
- hdrs_out[i].gx = NINT(1000*(f2+i*d2));
+ hdrs_out[i].gx = NINT(1000.0*(f2+i*d2));
hdrs_out[i].offset = (long)NINT((f2+i*d2) - xsyn[l]);
hdrs_out[i].scalel = -1000;
hdrs_out[i].selev = NINT(zsyn[l]*1000);