diff --git a/FFTlib/genfft.h b/FFTlib/genfft.h
index fbef0291f14755a9e871f03f4ed181bf3f2cd10d..4c7ccc6e9a5a5bdc2e2486e870d5a038a3bae778 100644
--- a/FFTlib/genfft.h
+++ b/FFTlib/genfft.h
@@ -79,6 +79,7 @@ extern "C" {
int optncc(int n);
int optncr(int n);
+long loptncr(long n);
void cc1fft(complex *data, int n, int sign);
void ccmfft(complex *data, int n1, int n2, int ld1, int sign);
diff --git a/FFTlib/optnumber.c b/FFTlib/optnumber.c
index bd903b695378a4c706c1a26ca593c8e2f378932e..510ae6ae031e2c30569537c9436b6ea850ed219b 100644
--- a/FFTlib/optnumber.c
+++ b/FFTlib/optnumber.c
@@ -2,6 +2,7 @@
#ifdef SGI
#include "sgintab.h"
int optnfft(int n);
+long loptnfft(long n);
#endif
#if defined(CRAY_MPP_64)
@@ -89,6 +90,45 @@ int optnfft(int n)
}
#endif
+long loptncr(long n)
+{
+
+#ifdef SGI
+ return loptnfft(n);
+#else
+ long n2, n3;
+
+ n2 = pow(2.0, 1.0*(long)(log((float)n)/log(2.0)+0.9999));
+ if (n2 != n) {
+ n3 = npfar(n);
+ if((n3-n) < (n2-n)) return npfar(n);
+ else return n2;
+ }
+ else return n;
+#endif
+}
+
+
+#ifdef SGI
+long loptnfft(long n)
+{
+ long i,j, nmax;
+
+ if (n > NTAB+3) {
+ i=13;
+ nmax=NTAB+3;
+ while (nmax < n) { nmax=(long)pow(2.0,(double)++i); }
+ return nmax;
+ }
+
+ nmax = NTAB;
+ for (i=0; i<NTAB-1 && ntab[i].n<n; ++i);
+ for (j=i+1; j<NTAB-1 && ntab[j].n<=n+nmax; ++j)
+ if (ntab[j].c<ntab[i].c) i = j;
+ return ntab[i].n;
+}
+#endif
+
#if defined(CRAY_MPP_64)
int factorized(int n)
{
diff --git a/include/genfft.h b/include/genfft.h
index fbef0291f14755a9e871f03f4ed181bf3f2cd10d..4c7ccc6e9a5a5bdc2e2486e870d5a038a3bae778 100644
--- a/include/genfft.h
+++ b/include/genfft.h
@@ -79,6 +79,7 @@ extern "C" {
int optncc(int n);
int optncr(int n);
+long loptncr(long n);
void cc1fft(complex *data, int n, int sign);
void ccmfft(complex *data, int n1, int n2, int ld1, int sign);
diff --git a/marchenko3D/Makefile b/marchenko3D/Makefile
index 6a8acbd37b6d7fbef7ae55138ff5e673656ca4eb..415254f5306f8a774072bb3d0eaf79cc84e867a5 100644
--- a/marchenko3D/Makefile
+++ b/marchenko3D/Makefile
@@ -7,18 +7,7 @@ LIBS += -L$L -lgenfft -lm $(LIBSM)
#ALL: fmute marchenko marchenko2
-ALL: fmute marchenko marchenko3D fmute3D
-
-SRCJ = fmute.c \
- getFileInfo.c \
- readData.c \
- applyMute.c \
- writeData.c \
- wallclock_time.c \
- verbosepkg.c \
- atopkge.c \
- docpkge.c \
- getpars.c
+ALL: marchenko3D fmute3D
SRCT = marchenko3D.c \
getFileInfo3D.c \
@@ -26,23 +15,9 @@ SRCT = marchenko3D.c \
readShotData3D.c \
readTinvData3D.c \
synthesis3D.c \
- applyMute.c \
- writeData.c \
- writeDataIter.c \
- wallclock_time.c \
- name_ext.c \
- verbosepkg.c \
- atopkge.c \
- docpkge.c \
- getpars.c
-
-SRCH = marchenko.c \
- getFileInfo.c \
- readData.c \
- readShotData.c \
- readTinvData.c \
- applyMute.c \
- writeData.c \
+ applyMute3D.c \
+ writeData3D.c \
+ ampest3D.c \
writeDataIter.c \
wallclock_time.c \
name_ext.c \
@@ -54,49 +29,30 @@ SRCH = marchenko.c \
SRCJ3 = fmute3D.c \
getFileInfo3D.c \
readData3D.c \
- applyMute.c \
- writeData.c \
+ applyMute3D.c \
+ writeData3D.c \
wallclock_time.c \
verbosepkg.c \
atopkge.c \
docpkge.c \
getpars.c
-OBJJ = $(SRCJ:%.c=%.o)
-
-fmute: $(OBJJ)
- $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o fmute $(OBJJ) $(LIBS)
-
OBJT = $(SRCT:%.c=%.o)
marchenko3D: $(OBJT)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko3D $(OBJT) $(LIBS)
-OBJH = $(SRCH:%.c=%.o)
-
-marchenko: $(OBJH)
- $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko $(OBJH) $(LIBS)
-
-OBJH2 = $(SRCH2:%.c=%.o)
-
-marchenko2: $(OBJH2)
- $(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko2 $(OBJH2) $(LIBS)
-
OBJJ3 = $(SRCJ3:%.c=%.o)
fmute3D: $(OBJJ3)
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o fmute3D $(OBJJ3) $(LIBS)
-install: fmute marchenko marchenko3D fmute3D
- cp fmute $B
- cp marchenko $B
+install: marchenko3D fmute3D
cp marchenko3D $B
cp fmute3D $B
-# cp marchenko2 $B
-
clean:
- rm -f core fmute $(OBJJ) marchenko $(OBJH) marchenko2 $(OBJH2) marchenko3D $(OBJT) fmute3D $(OBJJ3)
+ rm -f core marchenko3D $(OBJT) fmute3D $(OBJJ3)
realclean: clean
- rm -f $B/fmute $B/marchenko $B/marchenko2 $B/marchenko3D $B/fmute3D
+ rm -f $B/marchenko2 $B/marchenko3D $B/fmute3D
diff --git a/marchenko3D/applyMute3D.c b/marchenko3D/applyMute3D.c
index cc6f86d2800ac5671d601752c07efeffadab8776..82adead728e57f820e49dbf3ea14b337ae6f64dd 100644
--- a/marchenko3D/applyMute3D.c
+++ b/marchenko3D/applyMute3D.c
@@ -10,13 +10,13 @@
#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
-#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
-void applyMute3D( float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *ixpos, int npos, int shift)
+void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, long nxs, long nt, long *ixpos, long npos, long shift)
{
- int ix, iy, i, j, l, isyn;
+ long ix, iy, i, j, l, isyn;
float *costaper, scl;
- int imute, tmute;
+ long imute, tmute;
if (smooth) {
costaper = (float *)malloc(smooth*sizeof(float));
diff --git a/marchenko3D/demo/oneD/marchenko.scr b/marchenko3D/demo/oneD/marchenko.scr
index 8538681e22d6a73dbaef67b3c09efeebfbce2916..0c48a3c1379499e239eaa03debb8ffa5d937e3d7 100755
--- a/marchenko3D/demo/oneD/marchenko.scr
+++ b/marchenko3D/demo/oneD/marchenko.scr
@@ -5,12 +5,13 @@ export OMP_NUM_THREADS=1
#mute all events below the first arrival to get the intial focusing field
fmute file_shot=iniFocus_rp.su file_out=p0plus.su above=-1 shift=-8 verbose=1 check=0 hw=8
+suwind < p0plus.su itmax=1024 > p0plus2.su
#apply the Marchenko algorithm
-marchenko file_shot=shotsdx5_rp.su file_tinv=p0plus.su nshots=901 verbose=2 \
- tap=0 niter=8 hw=8 shift=12 smooth=3 scale=4 \
- file_green=pgreen3.su file_gplus=Gplus03.su file_gmin=Gmin03.su \
- file_f1plus=f1plus03.su file_f1min=f1min03.su file_f2=f23.su
+marchenko file_shot=shotsdx5_rp.su file_tinv=p0plus2.su nshots=901 verbose=2 \
+ tap=0 niter=8 hw=8 shift=12 smooth=3 scale=2 \
+ file_green=pgreen.su file_gplus=Gp.su file_gmin=Gm.su \
+ file_f1plus=f1p.su file_f1min=f1m.su file_f2=f2.su
exit
diff --git a/marchenko3D/demo/oneD/p5all.scr b/marchenko3D/demo/oneD/p5all.scr
index c749523043cbeb7f11ea3a1a86f72b397b805271..ca9877f16937fcf714309842e7a713a31c6d003c 100755
--- a/marchenko3D/demo/oneD/p5all.scr
+++ b/marchenko3D/demo/oneD/p5all.scr
@@ -4,9 +4,9 @@ export PATH=$HOME/src/OpenSource/bin:$PATH:
# Generate the full R matrix for a fixed spread geometry.
-dxshot=10000 # with scalco factor of 1000
+dxshot=5000 # with scalco factor of 1000
ishot=0
-nshots=451
+nshots=901
echo $1
@@ -22,7 +22,7 @@ do
(( ishot = $ishot + 1))
- suwind < shot_kxky.su key=tracl min=$tr1 max=$tr2 | \
+ suwind < shot5_rp.su key=tracl min=$tr1 max=$tr2 | \
sushw key=sx,gx,fldr,trwf \
a=$xsrc,-2250000,$ishot,$nshots b=0,$dxshot,0,0 j=0,$nshots,0,0 | \
suchw key1=offset key2=gx key3=sx c=-1 d=1000 >> shotsdx5_rp.su
diff --git a/marchenko3D/fmute3D.c b/marchenko3D/fmute3D.c
index 533c7eb5e1932106d48731f788172a7d49a3e863..fa03b817592628cb532ae04bc489f758ac704fbb 100644
--- a/marchenko3D/fmute3D.c
+++ b/marchenko3D/fmute3D.c
@@ -13,7 +13,7 @@
#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
-#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
#ifndef COMPLEX
typedef struct _complexStruct { /* complex number */
@@ -21,11 +21,12 @@ typedef struct _complexStruct { /* complex number */
} complex;
#endif/* complex */
-int getFileInfo3D(char *filename, int *n1, int *n2, int *n3, int *ngath, float *d1, float *d2, float *d3, float *f1, float *f2, float *f3, float *sclsxgxsygy, int *nxm);
-int readData3D(FILE *fp, float *data, segy *hdrs, int n1);
-int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
-int disp_fileinfo3D(char *file, int n1, int n2, int n3, float f1, float f2, float f3, float d1, float d2, float d3, segy *hdrs);
-void applyMute( float *data, int *mute, int smooth, int above, int Nfoc, int nxs, int nt, int *ixpos, int npos, int 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,
+ 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);
+long disp_fileinfo3D(char *file, long n1, long n2, long n3, float f1, float f2, float f3, float d1, float d2, float d3, segy *hdrs);
+void applyMute3D( float *data, long *mute, long smooth, long above, long Nfoc, long nxs, long nt, long *ixpos, long npos, long shift);
double wallclock_time(void);
/*********************** self documentation **********************/
@@ -60,13 +61,13 @@ char *sdoc[] = {
NULL};
/**************** end self doc ***********************************/
-int main (int argc, char **argv)
+long main (int argc, char **argv)
{
FILE *fp_in1, *fp_in2, *fp_out, *fp_chk, *fp_psline1, *fp_psline2;
- int verbose, shift, k, nx1, ny1, nt1, nx2, ny2, nt2, nxy;
- int ntmax, nxmax, nymax, ret, i, j, l, jmax, ixmax, iymax, above, check;
- int size, ntraces, ngath, *maxval, hw, smooth;
- int tstart, tend, scale, *xrcv;
+ long verbose, shift, k, nx1, ny1, nt1, nx2, ny2, nt2, nxy;
+ long ntmax, nxmax, nymax, ret, i, j, l, jmax, ixmax, iymax, above, check;
+ long size, ntraces, ngath, *maxval, hw, smooth;
+ long tstart, tend, scale, *xrcv;
float dt, dt1, dx1, dy1, ft1, fx1, fy1, t0, t1, dt2, dx2, dy2, ft2, fx2, fy2;
float w1, w2, dxrcv, dyrcv;
float *tmpdata, *tmpdata2, *costaper;
@@ -81,17 +82,17 @@ int main (int argc, char **argv)
if(!getparstring("file_mute", &file_mute)) file_mute=NULL;
if(!getparstring("file_shot", &file_shot)) file_shot=NULL;
if(!getparstring("file_out", &file_out)) file_out=NULL;
- if(!getparint("ntmax", &ntmax)) ntmax = 1024;
- if(!getparint("nxmax", &nxmax)) nxmax = 512;
- if(!getparint("above", &above)) above = 0;
- if(!getparint("check", &check)) check = 0;
- if(!getparint("scale", &scale)) scale = 0;
- if(!getparint("hw", &hw)) hw = 15;
- if(!getparint("smooth", &smooth)) smooth = 0;
+ if(!getparlong("ntmax", &ntmax)) ntmax = 1024;
+ if(!getparlong("nxmax", &nxmax)) nxmax = 512;
+ if(!getparlong("above", &above)) above = 0;
+ if(!getparlong("check", &check)) check = 0;
+ if(!getparlong("scale", &scale)) scale = 0;
+ if(!getparlong("hw", &hw)) hw = 15;
+ if(!getparlong("smooth", &smooth)) smooth = 0;
if(!getparfloat("w1", &w1)) w1=1.0;
if(!getparfloat("w2", &w2)) w2=1.0;
- if(!getparint("shift", &shift)) shift=0;
- if(!getparint("verbose", &verbose)) verbose=0;
+ if(!getparlong("shift", &shift)) shift=0;
+ if(!getparlong("verbose", &verbose)) verbose=0;
/* Reading input data for file_mute */
@@ -99,11 +100,11 @@ int main (int argc, char **argv)
ngath = 1;
ret = getFileInfo3D(file_mute, &nt1, &nx1, &ny1, &ngath, &dt1, &dx1, &dy1, &ft1, &fx1, &fy1, &sclsxgx, &ntraces);
- if (!getparint("ntmax", &ntmax)) ntmax = nt1;
- if (!getparint("nxmax", &nxmax)) nxmax = nx1;
- if (!getparint("nymax", &nymax)) nymax = ny1;
+ if (!getparlong("ntmax", &ntmax)) ntmax = nt1;
+ if (!getparlong("nxmax", &nxmax)) nxmax = nx1;
+ if (!getparlong("nymax", &nymax)) nymax = ny1;
if (verbose>=2 && (ntmax!=nt1 || nxmax!=nx1 || nymax!= ny1))
- vmess("dimensions overruled: %d x %d y %d",ntmax,nxmax,nymax);
+ vmess("dimensions overruled: %li x %li y %li",ntmax,nxmax,nymax);
if(!getparfloat("dt", &dt)) dt=dt1;
fp_in1 = fopen(file_mute, "r");
@@ -130,9 +131,9 @@ int main (int argc, char **argv)
ngath = 1;
ret = getFileInfo3D(file_shot, &nt2, &nx2, &ny2, &ngath, &dt2, &dx2, &dy2, &ft2, &fx2, &fy2, &sclshot, &ntraces);
- if (!getparint("ntmax", &ntmax)) ntmax = nt2;
- if (!getparint("nxmax", &nxmax)) nxmax = nx2;
- if (!getparint("nymax", &nymax)) nymax = ny2;
+ if (!getparlong("ntmax", &ntmax)) ntmax = nt2;
+ if (!getparlong("nxmax", &nxmax)) nxmax = nx2;
+ if (!getparlong("nymax", &nymax)) nymax = ny2;
size = ntmax * nxmax * nymax;
tmpdata2 = (float *)malloc(size*sizeof(float));
@@ -171,13 +172,13 @@ int main (int argc, char **argv)
sclsxgx = sclshot;
}
- if (verbose) vmess("sampling file_mute=%d, file_shot=%d", nt1, nt2);
+ if (verbose) vmess("sampling file_mute=%li, file_shot=%li", nt1, nt2);
/*================ initializations ================*/
nxy = nx1*ny1;
- maxval = (int *)calloc(nxy,sizeof(int));
- xrcv = (int *)calloc(nxy,sizeof(int));
+ maxval = (long *)calloc(nxy,sizeof(long));
+ xrcv = (long *)calloc(nxy,sizeof(long));
if (file_out==NULL) fp_out = stdout;
else {
@@ -206,7 +207,7 @@ int main (int argc, char **argv)
k=1;
while (nxy > 0) {
- if (verbose) vmess("processing input gather %d", k);
+ if (verbose) vmess("processing input gather %li", k);
/*================ loop over all shot records ================*/
@@ -267,7 +268,7 @@ int main (int argc, char **argv)
}
}
maxval[iymax*nx1+ixmax] = jmax;
- if (verbose >= 3) vmess("Mute max at src-trace x=%d y=%d is sample %d", ixmax, iymax, maxval[iymax*nx1+ixmax]);
+ if (verbose >= 3) vmess("Mute max at src-trace x=%li y=%li is sample %li", ixmax, iymax, maxval[iymax*nx1+ixmax]);
/* search forward in x-trace direction from maximum in file */
for (i = ixmax+1; i < nx1; i++) {
@@ -312,7 +313,7 @@ int main (int argc, char **argv)
}
}
maxval[i*nx1+ixmax] = jmax;
- if (verbose >= 8) vmess("Mute max at src-trace x=%d y=%d is sample %d", ixmax, i, maxval[i*nx1+ixmax]);
+ if (verbose >= 8) vmess("Mute max at src-trace x=%li y=%li is sample %li", ixmax, i, maxval[i*nx1+ixmax]);
/* search forward in x-trace direction from maximum in file */
for (l = ixmax+1; l < nx1; l++) {
tstart = MAX(0, (maxval[i*nx1+l-1]-hw));
@@ -357,7 +358,7 @@ int main (int argc, char **argv)
}
}
maxval[i*nx1+ixmax] = jmax;
- if (verbose >= 8) vmess("Mute max at src-trace x=%d y=%d is sample %d", ixmax, i, maxval[i*nx1+ixmax]);
+ if (verbose >= 8) vmess("Mute max at src-trace x=%li y=%li is sample %li", ixmax, i, maxval[i*nx1+ixmax]);
/* search forward in x-trace direction from maximum in file */
for (l = ixmax+1; l < nx1; l++) {
tstart = MAX(0, (maxval[i*nx1+l-1]-hw));
@@ -411,11 +412,11 @@ int main (int argc, char **argv)
/*================ apply mute window ================*/
- applyMute(tmpdata2, maxval, smooth, above, 1, nx2*ny2, nt2, xrcv, nx2*ny2, shift);
+ applyMute3D(tmpdata2, maxval, smooth, above, 1, nx2*ny2, nt2, xrcv, nx2*ny2, shift);
/*================ write result to output file ================*/
- ret = writeData(fp_out, tmpdata2, hdrs_in2, nt2, nx2*ny2);
+ ret = writeData3D(fp_out, tmpdata2, hdrs_in2, nt2, nx2*ny2);
if (ret < 0 ) verr("error on writing output file.");
/* put mute window in file to check correctness of mute */
@@ -434,7 +435,7 @@ int main (int argc, char **argv)
}
}
}
- ret = writeData(fp_chk, tmpdata, hdrs_in1, nt1, nx1*ny1);
+ ret = writeData3D(fp_chk, tmpdata, hdrs_in1, nt1, nx1*ny1);
if (ret < 0 ) verr("error on writing check file.");
for (l=0; l<ny1; l++) {
for (i=0; i<nx1; i++) {
@@ -463,10 +464,10 @@ int main (int argc, char **argv)
}
break;
}
- nt1 = (int)hdrs_in1[0].ns;
- if (nt1 > ntmax) verr("n_samples (%d) greater than ntmax", nt1);
- if (nx1 > nxmax) verr("n_traces (%d) greater than nxmax", nx1);
- if (ny1 > nymax) verr("n_traces (%d) greater than nymax", ny1);
+ nt1 = (long)hdrs_in1[0].ns;
+ if (nt1 > ntmax) verr("n_samples (%li) greater than ntmax", nt1);
+ if (nx1 > nxmax) verr("n_traces (%li) greater than nxmax", nx1);
+ if (ny1 > nymax) verr("n_traces (%li) greater than nymax", ny1);
if (verbose) {
disp_fileinfo3D(file_mute, nt1, nx1, ny1, ft1, fx1, fy1, dt, dx1, dy1, hdrs_in1);
}
@@ -481,10 +482,10 @@ int main (int argc, char **argv)
fclose(fp_in2);
break;
}
- nt2 = (int)hdrs_in2[0].ns;
- if (nt2 > ntmax) verr("n_samples (%d) greater than ntmax", nt2);
- if (nx2 > nxmax) verr("n_traces (%d) greater than nxmax", nx2);
- if (ny2 > nymax) verr("n_traces (%d) greater than nymax", ny2);
+ nt2 = (long)hdrs_in2[0].ns;
+ if (nt2 > ntmax) verr("n_samples (%li) greater than ntmax", nt2);
+ if (nx2 > nxmax) verr("n_traces (%li) greater than nxmax", nx2);
+ if (ny2 > nymax) verr("n_traces (%li) greater than nymax", ny2);
if (verbose) {
disp_fileinfo3D(file_shot, nt2, nx2, ny2, ft2, fx2, fy2, dt, dx2, dy2, hdrs_in2);
}
diff --git a/marchenko3D/getFileInfo3D.c b/marchenko3D/getFileInfo3D.c
index 1311901c2cdc36689b0b7ce0d8c9c4c77541f2ae..13553faf2f1c4d53f0cd995e6e766e238f3d0209 100644
--- a/marchenko3D/getFileInfo3D.c
+++ b/marchenko3D/getFileInfo3D.c
@@ -10,7 +10,7 @@
#define MAX(x,y) ((x) > (y) ? (x) : (y))
#define MIN(x,y) ((x) < (y) ? (x) : (y))
-#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
/**
* gets sizes, sampling and min/max values of a SU file
@@ -24,13 +24,13 @@ void vmess(char *fmt, ...);
void verr(char *fmt, ...);
int optncr(int n);
-int getFileInfo3D(char *filename, int *n1, int *n2, int *n3, int *ngath, float *d1, float *d2, float *d3, float *f1, float *f2, float *f3, float *sclsxgxsygy, int *nxm)
+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)
{
FILE *fp;
size_t nread, data_sz;
off_t bytes, ret, trace_sz, ntraces;
- int sx_shot, sy_shot, gx_start, gx_end, gy_start, gy_end, itrace, one_shot, igath, end_of_file, fldr_shot;
- int verbose=1, igy, nsx, nsy;
+ 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;
@@ -62,7 +62,7 @@ int getFileInfo3D(char *filename, int *n1, int *n2, int *n3, int *ngath, float *
data_sz = sizeof(float)*(*n1);
trace_sz = sizeof(float)*(*n1)+TRCBYTES;
- ntraces = (int) (bytes/trace_sz);
+ ntraces = (long) (bytes/trace_sz);
if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
else if (hdr.scalco == 0) scl = 1.0;
@@ -122,7 +122,7 @@ int getFileInfo3D(char *filename, int *n1, int *n2, int *n3, int *ngath, float *
/* expensive way to find out how many gathers there are */
-// fprintf(stderr, "ngath = %d dxrcv=%f d2=%f scl=%f \n", *ngath, dxrcv, *d2, scl);
+// fprintf(stderr, "ngath = %li dxrcv=%f d2=%f scl=%f \n", *ngath, dxrcv, *d2, scl);
if (*ngath == 0) {
*n2 = 0;
*n3 = 0;
@@ -175,7 +175,7 @@ int getFileInfo3D(char *filename, int *n1, int *n2, int *n3, int *ngath, float *
dysrc = (float)(hdr.sy - sy_shot)*scl;
}
if (verbose>1) {
- fprintf(stderr," . Scanning shot %d (%d) with %d traces dxrcv=%.2f dxsrc=%.2f %d %d dyrcv=%.2f dysrc=%.2f %d %d\n",sx_shot,igath,itrace,dxrcv*scl,dxsrc,gx_end,gx_start,dyrcv*scl,dysrc,gy_end,gy_start);
+ 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 );
@@ -197,7 +197,7 @@ int getFileInfo3D(char *filename, int *n1, int *n2, int *n3, int *ngath, float *
*ngath = ntraces/((*n2)*(*n3));
}
// *nxm = NINT((*xmax-*xmin)/dxrcv)+1;
- *nxm = (int)ntraces;
+ *nxm = (long)ntraces;
fclose( fp );
free(trace);
@@ -205,13 +205,13 @@ int getFileInfo3D(char *filename, int *n1, int *n2, int *n3, int *ngath, float *
return 0;
}
-int disp_fileinfo3D(char *file, int n1, int n2, int n3, float f1, float f2, float f3, float d1, float d2, float d3, segy *hdrs)
+long disp_fileinfo3D(char *file, long n1, long n2, long n3, float f1, float f2, float f3, float d1, float d2, float d3, segy *hdrs)
{
vmess("file %s contains", file);
- vmess("*** n1 = %d n2 = %d n3 = %d ntftt=%d", n1, n2, n3, optncr(n1));
+ vmess("*** n1 = %li n2 = %li n3 = %li ntftt=%li", n1, n2, n3, (long)optncr((int)n1));
vmess("*** d1 = %.5f d2 = %.5f d3 = %.5f", d1, d2, d3);
vmess("*** f1 = %.5f f2 = %.5f f3 = %.5f", f1, f2, f3);
- vmess("*** fldr = %d sx = %d sy = %d", hdrs[0].fldr, hdrs[0].sx, hdrs[0].sy);
+ vmess("*** fldr = %li sx = %li sy = %li", hdrs[0].fldr, hdrs[0].sx, hdrs[0].sy);
return 0;
}
diff --git a/marchenko3D/marchenko3D.c b/marchenko3D/marchenko3D.c
index 983aa15579a63adbc16529cf527ce878c1068e57..ace2b3179f7110a2d0884615c06fefa5ac7096bc 100644
--- a/marchenko3D/marchenko3D.c
+++ b/marchenko3D/marchenko3D.c
@@ -53,6 +53,9 @@ long writeData3D(FILE *fp, float *data, segy *hdrs, long n1, long n2);
long disp_fileinfo3D(char *file, long n1, long n2, long n3, float f1, float f2, float f3, float d1, float d2, float d3, segy *hdrs);
double wallclock_time(void);
+void AmpEst3D(float *f1d, float *Gd, float *ampest, long Nfoc, long nxs, long nys, long ntfft, long *ixpos, long npos,
+ char *file_wav, float dx, float dy, float dt);
+
void synthesisPositions3D(long nx, long ny, long nxs, long nys, long Nfoc, float *xrcv, float *yrcv, float *xsrc, float *ysrc,
long *xnx, float fxse, float fyse, float fxsb, float fysb, float dxs, float dys, long nshots, long nxsrc, long nysrc,
long *ixpos, long *npos, long reci, long verbose);
@@ -122,7 +125,7 @@ int main (int argc, char **argv)
long size, n1, n2, n3, ntap, tap, dxi, dyi, ntraces, pad;
long nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
long reci, countmin, mode, n2out, n3out, verbose, ntfft;
- long iter, niter, tracf, *muteW;
+ long iter, niter, tracf, *muteW, ampest;
long hw, smooth, above, shift, *ixpos, npos, ix;
long nshots_r, *isxcount, *reci_xsrc, *reci_xrcv;
float fmin, fmax, *tapersh, *tapersy, fxf, fyf, dxf, dyf, *xsrc, *ysrc, *xrcv, *yrcv, *zsyn, *zsrc, *xrcvsyn, *yrcvsyn;
@@ -131,9 +134,9 @@ int main (int argc, char **argv)
float *green, *f2p, *pmin, *G_d, dt, dx, dy, dxs, dys, scl, mem;
float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus;
float xmin, xmax, ymin, ymax, scale, tsq, Q, f0;
- float *ixmask, *iymask;
+ float *ixmask, *iymask, *ampscl, *Gd;
complex *Refl, *Fop;
- char *file_tinv, *file_shot, *file_green, *file_iter;
+ char *file_tinv, *file_shot, *file_green, *file_iter, *file_wav;
char *file_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_pmin;
segy *hdrs_out;
@@ -153,6 +156,7 @@ int main (int argc, char **argv)
if (!getparstring("file_f2", &file_f2)) file_f2 = NULL;
if (!getparstring("file_pmin", &file_pmin)) file_pmin = NULL;
if (!getparstring("file_iter", &file_iter)) file_iter = NULL;
+ if (!getparstring("file_wav", &file_wav)) file_wav = 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");
@@ -170,6 +174,7 @@ int main (int argc, char **argv)
if (!getparlong("tap", &tap)) tap = 0;
if (!getparlong("ntap", &ntap)) ntap = 0;
if (!getparlong("pad", &pad)) pad = 0;
+ if (!getparlong("ampest", &est)) ampest = 0;
if(!getparlong("niter", &niter)) niter = 10;
if(!getparlong("hw", &hw)) hw = 15;
@@ -199,7 +204,7 @@ int main (int argc, char **argv)
if (!getparfloat("dt", &dt)) dt = d1;
- ntfft = optncr(MAX(nt+pad, nts+pad));
+ ntfft = loptncr(MAX(nt+pad, nts+pad));
nfreq = ntfft/2+1;
nw_low = (long)MIN((fmin*ntfft*dt), nfreq-1);
nw_low = MAX(nw_low, 1);
@@ -236,7 +241,7 @@ int main (int argc, char **argv)
yrcv = (float *)calloc(nshots*ny*nx,sizeof(float)); // x-rcv postions of shots
xsrc = (float *)calloc(nshots,sizeof(float)); //x-src position of shots
ysrc = (float *)calloc(nshots,sizeof(float)); //x-src position of shots
- zsrc = (float *)calloc(nshots,sizeof(float)); // z-src position of shots
+ zsrc = (float *)calloc(nshots,sizeof(float)); //z-src position of shots
xnx = (long *)calloc(nshots,sizeof(long)); // number of traces per shot
if (reci!=0) {
@@ -268,7 +273,7 @@ int main (int argc, char **argv)
for (j = 0; j < nxs; j++) tapersy[j] = 1.0;
}
if (tap == 1 || tap == 3) {
- if (verbose) vmess("Taper for operator applied ntap=%d", ntap);
+ if (verbose) vmess("Taper for operator applied ntap=%li", ntap);
for (l = 0; l < Nfoc; l++) {
for (k = 0; k < nys; k++) {
for (i = 0; i < nxs; i++) {
@@ -335,7 +340,7 @@ int main (int argc, char **argv)
for (j = 0; j < nx; j++) tapersh[j] = 1.0;
}
if (tap == 2 || tap == 3) {
- if (verbose) vmess("Taper for shots applied ntap=%d", ntap);
+ if (verbose) vmess("Taper for shots applied ntap=%li", ntap);
for (l = 0; l < nshots; l++) {
for (j = 1; j < nw; j++) {
for (k = 0; k < ny; k++) {
@@ -400,20 +405,20 @@ int main (int argc, char **argv)
if ((NINT(dyi*dys) != NINT(dyf)) && verbose)
vwarn("dy in receiver (%.2f) and operator (%.2f) don't match",dy,dys);
if (nt != nts)
- vmess("Time samples in shot (%d) and focusing operator (%d) are not equal",nt, nts);
+ vmess("Time samples in shot (%li) and focusing operator (%li) are not equal",nt, nts);
if (verbose) {
- vmess("Number of focusing operators = %d", Nfoc);
- vmess("Number of receivers in focusop = x:%d y:%d total:%d", nxs, nys, nxs*nys);
- vmess("number of shots = %d", nshots);
- vmess("number of receiver/shot = x:%d y:%d total:%d", nx, ny, nx*ny);
+ vmess("Number of focusing operators = %li", Nfoc);
+ vmess("Number of receivers in focusop = x:%li y:%li total:%li", nxs, nys, nxs*nys);
+ vmess("number of shots = %li", nshots);
+ vmess("number of receiver/shot = x:%li y:%li total:%li", nx, ny, nx*ny);
vmess("first model position = x:%.2f y:%.2f", fxsb, fysb);
vmess("last model position = x:%.2f y:%.2f", fxse, fyse);
vmess("first source position = x:%.2f y:%.2f", fxf, fyf);
vmess("source distance = x:%.2f y:%.2f", dxsrc, dysrc);
vmess("last source position = x:%.2f y:%.2f", fxf+(nxshot-1)*dxsrc, fyf+(nyshot-1)*dysrc);
vmess("receiver distance = x:%.2f y:%.2f", dxf, dyf);
- vmess("direction of increasing traces = x:%d y:%d", dxi, dyi);
- vmess("number of time samples (nt,nts) = %d (%d,%d)", ntfft, nt, nts);
+ vmess("direction of increasing traces = x:%li y:%li", dxi, dyi);
+ vmess("number of time samples (nt,nts) = %li (%li,%li)", ntfft, nt, nts);
vmess("time sampling = %e ", dt);
if (file_green != NULL) vmess("Green output file = %s ", file_green);
if (file_gmin != NULL) vmess("Gmin output file = %s ", file_gmin);
@@ -437,8 +442,8 @@ int main (int argc, char **argv)
}
mem = Nfoc*n2out*n3out*ntfft*sizeof(float)/1048576.0;
if (verbose) {
- vmess("number of output traces = x:%d y:%d total:%d", n2out, n3out, n2out*n3out);
- vmess("number of output samples = %d", ntfft);
+ vmess("number of output traces = x:%li y:%li total:%li", n2out, n3out, n2out*n3out);
+ vmess("number of output samples = %li", ntfft);
vmess("Size of output data/file = %.1f MB", mem);
}
@@ -447,7 +452,7 @@ int main (int argc, char **argv)
synthesisPositions3D(nx, ny, nxs, nys, Nfoc, xrcv, yrcv, xsrc, ysrc, xnx,
fxse, fyse, fxsb, fysb, dxs, dys, nshots, nxshot, nyshot, ixpos, &npos, reci, verbose);
if (verbose) {
- vmess("synthesisPosistions: nxshot=%d nyshot=%d nshots=%d npos=%d", nxshot, nyshot, nshots, npos);
+ vmess("synthesisPosistions: nxshot=%li nyshot=%li nshots=%li npos=%li", nxshot, nyshot, nshots, npos);
}
/*================ set variables for output data ================*/
@@ -545,12 +550,12 @@ int main (int argc, char **argv)
}
}
if (iter==0) energyN0 = energyNi;
- if (verbose >=2) vmess(" - iSyn %d: Ni at iteration %d has energy %e; relative to N0 %e",
+ if (verbose >=2) vmess(" - iSyn %li: Ni at iteration %li has energy %e; relative to N0 %e",
l, iter, sqrt(energyNi), sqrt(energyNi/energyN0));
}
/* apply mute window based on times of direct arrival (in muteW) */
- applyMute(Ni, muteW, smooth, above, Nfoc, nxs*nys, nts, ixpos, npos, shift);
+ applyMute3D(Ni, muteW, smooth, above, Nfoc, nxs*nys, nts, ixpos, npos, shift);
/* update f2 */
for (l = 0; l < Nfoc; l++) {
@@ -589,12 +594,12 @@ int main (int argc, char **argv)
t2 = wallclock_time();
tcopy += t2 - t3;
- if (verbose) vmess("*** Iteration %d finished ***", iter);
+ if (verbose) vmess("*** Iteration %li finished ***", iter);
} /* end of iterations */
free(Ni);
- free(G_d);
+ if (ampest < 1) free(G_d);
/* compute full Green's function G = int R * f2(t) + f2(-t) = Pplus + Pmin */
for (l = 0; l < Nfoc; l++) {
@@ -611,7 +616,7 @@ int main (int argc, char **argv)
}
}
}
- applyMute(green, muteW, smooth, 4, Nfoc, nxs*nys, nts, ixpos, npos, shift);
+ applyMute3D(green, muteW, smooth, 4, Nfoc, nxs*nys, nts, ixpos, npos, shift);
/* compute upgoing Green's function G^+,- */
if (file_gmin != NULL) {
@@ -637,11 +642,11 @@ int main (int argc, char **argv)
}
}
/* Apply mute with window for Gmin */
- applyMute(Gmin, muteW, smooth, 4, Nfoc, nxs*nys, nts, ixpos, npos, shift);
+ applyMute3D(Gmin, muteW, smooth, 4, Nfoc, nxs*nys, nts, ixpos, npos, shift);
} /* end if Gmin */
/* compute downgoing Green's function G^+,+ */
- if (file_gplus != NULL) {
+ if (file_gplus != NULL || ampest > 0) {
Gplus = (float *)calloc(Nfoc*nys*nxs*ntfft,sizeof(float));
/* use f1-(*) as operator on R in frequency domain */
@@ -664,9 +669,33 @@ int main (int argc, char **argv)
}
}
/* Apply mute with window for Gplus */
- applyMute(Gplus, muteW, smooth, 4, Nfoc, nxs*nys, nts, ixpos, npos, shift);
+ applyMute3D(Gplus, muteW, smooth, 4, Nfoc, nxs*nys, nts, ixpos, npos, shift);
} /* end if Gplus */
+ /* Estimate the amplitude of the Marchenko Redatuming */
+ if (ampest>0) {
+ if (verbose>0) vmess("Estimating amplitude scaling");
+ Gd = (float *)calloc(Nfoc*nxs*nys*ntfft,sizeof(float));
+ memcpy(Gd,Gplus,sizeof(float)*Nfoc*nxs*nys*ntfft);
+ applyMute3D(Gd, muteW, smooth, 2, Nfoc, nxs*nys, nts, ixpos, npos, shift);
+ ampscl = (float *)calloc(Nfoc,sizeof(float));
+ AmpEst3D(G_d,Gd,ampscl,Nfoc,nxs,nys,ntfft,ixpos,npos,file_wav,dxs,dys,dt);
+ for (l=0; l<Nfoc; l++) {
+ for (j=0; j<nxs*nys*nts; j++) {
+ green[l*nxs*nts+j] *= ampscl[l];
+ if (file_gplus != NULL) Gplus[l*nxs*nys*nts+j] *= ampscl[l];
+ if (file_gmin != NULL) Gmin[l*nxs*nys*nts+j] *= ampscl[l];
+ if (file_f2 != NULL) f2p[l*nxs*nys*nts+j] *= ampscl[l];
+ if (file_pmin != NULL) pmin[l*nxs*nys*nts+j] *= ampscl[l];
+ if (file_f1plus != NULL) f1plus[l*nxs*nys*nts+j] *= ampscl[l];
+ if (file_f1min != NULL) f1min[l*nxs*nys*nts+j] *= ampscl[l];
+ }
+ if (verbose>1) vmess("Amplitude of focal position %li is equal to %.3e",l,ampscl[l]);
+ }
+ free(Gd);
+ if (file_gplus == NULL) free(Gplus);
+ }
+
t2 = wallclock_time();
if (verbose) {
vmess("Total CPU-time marchenko = %.3f", t2-t0);
@@ -731,23 +760,23 @@ int main (int argc, char **argv)
}
}
- ret = writeData(fp_out, (float *)&green[l*size], hdrs_out, n1, n2*n3);
+ ret = writeData3D(fp_out, (float *)&green[l*size], hdrs_out, n1, n2*n3);
if (ret < 0 ) verr("error on writing output file.");
if (file_gmin != NULL) {
- ret = writeData(fp_gmin, (float *)&Gmin[l*size], hdrs_out, n1, n2*n3);
+ ret = writeData3D(fp_gmin, (float *)&Gmin[l*size], hdrs_out, n1, n2*n3);
if (ret < 0 ) verr("error on writing output file.");
}
if (file_gplus != NULL) {
- ret = writeData(fp_gplus, (float *)&Gplus[l*size], hdrs_out, n1, n2*n3);
+ ret = writeData3D(fp_gplus, (float *)&Gplus[l*size], hdrs_out, n1, n2*n3);
if (ret < 0 ) verr("error on writing output file.");
}
if (file_f2 != NULL) {
- ret = writeData(fp_f2, (float *)&f2p[l*size], hdrs_out, n1, n2*n3);
+ ret = writeData3D(fp_f2, (float *)&f2p[l*size], hdrs_out, n1, n2*n3);
if (ret < 0 ) verr("error on writing output file.");
}
if (file_pmin != NULL) {
- ret = writeData(fp_pmin, (float *)&pmin[l*size], hdrs_out, n1, n2*n3);
+ ret = writeData3D(fp_pmin, (float *)&pmin[l*size], hdrs_out, n1, n2*n3);
if (ret < 0 ) verr("error on writing output file.");
}
if (file_f1plus != NULL) {
@@ -762,7 +791,7 @@ int main (int argc, char **argv)
f1plus[l*size+i*nts+j-n1/2] = trace[j];
}
}
- ret = writeData(fp_f1plus, (float *)&f1plus[l*size], hdrs_out, n1, n2*n3);
+ ret = writeData3D(fp_f1plus, (float *)&f1plus[l*size], hdrs_out, n1, n2*n3);
if (ret < 0 ) verr("error on writing output file.");
}
if (file_f1min != NULL) {
@@ -777,7 +806,7 @@ int main (int argc, char **argv)
f1min[l*size+i*nts+j-n1/2] = trace[j];
}
}
- ret = writeData(fp_f1min, (float *)&f1min[l*size], hdrs_out, n1, n2*n3);
+ ret = writeData3D(fp_f1min, (float *)&f1min[l*size], hdrs_out, n1, n2*n3);
if (ret < 0 ) verr("error on writing output file.");
}
}
@@ -788,7 +817,7 @@ int main (int argc, char **argv)
if (file_pmin != NULL) {ret += fclose(fp_pmin);}
if (file_f1plus != NULL) {ret += fclose(fp_f1plus);}
if (file_f1min != NULL) {ret += fclose(fp_f1min);}
- if (ret < 0) verr("err %d on closing output file",ret);
+ if (ret < 0) verr("err %li on closing output file",ret);
if (verbose) {
t1 = wallclock_time();
@@ -803,11 +832,11 @@ int main (int argc, char **argv)
exit(0);
}
-int unique_elements(float *arr, int len)
+long unique_elements(float *arr, long len)
{
if (len <= 0) return 0;
- int unique = 1;
- int outer, inner, is_unique;
+ long unique = 1;
+ long outer, inner, is_unique;
for (outer = 1; outer < len; ++outer)
{
diff --git a/marchenko3D/readData3D.c b/marchenko3D/readData3D.c
index 09b121a3c5f5ae311648fc55fa666d2b887bbdfc..8cef449aa345e68cc85413f822885e530d4a5f04 100644
--- a/marchenko3D/readData3D.c
+++ b/marchenko3D/readData3D.c
@@ -13,10 +13,10 @@
**/
-int readData3D(FILE *fp, float *data, segy *hdrs, int n1)
+int readData3D(FILE *fp, float *data, segy *hdrs, long n1)
{
size_t nread;
- int oneshot, itrace, sx, sy, fldr, gy, nx, ny;
+ long oneshot, itrace, sx, sy, fldr, gy, nx, ny;
segy hdr;
oneshot = 1;
diff --git a/marchenko3D/readShotData3D.c b/marchenko3D/readShotData3D.c
index f9c569ad4c44999bd820f57e59d77bcb3d52ce6f..c25327fd9a2b7be3e4683a19a6c5d7b61651a1d5 100644
--- a/marchenko3D/readShotData3D.c
+++ b/marchenko3D/readShotData3D.c
@@ -9,22 +9,22 @@ typedef struct { /* complex number */
float r,i;
} complex;
-#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
#define MAX(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);
-int readShotData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int nshots, int nx, int ny, int ntfft, int mode, float scale, int 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;
size_t nread;
- int fldr_shot, sx_shot, sy_shot, itrace, one_shot, igath, iw;
- int end_of_file, nt, nxy;
- int *isx, *igx, *isy, *igy, k, l, m, j;
- int samercv, samesrc, nxrk, nxrm, maxtraces, ixsrc;
+ long fldr_shot, sx_shot, sy_shot, itrace, one_shot, igath, iw;
+ long end_of_file, nt, nxy;
+ long *isx, *igx, *isy, *igy, k, l, m, j;
+ long samercv, samesrc, nxrk, nxrm, maxtraces, ixsrc;
float scl, scel, *trace, dt;
complex *ctrace;
@@ -58,10 +58,10 @@ int readShotData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
trace = (float *)calloc(ntfft,sizeof(float));
ctrace = (complex *)malloc(ntfft*sizeof(complex));
- isx = (int *)malloc((nxy*nshots)*sizeof(int));
- igx = (int *)malloc((nxy*nshots)*sizeof(int));
- isy = (int *)malloc((nxy*nshots)*sizeof(int));
- igy = (int *)malloc((nxy*nshots)*sizeof(int));
+ isx = (long *)malloc((nxy*nshots)*sizeof(long));
+ igx = (long *)malloc((nxy*nshots)*sizeof(long));
+ isy = (long *)malloc((nxy*nshots)*sizeof(long));
+ igy = (long *)malloc((nxy*nshots)*sizeof(long));
end_of_file = 0;
@@ -101,7 +101,7 @@ int readShotData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
if (ntfft > hdr.ns)
memset( &trace[nt-1], 0, sizeof(float)*(ntfft-nt) );
- rc1fft(trace,ctrace,ntfft,-1);
+ 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;
@@ -119,7 +119,7 @@ int readShotData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
if ((sx_shot != hdr.sx) || (sy_shot != hdr.sy) || (fldr_shot != hdr.fldr)) break;
}
if (verbose>2) {
- vmess("finished reading shot x=%d y=%d (%d) with %d traces",sx_shot,sy_shot,igath,itrace);
+ 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 */
diff --git a/marchenko3D/readTinvData3D.c b/marchenko3D/readTinvData3D.c
index 845d08076656e6987299f3fdbbbcf5d0aa55ec1f..203585d84f2a4584ad774d9dd499f6c22241d05d 100644
--- a/marchenko3D/readTinvData3D.c
+++ b/marchenko3D/readTinvData3D.c
@@ -15,18 +15,18 @@ typedef struct { /* complex number */
#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
-#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
-void findShotInMute(float *xrcvMute, float xrcvShot, int nxs, int *imute);
+void findShotInMute(float *xrcvMute, float xrcvShot, long nxs, long *imute);
-int readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc, int *xnx, int Nfoc, int nx, int ny, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose)
+long readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float *ysrc, float *zsrc, long *xnx, long Nfoc, long nx, long ny, long ntfft, long mode, long *maxval, float *tinv, long hw, long verbose)
{
FILE *fp;
segy hdr;
size_t nread;
- int fldr_shot, sx_shot, sy_shot, itrace, one_shot, ig, isyn, i, j, l;
- int end_of_file, nt, gx0, gx1, gy0, gy1;
- int nx1, ny1, jmax, imax, tstart, tend, nxy, ixmax, iymax;
+ long fldr_shot, sx_shot, sy_shot, itrace, one_shot, ig, isyn, i, j, l;
+ long end_of_file, nt, gx0, gx1, gy0, gy1;
+ long nx1, ny1, jmax, imax, tstart, tend, nxy, ixmax, iymax;
float xmax, tmax, lmax;
float scl, scel, *trace, dxrcv, dyrcv;
complex *ctrace;
@@ -113,7 +113,7 @@ int readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
}
}
if (verbose>2) {
- fprintf(stderr,"finished reading shot x=%d y=%d (%d) with %d traces\n",sx_shot,sy_shot,isyn,itrace);
+ fprintf(stderr,"finished reading shot x=%li y=%li (%li) with %li traces\n",sx_shot,sy_shot,isyn,itrace);
}
/* look for maximum in shot record to define mute window */
@@ -131,19 +131,19 @@ int readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
if (dyrcv==0.0) dyrcv=1.0;
iymax = NINT(((sy_shot-gy0)*scl)/dyrcv);
if (iymax > ny1-1) {
- vmess("source of y (%d) is past array, snapping to nearest y (%d)",iymax,ny1-1);
+ vmess("source of y (%li) is past array, snapping to nearest y (%li)",iymax,ny1-1);
iymax = ny1-1;
}
if (iymax < 0) {
- vmess("source of y (%d) is before array, snapping to nearest y (%d)",iymax,0);
+ vmess("source of y (%li) is before array, snapping to nearest y (%li)",iymax,0);
iymax = 0;
}
if (ixmax > nx1-1) {
- vmess("source of x (%d) is past array, snapping to nearest x (%d)",ixmax,nx1-1);
+ vmess("source of x (%li) is past array, snapping to nearest x (%li)",ixmax,nx1-1);
ixmax = nx1-1;
}
if (ixmax < 0) {
- vmess("source of x (%d) is before array, snapping to nearest x (%d)",ixmax,nx1-1);
+ vmess("source of x (%li) is before array, snapping to nearest x (%li)",ixmax,nx1-1);
ixmax = 0;
}
tmax=0.0;
@@ -159,7 +159,7 @@ int readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
}
}
maxval[isyn*nxy+iymax*nx+ixmax] = jmax;
- if (verbose >= 3) vmess("Mute max at src-trace x=%d y=%d is sample %d", ixmax, iymax, maxval[isyn*nxy+iymax*nx+ixmax]);
+ if (verbose >= 3) vmess("Mute max at src-trace x=%li y=%li is sample %li", ixmax, iymax, maxval[isyn*nxy+iymax*nx+ixmax]);
/* search forward in x-trace direction from maximum in file */
for (i = ixmax+1; i < nx1; i++) {
@@ -204,7 +204,7 @@ int readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
}
}
maxval[isyn*nxy+i*nx+ixmax] = jmax;
- if (verbose >= 8) vmess("Mute max at src-trace x=%d y=%d is sample %d", ixmax, i, maxval[isyn*nxy+i*nx+ixmax]);
+ if (verbose >= 8) vmess("Mute max at src-trace x=%li y=%li is sample %li", ixmax, i, maxval[isyn*nxy+i*nx+ixmax]);
/* search forward in x-trace direction from maximum in file */
for (l = ixmax+1; l < nx1; l++) {
tstart = MAX(0, (maxval[isyn*nxy+i*nx+(l-1)]-hw));
@@ -249,7 +249,7 @@ int readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
}
}
maxval[isyn*nxy+i*nx+ixmax] = jmax;
- if (verbose >= 8) vmess("Mute max at src-trace x=%d y=%d is sample %d", ixmax, i, maxval[isyn*nxy+i*nx+ixmax]);
+ if (verbose >= 8) vmess("Mute max at src-trace x=%li y=%li is sample %li", ixmax, i, maxval[isyn*nxy+i*nx+ixmax]);
/* search forward in x-trace direction from maximum in file */
for (l = ixmax+1; l < nx1; l++) {
tstart = MAX(0, (maxval[isyn*nxy+i*nx+(l-1)]-hw));
@@ -312,9 +312,9 @@ int readTinvData3D(char *filename, float *xrcv, float *yrcv, float *xsrc, float
/* simple sort algorithm */
-void findShotInMute(float *xrcvMute, float xrcvShot, int nxs, int *imute)
+void findShotInMute(float *xrcvMute, float xrcvShot, long nxs, long *imute)
{
- int i, sign;
+ long i, sign;
float diff1, diff2;
*imute=0;
diff --git a/marchenko3D/synthesis3D.c b/marchenko3D/synthesis3D.c
index 19d4c80b2ade811d3afb0b139a35d68cb34aafc5..d93200daee1113418bafc9dc823d8c253eca1879 100644
--- a/marchenko3D/synthesis3D.c
+++ b/marchenko3D/synthesis3D.c
@@ -17,9 +17,9 @@ void omp_set_num_threads(int num_threads);
#ifndef MIN
#define MIN(x,y) ((x) < (y) ? (x) : (y))
#endif
-#define NINT(x) ((int)((x)>0.0?(x)+0.5:(x)-0.5))
+#define NINT(x) ((long)((x)>0.0?(x)+0.5:(x)-0.5))
int compareInt(const void *a, const void *b)
-{ return (*(int *)a-*(int *)b); }
+{ return (*(long *)a-*(long *)b); }
#ifndef COMPLEX
typedef struct _complexStruct { /* complex number */
@@ -27,11 +27,13 @@ typedef struct _complexStruct { /* complex number */
} complex;
#endif/* complex */
-void synthesisPositions3D(int nx, int ny, int nxs, int nys, int Nfoc, float *xrcv, float *yrcv,
-float *xsrc, float *ysrc, int *xnx, float fxse, float fyse, float fxsb, float fysb, float dxs, float dys,
-int nshots, int nxsrc, int nysrc, int *ixpos, int *npos, int reci, int verbose)
+long linearsearch(long *array, size_t N, long value);
+
+void synthesisPositions3D(long nx, long ny, long nxs, long nys, long Nfoc, float *xrcv, float *yrcv,
+float *xsrc, float *ysrc, long *xnx, float fxse, float fyse, float fxsb, float fysb, float dxs, float dys,
+long nshots, long nxsrc, long nysrc, long *ixpos, long *npos, long reci, long verbose)
{
- int j, l, ixsrc, iysrc, isrc, k, *count, nxy;
+ long j, l, ixsrc, iysrc, isrc, k, *count, nxy;
float fxb, fxe, fyb, fye;
if (fxsb < 0) fxb = 1.001*fxsb;
@@ -45,7 +47,7 @@ int nshots, int nxsrc, int nysrc, int *ixpos, int *npos, int reci, int verbose)
nxy = nx*ny;
- count = (int *)calloc(nxs*nys,sizeof(int)); // number of traces that contribute to the integration over x
+ count = (long *)calloc(nxs*nys,sizeof(long)); // number of traces that contribute to the integration over x
/*================ SYNTHESIS ================*/
@@ -60,7 +62,7 @@ int nshots, int nxsrc, int nysrc, int *ixpos, int *npos, int reci, int verbose)
iysrc = NINT((ysrc[k] - fysb)/dys);
isrc = iysrc*nxs + ixsrc;
if (verbose>=3) {
- vmess("source position: x=%.2f y=%.2f in operator x=%d y=%d pos=%d", xsrc[k], ysrc[k], ixsrc, iysrc, isrc);
+ vmess("source position: x=%.2f y=%.2f in operator x=%li y=%li pos=%li", xsrc[k], ysrc[k], ixsrc, iysrc, isrc);
vmess("receiver positions: x:%.2f <--> %.2f y:%.2f <--> %.2f", xrcv[k*nxy+0], xrcv[k*nxy+nxy-1], yrcv[k*nxy+0], yrcv[k*nxy+nxy-1]);
vmess("focal point positions: x:%.2f <--> %.2f y:%.2f <--> %.2f", fxsb, fxse, fysb, fyse);
}
@@ -74,8 +76,8 @@ int nshots, int nxsrc, int nysrc, int *ixpos, int *npos, int reci, int verbose)
vwarn("source/receiver positions are outside synthesis aperture");
vmess("xsrc = %.2f xrcv_1 = %.2f xrvc_N = %.2f", xsrc[k], xrcv[k*nxy+0], xrcv[k*nxy+nxy-1]);
vmess("ysrc = %.2f yrcv_1 = %.2f yrvc_N = %.2f", ysrc[k], yrcv[k*nxy+0], yrcv[k*nxy+nxy-1]);
- vmess("source position x: %.2f in operator %d", xsrc[k], ixsrc);
- vmess("source position y: %.2f in operator %d", ysrc[k], iysrc);
+ vmess("source position x: %.2f in operator %li", xsrc[k], ixsrc);
+ vmess("source position y: %.2f in operator %li", ysrc[k], iysrc);
vmess("receiver positions x: %.2f <--> %.2f", xrcv[k*nxy+0], xrcv[k*nxy+nxy-1]);
vmess("receiver positions y: %.2f <--> %.2f", yrcv[k*nxy+0], yrcv[k*nxy+nxy-1]);
vmess("focal point positions x: %.2f <--> %.2f", fxsb, fxse);
@@ -94,7 +96,7 @@ int nshots, int nxsrc, int nysrc, int *ixpos, int *npos, int reci, int verbose)
count[*npos] += xnx[k];
*npos += 1;
}
- // vmess("source position %d is inside synthesis model %f *npos=%d count=%d", k, xsrc[k], *npos, count[*npos]);
+ // vmess("source position %li is inside synthesis model %f *npos=%li count=%li", k, xsrc[k], *npos, count[*npos]);
}
} /* end of nshots (k) loop */
@@ -103,21 +105,21 @@ int nshots, int nxsrc, int nysrc, int *ixpos, int *npos, int reci, int verbose)
if (verbose>=4) {
for (j=0; j < *npos; j++) {
- vmess("ixpos[%d] = %d count=%d", j, ixpos[j], count[j]);
+ vmess("ixpos[%li] = %li count=%li", j, ixpos[j], count[j]);
}
}
free(count);
/* sort ixpos into increasing values */
- qsort(ixpos, *npos, sizeof(int), compareInt);
+ qsort(ixpos, *npos, sizeof(long), compareInt);
return;
}
-int linearsearch(int *array, size_t N, int value)
+long linearsearch(long *array, size_t N, long value)
{
- int j;
+ long j;
/* Check is position is already in array */
j = 0;
while (j < N && value != array[j]) {
@@ -128,18 +130,18 @@ int linearsearch(int *array, size_t N, int value)
/*================ Convolution and Integration ================*/
-void synthesis3D(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int ny, int nt, int nxs, int nys, int nts, float dt, float *xsyn, float *ysyn,
-int Nfoc, float *xrcv, float *yrcv, float *xsrc, float *ysrc, int *xnx, float fxse, float fxsb, float fyse, float fysb, float dxs, float dys, float dxsrc,
-float dysrc, float dx, float dy, int ntfft, int nw, int nw_low, int nw_high, int mode, int reci, int nshots, int nxsrc, int nysrc,
-int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xrcv, float *ixmask, int verbose)
+void synthesis3D(complex *Refl, complex *Fop, float *Top, float *iRN, long nx, long ny, long nt, long nxs, long nys, long nts, float dt, float *xsyn, float *ysyn,
+long Nfoc, float *xrcv, float *yrcv, float *xsrc, float *ysrc, long *xnx, float fxse, float fxsb, float fyse, float fysb, float dxs, float dys, float dxsrc,
+float dysrc, float dx, float dy, long ntfft, long nw, long nw_low, long nw_high, long mode, long reci, long nshots, long nxsrc, long nysrc,
+long *ixpos, long npos, double *tfft, long *isxcount, long *reci_xsrc, long *reci_xrcv, float *ixmask, long verbose)
{
- int nfreq, size, inx;
+ long nfreq, size, inx;
float scl;
- int i, j, l, m, iw, ix, k, isrc, il, ik, nxy, nxys;
+ long i, j, l, m, iw, ix, k, isrc, il, ik, nxy, nxys;
float *rtrace, idxs, idys, fxb, fyb, fxe, fye;
complex *sum, *ctrace;
- int npe;
- static int first=1, *ircv;
+ long npe;
+ static long first=1, *ircv;
static double t0, t1, t;
if (fxsb < 0) fxb = 1.001*fxsb;
@@ -162,11 +164,11 @@ int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xr
scl = 1.0*dt/((float)ntfft);
#ifdef _OPENMP
- npe = omp_get_max_threads();
+ npe = (long)omp_get_max_threads();
/* parallelisation is over number of virtual source positions (Nfoc) */
if (npe > Nfoc) {
- vmess("Number of OpenMP threads set to %d (was %d)", Nfoc, npe);
- omp_set_num_threads(Nfoc);
+ vmess("Number of OpenMP threads set to %li (was %li)", Nfoc, npe);
+ omp_set_num_threads((int)Nfoc);
}
#endif
@@ -207,7 +209,7 @@ int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xr
}
idxs = 1.0/dxs;
idys = 1.0/dys;
- ircv = (int *)malloc(nshots*nxy*sizeof(int));
+ ircv = (long *)malloc(nshots*nxy*sizeof(long));
for (k=0; k<nshots; k++) {
for (i = 0; i < nxy; i++) {
ircv[k*nxy+i] = NINT((yrcv[k*nxy+i]-fysb)*idys)*nx+NINT((xrcv[k*nxy+i]-fxsb)*idxs);
@@ -266,18 +268,18 @@ int *ixpos, int npos, double *tfft, int *isxcount, int *reci_xsrc, int *reci_xr
#ifdef _OPENMP
#pragma omp single
- npe = omp_get_num_threads();
+ npe = (long)omp_get_num_threads();
#endif
} /* end of parallel region */
- if (verbose>4) vmess("*** Shot gather %d processed ***", k);
+ if (verbose>4) vmess("*** Shot gather %li processed ***", k);
} /* end of nshots (k) loop */
} /* end of if reci */
t = wallclock_time() - t0;
if (verbose) {
- vmess("OMP: parallel region = %f seconds (%d threads)", t, npe);
+ vmess("OMP: parallel region = %f seconds (%li threads)", t, npe);
}
return;