diff --git a/.gitignore b/.gitignore
index b67ba9f637168158a2e7faa8a9f36fbdffdb8fc1..3c9d60530844865ab966a8100bcac244d7f6d941 100644
--- a/.gitignore
+++ b/.gitignore
@@ -19,3 +19,4 @@ utils/mat2su
utils/syn2d
fdelmodc/fdelmodc
include/genfft.h
+Make_include
diff --git a/marchenko_full/AmpEst.c b/marchenko_full/AmpEst.c
index 2c95ed0833614e59aeda85ab989a5cb273ee6592..d8dbe82a51a7a6977e19bbacb59fd905aab147f5 100755
--- a/marchenko_full/AmpEst.c
+++ b/marchenko_full/AmpEst.c
@@ -62,6 +62,7 @@ int AmpEst(float *f1d, float *Gd, float *ampest, int Nsyn, int nxs, int ntfft, i
ampest[l] = Wmax/(Amax/((float)ntfft));
memset(&Af[0],0.0, sizeof(float)*2*nfreq);
}
+ free(Gdf);free(f1df);free(Af);free(At);free(wavelet);
return;
}
diff --git a/marchenko_full/Cost.c b/marchenko_full/Cost.c
index 57371da650ebda7ea69acd7d04d3dc2345629792..aa15b3a0e9ac95f245002ca271261c1c51e5884e 100755
--- a/marchenko_full/Cost.c
+++ b/marchenko_full/Cost.c
@@ -65,6 +65,7 @@ int Cost(float *f1p, float *f1d, float *Gm, float *Gm0, double *J, int Nsyn, int
R2L2 = 0.0;
R20L2 = 0.0;
}
+ free(f1pf);free(f1df);free(Gmf);free(Gm0f);free(R2f);free(R20f);free(R2);free(R20);
return;
}
diff --git a/marchenko_full/Makefile b/marchenko_full/Makefile
index 647e0734edea782194b070c50fd948322819cdd7..d738210371395a419f6d20399cdc427246cd10b8 100644
--- a/marchenko_full/Makefile
+++ b/marchenko_full/Makefile
@@ -37,7 +37,19 @@ SRCH = marchenko.c \
readSnapData.c \
Cost.c \
AmpEst.c \
- freqwave.c
+ freqwave.c \
+ getParameters.c \
+ getModelInfo.c \
+ recvPar.c \
+ raytime.c \
+ readModel.c \
+ JespersRayTracer.c \
+ getWaveletHeaders.c \
+ getWaveletInfo.c \
+ writeSrcRecPos.c \
+ writesufile.c \
+ gaussGen.c \
+ threadAffinity.c
OBJJ = $(SRCJ:%.c=%.o)
@@ -46,7 +58,7 @@ fmute: $(OBJJ)
OBJH = $(SRCH:%.c=%.o)
-marchenko: $(OBJH)
+marchenko: $(OBJH) raytime.h
$(CC) $(LDFLAGS) $(OPTC) $(CFLAGS) -o marchenko $(OBJH) $(LIBS)
install: fmute marchenko
diff --git a/marchenko_full/marchenko b/marchenko_full/marchenko
index e0daa07327fde0a9fae278ad67530b1f32d146d0..53bb50f04f67c4ed6e80dec058a734f62a9716df 100755
Binary files a/marchenko_full/marchenko and b/marchenko_full/marchenko differ
diff --git a/marchenko_full/marchenko.c b/marchenko_full/marchenko.c
index 18b81cf19815f8c8fac8fb99e6b0ca63cda81249..4e7ed45dba8cfafaae81e46d451c162c54695731 100644
--- a/marchenko_full/marchenko.c
+++ b/marchenko_full/marchenko.c
@@ -7,6 +7,7 @@
#include <assert.h>
#include <genfft.h>
#include "marchenko.h"
+#include "raytime.h"
int omp_get_max_threads(void);
int omp_get_num_threads(void);
@@ -119,6 +120,80 @@ char *sdoc[] = {
" file_iter= ............... output file with -Ni(-t) for each iteration",
" verbose=0 ................ silent option; >0 displays info",
" ",
+" RAYTIME PARAMETERS - Jesper Spetzler ray-trace modeling ",
+" ",
+" IO PARAMETERS:",
+" file_cp= .......... P (cp) velocity file",
+" file_src= ......... file with source signature",
+" file_rcv=recv.su .. base name for receiver files",
+" dx= ............... read from model file: if dx==0 then dx= can be used to set it",
+" dz= ............... read from model file: if dz==0 then dz= can be used to set it",
+" dt= ............... read from file_src: if dt==0 then dt= can be used to set it",
+"" ,
+" RAY TRACING PARAMETERS:",
+" smoothwindow=0 .... if set lenght of 2/3D smoothing window on slowness",
+" useT2=0 ........... 1: compute more accurate T2 pertubation correction",
+" geomspread=1 ...... 1: compute Geometrical Spreading Factor",
+" nraystep=5 ........ number of points on ray",
+" OPTIONAL PARAMETERS:",
+" ischeme=3 ......... 1=acoustic, 2=visco-acoustic 3=elastic, 4=visco-elastic",
+" sinkdepth=0 ....... receiver grid points below topography (defined bij cp=0.0)",
+" sinkdepth_src=0 ... source grid points below topography (defined bij cp=0.0)",
+" sinkvel=0 ......... use velocity of first receiver to sink through to next layer",
+" verbose=0 ......... silent mode; =1: display info",
+" ",
+" SHOT AND GENERAL SOURCE DEFINITION:",
+" xsrc=middle ....... x-position of (first) shot ",
+" zsrc=zmin ......... z-position of (first) shot ",
+" nshot=1 ........... number of shots to model",
+" dxshot=dx ......... if nshot > 1: x-shift in shot locations",
+" dzshot=0 .......... if nshot > 1: z-shift in shot locations",
+" xsrca= ............ defines source array x-positions",
+" zsrca= ............ defines source array z-positions",
+" wav_random=1 ...... 1 generates (band limited by fmax) noise signatures ",
+" src_multiwav=0 .... use traces in file_src as areal source",
+" src_at_rcv=1 ...... inject wavefield at receiver coordinates (1), inject at source (0)",
+"" ,
+" PLANE WAVE SOURCE DEFINITION:",
+" plane_wave=0 ...... model plane wave with nsrc= sources",
+" nsrc=1 ............ number of sources per (plane-wave) shot ",
+" src_angle=0 ....... angle of plane source array",
+" src_velo=1500 ..... velocity to use in src_angle definition",
+" src_window=0 ...... length of taper at edges of source array",
+"",
+" RANDOM SOURCE DEFINITION FOR SEISMIC INTERFEROMTERY:",
+" src_random=0 ...... 1 enables nsrc random sources positions in one modeling",
+" nsrc=1 ............ number of sources to use for one shot",
+" xsrc1=0 ........... left bound for x-position of sources",
+" xsrc2=0 ........... right bound for x-position of sources",
+" zsrc1=0 ........... left bound for z-position of sources",
+" zsrc2=0 ........... right bound for z-position of sources",
+" tsrc1=0.0 ......... begin time interval for random sources being triggered",
+" tsrc2=tmod ........ end time interval for random sources being triggered",
+" tactive=tsrc2 ..... end time for random sources being active",
+" tlength=tsrc2-tsrc1 average duration of random source signal",
+" length_random=1 ... duration of source is rand*tlength",
+" amplitude=0 ....... distribution of source amplitudes",
+" distribution=0 .... random function for amplitude and tlength 0=flat 1=Gaussian ",
+" seed=10 ........... seed for start of random sequence ",
+"" ,
+" RECEIVER SELECTION:",
+" xrcv1=xmin ........ first x-position of linear receiver array(s)",
+" xrcv2=xmax ........ last x-position of linear receiver array(s)",
+" dxrcv=dx .......... x-position increment of receivers in linear array(s)",
+" zrcv1=zmin ........ first z-position of linear receiver array(s)",
+" zrcv2=zrcv1 ....... last z-position of linear receiver array(s)",
+" dzrcv=0.0 ......... z-position increment of receivers in linear array(s)",
+" xrcva= ............ defines receiver array x-positions",
+" zrcva= ............ defines receiver array z-positions",
+" rrcv= ............. radius for receivers on a circle ",
+" arcv= ............. vertical arc-lenght for receivers on a ellipse (rrcv=horizontal)",
+" oxrcv=0.0 ......... x-center position of circle",
+" ozrcv=0.0 ......... z-center position of circle",
+" dphi=2 ............ angle between receivers on circle ",
+" rcv_txt=........... text file with receiver coordinates. Col 1: x, Col. 2: z",
+" rec_ntsam=nt ...... maximum number of time samples in file_rcv files",
+" ",
" ",
" author : Jan Thorbecke : 2016 (j.w.thorbecke@tudelft.nl)",
" ",
@@ -146,6 +221,14 @@ int main (int argc, char **argv)
char *file_f1plus, *file_f1min, *file_gmin, *file_gplus, *file_f2, *file_pmin, *wavtype;
segy *hdrs_out;
WavePar WP;
+ modPar mod;
+ recPar rec;
+ snaPar sna;
+ wavPar wav;
+ srcPar src;
+ bndPar bnd;
+ shotPar shot;
+ rayPar ray;
initargs(argc, argv);
requestdoc(1);
@@ -223,7 +306,34 @@ int main (int argc, char **argv)
/*================ Reading info about shot and initial operator sizes ================*/
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);
+ if (file_ray!=NULL && file_tinv==NULL) {
+ vmess("test ray");
+ ret = getFileInfo(file_ray, &n2, &n1, &ngath, &d2, &d1, &f2, &f1, &xmin, &xmax, &scl, &ntraces);
+ }
+ else if (file_ray==NULL && file_tinv==NULL) {
+ vmess("test raytime");
+ getParameters(&mod, &rec, &sna, &wav, &src, &shot, &bnd, &ray, verbose);
+ n1 = rec.nt;
+ n2 = rec.n;
+ ngath = shot.n;
+ d1 = mod.dt;
+ d2 = (rec.x[1]-rec.x[0])*mod.dx;
+ f1 = 0.0;
+ f2 = mod.x0+rec.x[0]*mod.dx;
+ xmin = mod.x0+rec.x[0]*mod.dx;
+ xmax = mod.x0+rec.x[rec.n-1]*mod.dx;
+ scl = 0.0010;
+ ntraces = n2*ngath;
+ }
+ else {
+ vmess("test tinv");
+ ret = getFileInfo(file_tinv, &n1, &n2, &ngath, &d1, &d2, &f1, &f2, &xmin, &xmax, &scl, &ntraces);
+ }
+
+ vmess("n1:%d n2:%d ngath:%d d1:%.4f d2:%.4f f1:%.4f f2:%.4f",n1,n2,ngath,d1,d2,f1,f2);
+ vmess("xmin:%.4f xmax:%.4f scl:%.4f ntraces:%d",xmin,xmax,scl,ntraces);
+
+
Nsyn = ngath;
nxs = n2;
nts = n1;
diff --git a/marchenko_full/raytime.c b/marchenko_full/raytime.c
index d4a8c839a511a30647c289c22c3fbd05c845979e..db62897b9990791a4a521216ea4eb5afa3489fcd 100644
--- a/marchenko_full/raytime.c
+++ b/marchenko_full/raytime.c
@@ -39,91 +39,7 @@ int defineSource(wavPar wav, srcPar src, modPar mod, float **src_nwav, int rever
int writeSrcRecPos(modPar *mod, recPar *rec, srcPar *src, shotPar *shot);
-/* Self documentation */
-char *sdoc[] = {
-" ",
-" raytime - Jesper Spetzler ray-trace modeling ",
-" ",
-" IO PARAMETERS:",
-" file_cp= .......... P (cp) velocity file",
-" file_src= ......... file with source signature",
-" file_rcv=recv.su .. base name for receiver files",
-" dx= ............... read from model file: if dx==0 then dx= can be used to set it",
-" dz= ............... read from model file: if dz==0 then dz= can be used to set it",
-" dt= ............... read from file_src: if dt==0 then dt= can be used to set it",
-"" ,
-" RAY TRACING PARAMETERS:",
-" smoothwindow=0 .... if set lenght of 2/3D smoothing window on slowness",
-" useT2=0 ........... 1: compute more accurate T2 pertubation correction",
-" geomspread=1 ...... 1: compute Geometrical Spreading Factor",
-" nraystep=5 ........ number of points on ray",
-" OPTIONAL PARAMETERS:",
-" ischeme=3 ......... 1=acoustic, 2=visco-acoustic 3=elastic, 4=visco-elastic",
-" sinkdepth=0 ....... receiver grid points below topography (defined bij cp=0.0)",
-" sinkdepth_src=0 ... source grid points below topography (defined bij cp=0.0)",
-" sinkvel=0 ......... use velocity of first receiver to sink through to next layer",
-" verbose=0 ......... silent mode; =1: display info",
-" ",
-" SHOT AND GENERAL SOURCE DEFINITION:",
-" xsrc=middle ....... x-position of (first) shot ",
-" zsrc=zmin ......... z-position of (first) shot ",
-" nshot=1 ........... number of shots to model",
-" dxshot=dx ......... if nshot > 1: x-shift in shot locations",
-" dzshot=0 .......... if nshot > 1: z-shift in shot locations",
-" xsrca= ............ defines source array x-positions",
-" zsrca= ............ defines source array z-positions",
-" wav_random=1 ...... 1 generates (band limited by fmax) noise signatures ",
-" src_multiwav=0 .... use traces in file_src as areal source",
-" src_at_rcv=1 ...... inject wavefield at receiver coordinates (1), inject at source (0)",
-"" ,
-" PLANE WAVE SOURCE DEFINITION:",
-" plane_wave=0 ...... model plane wave with nsrc= sources",
-" nsrc=1 ............ number of sources per (plane-wave) shot ",
-" src_angle=0 ....... angle of plane source array",
-" src_velo=1500 ..... velocity to use in src_angle definition",
-" src_window=0 ...... length of taper at edges of source array",
-"",
-" RANDOM SOURCE DEFINITION FOR SEISMIC INTERFEROMTERY:",
-" src_random=0 ...... 1 enables nsrc random sources positions in one modeling",
-" nsrc=1 ............ number of sources to use for one shot",
-" xsrc1=0 ........... left bound for x-position of sources",
-" xsrc2=0 ........... right bound for x-position of sources",
-" zsrc1=0 ........... left bound for z-position of sources",
-" zsrc2=0 ........... right bound for z-position of sources",
-" tsrc1=0.0 ......... begin time interval for random sources being triggered",
-" tsrc2=tmod ........ end time interval for random sources being triggered",
-" tactive=tsrc2 ..... end time for random sources being active",
-" tlength=tsrc2-tsrc1 average duration of random source signal",
-" length_random=1 ... duration of source is rand*tlength",
-" amplitude=0 ....... distribution of source amplitudes",
-" distribution=0 .... random function for amplitude and tlength 0=flat 1=Gaussian ",
-" seed=10 ........... seed for start of random sequence ",
-"" ,
-" RECEIVER SELECTION:",
-" xrcv1=xmin ........ first x-position of linear receiver array(s)",
-" xrcv2=xmax ........ last x-position of linear receiver array(s)",
-" dxrcv=dx .......... x-position increment of receivers in linear array(s)",
-" zrcv1=zmin ........ first z-position of linear receiver array(s)",
-" zrcv2=zrcv1 ....... last z-position of linear receiver array(s)",
-" dzrcv=0.0 ......... z-position increment of receivers in linear array(s)",
-" xrcva= ............ defines receiver array x-positions",
-" zrcva= ............ defines receiver array z-positions",
-" rrcv= ............. radius for receivers on a circle ",
-" arcv= ............. vertical arc-lenght for receivers on a ellipse (rrcv=horizontal)",
-" oxrcv=0.0 ......... x-center position of circle",
-" ozrcv=0.0 ......... z-center position of circle",
-" dphi=2 ............ angle between receivers on circle ",
-" rcv_txt=........... text file with receiver coordinates. Col 1: x, Col. 2: z",
-" rec_ntsam=nt ...... maximum number of time samples in file_rcv files",
-"",
-" Jan Thorbecke 2017",
-" TU Delft",
-" E-mail: janth@xs4all.nl ",
-"",
-NULL};
-
-
-int main(int argc, char **argv)
+int raytime(float *time, float *ampl, int *xnx, float *xrcv, float *xsrc, float *zsrc)
{
modPar mod;
recPar rec;
@@ -155,7 +71,7 @@ int main(int argc, char **argv)
int nRayStep;
fcoord coordsx, coordgx, Time;
icoord grid;
- float Jr, *ampl, *time;
+ float Jr;
segy hdr;
char filetime[1024], fileamp[1024];
size_t nwrite;
@@ -164,10 +80,10 @@ int main(int argc, char **argv)
double ddt;
t0= wallclock_time();
- initargs(argc,argv);
- requestdoc(0);
+ //initargs(argc,argv);
+ //requestdoc(0);
- if (!getparint("verbose",&verbose)) verbose=0;
+ //if (!getparint("verbose",&verbose)) verbose=0;
getParameters(&mod, &rec, &sna, &wav, &src, &shot, &bnd, &ray, verbose);
/* allocate arrays for model parameters: the different schemes use different arrays */
@@ -189,8 +105,8 @@ int main(int argc, char **argv)
/* allocate arrays for wavefield and receiver arrays */
size = shot.n*rec.n;
- time = (float *)calloc(size,sizeof(float));
- ampl = (float *)calloc(size,sizeof(float));
+ //time = (float *)calloc(size,sizeof(float));
+ //ampl = (float *)calloc(size,sizeof(float));
t1= wallclock_time();
if (verbose) {
@@ -233,7 +149,7 @@ int main(int argc, char **argv)
if (verbose>3) writeSrcRecPos(&mod, &rec, &src, &shot);
/* prepare output file and headers */
- strcpy(filetime, rec.file_rcv);
+ /*strcpy(filetime, rec.file_rcv);
name_ext(filetime, "_time");
fpt = fopen(filetime, "w");
assert(fpt != NULL);
@@ -243,15 +159,15 @@ int main(int argc, char **argv)
name_ext(fileamp, "_amp");
fpa = fopen(fileamp, "w");
assert(fpa != NULL);
- }
+ }*/
- ddt = (double)mod.dt;/* to avoid rounding in 32 bit precision */
+ /*ddt = (double)mod.dt;
hdr.dt = (unsigned short)lround((((double)1.0e6*ddt*rec.skipdt)));
hdr.scalco = -1000;
hdr.scalel = -1000;
hdr.trid = 1;
hdr.trwf = shot.n;
- hdr.ns = rec.n;
+ hdr.ns = rec.n;*/
/* Outer loop over number of shots */
for (ishot=0; ishot<shot.n; ishot++) {
@@ -272,6 +188,10 @@ int main(int argc, char **argv)
grid.z = mod.nz;
grid.y = 1;
+ xnx[ishot]=rec.n;
+ xsrc[ishot] = mod.x0+mod.dx*shot.x[ishot];
+ zsrc[ishot] = mod.z0+mod.dz*shot.z[ishot];
+
for (irec=0; irec<rec.n; irec++) {
coordgx.x=mod.x0+rec.xr[irec];
coordgx.z=mod.z0+rec.zr[irec];
@@ -279,12 +199,13 @@ int main(int argc, char **argv)
getWaveParameter(slowness, grid, mod.dx, coordsx, coordgx, ray, &Time, &Jr);
+ xrcv[ishot*rec.n+irec] = (mod.x0+rec.x[0]*mod.dx) + ((rec.x[1]-rec.x[0])*mod.dx*((float)irec));
time[ishot*rec.n + irec] = Time.x + Time.y + Time.z;
ampl[ishot*rec.n + irec] = Jr;
fprintf(stderr,"shot=%f,%f receiver at %f,%f T0=%f T1=%f T2=%f Jr=%f\n",coordsx.x, coordsx.z, coordgx.x, coordgx.z, Time.x, Time.y, Time.z, Jr);
}
- hdr.sx = 1000*(mod.x0+mod.dx*shot.x[ishot]);
+ /*hdr.sx = 1000*(mod.x0+mod.dx*shot.x[ishot]);
hdr.sdepth = 1000*(mod.z0+mod.dz*shot.z[ishot]);
hdr.selev = (int)(-1000.0*(mod.z0+mod.dz*shot.z[ishot]));
hdr.fldr = ishot+1;
@@ -307,11 +228,11 @@ int main(int argc, char **argv)
nwrite = fwrite( &l[ishot*rec.n], sizeof(float), rec.n, fpa);
assert(nwrite == rec.n);
fflush(fpa);
- }
+ }*/
} /* end of loop over number of shots */
- fclose(fpt);
- if (ray.geomspread) fclose(fpa);
+ //fclose(fpt);
+ //if (ray.geomspread) fclose(fpa);
t1= wallclock_time();
if (verbose) {
@@ -320,7 +241,7 @@ int main(int argc, char **argv)
/* free arrays */
- initargs(argc,argv); /* this will free the arg arrays declared */
+ //initargs(argc,argv); /* this will free the arg arrays declared */
free(velocity);
free(slowness);
diff --git a/marchenko_full/readTinvData.c b/marchenko_full/readTinvData.c
index 755f35ac91291ba8e12a626a75843fbed75949a0..bceab866f43e80bc7dc1550471dee92a8f1a83ca 100644
--- a/marchenko_full/readTinvData.c
+++ b/marchenko_full/readTinvData.c
@@ -23,21 +23,27 @@ typedef struct _complexStruct { /* complex number */
void findShotInMute(float *xrcvMute, float xrcvShot, int nxs, int *imute);
int readSnapData(char *filename, float *data, segy *hdrs, int nsnaps, int nx, int nz, int sx, int ex, int sz, int ez);
+int raytime(float *amp, float *time, int *xnx, float *xrcv, float *xsrc, float *zsrc);
int readTinvData(char *filename, WavePar WP, char *file_ray, char *file_amp, float dt, float *xrcv, float *xsrc, float *zsrc, int *xnx, int Nsyn, int nx, int ntfft, int mode, int *maxval, float *tinv, int hw, int verbose)
{
FILE *fp;
segy hdr, *hdrs_mute, *hdrs_amp, *hdrs_wav;
size_t nread;
+ char *file_cp;
int fldr_shot, sx_shot, itrace, one_shot, ig, isyn, i, j;
- int end_of_file, nt, gx0, gx1, nfreq;
+ int end_of_file, nt, gx0, gx1, nfreq, geosp;
int nx1, jmax, imax, tstart, tend, nwav;
float xmax, tmax, lmax, *wavelet, *wavelet2;
float scl, scel, *trace, dxrcv, *timeval, dw, *amp;
complex *cmute, *cwav;
+ if (!getparstring("file_cp", &file_cp)) file_cp=NULL;
+ if (!getparint("geomspread",&geosp)) geosp=1;
+ if (file_cp==NULL) geosp=0;
+
/*Check wheter the raytime is used or not*/
- if (file_ray!=NULL) {
+ if (file_ray!=NULL || file_cp!=NULL) {
/*Define parameters*/
nfreq = ntfft/2+1;
wavelet = (float *)calloc(ntfft,sizeof(float));
@@ -69,40 +75,48 @@ int readTinvData(char *filename, WavePar WP, char *file_ray, char *file_amp, flo
}
}
rc1fft(wavelet,cwav,ntfft,-1);
+ free(wavelet);
}
- /* Defining mute window using raytimes */
- vmess("Using raytime for mutewindow");
- hdrs_mute = (segy *) calloc(Nsyn,sizeof(segy));
- timeval = (float *)calloc(Nsyn*nx,sizeof(float));
- fp = fopen( file_ray, "r" );
- if ( fp == NULL ) {
- perror("Error opening file containing wavelet");
- }
- fclose(fp);
- readSnapData(file_ray, timeval, hdrs_mute, Nsyn, 1, nx, 0, 1, 0, nx);
-
- /*Check whether the amplitude is also used*/
- if (file_amp != NULL) {
- hdrs_amp = (segy *) calloc(Nsyn,sizeof(segy));
- amp = (float *)calloc(Nsyn*nx,sizeof(float));
- fp = fopen( file_amp, "r" );
- if ( fp == NULL ) {
+ timeval = (float *)calloc(Nsyn*nx,sizeof(float));
+ amp = (float *)calloc(Nsyn*nx,sizeof(float));
+
+ if (file_ray!=NULL) {
+
+ /* Defining mute window using raytimes */
+ vmess("Using raytime for mutewindow");
+ hdrs_mute = (segy *) calloc(Nsyn,sizeof(segy));
+ fp = fopen( file_ray, "r" );
+ if ( fp == NULL ) {
perror("Error opening file containing wavelet");
- }
- fclose(fp);
- readSnapData(file_amp, amp, hdrs_amp, Nsyn, 1, nx, 0, 1, 0, nx);
- }
+ }
+ fclose(fp);
+ readSnapData(file_ray, timeval, hdrs_mute, Nsyn, 1, nx, 0, 1, 0, nx);
+
+ /*Check whether the amplitude is also used*/
+ if (file_amp != NULL) {
+ hdrs_amp = (segy *) calloc(Nsyn,sizeof(segy));
+ fp = fopen( file_amp, "r" );
+ if ( fp == NULL ) {
+ perror("Error opening file containing wavelet");
+ }
+ fclose(fp);
+ readSnapData(file_amp, amp, hdrs_amp, Nsyn, 1, nx, 0, 1, 0, nx);
+ }
- /*Define source and receiver locations from the raytime*/
- for (isyn=0; isyn<Nsyn; isyn++) {
- for (itrace=0; itrace<nx; itrace++) {
- xrcv[isyn*nx+itrace] = (hdrs_mute[isyn].f1 + hdrs_mute[isyn].d1*((float)itrace));
- }
- xnx[isyn]=nx;
- xsrc[isyn] = hdrs_mute[isyn].sx;
- zsrc[isyn] = hdrs_mute[isyn].sdepth;
- }
+ /*Define source and receiver locations from the raytime*/
+ for (isyn=0; isyn<Nsyn; isyn++) {
+ for (itrace=0; itrace<nx; itrace++) {
+ xrcv[isyn*nx+itrace] = (hdrs_mute[isyn].f1 + hdrs_mute[isyn].d1*((float)itrace));
+ }
+ xnx[isyn]=nx;
+ xsrc[isyn] = hdrs_mute[isyn].sx;
+ zsrc[isyn] = hdrs_mute[isyn].sdepth;
+ }
+ }
+ else {
+ raytime(timeval,amp,xnx,xrcv,xsrc,zsrc);
+ }
/*Determine the mutewindow*/
for (j=0; j<Nsyn; j++) {
@@ -110,7 +124,7 @@ int readTinvData(char *filename, WavePar WP, char *file_ray, char *file_amp, flo
maxval[j*nx+i] = (int)roundf(timeval[j*nx+i]/dt);
if (maxval[j*nx+i] > ntfft-1) maxval[j*nx+i] = ntfft-1;
if (WP.wav) { /*Apply the wavelet to create a first arrival*/
- if (file_amp != NULL) {
+ if (file_amp != NULL || geosp==1) {
for (ig=0; ig<nfreq; ig++) {
cmute[ig].r = (cwav[ig].r*cos(ig*dw*timeval[j*nx+i]-M_PI/4.0)-cwav[ig].i*sin(ig*dw*timeval[j*nx+i]-M_PI/4.0))/amp[j*nx+i];
cmute[ig].i = (cwav[ig].i*cos(ig*dw*timeval[j*nx+i]-M_PI/4.0)+cwav[ig].r*sin(ig*dw*timeval[j*nx+i]-M_PI/4.0))/amp[j*nx+i];
@@ -129,7 +143,7 @@ int readTinvData(char *filename, WavePar WP, char *file_ray, char *file_amp, flo
}
}
- if (!WP.wav) {
+ if (WP.wav == 0 && file_ray==NULL && filename!=NULL) {
/* Reading first header */