diff --git a/Make_include_OSX b/Make_include_OSX
index ffcf0c7f947965e4fb6c9be8d77a80cd9fd22bd5..f84793d29350a7f7e0d20d47ca0f8bd42571af07 100644
--- a/Make_include_OSX
+++ b/Make_include_OSX
@@ -19,12 +19,12 @@ L = $(ROOT)/lib
#GNU
#CC = gcc-mp-5
-CC = gcc
-FC = gfortran
+#CC = gcc
+#FC = gfortran
# Linux gcc version 4.x
-OPTC = -O3 -ffast-math
+#OPTC = -O3 -ffast-math
#to include parallelisation with OpenMP
-OPTC += -fopenmp
+#OPTC += -fopenmp
# for better performing code you can try:
#OPTC = -O3 -fno-trapping-math -ffast-math -funroll-all-loops -mavx -fomit-frame-pointer -mfpmath=sse -ftree-vectorizer-verbose=1
@@ -83,7 +83,9 @@ LDFLAGS = -Ofast
#############################################################################
# BLAS and LAPACK libraries
+# and FFT LIBRARIES
MKLLIB=${MKLROOT}/lib
+OPTC += -DMKL -I$(MKLROOT)/include -I$(MKLROOT)/include/fftw
#for GNU compilers
#you might need to add intel64 to : ${MKLROOT}/lib/intel64
#BLAS = -L${MKLROOT}/lib/ -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
@@ -96,27 +98,14 @@ BLAS = -Wl,-rpath ${MKLLIB} -L${MKLROOT}/lib/ -lmkl_intel_lp64 -lmkl_sequential
#BLAS = -mkl
#############################################################################
-# FOR FFT LIBRARIES
-#AMD ACML 4.4.0
+# AMD ACML 4.4.0
#AMDROOT = /home/thorbcke/amdsdk/v1.0/acml/open64_64
#OPTC += -DACML440 -I$(AMDROOT)/include
#BLAS = -L$(AMDROOT)/lib -lacml -lfortran -lffio -lrt -lm
#Intel MKL
-MKLLIB=${MKLROOT}/lib
-OPTC += -DMKL -I$(MKLROOT)/include -I$(MKLROOT)/include/fftw
-#for GNU compilers
-#FFT = -Wl,-rpath ${MKLLIB} -L${MKLLIB} -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
-#for clang on OSX
-FFT = -L${MKLROOT}/lib/ -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
-#on linux you want to use groups and MKL is in lib/intel64
-#MKLLIB=${MKLROOT}/lib/intel64
-#FFT = -Wl,-rpath ${MKLLIB} -Wl,--start-group -L${MKLLIB} -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -Wl,--end-group -lpthread -lm -ldl
-#for Intel compilers
-#FFT = -Wl,-rpath ${MKLLIB} -L${MKLLIB} -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
#LIBARIES
-LIBS= -L$L -lgenfft $(FFT) $(BLAS)
-
+LIBS= -L$L -lgenfft $(BLAS)
########################################################################
# standard CFLAGS
diff --git a/Make_include_template b/Make_include_template
index ad3a45193b691d93f2508603b6bac3cd4f6b9ff6..002c45dc914840e33369eeefdd8cf28af97a87e4 100644
--- a/Make_include_template
+++ b/Make_include_template
@@ -83,13 +83,15 @@ OPTC += -fopenmp
#############################################################################
# BLAS and LAPACK libraries
-MKLROOT=/opt/intel/mkl/
+# and FFT LIBRARIES
+#MKLROOT=/opt/intel/mkl/
MKLLIB=${MKLROOT}/lib
+OPTC += -DMKL -I$(MKLROOT)/include -I$(MKLROOT)/include/fftw
#for GNU compilers
#you might need to add intel64 to : ${MKLROOT}/lib/intel64
BLAS = -L${MKLROOT}/lib/ -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
-#for GNU on OSX
-#BLAS = -L${MKLROOT}/lib/ -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
+#for GNU/clang on OSX
+#BLAS = -Wl,-rpath ${MKLLIB} -L${MKLROOT}/lib/ -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
#on linux you want to use groups and MKL is in lib/intel64
#MKLLIB=${MKLROOT}/lib/intel64
#BLAS = -Wl,-rpath ${MKLLIB} -Wl,--start-group -L${MKLLIB} -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -Wl,--end-group -lpthread -lm -ldl
@@ -97,27 +99,13 @@ BLAS = -L${MKLROOT}/lib/ -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm
#BLAS = -mkl
#############################################################################
-# FOR FFT LIBRARIES
-#AMD ACML 4.4.0
+# AMD ACML 4.4.0
#AMDROOT = /home/thorbcke/amdsdk/v1.0/acml/open64_64
#OPTC += -DACML440 -I$(AMDROOT)/include
#BLAS = -L$(AMDROOT)/lib -lacml -lfortran -lffio -lrt -lm
-#Intel MKL
-MKLROOT=/opt/intel/mkl/
-MKLLIB=${MKLROOT}/lib
-OPTC += -DMKL -I$(MKLROOT)/include
-#for GNU compilers
-FFT = -Wl,-rpath ${MKLLIB} -L${MKLLIB} -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
-#for GNU on OSX
-#FFT = -L${MKLROOT}/lib/ -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
-#on linux you want to use groups and MKL is in lib/intel64
-#MKLLIB=${MKLROOT}/lib/intel64
-#FFT = -Wl,-rpath ${MKLLIB} -Wl,--start-group -L${MKLLIB} -lmkl_gf_lp64 -lmkl_sequential -lmkl_core -Wl,--end-group -lpthread -lm -ldl
-#for Intel compilers
-#FFT = -Wl,-rpath ${MKLLIB} -L${MKLLIB} -lmkl_intel_lp64 -lmkl_sequential -lmkl_core -lpthread -lm -ldl
#LIBARIES
-LIBS= -L$L -lgenfft $(FFT) $(BLAS)
+LIBS= -L$L -lgenfft $(BLAS)
########################################################################
diff --git a/fdelmodc/getParameters.c b/fdelmodc/getParameters.c
index edf157883b26f007055e567fec5225bc541c9747..d882ba97c05a56c85c4314bc10ba87480ddcdb25 100644
--- a/fdelmodc/getParameters.c
+++ b/fdelmodc/getParameters.c
@@ -283,7 +283,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
/* Check stability and dispersion setting */
- if (cp_max > dx*stabfactor/dt) {
+ if (cp_max > (dx*stabfactor)/dt) {
vwarn("************ ! Stability ! ****************");
vwarn("From the input file maximum P-wave velocity");
vwarn("in the current model is %f !!", cp_max);
diff --git a/fdelmodc/getWaveletInfo.c b/fdelmodc/getWaveletInfo.c
index 2f3734aae6c38e54653fab909ec5e936a157d8ce..efeee9b4e9bfb8904f0b403f2881c6e0f4ed96ba 100644
--- a/fdelmodc/getWaveletInfo.c
+++ b/fdelmodc/getWaveletInfo.c
@@ -41,7 +41,7 @@ int getWaveletInfo(char *file_src, int *n1, int *n2, float *d1, float *d2, float
off_t bytes;
int ret, one_shot, ntraces;
int optn, nfreq, i, iwmax;
- float *trace;
+ float *trace, df;
float ampl, amplmax, tampl, tamplmax;
complex *ctrace;
segy hdr;
@@ -73,8 +73,9 @@ int getWaveletInfo(char *file_src, int *n1, int *n2, float *d1, float *d2, float
/* check to find out number of traces in shot gather */
- optn = optncr(*n1);
+ optn = optncr(*n1);
nfreq = optn/2 + 1;
+ df = 1.0/(optn*(*d1));
ctrace = (complex *)malloc(nfreq*sizeof(complex));
one_shot = 1;
trace = (float *)malloc(optn*sizeof(float));
@@ -107,12 +108,13 @@ int getWaveletInfo(char *file_src, int *n1, int *n2, float *d1, float *d2, float
iwmax = i;
}
}
+
/* from the maximum amplitude position look for the largest frequency
* which has an amplitude 400 times weaker than the maximum amplitude */
for (i=iwmax;i<nfreq;i++) {
ampl = sqrt(ctrace[i].r*ctrace[i].r+ctrace[i].i*ctrace[i].i);
if (400*ampl < amplmax) {
- *fmax = (i-1)*(1.0/(optn*(*d1)));
+ *fmax = (i-1)*df;
break;
}
}
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);