added time interpolation in Vz for up-down decomposition

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;

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

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");

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;

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));