diff --git a/fdelmodc/decomposition.c b/fdelmodc/decomposition.c
index 7f66c12162fa4da523a834300aedb7c3fa4128a6..c50c7888dd21cb500604b2d0dc21f4335d41d8a7 100644
--- a/fdelmodc/decomposition.c
+++ b/fdelmodc/decomposition.c
@@ -94,7 +94,7 @@ void kxwdecomp(complex *rp, complex *rvz, complex *up, complex *down,
}
avrp = avrp/89.0;
avrvz = avrvz/89.0;
- if (verbose>=5) {
+ if (verbose>=4) {
writesufile("anglerp.su", angle, 90, 1, 0, 0, 1, 1);
writesufile("anglervz.su", &angle[90], 90, 1, 0, 0, 1, 1);
}
@@ -102,7 +102,7 @@ void kxwdecomp(complex *rp, complex *rvz, complex *up, complex *down,
if (angle[av] < avrp) angle[av] = 0.0;
if (angle[90+av] < avrvz) angle[90+av] = 0.0;
}
- if (verbose>=5) {
+ if (verbose>=4) {
writesufile("anglerp0.su", angle, 90, 1, 0, 0, 1, 1);
writesufile("anglervz0.su", &angle[90], 90, 1, 0, 0, 1, 1);
}
@@ -166,7 +166,7 @@ void decud(float om, float rho, float cp, float dx, int nkx, float fangle, float
stab = eps*eps*kp*kp;
/* make kw filter at maximum angle alfa */
- perc = 0.10; /* percentage of band to use for smooth filter */
+ perc = 0.15; /* percentage of band to use for smooth filter */
filter = (float *)malloc(nkx*sizeof(float));
kpos = kp*sin(M_PI*fangle/180.0);
kneg = -kpos;
diff --git a/fdelmodc/demo/decomposition.scr b/fdelmodc/demo/decomposition.scr
index cfc27b23de6dc94189641d8da4d28ec161dd0f2c..abb782a8bc5597f91f7ebd0a53fd58b867fdb76f 100755
--- a/fdelmodc/demo/decomposition.scr
+++ b/fdelmodc/demo/decomposition.scr
@@ -53,7 +53,6 @@ supsgraph < anglerp.su hbox=2 wbox=4 style=normal \
labelsize=10 label2='energy' label1='angle in degrees' \
titlesize=-1 f1=0 d1=1 d1num=10.0 x1beg=1 > anglerp.eps
-
supsimage < shot_fd_Fz_zsrc1150_rvz.su \
wbox=2 hbox=4 titlesize=-1 labelsize=10 verbose=1 \
label1="time [s]" label2="lateral position [m]" > shot_fd_Fz_zsrc1150_rvz.eps
diff --git a/fdelmodc/getParameters.c b/fdelmodc/getParameters.c
index 8301167368556b41f59ed24ec8912efd08c2d684..addd2468b1c7528639b23d546fc09842c9dd46d9 100644
--- a/fdelmodc/getParameters.c
+++ b/fdelmodc/getParameters.c
@@ -90,7 +90,7 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
getModelInfo(mod->file_cp, &nz, &nx, &dz, &dx, &sub_z0, &sub_x0, &cp_min, &cp_max, &axis, 1, verbose);
getModelInfo(mod->file_ro, &n1, &n2, &d1, &d2, &zstart, &xstart, &ro_min, &ro_max, &axis, 0, verbose);
- assert(ro_min != 0.0);
+ assert( (ro_min != 0.0) );
if (NINT(100*(dx/d2)) != 100)
vwarn("dx differs for file_cp and file_den!");
if (NINT(100*(dz/d1)) != 100)
@@ -1146,14 +1146,14 @@ int getParameters(modPar *mod, recPar *rec, snaPar *sna, wavPar *wav, srcPar *sr
vmess("izmin = %d izmax = %d ", rec->z[0], rec->z[rec->n-1]);
vmess("ixmin = %d ixmax = %d ", rec->x[0], rec->x[rec->n-1]);
if (rec->type.vx) {
- fprintf(stderr," %s: Receiver interpolation for Vx:",xargv[0]);
+ fprintf(stderr," %s: Receiver interpolation for Vx: ",xargv[0]);
if(rec->int_vx==0) fprintf(stderr,"vx->vx\n");
if(rec->int_vx==1) fprintf(stderr,"vx->vz\n");
if(rec->int_vx==2) fprintf(stderr,"vx->txx/tzz\n");
if(rec->int_vx==3) fprintf(stderr,"interpolate to postion of receiver\n");
}
if (rec->type.vz) {
- fprintf(stderr," %s: Receiver interpolation for Vz:",xargv[0]);
+ fprintf(stderr," %s: Receiver interpolation for Vz: ",xargv[0]);
if(rec->int_vz==0) fprintf(stderr,"vz->vz\n");
if(rec->int_vz==1) fprintf(stderr,"vz->vx\n");
if(rec->int_vz==2) fprintf(stderr,"vz->txx/tzz\n");
diff --git a/fdelmodc/getRecTimes.c b/fdelmodc/getRecTimes.c
index 8582c70a4175495ef630347ecbab1325ec8d203c..ccfaeec672908d2778887d797431bb8cd480b592 100644
--- a/fdelmodc/getRecTimes.c
+++ b/fdelmodc/getRecTimes.c
@@ -21,6 +21,7 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
int n1, ibndx, ibndz;
int irec, ix, iz, ix2, iz2, ix1, iz1;
float rdz, rdx, C00, C10, C01, C11;
+ float *vz_t, c1, c2, roz;
ibndx = mod.ioPx;
ibndz = mod.ioPz;
@@ -38,6 +39,8 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
* vx velocities have one sample less in x-direction
* vz velocities have one sample less in z-direction
* txz stresses have one sample less in z-direction and x-direction
+*
+* Note, in the acoustic scheme P is stored in the Tzz array.
***********************************************************************/
for (irec=0; irec<rec.n; irec++) {
@@ -210,10 +213,26 @@ int getRecTimes(modPar mod, recPar rec, bndPar bnd, int itime, int isam, float *
if (rec.type.ud) {
iz = rec.z[0]+ibndz;
iz2 = iz+1;
+ vz_t = (float *)calloc(2*mod.nax,sizeof(float));
+ /* P and Vz are staggered in time and need to correct for this */
+ /* -1- compute Vz at next time step and average with current time step */
+ c1 = 9.0/8.0;
+ c2 = -1.0/24.0;
+ roz = mod.dt/(mod.dx*rec.rho);
+ for (ix=mod.ioZx; ix<mod.ieZx; ix++) {
+ vz_t[ix] = vz[ix*n1+iz] - roz*(
+ c1*(tzz[ix*n1+iz] - tzz[ix*n1+iz-1]) +
+ c2*(tzz[ix*n1+iz+1] - tzz[ix*n1+iz-2]));
+ vz_t[mod.nax+ix] = vz[ix*n1+iz2] - roz*(
+ c1*(tzz[ix*n1+iz2] - tzz[ix*n1+iz2-1]) +
+ c2*(tzz[ix*n1+iz2+1] - tzz[ix*n1+iz2-2]));
+ }
for (ix=0; ix<mod.nax; ix++) {
- rec_udvz[ix*rec.nt+isam] = 0.5*(vz[ix*n1+iz2]+vz[ix*n1+iz]);
+ /* -2- compute average in time and depth to get Vz at same depth and time as P */
+ rec_udvz[ix*rec.nt+isam] = 0.25*(vz[ix*n1+iz2]+vz[ix*n1+iz]+vz_t[mod.nax+ix]+vz_t[ix]);
rec_udp[ix*rec.nt+isam] = tzz[ix*n1+iz];
}
+ free(vz_t);
}
return 0;
diff --git a/fdelmodc/writeRec.c b/fdelmodc/writeRec.c
index d3af4fd77af30d2c0dcd3c1f83e0761fa633e052..8910ed3622a8f43321023e4923e78e0d87a0defb 100644
--- a/fdelmodc/writeRec.c
+++ b/fdelmodc/writeRec.c
@@ -110,7 +110,7 @@ int writeRec(recPar rec, modPar mod, bndPar bnd, int ixsrc, int izsrc, int nsam,
if (bnd.lef==4 || bnd.lef==2) ibndx += bnd.ntap;
cp = rec.cp;
rho = rec.rho;
- fprintf(stderr,"Decomposition array with cp=%f rho=%f\n", cp, rho);
+ if (verbose) vmess("Decomposition array at z=%.2f with cp=%.2f rho=%.2f", rec.zr[0]+mod.z0, cp, rho);
rec_up = (float *)calloc(ntfft*nkx,sizeof(float));
rec_down= (float *)calloc(ntfft*nkx,sizeof(float));
crec_vz = (complex *)malloc(nfreq*nkx*sizeof(complex));