To get smoother results I took the time axis to be continuous and the signal value at (t,x) to be distributed between the two points and .The two time points and the x-value are mapped to two slownesses .The signal from the (t,x)-pixel is sprayed into the horizontal line .To enable you to reproduce the result, I include the program below:
subroutine vspray( conj, nt,dt,t0, nx,dx,x0, tx, ns,ds,s0, zs) integer conj, it, nt, iz, nz, ix, nx, is, ns, isp, ism real tx(nt,nx), zs(nt,ns), scale real z,dz,z0, t,dt,t0, x,dx,x0, s,ds,s0, sm,sp, xm,xp, tm,tp nz=nt; dz=dt; z0=t0; if( conj == 0 ) { do is=1,ns; do iz=1,nz; zs(iz,is) = 0.} else { do ix=1,nx; do it=1,nt; tx(it,ix) = 0.} if( conj == 0) { do ix=1,nx; call halfdif ( 1, nt, tx(1,ix), tx(1,ix) )} do iz= 1, nz { z = z0 + dz*(iz-1) do ix= 1, nx { x = x0 + dx*(ix-1) do it= iz, nt { t = t0 + dt*(it-1) tm = t-dt/2; xm = x tp = t+dt/2; xp = x sm = (tm**2 -z**2)/xp**2; ism = 1.5+(sm-s0)/ds sp = (tp**2 -z**2)/xm**2; isp = 1.5+(sp-s0)/ds if( ism<2 ) next if( isp>ns) next scale = sqrt( t / (1.+isp-ism) ) / x do is= ism, isp { if( conj == 0) zs(iz ,is) = zs(iz ,is) + tx(it ,ix) * scale else tx(it ,ix) = tx(it ,ix) + zs(iz ,is) * scale } } } } if( conj != 0) { do ix=1,nx; call halfdif ( 0, nt, tx(1,ix), tx(1,ix) )} return; end
Figure 5 shows the result on the same inputs as in Figures 3 and 4.