diff --git a/fdelmodc/applySource.c b/fdelmodc/applySource.c
index d4d91fcc1f67827c3203b8086d3cc82bf642865c..30464ce09c4acd06f8d53acfaa118043b020c752 100644
--- a/fdelmodc/applySource.c
+++ b/fdelmodc/applySource.c
@@ -141,14 +141,16 @@ int applySource(modPar mod, srcPar src, wavPar wav, bndPar bnd, int itime, int i
/* Force source */
if (src.type == 6) {
-// vx[ix*n1+iz] += src_ampl*(dt/mod.dx)/(l2m[ix*n1+iz]);
- vx[ix*n1+iz] += src_ampl*rox[ix*n1+iz]/(l2m[ix*n1+iz]);
+ //vx[ix*n1+iz] += src_ampl*rox[ix*n1+iz]/(l2m[ix*n1+iz]);
+ /* stable implementation from "Numerical Techniques for Conservation Laws with Source Terms" by Justin Hudson */
+ vx[ix*n1+iz] = 0.5*(vx[(ix+1)*n1+iz]+vx[(ix-1)*n1+iz])+src_ampl*rox[ix*n1+iz]/(l2m[ix*n1+iz]);
}
else if (src.type == 7) {
- /* old amplitude setting does not obey reciprocity */
-// vz[ix*n1+iz] += src_ampl*(dt/mod.dx)/(l2m[ix*n1+iz]);
- vz[ix*n1+iz] += src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
-// fprintf(stderr,"ix=%d iz=%d l2m=%e rox=%e rho=%e dt=%e dx=%e\n", ix, iz, l2m[ix*n1+iz], roz[ix*n1+iz], dt/(mod.dx*roz[ix*n1+iz]), dt, mod.dx);
+ /* old amplitude setting does not obey reciprocity */
+ //vz[ix*n1+iz] += src_ampl*(dt/mod.dx)/(l2m[ix*n1+iz]);
+ //vz[ix*n1+iz] += src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
+ /* stable implementation from "Numerical Techniques for Conservation Laws with Source Terms" by Justin Hudson */
+ vz[ix*n1+iz] = 0.5*(vz[ix*n1+iz-1]+vz[ix*n1+iz+1])+src_ampl*roz[ix*n1+iz]/(l2m[ix*n1+iz]);
} /* src.type */
diff --git a/fdelmodc/defineSource.c b/fdelmodc/defineSource.c
index 0e10fe3e94c4f909511115286c100de208a5ac1b..456410a36851510ebecd93019893840e783f2047 100644
--- a/fdelmodc/defineSource.c
+++ b/fdelmodc/defineSource.c
@@ -113,15 +113,14 @@ int defineSource(wavPar wav, srcPar src, modPar mod, float **src_nwav, int rever
optnscale = optn;
nfreqscale = optnscale/2 + 1;
}
+// fprintf(stderr,"define S optn=%d ns=%d %e nt=%d %e\n", optn, wav.ns, wav.ds, optnscale, wav.dt);
ctrace = (complex *)calloc(nfreqscale,sizeof(complex));
trace = (float *)calloc(optnscale,sizeof(float));
- fprintf(stderr,"define S optn=%d ns=%d %e nt=%d %e\n", optn, wav.ns, wav.ds, optnscale, wav.dt);
df = 1.0/(optn*wav.ds);
deltom = 2.*M_PI*df;
scl = 1.0/optn;
- maxampl=0.0;
iwmax = nfreq;
for (i=0; i<wav.nx; i++) {
@@ -180,16 +179,29 @@ int defineSource(wavPar wav, srcPar src, modPar mod, float **src_nwav, int rever
/* avoid a (small) spike in the last sample
this is done to avoid diffraction from last wavelet sample
which will act as a pulse */
+ maxampl=0.0;
if (reverse) {
- for (j=0; j<wav.nt; j++) src_nwav[i][j] = scl*(trace[wav.nt-j-1]-trace[0]);
+ for (j=0; j<wav.nt; j++) {
+ src_nwav[i][j] = scl*(trace[wav.nt-j-1]-trace[0]);
+ maxampl = MAX(maxampl,fabs(src_nwav[i][j]));
+ }
}
else {
- for (j=0; j<wav.nt; j++) src_nwav[i][j] = scl*(trace[j]-trace[wav.nt-1]);
+ for (j=0; j<wav.nt; j++) {
+ src_nwav[i][j] = scl*(trace[j]-trace[wav.nt-1]);
+ maxampl = MAX(maxampl,fabs(src_nwav[i][j]));
+ }
}
- vmess("Wavelet sampling FFT-interpolated done for trace %d", i);
+ if (verbose > 3) vmess("Wavelet sampling (FFT-interpolated) done for trace %d", i);
}
-
}
+ /* set values smaller than 1e-5 maxampl to zero */
+ maxampl *= 1e-5;
+ for (i=0; i<wav.nx; i++) {
+ for (j=0; j<wav.nt; j++) {
+ if (fabs(src_nwav[i][j]) < maxampl) src_nwav[i][j] = 0.0;
+ }
+ }
free(ctrace);
free(trace);
diff --git a/fdelmodc/demo/fdelmodc_circ_medium.scr b/fdelmodc/demo/fdelmodc_circ_medium.scr
index ba7ebe9a7e42b4ec83c832e72bc4650be549fefc..390e4d53d42324a3cf123d26423e7e573d43472b 100755
--- a/fdelmodc/demo/fdelmodc_circ_medium.scr
+++ b/fdelmodc/demo/fdelmodc_circ_medium.scr
@@ -7,11 +7,11 @@
# receivers are placed on a circle
-export PATH=../../utils:$PATH:
+export PATH=../../bin:$PATH:
-dt=0.0001
-#makewave file_out=wavelet.su dt=0.0020 nt=1024 fp=13 shift=1 w=g2 verbose=1
-makewave w=fw fmin=0 flef=6 frig=89 fmax=95 dt=$dt file_out=wavefw.su nt=16384 t0=0.3 scale=0 scfft=1
+dt=0.0008
+makewave file_out=wave.su dt=$dt nt=1024 fp=25 shift=1 w=g2 verbose=1
+#makewave w=fw fmin=0 flef=6 frig=89 fmax=95 dt=$dt file_out=wavefw.su nt=16384 t0=0.3 scale=0 scfft=1
#sufft < wavefw.su | suamp | sugain scale=0.0001 | suxgraph style=normal
makemod file_base=model.su \
@@ -38,7 +38,7 @@ set -x
dxrcv=10 \
dtrcv=0.004 \
xsrc=0 zsrc=-2000 nshot=1 nsrc=1 \
- src_type=1 tmod=4.092 \
+ src_type=1 tmod=1.020 \
ntaper=100 \
left=2 right=2 bottom=2 top=2
@@ -46,17 +46,17 @@ export filecp=hom_cp.su
export filero=hom_ro.su
file_rcv=hom.su
-../fdelmodc \
+fdelmodc \
file_cp=$filecp file_den=$filero \
ischeme=1 \
- file_src=wavefw.su verbose=2 \
+ src_type=7 \
+ file_src=wave.su verbose=2 \
file_rcv=$file_rcv \
- rec_type_vz=0 rec_type_p=1 \
- xrcv1=-2000 xrcv2=2000 zrcv1=-2000 zrcv2=-2000 \
- dxrcv=10 \
+ rec_type_vz=1 rec_type_p=1 \
+ rrcv=500 dphi=1 \
dtrcv=0.004 \
- xsrc=0 zsrc=-2000 nshot=1 nsrc=1 \
- src_type=1 tmod=4.092 \
+ xsrc=0 zsrc=0 nshot=1 nsrc=1 \
+ tmod=1.020 \
ntaper=100 \
left=2 right=2 bottom=2 top=2
diff --git a/fdelmodc/demo/modelOilGas.scr b/fdelmodc/demo/modelOilGas.scr
index caa249863228abaaddb3c7150afafc7a46a9e19c..3ff23e18eb014ff7df15f09fe79fa0ae5d18e55c 100755
--- a/fdelmodc/demo/modelOilGas.scr
+++ b/fdelmodc/demo/modelOilGas.scr
@@ -7,47 +7,48 @@ rho=2500
dx=2.5
dt=0.0005
-
makemod sizex=5000 sizez=2500 dx=$dx dz=$dx cp0=$cp ro0=$rho \
orig=-2500,0 file_base=syncl.su \
intt=def x=-2500,0,2500 z=250,250,250 poly=0 cp=2300 ro=2000 \
intt=def x=-2500,-2000,-1000,-800,0,800,2500 z=650,650,700,750,900,750,600 poly=2 cp=2600 ro=2500 \
intt=def x=-2500,0,2500 z=1390,1390,1390 poly=0 cp=2000 ro=2000
-makewave fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
+makewave w=g1 fp=20 dt=$dt file_out=wave.su nt=4096 t0=0.1
export OMP_NUM_THREADS=2
-zsrc=0
zsrc=1100
+zsrc=15
-../fdelmodc \
+fdelmodc \
file_cp=syncl_cp.su ischeme=1 iorder=4 \
file_den=syncl_ro.su \
file_src=wave.su \
file_rcv=shot_fd.su \
- src_type=7 \
+ src_type=6 \
src_orient=1 \
- src_injectionrate=1 \
+ src_injectionrate=0 \
rec_type_vz=1 \
+ rec_type_vx=1 \
rec_type_p=1 \
rec_int_vz=2 \
dtrcv=0.004 \
rec_delay=0.1 \
- verbose=2 \
- tmod=2.55 \
+ verbose=4 \
+ tmod=2.0 \
dxrcv=10.0 \
xrcv1=-2250 xrcv2=2250 \
zrcv1=0 zrcv2=0 \
xsrc=0 zsrc=$zsrc \
- file_snap=snapF_$zsrc \
- tsnap1=0.1 tsnap2=4.0 dtsnap=0.05 dxsnap=$dx dzsnap=$dx \
ntaper=101 \
snapwithbnd=1 \
- left=1 right=1 top=1 bottom=1
+ file_snap=snapF_$zsrc \
+ tsnap1=0.1 tsnap2=4.0 dtsnap=0.05 dxsnap=$dx dzsnap=$dx \
+ left=2 right=2 top=2 bottom=2
exit
+
makemod sizex=5000 sizez=2500 dx=$dx dz=$dx cp0=$cp ro0=$rho \
orig=-2500,0 file_base=hom.su