From 4428129561aea6dfdd730c2700ebc934e391c93f Mon Sep 17 00:00:00 2001
From: JanThorbecke <janth@xs4all.nl>
Date: Thu, 25 Feb 2021 15:46:38 +0100
Subject: [PATCH] adding negative-time for dipping plane waves
---
fdelmodc/applySource.c | 3 ++-
fdelmodc/fdelmodc.c | 8 +++++++-
fdelmodc/fdelmodc.h | 1 +
fdelmodc/getParameters.c | 17 ++++++++++++++++-
4 files changed, 26 insertions(+), 3 deletions(-)
diff --git a/fdelmodc/applySource.c b/fdelmodc/applySource.c
index d29016e..ea1a613 100644
--- a/fdelmodc/applySource.c
+++ b/fdelmodc/applySource.c
@@ -97,7 +97,8 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
ix = ixsrc + ibndx + is0 + isrc;
iz = izsrc + ibndz;
}
- time = itime*dt - src.tbeg[isrc];
+ /* mod.t0 negative time correction for dipping plane waves */
+ time = itime*dt - src.tbeg[isrc] + mod.t0;
id1 = floor(time/dt);
id2 = id1+1;
diff --git a/fdelmodc/fdelmodc.c b/fdelmodc/fdelmodc.c
index 76229d4..694f71e 100644
--- a/fdelmodc/fdelmodc.c
+++ b/fdelmodc/fdelmodc.c
@@ -610,7 +610,11 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, izsrc, it, src_nwav, verbose)
writeToFile = ! ( (((it-rec.delay)/rec.skipdt)+1)%rec.nt );
itwritten = fileno*(rec.nt)*rec.skipdt;
/* Note that time step it=0 (t=0 for t**-fields t=-1/2 dt for v*-field) is not recorded */
- isam = (it-rec.delay-itwritten)/rec.skipdt+1;
+ //isam = (it-rec.delay-itwritten)/rec.skipdt+1;
+ /* negative time correction with mod.t0 for dipping plane waves modeling */
+ isam = (it-rec.delay-itwritten+NINT(mod.t0/mod.dt))/rec.skipdt+1;
+ if (isam < 0) isam = rec.nt+isam;
+
/* store time at receiver positions */
getRecTimes(mod, rec, bnd, it, isam, vx, vz, tzz, txx, txz,
l2m, rox, roz,
@@ -682,6 +686,8 @@ shared (shot, bnd, mod, src, wav, rec, ixsrc, izsrc, it, src_nwav, verbose)
}
}
}
+ /* negative time correction with mod.t0 for dipping plane waves modeling */
+ isam += NINT(-mod.t0/mod.dt)/rec.skipdt;
writeRec(rec, mod, bnd, wav, ixsrc, izsrc, isam+1, ishot, fileno,
rec_vx, rec_vz, rec_txx, rec_tzz, rec_txz,
rec_p, rec_pp, rec_ss, rec_udp, rec_udvz, verbose);
diff --git a/fdelmodc/fdelmodc.h b/fdelmodc/fdelmodc.h
index f68d00c..860a035 100644
--- a/fdelmodc/fdelmodc.h
+++ b/fdelmodc/fdelmodc.h
@@ -73,6 +73,7 @@ typedef struct _modelPar { /* Model Parameters */
int nt;
float z0;
float x0;
+ float t0;
/* medium max/min values */
float cp_min;
float cp_max;
diff --git a/fdelmodc/getParameters.c b/fdelmodc/getParameters.c
index 0894e7c..988db0e 100644
--- a/fdelmodc/getParameters.c
+++ b/fdelmodc/getParameters.c
@@ -174,6 +174,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
if (!getparfloat("rec_delay",&rdelay)) rdelay=0.0;
rec->delay=NINT(rdelay/mod->dt);
mod->nt = NINT(mod->tmod/mod->dt);
+ mod->t0 = 0.0;
dt = mod->dt;
if (!getparint("src_type",&src->type)) src->type=1;
@@ -878,6 +879,10 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
vwarn("number of gridpoints in X. Plane wave will be clipped to the edges of the model");
nsrc = mod->nx;
}
+ 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);
/* for a source defined on mutliple gridpoint calculate p delay factor */
@@ -887,6 +892,11 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
src->tend = (float *)malloc(nsrc*sizeof(float));
grad2rad = 17.453292e-3;
p = sin(src_angle*grad2rad)/src_velo;
+ for (is=0; is<nsrc; is++) {
+ src->tbeg[is] = (is-src_ix0)*dx*p;
+ }
+ mod->t0 = MIN(-src_ix0*dx*p,(nsrc-src_ix0)*dx*p);
+/*
if (p < 0.0) {
for (is=0; is<nsrc; is++) {
src->tbeg[is] = fabsf((nsrc-is-1)*dx*p);
@@ -899,11 +909,13 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
src->tbeg[is] = is*dx*p;
}
}
+*/
for (is=0; is<nsrc; is++) {
src->tend[is] = src->tbeg[is] + (wav->nt-1)*wav->dt;
}
is0 = -1*floor((nsrc-1)/2);
+ is0 = -src_ix0;
for (is=0; is<nsrc; is++) {
src->x[is] = is0 + is;
src->z[is] = 0;
@@ -986,7 +998,10 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
vmess("Areal source array is defined with %d sources.",nsrc);
/* vmess("Memory requirement for sources = %.2f MB.",sizeof(float)*(wav->nx*(wav->nt/(1024.0*1024.0))));*/
vmess("Memory requirement for sources = %.2f MB.",sizeof(float)*(nsamp/(1024.0*1024.0)));
- if (src->plane) vmess("Computed p-value = %f.",p);
+ if (src->plane) {
+ vmess("Computed p-value = %f.",p);
+ vmess("Maximum negative time delay is %f\n", mod->t0);
+ }
}
if (src->random) {
vmess("Sources are placed at random locations in domain: ");
--
GitLab