diff --git a/marchenko_full/getParameters.c b/marchenko_full/getParameters.c
index 8f82e3d264bca9e9007c8258088284b53c2d1d2e..265b823d1c1e37d2e9fd9c0252645ff545b6a6f6 100644
--- a/marchenko_full/getParameters.c
+++ b/marchenko_full/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/marchenko_full/marchenko b/marchenko_full/marchenko
index 6266b0c58a7930cb6bf1cab9d8c7c11aef5615c3..44eb6bcfd0a0210e25b80db9c080bce927619db1 100755
Binary files a/marchenko_full/marchenko and b/marchenko_full/marchenko differ
diff --git a/marchenko_full/raytime.c b/marchenko_full/raytime.c
index db62897b9990791a4a521216ea4eb5afa3489fcd..21932168e1e107830560bca868be2e4404a4c362 100644
--- a/marchenko_full/raytime.c
+++ b/marchenko_full/raytime.c
@@ -64,7 +64,7 @@ int raytime(float *time, float *ampl, int *xnx, float *xrcv, float *xsrc, float
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;
@@ -139,11 +139,13 @@ int raytime(float *time, float *ampl, int *xnx, float *xrcv, float *xsrc, float
*/
/* 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);
@@ -170,40 +172,41 @@ int raytime(float *time, float *ampl, int *xnx, float *xrcv, float *xsrc, float
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;
-
- 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];
- coordgx.y = 0;
-
- 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);
- }
+ 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;
+
+ xnx[(izshot*shot.nx)+ixshot] = rec.n;
+ xsrc[(izshot*shot.nx)+ixshot] = mod.x0+mod.dx*shot.x[ixshot];
+ zsrc[(izshot*shot.nx)+ixshot] = mod.z0+mod.dz*shot.z[izshot];
+
+ 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);
+
+ xrcv[((izshot*shot.nx)+ixshot)*rec.n + irec] = (mod.x0+rec.x[0]*mod.dx) + ((rec.x[1]-rec.x[0])*mod.dx*((float)irec));
+ 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[ishot]);
hdr.sdepth = 1000*(mod.z0+mod.dz*shot.z[ishot]);
@@ -229,7 +232,7 @@ int raytime(float *time, float *ampl, int *xnx, float *xrcv, float *xsrc, float
assert(nwrite == rec.n);
fflush(fpa);
}*/
-
+ }
} /* end of loop over number of shots */
//fclose(fpt);
//if (ray.geomspread) fclose(fpa);
diff --git a/marchenko_full/raytime.h b/marchenko_full/raytime.h
index dea8b1be838096a00a27998b005eea0983230a5f..464553958272a713004412b4cc1afdd8b7fd4d9e 100644
--- a/marchenko_full/raytime.h
+++ b/marchenko_full/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;