diff --git a/marchenko_full/marchenko b/marchenko_full/marchenko
index 53bb50f04f67c4ed6e80dec058a734f62a9716df..6266b0c58a7930cb6bf1cab9d8c7c11aef5615c3 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 4e7ed45dba8cfafaae81e46d451c162c54695731..71b0f76d2afdbadfb0b0cfda11b9705443692829 100644
--- a/marchenko_full/marchenko.c
+++ b/marchenko_full/marchenko.c
@@ -307,11 +307,9 @@ int main (int argc, char **argv)
ngath = 0; /* setting ngath=0 scans all traces; n2 contains maximum traces/gather */
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;
@@ -326,14 +324,9 @@ int main (int argc, char **argv)
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/raytime/JespersRayTracer.c b/raytime/JespersRayTracer.c
index e757defd634746ce353a2b28f4d8ad50a55bdaa9..67e590a6f979f10e416a9484fb331b214e3c2a02 100644
--- a/raytime/JespersRayTracer.c
+++ b/raytime/JespersRayTracer.c
@@ -110,7 +110,7 @@ int getWaveParameter(float *slowness, icoord size, float dgrid, fcoord s, fcoord
ds = lengthRefRay/(nRayTmp-1);
J = lengthRefRay;
dQdPhi = 0;
-
+
for (i = 0; i < nRayTmp-1; i++)
{
x = 0.5*(rayReference3D[i+1].x + rayReference3D[i].x);
diff --git a/raytime/getParameters.c b/raytime/getParameters.c
index 8f82e3d264bca9e9007c8258088284b53c2d1d2e..265b823d1c1e37d2e9fd9c0252645ff545b6a6f6 100644
--- a/raytime/getParameters.c
+++ b/raytime/getParameters.c
@@ -123,40 +123,50 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
if (!getparfloat("zsrc",&zsrc)) zsrc=sub_z0;
// if (!getparint("nsrc",&nsrc)) nsrc=1;
- if (!getparint("nshot",&shot->n)) shot->n=1;
+ //if (!getparint("nshot",&shot->n)) shot->n=1;
+ if (!getparint("nxshot",&shot->nx)) shot->nx=1;
+ if (!getparint("nzshot",&shot->nz)) shot->nz=1;
if (!getparfloat("dxshot",&dxshot)) dxshot=dx;
- if (!getparfloat("dzshot",&dzshot)) dzshot=0.0;
+ if (!getparfloat("dzshot",&dzshot)) dzshot=dz;
- if (shot->n>1) {
+ shot->n = (shot->nx)*(shot->nz);
+
+ if (shot->nx>1) {
idxshot=MAX(0,NINT(dxshot/dx));
- idzshot=MAX(0,NINT(dzshot/dz));
}
else {
idxshot=0.0;
- idzshot=0.0;
}
-
+ if (shot->nz>1) {
+ idzshot=MAX(0,NINT(dzshot/dz));
+ }
+ else {
+ idzshot=0.0;
+ }
+
/* calculate the shot positions */
src_ix0=MAX(0,NINT((xsrc-sub_x0)/dx));
src_ix0=MIN(src_ix0,nx);
src_iz0=MAX(0,NINT((zsrc-sub_z0)/dz));
src_iz0=MIN(src_iz0,nz);
- srcendx=(shot->n-1)*dxshot+xsrc;
- srcendz=(shot->n-1)*dzshot+zsrc;
+ srcendx=(shot->nx-1)*dxshot+xsrc;
+ srcendz=(shot->nz-1)*dzshot+zsrc;
src_ix1=MAX(0,NINT((srcendx-sub_x0)/dx));
src_ix1=MIN(src_ix1,nx);
src_iz1=MAX(0,NINT((srcendz-sub_z0)/dz));
src_iz1=MIN(src_iz1,nz);
- shot->x = (int *)calloc(shot->n,sizeof(int));
- shot->z = (int *)calloc(shot->n,sizeof(int));
- for (is=0; is<shot->n; is++) {
+ shot->x = (int *)calloc(shot->nx,sizeof(int));
+ shot->z = (int *)calloc(shot->nz,sizeof(int));
+ for (is=0; is<shot->nx; is++) {
shot->x[is] = src_ix0+is*idxshot;
- shot->z[is] = src_iz0+is*idzshot;
- if (shot->x[is] > nx-1) shot->n = is-1;
- if (shot->z[is] > nz-1) shot->n = is-1;
+ if (shot->x[is] > nx-1) shot->nx = is-1;
}
+ for (is=0; is<shot->nz; is++) {
+ shot->z[is] = src_iz0+is*idzshot;
+ if (shot->z[is] > nz-1) shot->nz = is-1;
+ }
/* check if source array is defined */
diff --git a/raytime/model.scr b/raytime/model.scr
index ad760e5ea328923107f6a18005aa79231ee3dcff..6cf97b84ae5d137521901c2fb2de5238e257d474 100755
--- a/raytime/model.scr
+++ b/raytime/model.scr
@@ -15,12 +15,38 @@ export OMP_NUM_THREADS=2
./raytime \
file_cp=syncl_cp.su \
- file_rcv=shot_fd.su \
+ file_rcv=z1.su \
smoothwindow=15 \
- verbose=2 \
+ verbose=3 \
dxrcv=10.0 \
xrcv1=-2500 xrcv2=2500 \
zrcv1=0 zrcv2=0 \
- xsrc=0 zsrc=1100
+ xsrc=0 zsrc=1100 \
+ nxshot=1 dxshot=10 \
+ nzshot=1 dzshot=10
+
+./raytime \
+ file_cp=syncl_cp.su \
+ file_rcv=z2.su \
+ smoothwindow=15 \
+ verbose=3 \
+ dxrcv=10.0 \
+ xrcv1=-2500 xrcv2=2500 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=1200 \
+ nxshot=1 dxshot=10 \
+ nzshot=1 dzshot=10
+
+./raytime \
+ file_cp=syncl_cp.su \
+ file_rcv=z0.su \
+ smoothwindow=15 \
+ verbose=3 \
+ dxrcv=10.0 \
+ xrcv1=-2500 xrcv2=2500 \
+ zrcv1=0 zrcv2=0 \
+ xsrc=0 zsrc=1100 \
+ nxshot=1 dxshot=10 \
+ nzshot=2 dzshot=100
echo "testJan"
diff --git a/raytime/raytime.c b/raytime/raytime.c
index d4a8c839a511a30647c289c22c3fbd05c845979e..be7fb12b2326875bdbd8b716c368b552412f3d82 100644
--- a/raytime/raytime.c
+++ b/raytime/raytime.c
@@ -148,7 +148,7 @@ int main(int argc, char **argv)
float sinkvel;
double t0, t1, t2, t3, tt, tinit;
size_t size, sizem, nsamp;
- int n1, ix, iz, ir, ishot, i;
+ int n1, ix, iz, ir, ixshot, izshot, i;
int ioPx, ioPz;
int it0, it1, its, it, fileno, isam;
int ixsrc, izsrc, irec;
@@ -223,11 +223,13 @@ int main(int argc, char **argv)
*/
/* sink sources to value different than zero */
- for (ishot=0; ishot<shot.n; ishot++) {
- iz = shot.z[ishot];
- ix = shot.x[ishot];
- while(velocity[(ix)*n1+iz] == 0.0) iz++;
- shot.z[ishot]=iz+src.sinkdepth;
+ for (izshot=0; izshot<shot.nz; izshot++) {
+ for (ixshot=0; ixshot<shot.nx; ixshot++) {
+ iz = shot.z[izshot];
+ ix = shot.x[ixshot];
+ while(velocity[(ix)*n1+iz] == 0.0) iz++;
+ shot.z[izshot]=iz+src.sinkdepth;
+ }
}
if (verbose>3) writeSrcRecPos(&mod, &rec, &src, &shot);
@@ -254,61 +256,62 @@ int main(int argc, char **argv)
hdr.ns = rec.n;
/* Outer loop over number of shots */
- for (ishot=0; ishot<shot.n; ishot++) {
-
- if (verbose) {
- vmess("Modeling source %d at gridpoints ix=%d iz=%d", ishot, shot.x[ishot], shot.z[ishot]);
- vmess(" which are actual positions x=%.2f z=%.2f", mod.x0+mod.dx*shot.x[ishot], mod.z0+mod.dz*shot.z[ishot]);
- vmess("Receivers at gridpoint x-range ix=%d - %d", rec.x[0], rec.x[rec.n-1]);
- vmess(" which are actual positions x=%.2f - %.2f", mod.x0+rec.xr[0], mod.x0+rec.xr[rec.n-1]);
- vmess("Receivers at gridpoint z-range iz=%d - %d", rec.z[0], rec.z[rec.n-1]);
- vmess(" which are actual positions z=%.2f - %.2f", mod.z0+rec.zr[0], mod.z0+rec.zr[rec.n-1]);
- }
-
- coordsx.x = mod.x0+shot.x[ishot]*mod.dx;
- coordsx.z = mod.z0+shot.z[ishot]*mod.dz;
- coordsx.y = 0;
- grid.x = mod.nx;
- grid.z = mod.nz;
- grid.y = 1;
-
- for (irec=0; irec<rec.n; irec++) {
- coordgx.x=mod.x0+rec.xr[irec];
- coordgx.z=mod.z0+rec.zr[irec];
- coordgx.y = 0;
-
- getWaveParameter(slowness, grid, mod.dx, coordsx, coordgx, ray, &Time, &Jr);
-
- 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.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;
- hdr.tracl = ishot+1;
- hdr.tracf = ishot+1;
- hdr.ntr = shot.n;
- hdr.d1 = (rec.x[1]-rec.x[0])*mod.dx;
- hdr.f1 = mod.x0+rec.x[0]*mod.dx;
- hdr.d2 = (shot.x[1]-shot.x[0])*mod.dx;
- hdr.f2 = mod.x0+shot.x[0]*mod.dx;
-
- nwrite = fwrite( &hdr, 1, TRCBYTES, fpt);
- assert(nwrite == TRCBYTES);
- nwrite = fwrite( &time[ishot*rec.n], sizeof(float), rec.n, fpt);
- assert(nwrite == rec.n);
- fflush(fpt);
- if (ray.geomspread) {
- nwrite = fwrite( &hdr, 1, TRCBYTES, fpa);
- assert(nwrite == TRCBYTES);
- nwrite = fwrite( &l[ishot*rec.n], sizeof(float), rec.n, fpa);
- assert(nwrite == rec.n);
- fflush(fpa);
- }
+ for (izshot=0; izshot<shot.nz; izshot++) {
+ for (ixshot=0; ixshot<shot.nx; ixshot++) {
+
+ if (verbose) {
+ vmess("Modeling source %d at gridpoints ix=%d iz=%d", (izshot*shot.n)+ixshot, shot.x[ixshot], shot.z[izshot]);
+ vmess(" which are actual positions x=%.2f z=%.2f", mod.x0+mod.dx*shot.x[ixshot], mod.z0+mod.dz*shot.z[izshot]);
+ vmess("Receivers at gridpoint x-range ix=%d - %d", rec.x[0], rec.x[rec.n-1]);
+ vmess(" which are actual positions x=%.2f - %.2f", mod.x0+rec.xr[0], mod.x0+rec.xr[rec.n-1]);
+ vmess("Receivers at gridpoint z-range iz=%d - %d", rec.z[0], rec.z[rec.n-1]);
+ vmess(" which are actual positions z=%.2f - %.2f", mod.z0+rec.zr[0], mod.z0+rec.zr[rec.n-1]);
+ }
+
+ coordsx.x = mod.x0+shot.x[ixshot]*mod.dx;
+ coordsx.z = mod.z0+shot.z[izshot]*mod.dz;
+ coordsx.y = 0;
+ grid.x = mod.nx;
+ grid.z = mod.nz;
+ grid.y = 1;
+
+ for (irec=0; irec<rec.n; irec++) {
+ coordgx.x=mod.x0+rec.xr[irec];
+ coordgx.z=mod.z0+rec.zr[irec];
+ coordgx.y = 0;
+
+ getWaveParameter(slowness, grid, mod.dx, coordsx, coordgx, ray, &Time, &Jr);
+
+ time[((izshot*shot.nx)+ixshot)*rec.n + irec] = Time.x + Time.y + Time.z;
+ ampl[((izshot*shot.nx)+ixshot)*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[ixshot]);
+ hdr.sdepth = 1000*(mod.z0+mod.dz*shot.z[izshot]);
+ hdr.selev = (int)(-1000.0*(mod.z0+mod.dz*shot.z[izshot]));
+ hdr.fldr = ((izshot*shot.nx)+ixshot)+1;
+ hdr.tracl = ((izshot*shot.nx)+ixshot)+1;
+ hdr.tracf = ((izshot*shot.nx)+ixshot)+1;
+ hdr.ntr = shot.n;
+ hdr.d1 = (rec.x[1]-rec.x[0])*mod.dx;
+ hdr.f1 = mod.x0+rec.x[0]*mod.dx;
+ hdr.d2 = (shot.x[1]-shot.x[0])*mod.dx;
+ hdr.f2 = mod.x0+shot.x[0]*mod.dx;
+ nwrite = fwrite( &hdr, 1, TRCBYTES, fpt);
+ assert(nwrite == TRCBYTES);
+ nwrite = fwrite( &time[((izshot*shot.nx)+ixshot)*rec.n], sizeof(float), rec.n, fpt);
+ assert(nwrite == rec.n);
+ fflush(fpt);
+ if (ray.geomspread) {
+ nwrite = fwrite( &hdr, 1, TRCBYTES, fpa);
+ assert(nwrite == TRCBYTES);
+ nwrite = fwrite( &l[((izshot*shot.nx)+ixshot)*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);
diff --git a/raytime/raytime.h b/raytime/raytime.h
index dea8b1be838096a00a27998b005eea0983230a5f..464553958272a713004412b4cc1afdd8b7fd4d9e 100644
--- a/raytime/raytime.h
+++ b/raytime/raytime.h
@@ -152,6 +152,8 @@ typedef struct _sourcePar { /* Source Array Parameters */
typedef struct _shotPar { /* Shot Parameters */
int n;
+ int nx;
+ int nz;
int *z;
int *x;
} shotPar;