diff --git a/README b/README
index 54dd9c38ccb2d75202b54329714629f8d67e4640..8dd89ff767db5b607815a3ab35059e82eb80ebef 100644
--- a/README
+++ b/README
@@ -13,6 +13,9 @@ For more information, go to http://software.seg.org/2017/00XX .
You must read and accept usage terms at:
http://software.seg.org/disclaimer.txt before use.
+DOI reference of release
+1.4.0 https://zenodo.org/badge/latestdoi/23060862
+
REFERENCES
---------
-1- If the Finite Difference code has helped you in your research please refer to our paper in your publications:
diff --git a/marchenko/demo/oneD/iterations.scr b/marchenko/demo/oneD/iterations.scr
new file mode 100755
index 0000000000000000000000000000000000000000..07839ed9456ba7c0b8d3284183b59f4a4194fff2
--- /dev/null
+++ b/marchenko/demo/oneD/iterations.scr
@@ -0,0 +1,24 @@
+#!/bin/bash
+
+#second reflector at time:
+# 800/1800+600/2300
+# .70531400966183574878 => sample 176
+#thrids reflector at time:
+# 800/1800+600/2300+800/2000
+# 1.10531400966183574878 sample 276
+
+select=451
+
+#for istart in 176 200 276 400
+for istart in 176
+do
+(( iend = istart + 1 ))
+
+../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
+ nshots=601 verbose=1 istart=$istart iend=$iend fmax=90 \
+ pad=1024 niter=45 smooth=10 niterskip=600 shift=20 file_rr=pred_rr_${istart}.su T=0 file_iter=iter.su
+
+suximage < f1min_${istart}043.su clip=1 title="${istart}" &
+
+done
+
diff --git a/marchenko/demo/oneD/primaries.scr b/marchenko/demo/oneD/primaries.scr
index 69cbf385a49b1ce148a7125ec07ef725c69a83b8..f5a82f2ebb581b3ad54f6bea1ec0e3e7f0550389 100755
--- a/marchenko/demo/oneD/primaries.scr
+++ b/marchenko/demo/oneD/primaries.scr
@@ -13,14 +13,14 @@ export OMP_NUM_THREADS=20
#shot record to remove internal multiples
select=451
-makewave fp=20 dt=0.004 file_out=wave.su nt=1024 t0=0.0 scale=0 scfft=1
+makewave fp=20 dt=0.004 file_out=wave.su nt=2048 t0=0.0 scale=0 scfft=1
../../marchenko_primaries file_shot=shotsdx5_rp.su ishot=$select file_src=wave.su \
- nshots=601 verbose=1 istart=40 fmax=90 \
- niter=15 niterskip=600 shift=20 file_rr=pred_rr.su T=0
+ nshots=601 verbose=2 istart=40 iend=500 fmax=90 pad=1024 \
+ niter=25 smooth=10 niterskip=600 shift=20 file_rr=pred_rr.su T=0
#for reference original shot record from Reflection matrix
-suwind key=fldr min=$select max=$select < shotsdx5_rp.su > shotsx0.su
+suwind key=fldr min=$select max=$select < shotsdx5_rp.su itmax=2047 > shotsx0.su
fconv file_in1=shotsx0.su file_in2=wave.su file_out=shotw.su
# for displaying results
diff --git a/marchenko/demo/oneD/primariesFrame.scr b/marchenko/demo/oneD/primariesFrame.scr
index c1a7bf4ee63cc85e3ecabad2eed17740a8eb664b..5017c4b908a441e7cec67ef68cdb183adeb513a5 100755
--- a/marchenko/demo/oneD/primariesFrame.scr
+++ b/marchenko/demo/oneD/primariesFrame.scr
@@ -23,13 +23,13 @@ fconv file_in1=shotsx0.su file_in2=wave.su file_out=shotw.su
basop file_in=shotw.su choice=5 | sugain scale=-1 > DD.su
nts=1024
-itime=300 # select time sample to work on
+itime=300 # select time sample to work on 1.2 seconds
shift=20
#mute time nts-ii+shift to compute G_d
(( itmax = nts-itime+shift ))
suzero itmax=$itmax < DD.su > G_d.su
#f1min = -DD(-t) = shotw.su
-cp shotw.su f1min.su
+cp shotw.su f1min0.su
#first iteration
cp G_d.su N0.su
@@ -39,6 +39,26 @@ fconv file_in1=shotsdx5_rp.su file_in2=N0.su file_out=fconvN1.su verbose=1 fmax=
suwind key=fldr min=451 max=451 < fconvN1.su | suximage
+# Ni(t) = -1*RN0(-t)
+sustack < fconvN1.su key=fldr > RN0.su
+basop file_in=RN0.su choice=5 | sugain scale=-1 > N1.su
+
+#apply mute window for samples above nts-ii
+(( itmax = nts-itime+shift ))
+suzero itmax=$itmax < N1.su > N1m.su
+
+
+exit
+#display
+file=fconvN1.su
+file_base=${file%.su}
+sumax < $file mode=abs outpar=nep
+clipref=`cat nep | awk '{print $1/2}'`
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < $file \
+ label1="time (s)" label2="distance (m)" \
+ n1tic=2 d2=5 x1beg=0 x1end=1.5 d1num=0.4 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_flip.eps
+
exit;
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index 3182ff4e7232ca2abdece4a31b87fd212e262f71..69a8977c90229b0fe5641c359bf7d8914a15c5ed 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -35,7 +35,7 @@ typedef struct _complexStruct { /* complex number */
int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int nshots,
int nx, int nxs, float fxsb, float dxs, int ntfft, int mode, float scale, float tsq, float Q, float f0, int reci, int *nshots_r, int *isxcount, int *reci_xsrc, int *reci_xrcv, float *ixmask, int verbose);
int readTinvData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nfoc, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
-int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int iter);
+int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int t0shift, int iter);
void name_ext(char *filename, char *extension);
@@ -464,7 +464,7 @@ int main (int argc, char **argv)
tsyn += t3 - t2;
if (file_iter != NULL) {
- writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, iter);
+ writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 1, iter);
}
/* N_k(x,t) = -N_(k-1)(x,-t) */
/* p0^-(x,t) += iRN = (R * T_d^inv)(t) */
diff --git a/marchenko/marchenko_primaries.c b/marchenko/marchenko_primaries.c
index 7942e75774d94f5eccf159397d45ca2dc093dbf9..80c23983eca7e283ff611ec2e8fa2e149f2f72dc 100644
--- a/marchenko/marchenko_primaries.c
+++ b/marchenko/marchenko_primaries.c
@@ -29,7 +29,7 @@ int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx
int readTinvData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nfoc, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose);
-int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int iter);
+int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int t0shift, int iter);
int getFileInfo(char *filename, int *n1, int *n2, int *ngath, float *d1, float *d2, float *f1, float *f2, float *xmin, float *xmax, float *sclsxgx, int *nxm);
int readData(FILE *fp, float *data, segy *hdrs, int n1);
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
@@ -102,7 +102,7 @@ int main (int argc, char **argv)
{
FILE *fp_out, *fp_rr, *fp_w;
size_t nread;
- int i, j, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath;
+ int i, j, k, l, ret, nshots, Nfoc, nt, nx, nts, nxs, ngath;
int size, n1, n2, ntap, tap, di, ntraces;
int nw, nw_low, nw_high, nfreq, *xnx, *xnxsyn;
int reci, countmin, mode, ixa, ixb, n2out, verbose, ntfft;
@@ -114,9 +114,9 @@ int main (int argc, char **argv)
double t0, t1, t2, t3, t4, tsyn, tread, tfft, tcopy, tii;
float d1, d2, f1, f2, fxsb, fxse, ft, fx, *xsyn, dxsrc;
float *G_d, *DD, *RR, dt, dx, dxs, scl, mem;
- float *rtrace, *tmpdata, *f1min, *iRN, *Ni, *trace;
+ float *rtrace, *tmpdata, *f1min, *f1plus, *iRN, *Ni, *trace;
float xmin, xmax, scale, tsq;
- float Q, f0, *ixmask;
+ float Q, f0, *ixmask, *costaper;
complex *Refl, *Fop, *ctrace, *cwave;
char *file_tinv, *file_shot, *file_rr, *file_src, *file_iter;
segy *hdrs_out, hdr;
@@ -157,9 +157,17 @@ int main (int argc, char **argv)
if(!getparint("shift", &shift)) shift=20;
if(!getparint("smooth", &smooth)) smooth = 5;
if(!getparint("ishot", &ishot)) ishot=300;
-
if (reci && ntap) vwarn("tapering influences the reciprocal result");
+ smooth = MIN(smooth, shift);
+ if (smooth) {
+ costaper = (float *)malloc(smooth*sizeof(float));
+ scl = M_PI/((float)smooth);
+ for (i=0; i<smooth; i++) {
+ costaper[i] = 0.5*(1.0+cos((i+1)*scl));
+ }
+ }
+
/*================ Reading info about shot and initial operator sizes ================*/
if (file_tinv != NULL) { /* G_d is read from file_tinv */
@@ -204,6 +212,7 @@ int main (int argc, char **argv)
Fop = (complex *)calloc(nxs*nw*Nfoc,sizeof(complex));
f1min = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
+ f1plus = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
iRN = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
Ni = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
G_d = (float *)calloc(Nfoc*nxs*ntfft,sizeof(float));
@@ -411,20 +420,21 @@ 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 shots = %d", nshots);
- vmess("number of receiver/shot = %d", nx);
- vmess("first model position = %.2f", fxsb);
- vmess("last model position = %.2f", fxse);
- vmess("first source position fxf = %.2f", fxf);
- vmess("source distance dxsrc = %.2f", dxsrc);
- vmess("last source position = %.2f", fxf+(nshots-1)*dxsrc);
- vmess("receiver distance dxf = %.2f", dxf);
+ vmess("Number of receivers in focusop = %d", nxs);
+ vmess("number of shots = %d", nshots);
+ vmess("number of receiver/shot = %d", nx);
+ vmess("first model position = %.2f", fxsb);
+ vmess("last model position = %.2f", fxse);
+ vmess("first source position fxf = %.2f", fxf);
+ vmess("source distance dxsrc = %.2f", dxsrc);
+ vmess("last source position = %.2f", fxf+(nshots-1)*dxsrc);
+ vmess("receiver distance dxf = %.2f", dxf);
vmess("direction of increasing traces = %d", di);
- vmess("number of time samples (nt,nts) = %d (%d,%d)", ntfft, nt, nts);
- vmess("time sampling = %e ", dt);
- if (file_rr != NULL) vmess("RR output file = %s ", file_rr);
- if (file_iter != NULL) vmess("Iterations output file = %s ", file_iter);
+ vmess("number of time samples fft nt nts = %d %d %d", ntfft, nt, nts);
+ vmess("time sampling = %e ", dt);
+ vmess("smoothing taper for time-window = %d ", smooth);
+ if (file_rr != NULL) vmess("RR output file = %s ", file_rr);
+ if (file_iter != NULL) vmess("Iterations output file = %s ", file_iter);
}
/*================ initializations ================*/
@@ -433,9 +443,10 @@ int main (int argc, char **argv)
else n2out = nshots;
mem = Nfoc*n2out*ntfft*sizeof(float)/1048576.0;
if (verbose) {
- vmess("number of output traces = %d", n2out);
- vmess("number of output samples = %d", ntfft);
- vmess("Size of output data/file = %.1f MB", mem);
+ vmess("Time-sample range processed = %d:%d", istart, iend);
+ vmess("number of output traces = %d", n2out);
+ vmess("number of output samples = %d", ntfft);
+ vmess("Size of output data/file = %.1f MB", mem);
}
ixa=0; ixb=0;
@@ -484,7 +495,7 @@ int main (int argc, char **argv)
vmess("*******************************************");
fprintf(stderr," %s: Progress: %3d%%",xargv[0],0);
}
- perc=(iend-istart)/10;if(!perc)perc=1;
+ perc=(iend-istart)/100;if(!perc)perc=1;
/*================ start loop over number of time-samples ================*/
@@ -503,9 +514,15 @@ int main (int argc, char **argv)
G_d[l*nxs*nts+i*nts+j] = DD[l*nxs*nts+i*nts+j];
}
/* apply mute window for samples above nts-ii */
- for (j = 0; j < nts-ii+T*shift; j++) {
+ for (j = 0; j < nts-ii+T*shift-smooth; j++) {
G_d[l*nxs*nts+i*nts+j] = 0.0;
}
+ for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) {
+ G_d[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
+ }
+ //for (j = 0; j < nts-ii+T*shift; j++) {
+ // G_d[l*nxs*nts+i*nts+j] = 0.0;
+ //}
}
for (i = 0; i < npos; i++) {
ix = ixpos[i];
@@ -514,6 +531,10 @@ int main (int argc, char **argv)
for (j = 1; j < nts; j++) {
f1min[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j];
}
+ /* apply mute window for samples above nts-ii */
+ //for (j = 0; j < nts-ii+T*shift; j++) {
+ // f1min[l*nxs*nts+i*nts+j] = 0.0;
+ //}
}
}
}
@@ -529,12 +550,21 @@ int main (int argc, char **argv)
for (j = 1; j < nts; j++) {
G_d[l*nxs*nts+i*nts+j] = -DD[l*nxs*nts+ix*nts+nts-j] - f1min[l*nxs*nts+i*nts+nts-j];
}
- for (j = 0; j < nts-ii+T*shift; j++) {
+ /* apply mute window for samples above nts-ii */
+ for (j = 0; j < nts-ii+T*shift-smooth; j++) {
G_d[l*nxs*nts+i*nts+j] = 0.0;
}
+ for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) {
+ G_d[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
+ }
+ //for (j = 0; j < nts-ii+T*shift; j++) {
+ // G_d[l*nxs*nts+i*nts+j] = 0.0;
+ //}
}
}
}
+/*================ initialization ================*/
+
memcpy(Ni, G_d, Nfoc*nxs*ntfft*sizeof(float));
/*================ number of Marchenko iterations ================*/
@@ -550,7 +580,7 @@ int main (int argc, char **argv)
reci, nshots, ixpos, npos, &tfft, isxcount, reci_xsrc, reci_xrcv, ixmask, verbose);
if (file_iter != NULL) {
- writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 1000*ii+iter);
+ writeDataIter(file_iter, iRN, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 1, 1000*ii+iter);
}
t3 = wallclock_time();
@@ -568,20 +598,35 @@ int main (int argc, char **argv)
}
}
- if (iter % 2 == 0) { /* even iterations: => f_1^-(t) */
+ if (iter % 2 == 0) { /* even iterations: => f_1^+(t) */
/* apply muting for the acausal part */
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
+ /* apply mute window for samples after ii */
for (j = ii-T*shift; j < nts; j++) {
Ni[l*nxs*nts+i*nts+j] = 0.0;
}
- for (j = 0; j < shift; j++) {
+ for (j = ii-T*shift+smooth, k=0; j < ii-T*shift; j++, k++) {
+ Ni[l*nxs*nts+i*nts+j] *= costaper[k];
+ }
+ /* apply mute window for delta function at t=0*/
+ for (j = 0; j < shift-smooth; j++) {
Ni[l*nxs*nts+i*nts+j] = 0.0;
}
+ for (j = shift-smooth, k=1; j < shift; j++, k++) {
+ Ni[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
+ }
+ f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ f1plus[l*nxs*nts+i*nts+j] += Ni[l*nxs*nts+i*nts+j];
+ }
}
}
+ if (file_iter != NULL) {
+ writeDataIter("f1plus.su", f1plus, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter);
+ }
}
- else {/* odd iterations: => f_1^+(t) */
+ else {/* odd iterations: => f_1^-(t) */
for (l = 0; l < Nfoc; l++) {
for (i = 0; i < npos; i++) {
ix = ixpos[i];
@@ -602,17 +647,30 @@ int main (int argc, char **argv)
f1min[l*nxs*nts+i*nts+j] -= Ni[l*nxs*nts+i*nts+nts-j];
}
}
- /* apply mute window */
- for (j = nts-shift; j < nts; j++) {
- Ni[l*nxs*nts+i*nts+j] = 0.0;
- }
+ /* apply mute window for delta function at t=0*/
+ for (j = nts-shift+smooth; j < nts; j++) {
+ Ni[l*nxs*nts+i*nts+j] = 0.0;
+ }
+ for (j = nts-shift, k=0; j < nts-shift+smooth; j++, k++) {
+ Ni[l*nxs*nts+i*nts+j] *= costaper[k];
+ }
/* apply mute window for samples above nts-ii */
- for (j = 0; j < nts-ii+T*shift; j++) {
- Ni[l*nxs*nts+i*nts+j] = 0.0;
- }
+ /* important to apply this mute after updating f1min */
+ //for (j = 0; j < nts-ii+T*shift; j++) {
+ // Ni[l*nxs*nts+i*nts+j] = 0.0;
+ //}
+ /* apply mute window for samples above nts-ii */
+ for (j = 0; j < nts-ii+T*shift-smooth; j++) {
+ Ni[l*nxs*nts+i*nts+j] = 0.0;
+ }
+ for (j = nts-ii+T*shift-smooth, k=1; j < nts-ii+T*shift; j++, k++) {
+ Ni[l*nxs*nts+i*nts+j] *= costaper[smooth-k];
+ }
}
}
- //fprintf(stderr,"f1min at %d = %f NI=%f\n", ii, f1min[(nxs/2)*nts+ii], Ni[(nxs/2)*nts+ii]);
+ if (file_iter != NULL) {
+ writeDataIter("f1min.su", f1min, hdrs_out, ntfft, nxs, d2, f2, n2out, Nfoc, xsyn, zsyn, ixpos, npos, 0, 1000*ii+iter);
+ }
} /* end else (iter) branch */
t2 = wallclock_time();
@@ -627,13 +685,13 @@ int main (int argc, char **argv)
}
}
- /* To Do optional write intermediate RR results to file */
+ /* To Do? optional write intermediate RR results to file */
if (verbose) {
if(!((iend-ii-istart)%perc)) fprintf(stderr,"\b\b\b\b%3d%%",(ii-istart)*100/(iend-istart));
if((ii-istart)==10)t4=wallclock_time();
- if((ii-istart)==50){
- t4=(wallclock_time()-t4)*((iend-istart)/40.0);
+ if((ii-istart)==20){
+ t4=(wallclock_time()-t4)*((iend-istart)/10.0);
fprintf(stderr,"\r %s: Estimated total compute time = %.2fs.\n %s: Progress: %3d%%",xargv[0],t4,xargv[0],(ii-istart)/((iend-istart)/100));
}
//t4=wallclock_time();
@@ -645,6 +703,8 @@ int main (int argc, char **argv)
free(Ni);
free(G_d);
+ free(f1min);
+ free(f1plus);
t2 = wallclock_time();
if (verbose) {
diff --git a/marchenko/writeDataIter.c b/marchenko/writeDataIter.c
index 5de0b5d49aba95acbd476a323c8f825c6fc1d740..b1a51f3d0aa776be14f3a765d1e23e103eb42eab 100644
--- a/marchenko/writeDataIter.c
+++ b/marchenko/writeDataIter.c
@@ -17,7 +17,7 @@
void name_ext(char *filename, char *extension);
int writeData(FILE *fp, float *data, segy *hdrs, int n1, int n2);
-int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int iter)
+int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, float d2, float f2, int n2out, int Nfoc, float *xsyn, float *zsyn, int *ixpos, int npos, int t0shift, int iter)
{
FILE *fp_iter;
size_t nwrite;
@@ -43,14 +43,19 @@ int writeDataIter(char *file_iter, float *data, segy *hdrs, int n1, int n2, floa
hdrs[i].selev = NINT(zsyn[l]*1000);
hdrs[i].sdepth = NINT(-zsyn[l]*1000);
/* rotate to get t=0 in the middle */
- hdrs[i].f1 = -n1*0.5*hdrs[i].d1;
- memcpy(&trace[0],&data[l*size+ix*n1],n1*sizeof(float));
- for (j = 0; j < n1/2; j++) {
- trace[n1/2+j] = data[l*size+ix*n1+j];
- }
- for (j = n1/2; j < n1; j++) {
- trace[j-n1/2] = data[l*size+ix*n1+j];
- }
+ if (t0shift) {
+ hdrs[i].f1 = -n1*0.5*hdrs[i].d1;
+ for (j = 0; j < n1/2; j++) {
+ trace[n1/2+j] = data[l*size+ix*n1+j];
+ }
+ for (j = n1/2; j < n1; j++) {
+ trace[j-n1/2] = data[l*size+ix*n1+j];
+ }
+ }
+ else {
+ hdrs[i].f1 = 0.0;
+ memcpy(&trace[0],&data[l*size+ix*n1],n1*sizeof(float));
+ }
nwrite = fwrite(&hdrs[i], 1, TRCBYTES, fp_iter);
assert(nwrite == TRCBYTES);
nwrite = fwrite(trace, sizeof(float), n1, fp_iter);