From cdd70d64153b8c4ff62eae0e5bb49f6da3b0a8ca Mon Sep 17 00:00:00 2001
From: JanThorbecke <janth@xs4all.nl>
Date: Mon, 16 Mar 2020 09:44:02 +0100
Subject: [PATCH] MME algorithm changes for paper
---
fdelmodc/fdelmodc.c | 4 +-
marchenko/Makefile | 2 +-
marchenko/demo/oneD/README_PRIMARIES | 7 +-
marchenko/demo/oneD/epsPrimaries.scr | 41 ++++-
marchenko/demo/oneD/iterations.scr | 23 ++-
marchenko/demo/oneD/model.scr | 24 +++
marchenko/demo/oneD/primaries.scr | 7 +-
marchenko/marchenko_primaries.c | 215 +++++++++++++++------------
marchenko/synthesis.c | 11 +-
utils/getFileInfo.c | 28 ++--
10 files changed, 236 insertions(+), 126 deletions(-)
diff --git a/fdelmodc/fdelmodc.c b/fdelmodc/fdelmodc.c
index 878f6b6..f780b78 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 cd1c541..704ba08 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 5eda719..7b97347 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 00949d8..2d05f4c 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 3cc1ec7..f6f774f 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 9165e87..1a5e0de 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 b5398ef..ca454a6 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 d893a33..b6775fd 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 e975d87..237a016 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 5aec1a4..d93010a 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);
--
GitLab