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
#
subroutine kfastaax(conj,add,anti,vrms,    t0,dt,dx,nt,nx,modl,      data)
integer 	    conj,add,                       nt,nx
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