raytime FD-option added amplitude estimation

fdelmodc/boundaries.c

@@ -1433,6 +1433,9 @@ int boundariesV(modPar mod, bndPar bnd, float *vx, float *vz, float *tzz, float
 			for (ix=mod.ioPx; ix<mod.iePx; ix++) {
 				iz = bnd.surface[ix];
 				tzz[ix*n1+iz] = 0.0;
+                //vz[ix*n1+iz] = -vz[ix*n1+iz+1];
+                //vz[ix*n1+iz-1] = -vz[ix*n1+iz+2];
+
 			}
 		}
 		if (bnd.rig==1) { /* free surface at right */

raytime/SrcPositions1.txt

-0.000000 0.000000
+1100.000000 0.000000

raytime/model.scr

@@ -17,13 +17,13 @@ export OMP_NUM_THREADS=4
     file_cp=syncl_cp.su  \
     file_rcv=z1.su \
     smoothwindow=15 \
-    verbose=3 \
+    verbose=5 \
     dxrcv=10.0 \
     xrcv1=-2500 xrcv2=2500 \
     zrcv1=0 zrcv2=0 \
     xsrc=0 zsrc=1100 \
 	nxshot=1 dxshot=10 \
-	nzshot=1 dzshot=10 useT2=1
+	nzshot=1 dzshot=10 useT2=1 method=fd
 
 exit
 

raytime/raytime.c

@@ -117,10 +117,10 @@ int main(int argc, char **argv)
 	double t0, t1, t2, tinit, tray, tio;
 	size_t size;
 	int nw, n1, ix, iz, ir, ixshot, izshot;
-	int irec, sbox, ipos;
+	int irec, sbox, ipos, nrx, nrz, nr;
     fcoord coordsx, coordgx, Time;
     icoord grid, isrc; 
-    float Jr, *ampl, *time, *ttime, *ttime_p;
+    float Jr, *ampl, *time, *ttime, *ttime_p, cp_average;
 	float dxrcv, dzrcv;
     segy hdr;
     char filetime[1024], fileamp[1024], *method;
@@ -285,8 +285,20 @@ private (coordgx,irec,Time,Jr)
 					ipos = ((izshot*shot.nx)+ixshot)*rec.n + irec;
 	
             		time[ipos] = ttime[rec.z[irec]*mod.nx+rec.x[irec]];
-            		ampl[ipos] = sqrt(time[ipos]);
-            		if (verbose>4) vmess("FD: shot=%f,%f receiver at %f(%d),%f(%d) T=%e Ampl=%f",coordsx.x, coordsx.z, coordgx.x, rec.x[irec], coordgx.z, rec.z[irec], time[ipos], ampl[ipos]); 
+					/* compute average velocity between source and receiver */
+ 					nrx = (rec.x[irec]-isrc.x);
+ 					nrz = (rec.z[irec]-isrc.z);
+ 					nr = abs(nrx) + abs(nrz);
+					cp_average = 0.0;
+        			for (ir=0; ir<nr; ir++) {
+						ix = isrc.x + floor((ir*nrx)/nr);
+						iz = isrc.z + floor((ir*nrz)/nr);
+						//fprintf(stderr,"ir=%d ix=%d iz=%d velocity=%f\n", ir, ix, iz, velocity[ix*mod.nz+iz]);
+						cp_average += velocity[ix*mod.nz+iz];
+					}
+					cp_average = cp_average/((float)nr);
+            		ampl[ipos] = sqrt(time[ipos]*cp_average);
+            		if (verbose>4) vmess("FD: shot=%f,%f receiver at %f(%d),%f(%d) T=%e V=%f Ampl=%f",coordsx.x, coordsx.z, coordgx.x, rec.x[irec], coordgx.z, rec.z[irec], time[ipos], cp_average, ampl[ipos]); 
         		}
             }
 	        t2=wallclock_time();