Next: TIMING TESTS Up: Bevc & Claerbout: Fast Previous: Introduction

IMPLEMENTATION

The fast anti-aliasing routine kfastaax() achieves most of its speed by re-ordering the do-loops of aamig() and trimo(). Re-ordering is done in the same way as for kirchslow() and kirchfast() so that the outer loop runs over offset. This allows overhead calculations to be reused many times. Additional speed is obtained by transposing the input and output so that the inner loop over ix is along stride instead of across stride. This inner loop over ix is performed in the linear interpolation subroutine spot3x(). The speed gained by transposing the input is achieved at the cost of assuming that velocity does not vary with offset.

Program listings are provided for kfastaax() and spot3x(). The subroutine doubint() performs double integration. For conj=0 it integrates once causally and then once anticausally.

For conj=1 kfastaax() migrates the input by summing along hyperbolic trajectories with triangular smoothing in time. For conj=0 kfastaax() models the data by spraying the input along hyperbolic trajectories.

For data modeling the program proceeds as follows:

• Transpose the input so that the inner loop is along stride.
• Loop over offset ih
• Within the iz loop the slope of the hyperbola is calculated and used to calculate the end points of the smoothing triangles. An amplitude factor which incorporates spherical divergence, obliquity, and triangle width is calculated.
• The triangle endpoints and amplitudes are passed to the interpolation routine spot3x() along with a column of the transposed input data. spot3x() returns six column vectors of linearly interpolated output.
• The output is transposed and integrated twice. The causal integration followed by anti-causal integration completes the triangular smoothing.

```# Kirchhoff migration and diffraction.  (fast anti-aliased)
# combination of kirchfast() and trimo().  Written by Dimitri Bevc 082492
#
real   			          vrms(nt),t0,dt,dx,     modl(nt,nx),data(nt,nx)
real        z, t, h,       tp, tm,   amp, anti, slope
integer ix,iz,it,ih,      itp,itm,  xstart, xend
temporary real ss(nt),trdata(nx,nt),trmodl(nx,nt)

call conjnull(      conj, add,                         modl,nt*nx,  data,nt*nx)
call null(     ss,nt) ; call null(trmodl,nx*nt) ; call null(trdata,nx*nt)

if(conj!=0)
do ix=1,nx {
call doubint( 1,  0, nt, ss, data(1,ix))
do it=1,nt
trdata(ix,it)=ss(it)
}
if(conj==0)
do ix=1,nx
do it=1,nt
trmodl(ix,it)=modl(it,ix)
do ih= -nx, nx {	h = dx * ih 			# h = offset
do iz= 1, nt {	z = t0 + dt * (iz-1)		# z = travel-time depth
t = sqrt( z**2  +  (h/vrms(iz))**2 )
slope = anti * ( (1./vrms(iz))**2 ) * h / (t+dt)
it  = 1.0 + (t-t0 ) / dt
tp = t + abs(slope * dx) + dt;	 itp = 1.0 + (tp-t0) / dt
tm = t - abs(slope * dx) - dt;	 itm = 1.0 + (tm-t0) / dt
if( itm <= 1  )  next
if( itp >= nt )  break
amp = 1.0 * sqrt( nt*dt/(t+dt)) * (z+dt)/(t+dt) * (dt/(dt+tp-tm)) ** 2
xstart = 1-ih;		xstart = max0(  1, xstart)
xend   = nx-ih;		xend   = min0( nx, xend)
call spot3x(conj, -amp,2*amp,-amp, t0,dt, tm,t,tp,xstart,xend,
ih,iz,nt,nx,trmodl(1,iz),trdata(1,itm  ),trdata(1,it  ),trdata(1,itp  ),
trdata(1,itm+1),trdata(1,it+1),trdata(1,itp+1))
}
}
if(conj!=0)
do ix=1,nx
do it=1,nt
modl(it,ix)=trmodl(ix,it)
if(conj==0)
do ix=1,nx {
do it=1,nt
ss(it)=trdata(ix,it)
call doubint( 0, add, nt, ss, data(1,ix))
}
return; end
```

```#  Scaled linear interpolation.
#  Modified for fast triangular weighted anti-aliased migration
#  by Dimitri Bevc 8-14-92
#
subroutine spot3x(conj, scm,sc,scp, t0,dt, tm,t,tp, xstart,xend,ih,iz,
nt,nx,  val,vecm ,vec, vecp,
vecm1,vec1,vecp1)
integer 	  conj,                             xstart,xend,ih,iz,nt,nx
real                        scm,sc,scp, t0,dt, tm,t,tp
real                 val(nx),vecm(nx) ,vec(nx), vecp(nx)
real                         vecm1(nx),vec1(nx),vecp1(nx)
integer itcm,itc,itcp, ix
real    tcm,tcp,tc,  frm,fr,frp

tcm = (tm-t0) / dt;	itcm = tcm;	 frm = tcm - itcm
tc  = (t -t0) / dt;	itc  = tc;	 fr  = tc  - itc
tcp = (tp-t0) / dt;	itcp = tcp;	 frp = tcp - itcp
if( conj == 0) {
do ix=xstart,xend
vecm(ix+ih) = vecm(ix+ih)  + (1.-frm) * val(ix) * scm
do ix=xstart,xend
vec(ix+ih) = vec(ix+ih)   + (1.-fr ) * val(ix) * sc
do ix=xstart,xend
vecp(ix+ih) = vecp(ix+ih)  + (1.-frp) * val(ix) * scp
do ix=xstart,xend
vecm1(ix+ih) = vecm1(ix+ih) +  frm     * val(ix) * scm
do ix=xstart,xend
vec1(ix+ih) = vec1(ix+ih)  +  fr      * val(ix) * sc
do ix=xstart,xend
vecp1(ix+ih) = vecp1(ix+ih) +  frp     * val(ix) * scp
}
else	{
do ix= xstart, xend
val(ix) = val(ix)+
((1.-fr )*vec (ix+ih) + fr *vec1 (ix+ih))*sc   +
((1.-frm)*vecm(ix+ih) + frm*vecm1(ix+ih))*scm  +
((1.-frp)*vecp(ix+ih) + frp*vecp1(ix+ih))*scp
}
return;	end
```

Next: TIMING TESTS Up: Bevc & Claerbout: Fast Previous: Introduction
Stanford Exploration Project
11/17/1997