REFERENCES

Bevc, D., and Claerbout, J., 1992, Fast anti-aliased Kirchhoff migration and modeling: SEP-75, 91-96.

Blondel, P., 1993, Constant-velocity anti-aliasing three-dimensional integral dip moveout: SEP-77, 49-58.

Claerbout, J. F., 1992, Anti aliasing: SEP-73, 371-390.

Lumley, D. E., 1993, Anti-aliased Kirchhoff 3-D migration: A salt intrusion example: SEP-77, 1-18.

Nichols, D., 1993, Integration along a line in a sampled space: SEP-77, 283-294.

APPENDIX A

The code used to generate the figures in this report is presented here.

# Kirchhoff migration and diffraction.  (fast anti-aliased)
#
subroutine kfastaax(conj,add,anti,velocity,t0,dt,dx,nt,nx, modl,       data)
integer 	    conj,add,                       nt,nx
real   		  	     anti,velocity,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,nx)

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

if(conj!=0) 
	do ix=1,nx
		call doubint( 1, 0, nt, ss(1,ix), data(1,ix))
do ih= -nx, nx {	h = dx * ih 			# h = offset
	do iz= 2, nt {	z = t0 + dt * (iz-1)		# z = travel-time depth
			t = sqrt( z**2  +  (h/velocity)**2 )
	slope = anti * ( (1./velocity)**2 ) * h / (t+dt)
	tp = t + abs(slope * dx) + dt;	 itp = 1.5 + (tp-t0) / dt
	tm = t - abs(slope * dx) - dt;	 itm = 1.5 + (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,1, -amp,2*amp,-amp, t0,dt, tm,t,tp,xstart,xend,ih,iz,nt,nx,modl,ss)
		}
	}
if(conj==0)
	do ix=1,nx
		call doubint( 0, 0, nt, ss(1,ix), data(1,ix))
return; end

# Double integration, first causal, then anticausal.
#
subroutine doubint( conj, add, n,             pp   , qq  )
integer             conj, add, n;	real  pp(n), qq(n)
temporary real tt(n)
call conjnull(      conj, add,                pp,n,  qq,n)
if( conj == 0 ) {	call causint( 0,   0, n,pp,  tt  )
			call causint( 1, add, n,qq,  tt  )
		}
else {			call causint( 1,   0, n,tt,  qq  )
			call causint( 0, add, n,tt,  pp  )
		}
return; end

# causal integral the old way with 1's on diagonal (SEP73 p375)
#
subroutine causintold( conj, add,        n,pp,  qq   )
integer i,          n, conj, add;   real pp(n), qq(n )
temporary real tt( n)
call conjnull(      conj, add,        pp,n,  qq,n )
if( conj == 0){				tt(1) = pp(1) 
				do i= 2, n
					tt(i) = tt(i-1) + pp(i)
				do i= 1, n
		  			qq(i) = qq(i) + tt(i)
		}

else {					tt(n) = qq(n)
				do i= n, 2, -1 
					tt(i-1) = tt(i) + qq(i-1)
				do i= 1, n
		  			pp(i) = pp(i) + tt(i)
		}
return; end

# causal integral with 1/2's on the diagonal
#
subroutine causintnew( conj, add,        n,pp,  qq   )
integer i,          n, conj, add;   real pp(n), qq(n )
temporary real tt( n)
call conjnull(      conj, add,        pp,n,  qq,n )
if( conj == 0){				tt(1) = pp(1) / 2.
				do i= 2, n
					tt(i) = tt(i-1) + (pp(i) + pp(i-1)) / 2.
		do i= 1, n
		  qq(i) = qq(i) + tt(i)
		}

else {					tt(n) = qq(n) / 2.
				do i= n, 2, -1 
					tt(i-1) = tt(i) + (qq(i) + qq(i-1)) / 2.
		do i= 1, n
		  pp(i) = pp(i) + tt(i)
		}
return; end

#  Scaled linear interpolation.
#  Modified for fast triangular weighted anti-aliased migration. DB 8-14-92
#  This routine loops over ix.
#
subroutine spot3x(conj,add, scm,sc,scp, t0,dt, tm,t,tp, xstart,xend,ih,iz,nt,nx,val,ss)
integer 	  conj,add,                             xstart,xend,ih,iz,nt,nx
real                        scm,sc,scp, t0,dt, tm,t,tp
real val(nt,nx),ss(nt,nx)
integer itm,it, itp, itcm,itc,itcp, ix
real    tcm,tcp,tc,  frm,fr,frp
temporary real vecm(nx),   vec(nx),  vecp(nx)
temporary real vecm1(nx), vec1(nx), vecp1(nx)

call conjnull(    conj, add,                     val,nt*nx, ss,nt*nx)
if (conj == 0){
	call null(vecm ,nx)
	call null(vec  ,nx)
	call null(vecp ,nx)
	call null(vecm1,nx)
	call null(vec1 ,nx)
	call null(vecp1,nx)
	}
tcm = (tm-t0) / dt;	itcm = tcm;	 itm = 1 + itcm;	frm = tcm - itcm
tc  = (t -t0) / dt;	itc  = tc;	 it  = 1 + itc;   	fr  = tc  - itc
tcp = (tp-t0) / dt;	itcp = tcp;	 itp = 1 + itcp;	frp = tcp - itcp

if(1 <= it & it < nt) { 
if( conj != 0) {
	do ix=1,nx {
		 vecm(ix) = ss(itm  ,ix)
		  vec(ix) = ss(it   ,ix)
		 vecp(ix) = ss(itp  ,ix)
		vecm1(ix) = ss(itm+1,ix)
		 vec1(ix) = ss(it+1 ,ix)
		vecp1(ix) = ss(itp+1,ix)
		}
	}
do ix=xstart,xend {
if( conj == 0) {
	 vecm(ix+ih) = vecm(ix+ih)  + (1.-frm) * val(iz,ix) * scm
 	  vec(ix+ih) = vec(ix+ih)   + (1.-fr ) * val(iz,ix) * sc
	 vecp(ix+ih) = vecp(ix+ih)  + (1.-frp) * val(iz,ix) * scp
	vecm1(ix+ih) = vecm1(ix+ih) +  frm     * val(iz,ix) * scm
	 vec1(ix+ih) = vec1(ix+ih)  +  fr      * val(iz,ix) * sc
	vecp1(ix+ih) = vecp1(ix+ih) +  frp     * val(iz,ix) * scp
	}
else	{
	val(iz,ix)=val(iz,ix)+((1.-fr )*vec(ix+ih )+fr *vec1(ix+ih ))*sc 
	val(iz,ix)=val(iz,ix)+((1.-frm)*vecm(ix+ih)+frm*vecm1(ix+ih))*scm
	val(iz,ix)=val(iz,ix)+((1.-frp)*vecp(ix+ih)+frp*vecp1(ix+ih))*scp
	}
}
if( conj == 0) {
	do ix=1,nx {
		ss(itm  ,ix) = vecm(ix)  + ss(itm  ,ix)
		ss(it   ,ix) = vec(ix)   + ss(it   ,ix)
		ss(itp  ,ix) = vecp(ix)  + ss(itp  ,ix)
		ss(itm+1,ix) = vecm1(ix) + ss(itm+1,ix)
		ss(it+1 ,ix) = vec1(ix)  + ss(it+1 ,ix)
		ss(itp+1,ix) = vecp1(ix) + ss(itp+1,ix)
		}
	}
}
else
	call erexit('spot3: at boundary')
return;	end