diff --git a/.gitignore b/.gitignore
index 11ccdb196da8efd589de46b58e2f84b4b677156b..34223bc79142f78e7b583d9d61d62916c0fabf66 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,4 @@
bin/*
*.su
*.bin
+.DS*
diff --git a/FFTlib/Makefile b/FFTlib/Makefile
index 7585a210fad86ae97fbb530dab5f04584b45c3af..a62691e8dbd8ada7b4b9386dbb04623db0cb7ef3 100644
--- a/FFTlib/Makefile
+++ b/FFTlib/Makefile
@@ -43,7 +43,7 @@ $(LIB) : $(ARCH)
mkdirs:
-mkdir -p $(ROOT)/lib
-mkdir -p $(ROOT)/include
- -ln -sf genfft.h ../include/genfft.h
+ -cp -rp genfft.h ../include/genfft.h
bins:
cd test ; $(MAKE)
diff --git a/fdelmodc/demo/benchmark.scr b/fdelmodc/demo/benchmark.scr
new file mode 100755
index 0000000000000000000000000000000000000000..e689f02b8421bfb6c90b7a90e46df2567eb7a35f
--- /dev/null
+++ b/fdelmodc/demo/benchmark.scr
@@ -0,0 +1,39 @@
+#!/bin/bash
+
+cp=2000
+rho=1000
+dx=2.5
+dt=0.0005
+
+makemod sizex=6000 sizez=4000 dx=$dx dz=$dx cp0=$cp cs0=$cs ro0=$rho \
+ orig=-3000,-1000 file_base=syncl.su \
+ intt=def x=-3000,0,3000 z=250,250,250 poly=0 cp=2300 ro=5000 \
+ intt=def x=-3000,-2000,-1000,-800,0,800,3000 z=650,650,700,750,900,750,600 poly=2 cp=2600 ro=1000 \
+ intt=def x=-3000,0,3000 z=1390,1390,1390 poly=0 cp=2000 ro=5000
+
+makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
+
+xsrc=-200
+
+time ../fdelmodc \
+ file_cp=syncl_cp.su ischeme=1 \
+ file_den=syncl_ro.su \
+ file_src=wave.su \
+ file_rcv=shot_xsrc${xsrc}.su \
+ src_type=7 \
+ src_orient=1 \
+ src_injectionrate=1 \
+ rec_type_vz=0 \
+ rec_type_p=1 \
+ rec_int_vz=2 \
+ dtrcv=0.004 \
+ rec_delay=0.1 \
+ verbose=2 \
+ tmod=4.10 \
+ dxrcv=10.0 \
+ xrcv1=-2250 xrcv2=2250 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=1150 \
+ ntaper=200 \
+ left=4 right=4 top=1 bottom=4
+
diff --git a/fdelmodc/getWaveletHeaders.c b/fdelmodc/getWaveletHeaders.c
index 0f291e3c32fc8bf40a22ad4972922d2fa8e62f1a..5bff37528015722251741fcdb434db218e06ed90 100644
--- a/fdelmodc/getWaveletHeaders.c
+++ b/fdelmodc/getWaveletHeaders.c
@@ -18,7 +18,8 @@ int getWaveletHeaders(char *file_src, int n1, int n2, float *gx, float *sx, floa
FILE *fp;
size_t nread;
int ix;
- off_t offset, trace_sz;
+ size_t trace_sz;
+ off_t offset;
float scl, scll;
segy hdr;
@@ -33,7 +34,7 @@ int getWaveletHeaders(char *file_src, int n1, int n2, float *gx, float *sx, floa
if (hdr.scalel < 0) scll = 1.0/fabs(hdr.scalel);
else if (hdr.scalel == 0) scll = 1.0;
else scll = hdr.scalel;
- trace_sz = sizeof(float)*(n1)+TRCBYTES;
+ trace_sz = (size_t)sizeof(float)*(n1)+TRCBYTES;
for (ix=0; ix<n2; ix++) {
offset = ix*trace_sz;
diff --git a/marchenko/applyMute.c b/marchenko/applyMute.c
index 4ba1798d4caa084cfeba66a8f93c8e4abd3d8fdb..4dd77b58dbbfaffba021d80f9fc56d82a6a21594 100644
--- a/marchenko/applyMute.c
+++ b/marchenko/applyMute.c
@@ -14,11 +14,11 @@
void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs, int nt, int *xrcvsyn, int npossyn, int shift)
{
- int i, j, l, isyn;
- float *costaper, scl;
- int imute, tmute;
+ int i, j, l, isyn;
+ float *costaper, scl;
+ int imute, tmute;
- if (smooth) {
+ if (smooth) {
costaper = (float *)malloc(smooth*sizeof(float));
scl = M_PI/((float)smooth);
for (i=0; i<smooth; i++) {
@@ -26,11 +26,11 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs
}
}
- for (isyn = 0; isyn < Nsyn; isyn++) {
- if (above==1) {
+ for (isyn = 0; isyn < Nsyn; isyn++) {
+ if (above==1) {
for (i = 0; i < npossyn; i++) {
- imute = xrcvsyn[i];
- tmute = mute[isyn*nxs+imute];
+ imute = xrcvsyn[i];
+ tmute = mute[isyn*nxs+imute];
for (j = 0; j < tmute-shift-smooth; j++) {
data[isyn*nxs*nt+i*nt+j] = 0.0;
}
@@ -41,8 +41,12 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs
}
else if (above==0){
for (i = 0; i < npossyn; i++) {
- imute = xrcvsyn[i];
- tmute = mute[isyn*nxs+imute];
+ imute = xrcvsyn[i];
+ tmute = mute[isyn*nxs+imute];
+ if (tmute >= nt/2) {
+ memset(&data[isyn*nxs*nt+i*nt],0, sizeof(float)*nt);
+ continue;
+ }
for (j = tmute-shift,l=0; j < tmute-shift+smooth; j++,l++) {
data[isyn*nxs*nt+i*nt+j] *= costaper[l];
}
@@ -56,8 +60,8 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs
}
else if (above==-1){
for (i = 0; i < npossyn; i++) {
- imute = xrcvsyn[i];
- tmute = mute[isyn*nxs+imute];
+ imute = xrcvsyn[i];
+ tmute = mute[isyn*nxs+imute];
for (j = tmute-shift,l=0; j < tmute-shift+smooth; j++,l++) {
data[isyn*nxs*nt+i*nt+j] *= costaper[l];
}
@@ -66,7 +70,7 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs
}
}
}
- else if (above==4) { //Psi gate which is the inverse of the Theta gate (above=0)
+ else if (above==4) { //Psi gate which is the inverse of the Theta gate (above=0)
for (i = 0; i < npossyn; i++) {
imute = xrcvsyn[i];
tmute = mute[isyn*nxs+imute];
@@ -104,8 +108,8 @@ void applyMute( float *data, int *mute, int smooth, int above, int Nsyn, int nxs
}
}
- if (smooth) free(costaper);
+ if (smooth) free(costaper);
- return;
+ return;
}
diff --git a/marchenko/demo/oneD/README b/marchenko/demo/oneD/README
index 22d2b75a71abf822c42eea4bc017068ed88ccf51..16b05691a7a83ab74c67d86e852c2f6d591a67f8 100644
--- a/marchenko/demo/oneD/README
+++ b/marchenko/demo/oneD/README
@@ -159,7 +159,6 @@ The postscript files of Figure 8 are generated with
==> run epsBackprop.scr
-
-- Figure 8 column 1
backpropf2_-0.30.eps
backpropf2_-0.15.eps
@@ -183,5 +182,11 @@ The figures in the appendix, to explain the different options in the programs, a
==> run figAppendi.scr
-- Figure A-1
+noise_above0.eps
+noise_above1.eps
+noise_above-1.eps
+noise_above2.eps
+noise_above4.eps
-- Figure A-2
+iniFocus_shifts.eps
diff --git a/marchenko/demo/oneD/figAppendi.scr b/marchenko/demo/oneD/figAppendi.scr
new file mode 100755
index 0000000000000000000000000000000000000000..22475c175d060aeab2dee98e8aad5fcc5b4dd63f
--- /dev/null
+++ b/marchenko/demo/oneD/figAppendi.scr
@@ -0,0 +1,43 @@
+#!/bin/bash
+
+file=iter_002.su
+file_base=${file%.su}
+
+ns=`surange < $file | grep ns | awk '{print $2}'`
+dtrcv=`surange < $file | grep dt | awk '{print $2/1000000.0}'`
+shift=$(echo "scale=4; ($dtrcv*($ns/2.0-1))" | bc -l)
+suzero < $file itmax=$ns | suaddnoise | sushw key=f1 a=0 > noise.su
+file_base=noise
+sumax < ${file_base}.su mode=abs outpar=nep
+clipiter=`cat nep | awk '{print $1/6}'`
+clipref=$clipiter
+
+#basop choice=shift shift=$shift file_in=$file file_out=${file_base}_t0.su
+
+for above in 0 1 -1 2 4
+do
+fmute file_mute=iniFocus_rp.su file_shot=${file_base}.su file_out=nep.su above=${above} shift=8 verbose=1 check=1 hw=4
+
+basop choice=shift shift=-$shift file_in=nep.su file_out=nep_t0.su
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < nep.su \
+ n1tic=2 d2=5 x1beg=0 d1num=0.5 \
+ curve=pslinepos.asci,pslineneg.asci npair=901,901 curvewidth=2,2 curvecolor=black,black curvedash=3,3 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > ${file_base}_above${above}.eps
+done
+
+for shift in 0 20 -20
+do
+fmute file_mute=iniFocus_rp.su file_shot=${file_base}.su file_out=nep.su above=${above} shift=$shift verbose=1 check=1 hw=4
+mv pslinepos.asci pslinepos${shift}.asci
+done
+
+suzero < $file itmax=$ns | sushw key=f1 a=0 > zero.su
+sumax < iniFocus_rp.su mode=abs outpar=nep
+clipiter=`cat nep | awk '{print $1/6}'`
+clipref=$clipiter
+supsimage hbox=6 wbox=4 labelsize=10 linewidth=0.0 < iniFocus_rp.su \
+ n1tic=2 d2=5 x1beg=0 d1num=0.5 \
+ curve=pslinepos0.asci,pslinepos15.asci,pslinepos-15.asci npair=901,901,901 \
+ curvewidth=1,1,1 curvecolor=white,black,black curvedash=3,3,3 \
+ f2=-2250 f2num=-2000 d2num=1000 clip=$clipref > iniFocus_shifts.eps
+
diff --git a/marchenko/fmute.c b/marchenko/fmute.c
index 5abbae6d0644bb9206b387bdd76c583229b51fe3..6554a39d72585a9caad17119e0bfadf6f26599c4 100644
--- a/marchenko/fmute.c
+++ b/marchenko/fmute.c
@@ -60,36 +60,36 @@ NULL};
int main (int argc, char **argv)
{
- FILE *fp_in1, *fp_in2, *fp_out, *fp_chk, *fp_psline1, *fp_psline2;
- int verbose, shift, k, nx1, nt1, nx2, nt2;
- int ntmax, nxmax, ret, i, j, jmax, imax, above, check;
- int size, ntraces, ngath, *maxval, hw, smooth;
+ FILE *fp_in1, *fp_in2, *fp_out, *fp_chk, *fp_psline1, *fp_psline2;
+ int verbose, shift, k, nx1, nt1, nx2, nt2;
+ int ntmax, nxmax, ret, i, j, jmax, imax, above, check;
+ int size, ntraces, ngath, *maxval, hw, smooth;
int tstart, tend, scale, *xrcv;
- float dt, d2, f1, f2, t0, t1, f1b, f2b, d1, d1b, d2b;
- float w1, w2, dxrcv;
- float *tmpdata, *tmpdata2, *costaper;
- char *file_mute, *file_shot, *file_out;
- float scl, sclsxgx, sclshot, xmin, xmax, tmax, lmax;
- segy *hdrs_in1, *hdrs_in2;
-
- t0 = wallclock_time();
- initargs(argc, argv);
- requestdoc(1);
-
- 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;
+ float dt, d2, f1, f2, t0, t1, f1b, f2b, d1, d1b, d2b;
+ float w1, w2, dxrcv;
+ float *tmpdata, *tmpdata2, *costaper;
+ char *file_mute, *file_shot, *file_out;
+ float scl, sclsxgx, sclshot, xmin, xmax, tmax, lmax;
+ segy *hdrs_in1, *hdrs_in2;
+
+ t0 = wallclock_time();
+ initargs(argc, argv);
+ requestdoc(1);
+
+ 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(!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(!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;
/* Reading input data for file_mute */
@@ -101,7 +101,7 @@ int main (int argc, char **argv)
if (!getparint("nxmax", &nxmax)) nxmax = nx1;
if (verbose>=2 && (ntmax!=nt1 || nxmax!=nx1))
vmess("dimensions overruled: %d x %d",ntmax,nxmax);
- if(!getparfloat("dt", &dt)) dt=d1;
+ if(!getparfloat("dt", &dt)) dt=d1;
fp_in1 = fopen(file_mute, "r");
if (fp_in1 == NULL) verr("error on opening input file_mute=%s", file_mute);
@@ -123,33 +123,33 @@ int main (int argc, char **argv)
/* Reading input data for file_shot */
- ngath = 1;
- getFileInfo(file_shot, &nt2, &nx2, &ngath, &d1b, &d2b, &f1b, &f2b, &xmin, &xmax, &sclshot, &ntraces);
-
- if (!getparint("ntmax", &ntmax)) ntmax = nt2;
- if (!getparint("nxmax", &nxmax)) nxmax = nx2;
-
- size = ntmax * nxmax;
- tmpdata2 = (float *)malloc(size*sizeof(float));
- hdrs_in2 = (segy *) calloc(nxmax,sizeof(segy));
-
- if (file_shot != NULL) fp_in2 = fopen(file_shot, "r");
- else fp_in2=stdin;
- if (fp_in2 == NULL) verr("error on opening input file_shot=%s", file_shot);
-
- nx2 = readData(fp_in2, tmpdata2, hdrs_in2, nt2);
- if (nx2 == 0) {
- fclose(fp_in2);
- if (verbose) vmess("end of file_shot data reached");
- }
- nt2 = hdrs_in2[0].ns;
- f1b = hdrs_in2[0].f1;
- f2b = hdrs_in2[0].f2;
- d1b = (float)hdrs_in2[0].dt*1e-6;
-
- if (verbose) {
- disp_fileinfo(file_shot, nt2, nx2, f1b, f2b, d1b, d2b, hdrs_in2);
- }
+ ngath = 1;
+ getFileInfo(file_shot, &nt2, &nx2, &ngath, &d1b, &d2b, &f1b, &f2b, &xmin, &xmax, &sclshot, &ntraces);
+
+ if (!getparint("ntmax", &ntmax)) ntmax = nt2;
+ if (!getparint("nxmax", &nxmax)) nxmax = nx2;
+
+ size = ntmax * nxmax;
+ tmpdata2 = (float *)malloc(size*sizeof(float));
+ hdrs_in2 = (segy *) calloc(nxmax,sizeof(segy));
+
+ if (file_shot != NULL) fp_in2 = fopen(file_shot, "r");
+ else fp_in2=stdin;
+ if (fp_in2 == NULL) verr("error on opening input file_shot=%s", file_shot);
+
+ nx2 = readData(fp_in2, tmpdata2, hdrs_in2, nt2);
+ if (nx2 == 0) {
+ fclose(fp_in2);
+ if (verbose) vmess("end of file_shot data reached");
+ }
+ nt2 = hdrs_in2[0].ns;
+ f1b = hdrs_in2[0].f1;
+ f2b = hdrs_in2[0].f2;
+ d1b = (float)hdrs_in2[0].dt*1e-6;
+
+ if (verbose) {
+ disp_fileinfo(file_shot, nt2, nx2, f1b, f2b, d1b, d2b, hdrs_in2);
+ }
/* file_shot will be used as well to define the mute window */
if (file_mute == NULL) {
@@ -163,41 +163,41 @@ 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=%d, file_shot=%d", nt1, nt2);
/*================ initializations ================*/
- maxval = (int *)calloc(nx1,sizeof(int));
- xrcv = (int *)calloc(nx1,sizeof(int));
-
- if (file_out==NULL) fp_out = stdout;
- else {
- fp_out = fopen(file_out, "w+");
- if (fp_out==NULL) verr("error on ceating output file");
- }
+ maxval = (int *)calloc(nx1,sizeof(int));
+ xrcv = (int *)calloc(nx1,sizeof(int));
+
+ if (file_out==NULL) fp_out = stdout;
+ else {
+ fp_out = fopen(file_out, "w+");
+ if (fp_out==NULL) verr("error on ceating output file");
+ }
if (check!=0){
fp_chk = fopen("check.su", "w+");
- if (fp_chk==NULL) verr("error on ceating output file");
+ if (fp_chk==NULL) verr("error on ceating output file");
fp_psline1 = fopen("pslinepos.asci", "w+");
- if (fp_psline1==NULL) verr("error on ceating output file");
+ if (fp_psline1==NULL) verr("error on ceating output file");
fp_psline2 = fopen("pslineneg.asci", "w+");
- if (fp_psline2==NULL) verr("error on ceating output file");
+ if (fp_psline2==NULL) verr("error on ceating output file");
}
- 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));
-/* fprintf(stderr,"costaper[%d]=%f\n",i,costaper[i]);*/
- }
- }
+ 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));
+/* fprintf(stderr,"costaper[%d]=%f\n",i,costaper[i]);*/
+ }
+ }
/*================ loop over all shot records ================*/
- k=1;
- while (nx1 > 0) {
- if (verbose) vmess("processing input gather %d", k);
+ k=1;
+ while (nx1 > 0) {
+ if (verbose) vmess("processing input gather %d", k);
/*================ loop over all shot records ================*/
@@ -205,42 +205,42 @@ int main (int argc, char **argv)
/* find global maximum
xmax=0.0;
- for (i = 0; i < nx1; i++) {
- tmax=0.0;
- jmax = 0;
- for (j = 0; j < nt1; j++) {
+ for (i = 0; i < nx1; i++) {
+ tmax=0.0;
+ jmax = 0;
+ for (j = 0; j < nt1; j++) {
lmax = fabs(tmpdata[i*nt1+j]);
- if (lmax > tmax) {
- jmax = j;
- tmax = lmax;
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
if (lmax > xmax) {
imax = i;
xmax=lmax;
}
- }
- }
- maxval[i] = jmax;
- }
- */
+ }
+ }
+ maxval[i] = jmax;
+ }
+ */
- /* alternative find maximum at source position */
+ /* alternative find maximum at source position */
dxrcv = (hdrs_in1[nx1-1].gx - hdrs_in1[0].gx)*sclsxgx/(float)(nx1-1);
- imax = NINT(((hdrs_in1[0].sx-hdrs_in1[0].gx)*sclsxgx)/dxrcv);
- tmax=0.0;
- jmax = 0;
- xmax=0.0;
- for (j = 0; j < nt1; j++) {
+ imax = NINT(((hdrs_in1[0].sx-hdrs_in1[0].gx)*sclsxgx)/dxrcv);
+ tmax=0.0;
+ jmax = 0;
+ xmax=0.0;
+ for (j = 0; j < nt1; j++) {
lmax = fabs(tmpdata[imax*nt1+j]);
- if (lmax > tmax) {
- jmax = j;
- tmax = lmax;
+ if (lmax > tmax) {
+ jmax = j;
+ tmax = lmax;
if (lmax > xmax) {
xmax=lmax;
}
- }
- }
- maxval[imax] = jmax;
- if (verbose >= 3) vmess("Mute max at src-trace %d is sample %d", imax, maxval[imax]);
+ }
+ }
+ maxval[imax] = jmax;
+ if (verbose >= 3) vmess("Mute max at src-trace %d is sample %d", imax, maxval[imax]);
/* search forward */
for (i = imax+1; i < nx1; i++) {
@@ -275,46 +275,46 @@ int main (int argc, char **argv)
/* scale with maximum ampltiude */
- if (scale==1) {
- for (i = 0; i < nx2; i++) {
- lmax = fabs(tmpdata2[i*nt2+maxval[i]]);
- for (j = 0; j < nt2; j++) {
- tmpdata2[i*nt2+j] = tmpdata2[i*nt2+j]/lmax;
- }
- }
- }
+ if (scale==1) {
+ for (i = 0; i < nx2; i++) {
+ lmax = fabs(tmpdata2[i*nt2+maxval[i]]);
+ for (j = 0; j < nt2; j++) {
+ tmpdata2[i*nt2+j] = tmpdata2[i*nt2+j]/lmax;
+ }
+ }
+ }
for (i = 0; i < nx2; i++) xrcv[i] = i;
/*================ apply mute window ================*/
- applyMute(tmpdata2, maxval, smooth, above, 1, nx2, nt2, xrcv, nx2, shift);
+ applyMute(tmpdata2, maxval, smooth, above, 1, nx2, nt2, xrcv, nx2, shift);
/*================ write result to output file ================*/
- ret = writeData(fp_out, tmpdata2, hdrs_in2, nt2, nx2);
- if (ret < 0 ) verr("error on writing output file.");
+ ret = writeData(fp_out, tmpdata2, hdrs_in2, nt2, nx2);
+ if (ret < 0 ) verr("error on writing output file.");
/* put mute window in file to check correctness of mute */
if (check !=0) {
for (i = 0; i < nx1; i++) {
jmax = maxval[i]-shift;
tmpdata[i*nt1+jmax] = 2*xmax;
- }
- if (above==0){
- for (i = 0; i < nx1; i++) {
- jmax = nt2-maxval[i]+shift;
- tmpdata[i*nt1+jmax] = 2*xmax;
- }
- }
+ }
+ if (above==0){
+ for (i = 0; i < nx1; i++) {
+ jmax = nt2-maxval[i]+shift;
+ tmpdata[i*nt1+jmax] = 2*xmax;
+ }
+ }
ret = writeData(fp_chk, tmpdata, hdrs_in1, nt1, nx1);
if (ret < 0 ) verr("error on writing check file.");
- for (i=0; i<nx1; i++) {
- jmax = maxval[i]-shift;
- ret = fprintf(fp_psline1, "%.5f %.5f \n",jmax*dt,hdrs_in1[i].gx*sclshot);
- jmax =-maxval[i]+shift;
- ret = fprintf(fp_psline2, "%.5f %.5f \n",jmax*dt,hdrs_in1[i].gx*sclshot);
- }
+ for (i=0; i<nx1; i++) {
+ jmax = maxval[i]-shift;
+ ret = fprintf(fp_psline1, "%.5f %.5f \n",jmax*dt,hdrs_in1[i].gx*sclshot);
+ jmax =-maxval[i]+shift;
+ ret = fprintf(fp_psline2, "%.5f %.5f \n",jmax*dt,hdrs_in1[i].gx*sclshot);
+ }
}
/*================ Read next record for muting ================*/
@@ -328,9 +328,9 @@ int main (int argc, char **argv)
if (fp_out!=stdout) fclose(fp_out);
if (check!=0) fclose(fp_chk);
if (check!=0) {
- fclose(fp_psline1);
- fclose(fp_psline2);
- }
+ fclose(fp_psline1);
+ fclose(fp_psline2);
+ }
break;
}
nt1 = (int)hdrs_in1[0].ns;
@@ -353,17 +353,17 @@ int main (int argc, char **argv)
if (file_mute == NULL) {
nx1=nx2;
nt1=nt2;
- hdrs_in1 = hdrs_in2;
+ hdrs_in1 = hdrs_in2;
tmpdata = tmpdata2;
}
- k++;
- }
+ k++;
+ }
- t1 = wallclock_time();
- if (verbose) vmess("Total CPU-time = %f",t1-t0);
-
+ t1 = wallclock_time();
+ if (verbose) vmess("Total CPU-time = %f",t1-t0);
+
- return 0;
+ return 0;
}
diff --git a/marchenko/marchenko.c b/marchenko/marchenko.c
index d0af00dd3776fbe57d046ac49504e86c0d576042..dc31009bcdc5494e9885111662e29dd84313bfec 100644
--- a/marchenko/marchenko.c
+++ b/marchenko/marchenko.c
@@ -25,7 +25,7 @@ typedef struct _complexStruct { /* complex number */
} complex;
#endif/* complex */
-int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int ngath, int nx, int nxm, int ntfft, int mode, float weight, int verbose);
+int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int ngath, int nx, int nxm, int ntfft, int mode, float weight, float tsq, int verbose);
int readTinvData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nsyn, 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 Nsyn, float *xsyn, float *zsyn, int iter);
void name_ext(char *filename, char *extension);
@@ -63,11 +63,13 @@ char *sdoc[] = {
" fmax=70 .................. maximum frequency",
" MARCHENKO ITERATIONS ",
" niter=10 ................. number of iterations",
-" MUTE WINDOW ",
+" MUTE-WINDOW ",
" above=0 .................. mute above(1), around(0) or below(-1) the first travel times of file_tinv",
" shift=12 ................. number of points above(positive) / below(negative) travel time for mute",
" hw=8 ..................... window in time samples to look for maximum in next trace",
" smooth=5 ................. number of points to smooth mute with cosine window",
+" REFLECTION RESPONSE CORRECTION ",
+" tsq=0.0 .................. weight factor n for t^n for true amplitude recovery",
" weight=1 ................. weight factor of R for summation of Ni with G_d",
" OUTPUT DEFINITION ",
" file_green= .............. output file with full Green function(s)",
@@ -102,7 +104,7 @@ int main (int argc, char **argv)
float d1, d2, f1, f2, fxs, ft, fx, *xsyn, dxsrc;
float *green, *f2p, *pmin, *G_d, dt, dx, dxs, scl, mem;
float *f1plus, *f1min, *iRN, *Ni, *trace, *Gmin, *Gplus;
- float xmin, xmax, weight;
+ float xmin, xmax, weight, tsq;
complex *Refl, *Fop;
char *file_tinv, *file_shot, *file_green, *file_iter;
char *file_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_pmin;
@@ -136,8 +138,9 @@ int main (int argc, char **argv)
if (!getparint("ixa", &ixa)) ixa = 0;
if (!getparint("ixb", &ixb)) ixb = ixa;
// if (!getparint("reci", &reci)) reci = 0;
- reci=0; // source-receiver reciprocity is not yet fully build into the code
+ reci=0; // source-receiver reciprocity is not yet fully build into the code
if (!getparfloat("weight", &weight)) weight = 1.0;
+ if (!getparfloat("tsq", &tsq)) tsq = 0.0;
if (!getparint("tap", &tap)) tap = 0;
if (!getparint("ntap", &ntap)) ntap = 0;
@@ -162,7 +165,7 @@ int main (int argc, char **argv)
ngath = 0; /* setting ngath=0 scans all traces; nx contains maximum traces/gather */
ret = getFileInfo(file_shot, &nt, &nx, &ngath, &d1, &dx, &ft, &fx, &xmin, &xmax, &scl, &ntraces);
nshots = ngath;
- assert (nxs >= nshots);
+ assert (nxs >= nshots);
if (!getparfloat("dt", &dt)) dt = d1;
@@ -207,10 +210,10 @@ int main (int argc, char **argv)
mode=-1; /* apply complex conjugate to read in data */
readTinvData(file_tinv, xrcvsyn, xsyn, zsyn, xnxsyn, Nsyn, nxs, ntfft,
mode, muteW, G_d, hw, verbose);
- /* reading data added zero's to the number of time samples to be the same as ntfft */
+ /* reading data added zero's to the number of time samples to be the same as ntfft */
nts = ntfft;
- /* define tapers to taper edges of acquisition */
+ /* define tapers to taper edges of acquisition */
if (tap == 1 || tap == 3) {
for (j = 0; j < ntap; j++)
tapersy[j] = (cos(PI*(j-ntap)/ntap)+1)/2.0;
@@ -233,7 +236,7 @@ int main (int argc, char **argv)
}
}
- /* check consistency of header values */
+ /* check consistency of header values */
if (xrcvsyn[0] != 0 || xrcvsyn[1] != 0 ) fxs = xrcvsyn[0];
fxs2 = fxs + (float)(nxs-1)*dxs;
dxf = (xrcvsyn[nxs-1] - xrcvsyn[0])/(float)(nxs-1);
@@ -247,7 +250,7 @@ int main (int argc, char **argv)
mode=1;
readShotData(file_shot, xrcv, xsrc, zsrc, xnx, Refl, nw, nw_low, ngath, nx, nx, ntfft,
- mode, weight, verbose);
+ mode, weight, tsq, verbose);
tapersh = (float *)malloc(nx*sizeof(float));
if (tap == 2 || tap == 3) {
@@ -274,7 +277,7 @@ int main (int argc, char **argv)
}
free(tapersh);
- /* check consistency of header values */
+ /* check consistency of header values */
fxf = xsrc[0];
if (nx > 1) dxf = (xrcv[0] - xrcv[nx-1])/(float)(nx-1);
else dxf = d2;
@@ -409,7 +412,7 @@ int main (int argc, char **argv)
}
}
- /* apply mute window based on times of direct arrival (in muteW) */
+ /* apply mute window based on times of direct arrival (in muteW) */
applyMute(Ni, muteW, smooth, above, Nsyn, nxs, nts, ixpossyn, npossyn, shift);
/* initialization */
@@ -519,6 +522,11 @@ int main (int argc, char **argv)
for (l = 0; l < Nsyn; l++) {
for (i = 0; i < npossyn; i++) {
j = 0;
+ /* set green to zero if mute-window exceeds nt/2 */
+ if (muteW[l*nxs+ixpossyn[i]] >= nts/2) {
+ memset(&green[l*nxs*nts+i*nts],0, sizeof(float)*nt);
+ continue;
+ }
green[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+j] + pmin[l*nxs*nts+i*nts+j];
for (j = 1; j < nts; j++) {
green[l*nxs*nts+i*nts+j] = f2p[l*nxs*nts+i*nts+nts-j] + pmin[l*nxs*nts+i*nts+j];
@@ -531,7 +539,7 @@ int main (int argc, char **argv)
Gmin = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
/* use f1+ as operator on R in frequency domain */
- mode=1;
+ mode=1;
synthesis(Refl, Fop, f1plus, iRN, nx, nt, nxs, nts, dt, xsyn, Nsyn,
xrcv, xsrc, fxs2, fxs, dxs, dxsrc, dx, ixa, ixb, ntfft, nw, nw_low, nw_high, mode,
reci, nshots, ixpossyn, npossyn, &tfft, verbose);
@@ -555,7 +563,7 @@ int main (int argc, char **argv)
Gplus = (float *)calloc(Nsyn*nxs*ntfft,sizeof(float));
/* use f1-(*) as operator on R in frequency domain */
- mode=-1;
+ mode=-1;
synthesis(Refl, Fop, f1min, iRN, nx, nt, nxs, nts, dt, xsyn, Nsyn,
xrcv, xsrc, fxs2, fxs, dxs, dxsrc, dx, ixa, ixb, ntfft, nw, nw_low, nw_high, mode,
reci, nshots, ixpossyn, npossyn, &tfft, verbose);
@@ -639,7 +647,7 @@ int main (int argc, char **argv)
hdrs_out[i].tracf = tracf++;
hdrs_out[i].selev = NINT(zsyn[l]*1000);
hdrs_out[i].sdepth = NINT(-zsyn[l]*1000);
- hdrs_out[i].f1 = f1;
+ hdrs_out[i].f1 = f1;
}
ret = writeData(fp_out, (float *)&green[l*size], hdrs_out, n1, n2);
@@ -725,7 +733,7 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
float *rtrace, idxs;
complex *sum, *ctrace;
int npe;
- static int first=1, *ixrcv;
+ static int first=1, *ixrcv;
static double t0, t1, t;
size = nxs*nts;
@@ -749,47 +757,47 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
/* reset output data to zero */
memset(&iRN[0], 0, Nsyn*nxs*nts*sizeof(float));
- ctrace = (complex *)calloc(ntfft,sizeof(complex));
- if (!first) {
+ ctrace = (complex *)calloc(ntfft,sizeof(complex));
+ if (!first) {
/* transform muted Ni (Top) to frequency domain, input for next iteration */
- for (l = 0; l < Nsyn; l++) {
- /* set Fop to zero, so new operator can be defined within ixpossyn points */
+ for (l = 0; l < Nsyn; l++) {
+ /* set Fop to zero, so new operator can be defined within ixpossyn points */
memset(&Fop[l*nxs*nw].r, 0, nxs*nw*2*sizeof(float));
for (i = 0; i < npossyn; i++) {
- rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
- ix = ixpossyn[i];
- 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;
- }
+ rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
+ ix = ixpossyn[i];
+ 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;
+ }
}
- }
- }
- else { /* only for first call to synthesis */
+ }
+ }
+ else { /* only for first call to synthesis */
/* transform G_d to frequency domain, over all nxs traces */
- first=0;
- for (l = 0; l < Nsyn; l++) {
- /* set Fop to zero, so new operator can be defined within all ix points */
+ first=0;
+ for (l = 0; l < Nsyn; l++) {
+ /* set Fop to zero, so new operator can be defined within all ix points */
memset(&Fop[l*nxs*nw].r, 0, nxs*nw*2*sizeof(float));
for (i = 0; i < nxs; i++) {
- rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
- for (iw=0; iw<nw; iw++) {
- Fop[l*nxs*nw+iw*nxs+i].r = ctrace[nw_low+iw].r;
- Fop[l*nxs*nw+iw*nxs+i].i = mode*ctrace[nw_low+iw].i;
- }
+ rc1fft(&Top[l*size+i*nts],ctrace,ntfft,-1);
+ for (iw=0; iw<nw; iw++) {
+ Fop[l*nxs*nw+iw*nxs+i].r = ctrace[nw_low+iw].r;
+ Fop[l*nxs*nw+iw*nxs+i].i = mode*ctrace[nw_low+iw].i;
+ }
}
- }
- idxs = 1.0/dxs;
- ixrcv = (int *)malloc(nshots*nx*sizeof(int));
- for (k=0; k<nshots; k++) {
+ }
+ idxs = 1.0/dxs;
+ ixrcv = (int *)malloc(nshots*nx*sizeof(int));
+ for (k=0; k<nshots; k++) {
for (i = 0; i < nx; i++) {
ixrcv[k*nx+i] = NINT((xrcv[k*nx+i]-fxs)*idxs);
- }
- }
- }
- free(ctrace);
+ }
+ }
+ }
+ free(ctrace);
t1 = wallclock_time();
- *tfft += t1 - t0;
+ *tfft += t1 - t0;
for (k=0; k<nshots; k++) {
@@ -808,7 +816,7 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
}
- iox = 0; inx = nx;
+ iox = 0; inx = nx;
/*================ SYNTHESIS ================*/
@@ -828,11 +836,11 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
ix = k;
- /* multiply R with Fop and sum over nx */
- memset(&sum[0].r,0,nfreq*2*sizeof(float));
+ /* multiply R with Fop and sum over nx */
+ memset(&sum[0].r,0,nfreq*2*sizeof(float));
//for (j = 0; j < nfreq; j++) sum[j].r = sum[j].i = 0.0;
- for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
- for (i = iox; i < inx; i++) {
+ for (j = nw_low, m = 0; j <= nw_high; j++, m++) {
+ for (i = iox; i < inx; i++) {
sum[j].r += Refl[k*nw*nx+m*nx+i].r*Fop[l*nw*nxs+m*nxs+ixrcv[k*nx+i]].r -
Refl[k*nw*nx+m*nx+i].i*Fop[l*nw*nxs+m*nxs+ixrcv[k*nx+i]].i;
sum[j].i += Refl[k*nw*nx+m*nx+i].i*Fop[l*nw*nxs+m*nxs+ixrcv[k*nx+i]].r +
@@ -840,7 +848,7 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
}
}
- /* transfrom result back to time domain */
+ /* transfrom result back to time domain */
cr1fft(sum, rtrace, ntfft, 1);
/* dx = receiver distance */
@@ -860,7 +868,7 @@ void synthesis(complex *Refl, complex *Fop, float *Top, float *iRN, int nx, int
}
} /* end of parallel region */
- if (verbose>3) vmess("*** Shot gather %d processed ***", k);
+ if (verbose>3) vmess("*** Shot gather %d processed ***", k);
} /* end of nshots (k) loop */
@@ -902,7 +910,7 @@ void synthesisPosistions(int nx, int nt, int nxs, int nts, float dt, float *xsyn
break;
}
- iox = 0; inx = nx;
+ iox = 0; inx = nx;
if (ixa || ixb) {
if (reci == 0) {
@@ -944,7 +952,7 @@ void synthesisPosistions(int nx, int nt, int nxs, int nts, float dt, float *xsyn
*npossyn += 1;
}
if (verbose>=3) {
- vmess("ixpossyn[%d] = %d ixsrc=%d ix=%d", *npossyn-1, ixpossyn[*npossyn-1], ixsrc, ix);
+ vmess("ixpossyn[%d] = %d ixsrc=%d ix=%d", *npossyn-1, ixpossyn[*npossyn-1], ixsrc, ix);
}
if (reci == 1 || reci == 2) {
@@ -972,28 +980,28 @@ void update(float *field, float *term, int Nsyn, int nx, int nt, int reverse, in
{
int i, j, l, ix;
- if (reverse) {
- for (l = 0; l < Nsyn; l++) {
- for (i = 0; i < npossyn; i++) {
- j = 0;
- Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+nts-j];
- }
- }
- }
- }
- else {
- for (l = 0; l < Nsyn; l++) {
- for (i = 0; i < npossyn; i++) {
- j = 0;
- Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+j];
- for (j = 1; j < nts; j++) {
- Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+nts-j];
- }
- }
- }
- }
- return;
+ if (reverse) {
+ for (l = 0; l < Nsyn; l++) {
+ for (i = 0; i < npossyn; i++) {
+ j = 0;
+ Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+nts-j];
+ }
+ }
+ }
+ }
+ else {
+ for (l = 0; l < Nsyn; l++) {
+ for (i = 0; i < npossyn; i++) {
+ j = 0;
+ Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+j];
+ for (j = 1; j < nts; j++) {
+ Ni[l*nxs*nts+i*nts+j] = -iRN[l*nxs*nts+i*nts+nts-j];
+ }
+ }
+ }
+ }
+ return;
}
*/
diff --git a/marchenko/readShotData.c b/marchenko/readShotData.c
index 83f83e08a13e0d633487f575187cde9e8cb09606..5b99e601edefe8783ff0ae6d005a85151c25015b 100644
--- a/marchenko/readShotData.c
+++ b/marchenko/readShotData.c
@@ -18,108 +18,116 @@ void rc1fft(float *rdata, complex *cdata, int n, int sign);
int compare(const void *a, const void *b)
{ return (*(float *)b-*(float *)a); }
-int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int ngath, int nx, int nxm, int ntfft, int mode, float weight, int verbose)
+int readShotData(char *filename, float *xrcv, float *xsrc, float *zsrc, int *xnx, complex *cdata, int nw, int nw_low, int ngath, int nx, int nxm, int ntfft, int mode, float weight, float tsq, int verbose)
{
- FILE *fp;
- segy hdr;
- size_t nread;
- int fldr_shot, sx_shot, itrace, one_shot, igath, iw;
- int end_of_file, nt;
- float scl, scel, *trace;
- complex *ctrace;
-
- /* Reading first header */
-
- if (filename == NULL) fp = stdin;
- else fp = fopen( filename, "r" );
- if ( fp == NULL ) {
- fprintf(stderr,"input file %s has an error\n", filename);
- perror("error in opening file: ");
- fflush(stderr);
- return -1;
- }
-
- fseek(fp, 0, SEEK_SET);
- nread = fread( &hdr, 1, TRCBYTES, fp );
- assert(nread == TRCBYTES);
- if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
- else if (hdr.scalco == 0) scl = 1.0;
- else scl = hdr.scalco;
- if (hdr.scalel < 0) scel = 1.0/fabs(hdr.scalel);
- else if (hdr.scalel == 0) scel = 1.0;
- else scel = hdr.scalel;
-
- fseek(fp, 0, SEEK_SET);
-
- nt = hdr.ns;
-
- trace = (float *)calloc(ntfft,sizeof(float));
- ctrace = (complex *)malloc(ntfft*sizeof(complex));
-
- end_of_file = 0;
- one_shot = 1;
- igath = 0;
-
- /* Read shots in file */
-
- while (!end_of_file) {
-
- /* start reading data (shot records) */
- itrace = 0;
- nread = fread( &hdr, 1, TRCBYTES, fp );
- if (nread != TRCBYTES) { /* no more data in file */
- break;
- }
-
- sx_shot = hdr.sx;
- fldr_shot = hdr.fldr;
- xsrc[igath] = sx_shot*scl;
- zsrc[igath] = hdr.selev*scel;
- xnx[igath]=0;
- while (one_shot) {
- xrcv[igath*nxm+itrace] = hdr.gx*scl;
- nread = fread( trace, sizeof(float), nt, fp );
- assert (nread == hdr.ns);
-
- /* transform to frequency domain */
- if (ntfft > hdr.ns)
- memset( &trace[nt-1], 0, sizeof(float)*(ntfft-nt) );
-
- rc1fft(trace,ctrace,ntfft,-1);
- for (iw=0; iw<nw; iw++) {
- cdata[igath*nx*nw+iw*nx+itrace].r = weight*ctrace[nw_low+iw].r;
- cdata[igath*nx*nw+iw*nx+itrace].i = weight*mode*ctrace[nw_low+iw].i;
- }
- itrace++;
- xnx[igath]+=1;
-
- /* read next hdr of next trace */
- nread = fread( &hdr, 1, TRCBYTES, fp );
- if (nread != TRCBYTES) {
- one_shot = 0;
- end_of_file = 1;
- break;
- }
- if ((sx_shot != hdr.sx) || (fldr_shot != hdr.fldr) ) break;
- }
- if (verbose>2) {
- fprintf(stderr,"finished reading shot %d (%d) with %d traces\n",sx_shot,igath,itrace);
- //disp_fileinfo(filename, nt, xnx[igath], hdr.f1, xrcv[igath*nxm], d1, d2, &hdr);
- }
-
- if (itrace != 0) { /* end of shot record */
- fseek( fp, -TRCBYTES, SEEK_CUR );
- igath++;
- }
- else {
- end_of_file = 1;
- }
- }
-
- free(ctrace);
- free(trace);
-
- return 0;
+ FILE *fp;
+ segy hdr;
+ size_t nread;
+ int fldr_shot, sx_shot, itrace, one_shot, igath, iw;
+ int end_of_file, nt;
+ float scl, scel, *trace, dt;
+ complex *ctrace;
+
+ /* Reading first header */
+
+ if (filename == NULL) fp = stdin;
+ else fp = fopen( filename, "r" );
+ if ( fp == NULL ) {
+ fprintf(stderr,"input file %s has an error\n", filename);
+ perror("error in opening file: ");
+ fflush(stderr);
+ return -1;
+ }
+
+ fseek(fp, 0, SEEK_SET);
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ assert(nread == TRCBYTES);
+ if (hdr.scalco < 0) scl = 1.0/fabs(hdr.scalco);
+ else if (hdr.scalco == 0) scl = 1.0;
+ else scl = hdr.scalco;
+ if (hdr.scalel < 0) scel = 1.0/fabs(hdr.scalel);
+ else if (hdr.scalel == 0) scel = 1.0;
+ else scel = hdr.scalel;
+
+ fseek(fp, 0, SEEK_SET);
+
+ nt = hdr.ns;
+ dt = hdr.dt/(1E6);
+
+ trace = (float *)calloc(ntfft,sizeof(float));
+ ctrace = (complex *)malloc(ntfft*sizeof(complex));
+
+ end_of_file = 0;
+ one_shot = 1;
+ igath = 0;
+
+ /* Read shots in file */
+
+ while (!end_of_file) {
+
+ /* start reading data (shot records) */
+ itrace = 0;
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) { /* no more data in file */
+ break;
+ }
+
+ sx_shot = hdr.sx;
+ fldr_shot = hdr.fldr;
+ xsrc[igath] = sx_shot*scl;
+ zsrc[igath] = hdr.selev*scel;
+ xnx[igath]=0;
+ while (one_shot) {
+ xrcv[igath*nxm+itrace] = hdr.gx*scl;
+ nread = fread( trace, sizeof(float), nt, fp );
+ assert (nread == hdr.ns);
+
+ /* True Amplitude Recovery */
+ if (tsq != 0.0) {
+ for (iw=0; iw<nt; iw++) {
+ trace[iw] *= powf(dt*iw,tsq);
+ }
+ }
+
+ /* transform to frequency domain */
+ if (ntfft > hdr.ns)
+ memset( &trace[nt-1], 0, sizeof(float)*(ntfft-nt) );
+
+ rc1fft(trace,ctrace,ntfft,-1);
+ for (iw=0; iw<nw; iw++) {
+ cdata[igath*nx*nw+iw*nx+itrace].r = weight*ctrace[nw_low+iw].r;
+ cdata[igath*nx*nw+iw*nx+itrace].i = weight*mode*ctrace[nw_low+iw].i;
+ }
+ itrace++;
+ xnx[igath]+=1;
+
+ /* read next hdr of next trace */
+ nread = fread( &hdr, 1, TRCBYTES, fp );
+ if (nread != TRCBYTES) {
+ one_shot = 0;
+ end_of_file = 1;
+ break;
+ }
+ if ((sx_shot != hdr.sx) || (fldr_shot != hdr.fldr) ) break;
+ }
+ if (verbose>2) {
+ fprintf(stderr,"finished reading shot %d (%d) with %d traces\n",sx_shot,igath,itrace);
+ //disp_fileinfo(filename, nt, xnx[igath], hdr.f1, xrcv[igath*nxm], d1, d2, &hdr);
+ }
+
+ if (itrace != 0) { /* end of shot record */
+ fseek( fp, -TRCBYTES, SEEK_CUR );
+ igath++;
+ }
+ else {
+ end_of_file = 1;
+ }
+ }
+
+ free(ctrace);
+ free(trace);
+
+ return 0;
}