diff --git a/fdemmodc/demo/em.scr b/fdemmodc/demo/em.scr
index 3565306571d3fd078c01705320dd17ac9cfb63f9..441956ffa2e97a4faf858e7294d42ed2722d7e41 100755
--- a/fdemmodc/demo/em.scr
+++ b/fdemmodc/demo/em.scr
@@ -40,12 +40,12 @@ makewave fp=250e+6 dt=$dt file_out=wave.su nt=4096 t0=6e-9 verbose=1
file_ks=syncl_ro.su \
file_src=wave.su \
file_rcv=shot_zsrc1150.su \
- src_type=1 \
+ src_type=7 \
src_orient=1 \
src_injectionrate=1 \
dtrcv=1e-10 \
rec_delay=6e-9 \
- verbose=3 \
+ verbose=6 \
tmod=2e-7 \
dt=$dt \
dxrcv=0.02 \
diff --git a/fdemmodc/em4.c b/fdemmodc/em4.c
index 760399c1ea12f990e7959a404c40c8065efafaf5..98673278a357ea891e0035e204c0008cc711ee57 100644
--- a/fdemmodc/em4.c
+++ b/fdemmodc/em4.c
@@ -5,16 +5,15 @@
#include<assert.h>
#include"fdelmodc.h"
-int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float *vx, float *vz, float *tzz,
-float *txx, float *txz, float *rox, float *roz, float *l2m, float **src_nwav, int verbose);
+int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float **src_nwav, int verbose);
int storeSourceOnSurface(modPar mod, srcPar src, bndPar bnd, int ixsrc, int izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, int verbose);
int reStoreSourceOnSurface(modPar mod, srcPar src, bndPar bnd, int ixsrc, int izsrc, float *vx, float *vz, float *tzz, float *txx, float *txz, int verbose);
-int boundariesP(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int verbose);
+int boundariesP(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int itime, int verbose);
-int boundariesV(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int verbose);
+int boundariesV(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float *txx, float *txz, float *rox, float *roz, float *l2m, float *lam, float *mul, int itime, int verbose);
int em4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, int izsrc, float **src_nwav, float *hz, float *hx, float *Ey, float *eprs, float *ksigma, float *mu, int verbose)
{
@@ -93,6 +92,7 @@ int em4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, in
c1*(Ey[ix*n1+iz] - Ey[(ix-1)*n1+iz]) +
c2*(Ey[(ix+1)*n1+iz] - Ey[(ix-2)*n1+iz]));
//if (hz[ix*n1+iz] > 0.1*FLT_MAX) fprintf(stderr,"%d: hz[%d %d] = %e\n", itime, ix, iz, hz[ix*n1+iz]);
+ if (hz[ix*n1+iz] != 0.0) fprintf(stderr,"%d: hz[%d %d] = %e\n", itime, ix, iz, hz[ix*n1+iz]);
}
}
@@ -105,16 +105,19 @@ int em4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, in
c1*(Ey[ix*n1+iz] - Ey[ix*n1+iz-1]) +
c2*(Ey[ix*n1+iz+1] - Ey[ix*n1+iz-2]));
//if (hx[ix*n1+iz] > 0.1*FLT_MAX) fprintf(stderr,"%d: hx[%d %d] = %e\n", itime, ix, iz, hx[ix*n1+iz]);
+ if (fabs(hx[ix*n1+iz]) != 0.0) fprintf(stderr,"%d: hx[%d %d] = %e\n", itime, ix, iz, hx[ix*n1+iz]);
}
}
/* Add force source */
if (src.type > 5) {
applySource(mod, src, wav, bnd, itime, ixsrc, izsrc, hz, hx, Ey, NULL, NULL, mu, mu, eprs, src_nwav, verbose);
+// fprintf(stderr,"%d: hx[%d %d] = %e hz=%e \n", itime, ixsrc, izsrc, hx[ixsrc*n1+izsrc], hz[ixsrc*n1+izsrc]);
+//fprintf(stderr,"src.n = %d\n", src.n);
}
/* boundary condition clears velocities on boundaries */
- boundariesP(mod, bnd, hz, hx, Ey, NULL, NULL, mu, mu, eprs, NULL, NULL, verbose);
+ boundariesP(mod, bnd, hz, hx, Ey, NULL, NULL, mu, mu, eprs, NULL, NULL, itime, verbose);
/* calculate p/tzz for all grid points except on the virtual boundary */
for (ix=mod.ioPx; ix<mod.iePx; ix++) {
@@ -137,6 +140,7 @@ int em4(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int ixsrc, in
for (iz=mod.ioPz; iz<mod.iePz; iz++) {
Ey[ix*n1+iz] -= eprs[ix*n1+iz]*(dzhx[iz]+dxhz[iz]) + ksigma[ix*n1+iz]*Ey[ix*n1+iz];
//if (Ey[ix*n1+iz] > 0.1*FLT_MAX) fprintf(stderr,"%d: Ey[%d %d] = %e\n", itime, ix, iz, Ey[ix*n1+iz]);
+ if (Ey[ix*n1+iz] != 0.0) fprintf(stderr,"%d: Ey[%d %d] = %e\n", itime, ix, iz, Ey[ix*n1+iz]);
}
}
diff --git a/fdemmodc/fdemmodc.c b/fdemmodc/fdemmodc.c
index fbd2ffb3f9a29ec3e8e275112161a6bbb525e9da..85d876c2f8dbd1aa164f4df4186b3c465dbdcff3 100644
--- a/fdemmodc/fdemmodc.c
+++ b/fdemmodc/fdemmodc.c
@@ -16,7 +16,7 @@ int getEmParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *
int readEmModel(modPar mod, bndPar bnd, float *eprs, float *ksigma, float *mu);
-int defineSource(wavPar wav, srcPar src, float **src_nwav, int reverse, int verbose);
+int defineSource(wavPar wav, srcPar src, modPar mod, float **src_nwav, int reverse, int verbose);
int writeSrcRecPos(modPar *mod, recPar *rec, srcPar *src, shotPar *shot);
@@ -273,7 +273,7 @@ int main(int argc, char **argv)
}
}
- defineSource(wav, src, src_nwav, mod.grid_dir, verbose);
+ defineSource(wav, src, mod, src_nwav, mod.grid_dir, verbose);
/* allocate arrays for wavefield and receiver arrays */
diff --git a/fdemmodc/getEmParameters.c b/fdemmodc/getEmParameters.c
index 6e49d7eae143b782305953654b149198d8c3394b..95bb85f65b0034a34e7733e66568a95244c097c7 100644
--- a/fdemmodc/getEmParameters.c
+++ b/fdemmodc/getEmParameters.c
@@ -20,6 +20,8 @@
* The Netherlands
**/
+int optncr(int n);
+
float gaussGen();
int getModelInfo(char *file_name, int *n1, int *n2, float *d1, float *d2, float *f1, float *f2, float *min, float *max, int *axis, int zeroch, int verbose);
@@ -36,7 +38,7 @@ int writesufile(char *filename, float *data, int n1, int n2, float f1, float f2,
int getEmParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *src, shotPar *shot, bndPar *bnd, int verbose)
{
int isnapmax1, isnapmax2, isnapmax, sna_nrsna;
- int n1, n2, nx, nz, nsrc, ix, axis, ioPz, is0;
+ int n1, n2, nx, nz, nsrc, ix, axis, ioPz, is0, optn;
int idzshot, idxshot;
int src_ix0, src_iz0, src_ix1, src_iz1;
int disable_check;
@@ -49,7 +51,7 @@ int getEmParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *
float tsnap1, tsnap2, dtsnap, dxsnap, dzsnap, dtrcv;
float xsnap1, xsnap2, zsnap1, zsnap2, xmax, zmax;
float xsrc1, xsrc2, zsrc1, zsrc2, tsrc1, tsrc2, tlength, tactive;
- float src_angle, src_velo, p, grad2rad, rdelay;
+ float src_angle, src_velo, p, grad2rad, rdelay, scaledt;
float *xsrca, *zsrca, rrcv;
float rsrc, oxsrc, ozsrc, dphisrc, ncsrc;
size_t nsamp;
@@ -102,13 +104,28 @@ int getEmParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *
/* define wavelet(s), modeling time and wavelet maximum frequency */
if (wav->file_src!=NULL) {
- getWaveletInfo(wav->file_src, &wav->nt, &wav->nx, &wav->dt, &d2, &f1, &f2, &fmax, &ntraces, verbose);
- if (wav->dt <= 0.0) {
+ getWaveletInfo(wav->file_src, &wav->ns, &wav->nx, &wav->ds, &d2, &f1, &f2, &fmax, &ntraces, verbose);
+ if (wav->ds <= 0.0) {
vwarn("dt in wavelet (file_src) equal to 0.0 or negative.");
vwarn("Use parameter dt= to overule dt from file_src.");
}
+ wav->nt = wav->ns;
+ wav->dt = wav->ds;
if(!getparfloat("tmod",&mod->tmod)) mod->tmod = (wav->nt-1)*wav->dt;
if(!getparfloat("dt",&mod->dt)) mod->dt=wav->dt;
+ if (NINT(wav->ds*1000000) != NINT(mod->dt*1000000)) {
+ if (wav->dt > mod->dt) {
+ scaledt = wav->dt/mod->dt;
+ scaledt = floorf(wav->dt/mod->dt);
+ optn = optncr(wav->ns);
+ wav->nt = floorf(scaledt*optn);
+ vmess("file_src dt-scalefactor=%f : wav.dt=%e ==interpolated==> mod.dt=%e", scaledt, wav->dt, mod->dt);
+ wav->dt = mod->dt;
+ }
+ else {
+ wav->dt = mod->dt; /* in case if wav.dt is smaller than 1e-7 and can not be read by SU-getpar */
+ }
+ }
if(!getparfloat("fmax",&wav->fmax)) wav->fmax=fmax;
}
else {
@@ -714,7 +731,7 @@ int getEmParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *
vmess("*******************************************");
vmess("************* wavelet info ****************");
vmess("*******************************************");
- vmess("wav_nt = %6d wav_nx = %d", wav->nt, wav->nx);
+ vmess("wav_nt = %6d wav_nx = %d", wav->ns, wav->nx);
vmess("src_type = %6d src_orient = %d", src->type, src->orient);
vmess("fmax = %10.3e", fmax);
fprintf(stderr," %s: Source type : ",xargv[0]);