diff --git a/fdelmodc/fdelmodc.c b/fdelmodc/fdelmodc.c
index 878f6b61edb35ca53052e38db12355786b7a8e06..f780b78fadff0c9c9a1937dc45663e2df6a2739a 100644
--- a/fdelmodc/fdelmodc.c
+++ b/fdelmodc/fdelmodc.c
@@ -132,8 +132,8 @@ char *sdoc[] = {
" - 1=monopole",
" - 2=dipole +/- vertical oriented",
" - 3=dipole - + horizontal oriented",
-//" - 4=dipole +/0/-",
-//" - 5=dipole + -",
+" - 4=dipole +/0/-",
+" - 5=dipole + -",
" dip=0 ............. dip for double-couple source",
" strike=0 .......... strike for double-couple source",
" xsrc=middle ....... x-position of (first) shot ",
diff --git a/marchenko/Makefile b/marchenko/Makefile
index cd1c5416851cc624e7455be23ff1bf08c6a6a00a..704ba08d7723f9898e33a63a075ec70ce5d270c3 100644
--- a/marchenko/Makefile
+++ b/marchenko/Makefile
@@ -37,7 +37,7 @@ SRCH = marchenko.c \
getpars.c
SRCP = marchenko_primaries.c \
- synthesis.c \
+ synthesisp.c \
getFileInfo.c \
readData.c \
readShotData.c \
diff --git a/marchenko/demo/oneD/README_PRIMARIES b/marchenko/demo/oneD/README_PRIMARIES
index 5eda7197ef806033b2f1286d28a433d9d53b2c64..7b97347bbbf32380e50763d094dd84853cf14c60 100644
--- a/marchenko/demo/oneD/README_PRIMARIES
+++ b/marchenko/demo/oneD/README_PRIMARIES
@@ -2,10 +2,9 @@ Description of files:
1) model.scr computes the model and the 'basis' shot of R => shot5_rp.su
-2) p5all.scr create from basis shot full Reflection response matrix => shotsdx5_rp.su (3.3 GB)
-3) primaries.scr computes the internal multiple attenuated (middle) shot record.
-4) itertions.scr computes the intermediate results of the multiple attenutation scheme and produces all output files that are used in the manuscript.
-5) epsPrimaries.scr selected output form 4) are converted to .eps pictures that are used in the Figures to explain the method.
+2) primaries.scr computes the internal multiple attenuated (middle) shot record.
+3) itertions.scr computes the intermediate results of the multiple attenutation scheme and produces all output files that are used in the manuscript.
+4) epsPrimaries.scr selected output from 3) are converted to .eps pictures that are used in the Figures to explain the method.
reproduce the postscript files of the manuscript using SU postscript plotting programs.
extra scripts
diff --git a/marchenko/demo/oneD/epsPrimaries.scr b/marchenko/demo/oneD/epsPrimaries.scr
index 00949d8fe8d62007ba0bae6ec5dda71b5bf1ef44..2d05f4cb0f9dc81d988f09dbacfaf6f916ada837 100755
--- a/marchenko/demo/oneD/epsPrimaries.scr
+++ b/marchenko/demo/oneD/epsPrimaries.scr
@@ -133,6 +133,40 @@ supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
done
fi
+if [[ "$1" == "inaction" ]];
+then
+rm vplus.su
+istart=200
+for (( iter=1; iter<=29; iter+=2 ))
+do
+piter=$(printf %03d $iter)
+echo $piter
+file=f1plus_${istart}${piter}.su
+file_base=${file%.su}
+suwind < $file key=offset min=0 max=0 >> vplus.su
+done
+
+basop file_in=vplus.su file_out=vplus_env.su choice=4
+sumax < vplus_env.su verbose=1 mode=max
+supsmax < vplus_env.su n=15 verbose=1 mode=max \
+ style=normal linewidth=2.0 f1=1 labelsize=14 label1="iteration number" label2="amplitude" \
+ d1num=2 d2num=0.2 wbox=6 hbox=3 grid1=dot grid2=dot n1tic=2 n2tic=2 x2end=2.9 x2beg=1.8 > v1plus_max.eps
+
+file=f1plus_${istart}001.su
+file_base=${file%.su}
+supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
+ n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
+ label1="time sample number" label2="lateral distance" \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipiter > $file_base.eps
+
+file=f1min_${istart}030.su
+file_base=${file%.su}
+supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 < $file\
+ n1tic=2 d2=5 x1beg=0 x1end=$nsd d1=1 d1num=100 \
+ label1="time sample number" label2="lateral distance" \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > $file_base.eps
+fi
+
if [[ "$1" == "Figure5" ]];
then
for (( iter=1; iter<=21; iter+=2 ))
@@ -148,7 +182,8 @@ supsimage < $file hbox=6 wbox=4 labelsize=14 linewidth=0.0 \
done
fi
-exit;
+if [[ "$1" == "Figure7" ]];
+then
#iterations
for (( iter=0; iter<=31; iter+=2 ))
do
@@ -174,7 +209,10 @@ supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 < $file\
f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > $file_base.eps
done
+fi
+if [[ "$1" == "Figure8" ]];
+then
iter=32
piter=$(printf %03d $iter)
echo $piter
@@ -187,6 +225,7 @@ supsimage hbox=6 wbox=4 labelsize=14 linewidth=0.0 < $file\
label1="time sample number" label2="lateral distance" \
f2=-2250 f2num=-2000 d2num=1000 clip=$clipf1 > $file_base.eps
done
+fi
exit;
diff --git a/marchenko/demo/oneD/iterations.scr b/marchenko/demo/oneD/iterations.scr
index 3cc1ec7981379ceab761f5aed388890c6c2246fa..f6f774f27e0b1c55824ab5665dedad77b22e5156 100755
--- a/marchenko/demo/oneD/iterations.scr
+++ b/marchenko/demo/oneD/iterations.scr
@@ -3,10 +3,16 @@
#second reflector at time:
# 800/1800+600/2300
# .70531400966183574878 => sample 176
-#third reflector at time:
+#third reflector at time model.scr:
# 800/1800+600/2300+800/2000
# 1.10531400966183574878 sample 276
+#third reflector at time modepm.scr
+#800/1800+600/3200+520/2000
+#.96531400966183574878 sample 241
+
+R=shotsdx5_rp.su
+makewave fp=20 dt=0.004 file_out=wave.su nt=2048 t0=0.0 scale=0 scfft=1
select=451
if [[ "$1" == "Figure3467" ]];
@@ -15,7 +21,7 @@ for istart in 276
do
(( iend = istart + 1 ))
-../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
+../../marchenko_primaries file_shot=$R ishot=$select file_src=wave.su \
nshots=901 verbose=2 istart=$istart iend=$iend fmax=90 \
pad=512 niter=33 smooth=8 niterskip=600 shift=16 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su
done
@@ -27,17 +33,26 @@ for istart in 246 256 266 276 286 296 306 316
do
(( iend = istart + 1 ))
-../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
+../../marchenko_primaries file_shot=$R ishot=$select file_src=wave.su \
nshots=901 verbose=2 istart=$istart iend=$iend fmax=90 \
pad=512 niter=33 smooth=8 niterskip=600 shift=16 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su
done
fi
+if [[ "$1" == "inaction" ]];
+then
+istart=200
+(( iend = istart + 1 ))
+../../marchenko_primaries file_shot=$R ishot=$select file_src=wave.su \
+ nshots=901 verbose=2 istart=$istart iend=$iend fmax=90 \
+ pad=512 niter=30 smooth=8 niterskip=600 shift=16 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su
+fi
+
if [[ "$1" == "Figure10" ]];
then
istart=276
(( iend = istart + 1 ))
-../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
+../../marchenko_primaries file_shot=$R ishot=$select file_src=wave.su \
nshots=901 verbose=2 istart=$istart iend=$iend fmax=90 \
pad=512 niter=33 smooth=8 niterskip=600 shift=16 file_rr=pred_rr_${istart}.su T=1 file_iter=iterT.su
fi
diff --git a/marchenko/demo/oneD/model.scr b/marchenko/demo/oneD/model.scr
index 9165e8758862d4cef0b497ca585bea46336e9d13..1a5e0de133ac37005fe03cef01bea23b609ae04d 100755
--- a/marchenko/demo/oneD/model.scr
+++ b/marchenko/demo/oneD/model.scr
@@ -75,3 +75,27 @@ fdelmodc \
sudiff shot5_fd_rp.su shot5_hom_fd_rp.su > shot5_rp.su
+# Generate the full R matrix for a fixed spread geometry.
+
+dxshot=5000 # with scalco factor of 1000
+ishot=0
+nshots=901
+
+rm shotsdx5_rp.su
+
+while (( ishot < nshots ))
+do
+
+ (( xsrc = -2250000 + ${ishot}*${dxshot} ))
+ (( tr1 = $nshots - ${ishot} ))
+ (( tr2 = ${tr1} + $nshots - 1 ))
+ echo xsrc=$xsrc tr1=$tr1 tr2=$tr2
+
+ (( ishot = $ishot + 1))
+
+ 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
+
+done
diff --git a/marchenko/demo/oneD/primaries.scr b/marchenko/demo/oneD/primaries.scr
index b5398efee445d97e3c0a7ca5e7a326c9b3fb01c5..ca454a6c11f8a178e95ddb812bd6f2504ebbbc42 100755
--- a/marchenko/demo/oneD/primaries.scr
+++ b/marchenko/demo/oneD/primaries.scr
@@ -4,6 +4,8 @@
#SBATCH --ntasks=1
#SBATCH --time=2:40:00
+# Generate the model and Reflection data
+#modelpm.scr
#cd $SLURM_SUBMIT_DIR
@@ -11,18 +13,19 @@ export PATH=$HOME/src/OpenSource/bin:$PATH:
export OMP_NUM_THREADS=10
#shot record to remove internal multiples
+R=shotsdx5_rp.su
select=451
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 \
+../../marchenko_primaries file_shot=$R ishot=$select file_src=wave.su \
nshots=901 verbose=2 istart=40 iend=500 fmax=90 pad=1024 \
niter=31 smooth=10 niterec=2 niterskip=50 shift=20 file_rr=pred_rr.su T=0 file_update=update.su
exit;
#for reference original shot record from Reflection matrix
-suwind key=fldr min=$select max=$select < shotsdx5_rp.su itmax=2047 > shotsx0.su
+suwind key=fldr min=$select max=$select < $R itmax=2047 > shotsx0.su
fconv file_in1=shotsx0.su file_in2=wave.su file_out=shotw.su
# for displaying results
diff --git a/marchenko/marchenko_primaries.c b/marchenko/marchenko_primaries.c
index d893a330715027cf0e5069297341b693a8a82bc0..b6775fdcf8c988aa26d3783b6b67f92e610c76d9 100644
--- a/marchenko/marchenko_primaries.c
+++ b/marchenko/marchenko_primaries.c
@@ -37,7 +37,7 @@ 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
+void synthesisp(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);
@@ -72,7 +72,6 @@ char *sdoc[] = {
" src_angle=0 .............. angle with horizontal of plane source array",
" src_velo=1500 ............ velocity to use in src_angle definition",
" t0=0.1 ................... time shift in plane-wave source wavelet for migration",
-//" xorig=nrecv/2 ............ center of plane-wave",
" MARCHENKO ITERATIONS ",
" niter=22 ................. number of iterations to initialize and restart",
" niterec=2 ................ number of iterations in recursive part of the time-samples",
@@ -92,16 +91,17 @@ char *sdoc[] = {
//" countmin=0 ............... 0.3*nxrcv; minimum number of reciprocal traces for a contribution",
" OUTPUT DEFINITION ",
" file_rr= ................. output file with primary only shot record",
+" file_dd= ................. output file with input of the algorithm ",
" file_iter= ............... output file with -Mi(-t) for each iteration: writes",
" ............... M0.su=M0 : initialisation of algorithm",
" ............... RMi: iterative terms ",
-" ............... f1min.su: f1min terms ",
+" ............... u1min.su: u1min terms ",
" file_update= ............. output file with updates only => removed internal multiples",
" T=0 ...................... :1 compute transmission-losses compensated primaries ",
" verbose=0 ................ silent option; >0 displays info",
" ",
" ",
-" author : Lele Zhang & Jan Thorbecke : 2019 ",
+" author : Lele Zhang & Jan Thorbecke : 2020 ",
" ",
NULL};
/**************** end self doc ***********************************/
@@ -109,27 +109,28 @@ NULL};
int main (int argc, char **argv)
{
FILE *fp_out, *fp_rr, *fp_w, *fp_up;
- size_t nread;
- int i, j, k, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath;
- int size, n1, n2, ntap, tap, di, ntraces, tr;
+ size_t nread, sizea, size;
+ int i, j, k, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath, nacq;
+ int n1, n2, ntap, tap, di, ntraces, tr;
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
- int reci, countmin, mode, ixa, ixb, n2out, verbose, ntfft;
+ int reci, countmin, mode, n2out, verbose, ntfft;
int iter, niter, niterec, recur, niterskip, niterrun, tracf, *muteW;
int hw, ii, iw, ishot, istart, iend;
- int smooth, *ixpos, npos, ix, m, pad, T, isms, isme, perc;
- int nshots_r, *isxcount, *reci_xsrc, *reci_xrcv, shift, plane_wave, xorig;
+ int smooth, *ixpos, *ixp, npos, ix, ixrcv, m, pad, T, isms, isme, perc;
+ int nshots_r, *isxcount, *reci_xsrc, *reci_xrcv, shift, plane_wave;
float fmin, fmax, tom, deltom, *tapersh, *tapersy, fxf, dxf, *xsrc, *xrcv, *zsyn, *zsrc, *xrcvsyn;
double t0, t1, t2, t3, t4, tsyn, tread, tfft, tcopy, tii;
double energyMi, *energyM0;
float tt0, d1, d2, f1, f2, fxsb, fxse, ft, fx, *xsyn, dxsrc;
float *M0, *DD, *RR, *SRC, dt, dx, dxs, scl, mem, scltap;
- float *rtrace, *tmpdata, *f1min, *f1plus, *RMi, *Mi, *trace;
+ float *rtrace, *tmpdata, *u1min, *v1plus, *RMi, *Mi, *trace;
float *Mup, *Msp, *maxval;
float xmin, xmax, scale, tsq;
float Q, f0, *ixmask, *costaper;
float src_velo, src_angle, grad2rad, p, *twplane;
complex *Refl, *Fop, *ctrace, *cwave, csum, cwav;
char *file_tinv, *file_shot, *file_rr, *file_src, *file_iter, *file_update;
+ char *file_dd;
segy *hdrs_out, hdr;
initargs(argc, argv);
@@ -142,6 +143,7 @@ int main (int argc, char **argv)
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_dd", &file_rr)) file_dd = NULL;
if (!getparstring("file_iter", &file_iter)) file_iter = NULL;
if (!getparstring("file_update", &file_update)) file_update = NULL;
@@ -188,7 +190,7 @@ int main (int argc, char **argv)
smooth = MIN(smooth, shift);
if (smooth) {
costaper = (float *)malloc(smooth*sizeof(float));
- scl = M_PI/((float)smooth);
+ scl = M_PI/((float)smooth+1);
for (i=0; i<smooth; i++) {
costaper[i] = 0.5*(1.0+cos((i+1)*scl));
}
@@ -200,6 +202,7 @@ int main (int argc, char **argv)
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;
+ nacq = n2;
nxs = n2;
nts = n1;
dxs = d2;
@@ -220,27 +223,28 @@ int main (int argc, char **argv)
if(!getparint("ishot", &ishot)) ishot=1+(nshots-1)/2;
ishot -= 1; /* shot numbering starts at 0 */
Nfoc = 1;
- nxs = nx;
+ nacq = NINT((xmax-xmin)/dx)+1;
+ nxs = nshots;
nts = nt;
dxs = dx;
- fxsb = fx;
+ fxsb = xmin;
+ fxse = xmax;
}
- assert (nxs >= nshots); /* ToDo allow other geometries */
- //if(!getparint("xorig", &xorig)) xorig=-(nxs-1)/2;
+ assert (nacq >= nshots); /* ToDo allow other geometries */
/* compute time delay for plane-wave responses */
- twplane = (float *) calloc(nxs,sizeof(float)); /* initialize with zeros */
+ twplane = (float *) calloc(nacq,sizeof(float)); /* initialize with zeros */
if (plane_wave==1) {
grad2rad = 17.453292e-3;
p = sin(src_angle*grad2rad)/src_velo;
if (p < 0.0) {
- for (i=0; i<nxs; i++) {
- twplane[i] = fabsf((nxs-i-1)*dxs*p)+tt0;
- if (verbose >=3) vmess("plane-wave delay-time i=%d x=%f t=%f", i, dxs*(nxs-i-1), twplane[i]);
+ for (i=0; i<nacq; i++) {
+ twplane[i] = fabsf((nacq-i-1)*dxs*p)+tt0;
+ if (verbose >=3) vmess("plane-wave delay-time i=%d x=%f t=%f", i, dxs*(nacq-i-1), twplane[i]);
}
}
else {
- for (i=0; i<nxs; i++) {
+ for (i=0; i<nacq; i++) {
twplane[i] = (i)*dxs*p+tt0;
if (verbose >=3) vmess("plane-wave delay-time i=%d x=%f t=%f", i, dxs*(i), twplane[i]);
}
@@ -260,19 +264,19 @@ int main (int argc, char **argv)
/*================ 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));
- RMi = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ Fop = (complex *)calloc(nacq*nw*Nfoc,sizeof(complex));
+ RMi = (float *)calloc(Nfoc*nacq*ntfft,sizeof(float));
+ DD = (float *)calloc(Nfoc*nacq*ntfft,sizeof(float));
Mi = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
- M0 = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
- DD = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ M0 = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ u1min = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ v1plus = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
SRC = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
RR = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
Mup = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
Msp = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
trace = (float *)malloc(ntfft*sizeof(float));
- ixpos = (int *)malloc(nxs*sizeof(int));
+ ixpos = (int *)malloc(nacq*sizeof(int));
energyM0= (double *)malloc(Nfoc*sizeof(double));
xrcvsyn = (float *)calloc(Nfoc*nxs,sizeof(float));
xsyn = (float *)malloc(Nfoc*sizeof(float));
@@ -363,8 +367,8 @@ int main (int argc, char **argv)
/*================ 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);
+ readShotData(file_shot, xrcv, xsrc, zsrc, xnx, Refl, nw, nw_low, nshots, 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));
@@ -402,19 +406,30 @@ int main (int argc, char **argv)
dxsrc = dx;
}
+/*================ set focusing postions in R ================*/
+
+ synthesisPositions(nx, nt, nacq, 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);
+ }
+
/*================ Defining focusing operator(s) from R ================*/
/* M0 = -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 M0 from Refl of %s", file_shot);
- nts = ntfft;
+ nts = ntfft;
+ //nxs = xnx[ishot];
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++) {
+ ixrcv = NINT((xrcv[ishot*nx+i]-fxsb)/dxs);
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;;
@@ -422,7 +437,7 @@ int main (int argc, char **argv)
/* 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];
+ DD[0*nxs*nts+ixrcv*nts+j] = -1.0*scl*rtrace[j];
}
/* remove taper at edges for selected shot record */
//if (tap == 2 || tap == 3) scltap= scl/tapersh[i];
@@ -439,15 +454,16 @@ int main (int argc, char **argv)
//for (i=0; i<nxs; i++) twplane[i] = dt*(maxval[i]-maxval[ishot]);
//for (i = 0; i < xnx[ishot]; i++) fprintf(stderr,"maxval[%d] = %f tw=%f\n", i, dt*maxval[i], twplane[i]);
- /* construct plane wave (time reversed and multiplid with -1) from all shot records */
+ /* construct plane wave (time reversed and multiplied with -1) from all shot records */
if (plane_wave==1) {
- for (l=0; l<nshots; l++) {
+ nxs = nshots;
+ for (l=0; l<nxs; l++) {
+ ix = ixpos[l];
memset(ctrace, 0, sizeof(complex)*(nfreq+1));
for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
- tom = j*deltom*twplane[l];
+ tom = j*deltom*twplane[ix];
csum.r=0.0; csum.i=0.0;
for (i = 0; i < xnx[l]; i++) {
- tom = j*deltom*twplane[i];
csum.r += Refl[l*nw*nx+m*nx+i].r*cos(-tom) - Refl[l*nw*nx+m*nx+i].i*sin(-tom);
csum.i += Refl[l*nw*nx+m*nx+i].i*cos(-tom) + Refl[l*nw*nx+m*nx+i].r*sin(-tom);
}
@@ -458,11 +474,11 @@ int main (int argc, char **argv)
/* transfrom result back to time domain */
cr1fft(ctrace, rtrace, ntfft, 1);
for (j = 0; j < nts; j++) {
- DD[0*nxs*nts+l*nts+j] = -1.0*scl*rtrace[j];
+ DD[0*nxs*nts+ix*nts+j] = -1.0*scl*rtrace[j];
}
/* compute Source wavelet for plane-wave imaging */
for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
- tom = j*deltom*twplane[l];
+ tom = j*deltom*twplane[ix];
csum.r = cos(-tom);
csum.i = sin(-tom);
/* Optional add wavelet to SRC-field */
@@ -483,14 +499,17 @@ int main (int argc, char **argv)
xsyn[0] = xsrc[ishot];
zsyn[0] = zsrc[ishot];
xnxsyn[0] = xnx[ishot];
+ /* find minimum and maximum postions in header values */
+/*
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);
}
}
+*/
+/* not needed xmin and xmax are the same */
dxf = dx;
dxs = dx;
}
@@ -507,13 +526,16 @@ int main (int argc, char **argv)
if (verbose) vmess("Taper for operator applied ntap=%d", ntap);
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < nxs; i++) {
+ ix = ixpos[i];
for (j = 0; j < nts; j++) {
- DD[l*nxs*nts+i*nts+j] *= tapersy[i];
+ DD[l*nxs*nts+ix*nts+j] *= tapersy[i];
}
}
}
free(tapersy);
}
+ ixp = (int *)malloc(nxs*sizeof(int));
+ for (i=0; i<nxs; i++) ixp[i] = i;
/*================ Check the size of the files ================*/
@@ -525,9 +547,10 @@ int main (int argc, char **argv)
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 receivers in tinv = %d", nxs);
vmess("number of shots = %d", nshots);
vmess("number of receiver/shot = %d", nx);
+ vmess("number of receiver acquisition = %d", nacq);
vmess("first model position = %.2f", fxsb);
vmess("last model position = %.2f", fxse);
vmess("first source position fxf = %.2f", fxf);
@@ -550,14 +573,14 @@ int main (int argc, char **argv)
if (file_iter != NULL) {
vmess("Iterations output file = %s ", file_iter);
vmess("Initialisation input data = M0_%06d.su ", 1000*istart);
- vmess("f1min intermediate array = f1min_%03d'iter'.su ", istart);
+ vmess("u1min intermediate array = u1min_%03d'iter'.su ", istart);
vmess("Mi intermediate array = Mi_%03d'iter'.su ", istart);
}
}
/*================ initializations ================*/
- if (reci) n2out = nxs;
+ if (reci) n2out = nacq;
else n2out = nshots;
mem = Nfoc*n2out*ntfft*sizeof(float)/1048576.0;
if (verbose) {
@@ -567,15 +590,6 @@ int main (int argc, char **argv)
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;
@@ -616,16 +630,22 @@ int main (int argc, char **argv)
perc=(iend-istart)/100;if(!perc)perc=1;
if (plane_wave) {
- writeDataIter("SRCplane.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, NINT(src_angle));
- /* make DD causal again and undo the -1 multplication */
- for (l=0; l<nshots; l++) {
- j=0;
- SRC[l*nts+j] = -1.0*DD[l*nts+j];
- for (j = 1; j < nts; j++) {
- SRC[l*nts+j] = -1.0*DD[l*nts+nts-j];
- }
+ writeDataIter("SRCplane.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, nxs, 0, NINT(src_angle));
+ }
+ /* make DD causal again and undo the -1 multiplication */
+ for (l=0; l<npos; l++) {
+ ix = ixpos[l];
+ j=0;
+ SRC[l*nts+j] = -1.0*DD[ix*nts+j];
+ for (j = 1; j < nts; j++) {
+ SRC[l*nts+j] = -1.0*DD[ix*nts+nts-j];
}
- writeDataIter("DDplane.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, NINT(src_angle));
+ }
+ if (plane_wave) {
+ writeDataIter("DDplane.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, nxs, 0, NINT(src_angle));
+ }
+ else {
+ writeDataIter("DDshot.su", SRC, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, nxs, 0, ishot);
}
free(SRC);
@@ -635,40 +655,43 @@ int main (int argc, char **argv)
/*================ initialization ================*/
+ memset(M0, 0, Nfoc*nxs*ntfft*sizeof(float));
+ memset(v1plus, 0, Nfoc*nxs*ntfft*sizeof(float));
/* once every 'niterskip' time-steps start from fresh M0 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);
+ if (verbose>2) 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++) {
- iw = NINT((ii*dt+twplane[i])/dt);
- //iw = ii;
+ for (i = 0; i < npos; i++) {
+ ix = ixpos[i];
+ iw = NINT((ii*dt+twplane[ix])/dt);
for (j = 0; j < nts; j++) {
- M0[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+i*nts+j];
+ M0[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+ix*nts+j];
}
/* apply mute window for samples above nts-ii */
for (j = 0; j < MIN(nts, nts-iw+isms); j++) {
M0[l*nxs*nts+i*nts+j] = 0.0;
}
for (j = nts-iw+isms, k=1; j < MIN(nts, nts-iw+isme); j++, k++) {
+ //fprintf(stderr,"j=%d smooth-k=%d taper=%f\n",j,smooth-k,costaper[smooth-k]);
M0[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];
+ u1min[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];
+ u1min[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j];
}
}
}
if (file_iter != NULL) {
- writeDataIter("M0.su", M0, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii);
+ writeDataIter("M0.su", M0, hdrs_out, ntfft, nacq, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii);
}
}
- else { /* use f1min from previous iteration as starting point and do niterec iterations */
+ else { /* use u1min 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);
@@ -678,9 +701,9 @@ int main (int argc, char **argv)
ix = ixpos[i];
iw = NINT((ii*dt+twplane[ix])/dt);
//iw = ii;
- M0[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts] - f1min[l*nxs*nts+i*nts+j];
+ M0[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts] - u1min[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- M0[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j] - f1min[l*nxs*nts+i*nts+nts-j];
+ M0[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j] - u1min[l*nxs*nts+i*nts+nts-j];
}
/* apply mute window for samples above nts-ii */
for (j = 0; j < MIN(nts,nts-iw+isms); j++) {
@@ -705,25 +728,26 @@ int main (int argc, char **argv)
/*================ construction of Mi(-t) = - \int R(x,t) Mi(t) ================*/
- synthesis(Refl, Fop, Mi, RMi, nx, nt, nxs, nts, dt, xsyn, Nfoc,
+ synthesisp(Refl, Fop, Mi, RMi, nx, nt, nacq, 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, RMi, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter+1);
+ writeDataIter(file_iter, RMi, hdrs_out, ntfft, nacq, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter+1);
}
if (verbose >=2) {
for (l = 0; l < Nfoc; l++) {
energyMi = 0.0;
for (i = 0; i < npos; i++) {
+ ix = ixpos[i];
for (j = 0; j < nts; j++) {
- energyMi += RMi[l*nxs*nts+ix*nts+j]*RMi[l*nxs*nts+ix*nts+j];
+ energyMi += RMi[l*nacq*nts+ix*nts+j]*RMi[l*nacq*nts+ix*nts+j];
}
}
if (iter % 2 == 0) {
- if (iter==0) energyM0[Nfoc] = energyMi;
- vmess(" - iSyn %d: Mi at iteration %d has energy %e; relative to M0 %e", l, iter, sqrt(energyMi), sqrt(energyMi/energyM0[Nfoc]));
+ if (iter==0) energyM0[l] = energyMi;
+ vmess(" - ii %d: Mi at iteration %d has energy %e; relative to M0 %e", ii, iter, sqrt(energyMi), sqrt(energyMi/energyM0[l]));
}
}
}
@@ -736,9 +760,9 @@ int main (int argc, char **argv)
for (i = 0; i < npos; i++) {
j = 0;
ix = ixpos[i];
- Mi[l*nxs*nts+i*nts+j] = -RMi[l*nxs*nts+ix*nts+j];
+ Mi[l*nxs*nts+i*nts+j] = -RMi[l*nacq*nts+ix*nts+j];
for (j = 1; j < nts; j++) {
- Mi[l*nxs*nts+i*nts+j] = -RMi[l*nxs*nts+ix*nts+nts-j];
+ Mi[l*nxs*nts+i*nts+j] = -RMi[l*nacq*nts+ix*nts+nts-j];
}
}
}
@@ -765,44 +789,41 @@ int main (int argc, char **argv)
for (j = MAX(0,iw+shift-smooth), k=1; j < MAX(0,iw+shift); j++, k++) {
Mi[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
}
- f1plus[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- f1plus[l*nxs*nts+i*nts+j] += Mi[l*nxs*nts+i*nts+j];
+ for (j = 0; j < nts; j++) {
+ v1plus[l*nxs*nts+i*nts+j] += Mi[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);
+ writeDataIter("v1plus.su", v1plus, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii+iter+1);
}
-*/
}
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 */
+ if (recur==1) { /* use u1min from previous iteration */
for (j = 0; j < nts; j++) {
Mi[l*nxs*nts+i*nts+j] += DD[l*nxs*nts+ix*nts+j];
}
j = 0;
- f1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+j];
+ u1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- f1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+nts-j];
+ u1min[l*nxs*nts+i*nts+j] = -Mi[l*nxs*nts+i*nts+nts-j];
}
if (file_update != NULL) {
j=0;
- Mup[l*nxs*nts+i*nts+j] += f1min[l*nxs*nts+i*nts+j]+DD[l*nxs*nts+i*nts+j];
+ Mup[l*nxs*nts+i*nts+j] += u1min[l*nxs*nts+i*nts+j]+DD[l*nxs*nts+ix*nts+j];
for (j = 1; j < nts; j++) {
- Mup[l*nxs*nts+i*nts+j] += f1min[l*nxs*nts+i*nts+j]+DD[l*nxs*nts+i*nts+nts-j];
+ Mup[l*nxs*nts+i*nts+j] += u1min[l*nxs*nts+i*nts+j]+DD[l*nxs*nts+ix*nts+nts-j];
}
}
}
else {
j = 0;
- f1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+j];
+ u1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
- f1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+nts-j];
+ u1min[l*nxs*nts+i*nts+j] -= Mi[l*nxs*nts+i*nts+nts-j];
}
if (file_update != NULL) {
j=0;
@@ -832,11 +853,11 @@ int main (int argc, char **argv)
}
}
if (file_iter != NULL) {
- writeDataIter("f1min.su", f1min, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter+1);
+ writeDataIter("u1min.su", u1min, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii+iter+1);
}
} /* end else (iter) branch */
if (file_iter != NULL) {
- writeDataIter("Mi.su", Mi, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter+1);
+ writeDataIter("Mi.su", Mi, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixp, npos, 0, 1000*ii+iter+1);
}
t2 = wallclock_time();
@@ -850,7 +871,7 @@ int main (int argc, char **argv)
ix = ixpos[i];
iw = NINT((ii*dt+twplane[ix])/dt);
if ( iw<nts && iw>=0 ) {
- RR[l*nxs*nts+i*nts+iw] = f1min[l*nxs*nts+i*nts+iw];
+ RR[l*nxs*nts+i*nts+iw] = u1min[l*nxs*nts+i*nts+iw];
if (file_update != NULL) Msp[l*nxs*nts+i*nts+iw] = Mup[l*nxs*nts+i*nts+iw];
}
}
@@ -875,8 +896,8 @@ int main (int argc, char **argv)
free(Mi);
free(energyM0);
free(M0);
- free(f1min);
- free(f1plus);
+ free(u1min);
+ free(v1plus);
t2 = wallclock_time();
if (verbose) {
diff --git a/marchenko/synthesis.c b/marchenko/synthesis.c
index e975d8746ffa1b44fcf34a15f9b3fa0db03772fb..237a0164aba1fcc1db9606275068a0dcbc75d983 100644
--- a/marchenko/synthesis.c
+++ b/marchenko/synthesis.c
@@ -89,16 +89,16 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
ctrace = (complex *)calloc(ntfft,sizeof(complex));
/* this first check is done to support an acquisition geometry that has more receiver than source
- * postions. In the first iteration the int R(x_r,x_s) Fop(x_r) d x_r results in a grid on x_s.
- * so for the next interations onlt x_s traces have to be computed on Fop */
+ * positions. In the first iteration the int R(x_r,x_s) Fop(x_r) d x_r results in a grid on x_s.
+ * so for the next interations only x_s traces have to be computed on Fop */
if (!first) {
- /* transform muted Ni (Top) to frequency domain, input for next iteration */
+ /* 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], 0, nxs*nw*sizeof(complex));
for (i = 0; i < npos; i++) {
- rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
ix = ixpos[i];
+ rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
for (iw=0; iw<nw; iw++) {
Fop[l*nxs*nw+iw*nxs+ix].r = ctrace[nw_low+iw].r;
Fop[l*nxs*nw+iw*nxs+ix].i = mode*ctrace[nw_low+iw].i;
@@ -107,7 +107,7 @@ nw, int nw_low, int nw_high, int mode, int reci, int nshots, int *ixpos, int np
}
}
else { /* only for first call to synthesis using all nxs traces in G_d */
- /* transform G_d to frequency domain, over all nxs traces */
+ /* transform G_d to frequency domain, over all nxs traces */
first=0;
for (l = 0; l < Nfoc; l++) {
/* set Fop to zero, so new operator can be defined within all ix points */
@@ -271,6 +271,7 @@ float fxse, float fxsb, float dxs, float dxsrc, float dx, int nshots, int *ixpos
for (k=0; k<nshots; k++) {
ixsrc = NINT((xsrc[k] - fxsb)/dxs);
+ vmess("source position: %.2f in operator %d", xsrc[k], ixsrc);
if (verbose>=3) {
vmess("source position: %.2f in operator %d", xsrc[k], ixsrc);
vmess("receiver positions: %.2f <--> %.2f", xrcv[k*nx+0], xrcv[k*nx+nx-1]);
diff --git a/utils/getFileInfo.c b/utils/getFileInfo.c
index 5aec1a4def8a94cfdfe14f1d60a03ea434e19cdd..d93010a4f54f1023a0aafc64a57ea7c8d9731860 100644
--- a/utils/getFileInfo.c
+++ b/utils/getFileInfo.c
@@ -43,6 +43,8 @@ int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *
ret = fseeko( fp, 0, SEEK_END );
if (ret<0) perror("fseeko");
bytes = ftello( fp );
+ *xmax = hdr.gx;
+ *xmin = hdr.gx;
*n1 = hdr.ns;
if ( (hdr.trid == 1) && (hdr.dt != 0) ) {
@@ -112,9 +114,7 @@ int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *
one_shot = 1;
igath = 0;
fseeko( fp, 0, SEEK_SET );
- offset = (NINT((hdr.gx-hdr.sx)*scl*100)/100.0);
- *xmax = offset;
- *xmin = offset;
+ //offset = (NINT((hdr.gx-hdr.sx)*scl*100)/100.0);
dxrcv = *d2;
while (!end_of_file) {
@@ -125,6 +125,8 @@ int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *
sx_shot = hdr.sx;
gx_start = hdr.gx;
gx_end = hdr.gx;
+ *xmax = MAX(MAX(hdr.gx,*xmax),hdr.sx);
+ *xmin = MIN(MIN(hdr.gx,*xmin),hdr.sx);
itrace = 0;
while (one_shot) {
@@ -132,9 +134,11 @@ int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *
itrace++;
if (hdr.gx != gx_end) dxrcv = MIN(dxrcv,abs(hdr.gx-gx_end));
gx_end = hdr.gx;
- offset = (NINT((hdr.gx-hdr.sx)*scl*100)/100.0);
- *xmax = MAX(*xmax,offset);
- *xmin = MIN(*xmin,offset);
+ //offset = (NINT((hdr.gx-hdr.sx)*scl*100)/100.0);
+ //*xmax = MAX(*xmax,offset);
+ //*xmin = MIN(*xmin,offset);
+ *xmax = MAX(MAX(hdr.gx,*xmax),hdr.sx);
+ *xmin = MIN(MIN(hdr.gx,*xmin),hdr.sx);
nread = fread( &hdr, 1, TRCBYTES, fp );
if (nread != TRCBYTES) {
one_shot = 0;
@@ -167,14 +171,18 @@ int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *
fseeko( fp, -trace_sz, SEEK_END );
nread = fread( &hdr, 1, TRCBYTES, fp );
- offset = (NINT((hdr.gx-hdr.sx)*scl*100)/100.0);
- *xmax = MAX(*xmax,offset);
- *xmin = MIN(*xmin,offset);
- *xmin = MIN(sx_shot,hdr.sx*scl);
+ //offset = (NINT((hdr.gx-hdr.sx)*scl*100)/100.0);
+ //*xmax = MAX(*xmax,offset);
+ //*xmin = MIN(*xmin,offset);
+ //*xmin = MIN(sx_shot,hdr.sx*scl);
+ *xmax = MAX(MAX(hdr.gx,*xmax),hdr.sx);
+ *xmin = MIN(MIN(hdr.gx,*xmin),hdr.sx);
*ngath = ntraces/(*n2);
}
// *nxm = NINT((*xmax-*xmin)/dxrcv)+1;
*nxm = (int)ntraces;
+ *xmax *= scl;
+ *xmin *= scl;
fclose( fp );
free(trace);