POST-STACK TIME MIGRATION (FAST, ANTIALIASED)

module kirchnew {
  integer                      :: nt, nx, sw
  real                         :: t0, dt, dx
  real, dimension (:), pointer :: vrms
#% _init (vrms, t0,dt,dx, nt,nx, sw)
#% _lop (modl(nt,nx), data(nt,nx))
  integer :: ix,iz,it,ib,iy, minx(2),maxx(2), is,i
  real	:: amp,t,z,b,db,f,g
	maxx(1) = nx;	minx(2) = 1
  do iz= 1,nt-1 {     z = t0 + dt * (iz-1)		# vertical traveltime
  do it= nt,iz+1,-1 { t = t0 + dt * (it-1)		# time shift
		    b = sqrt(t*t - z*z); db = dx*b*2./(vrms(iz)*t)
						if(db < dt .or. sw == 1) exit
		    f = 0.5*vrms(iz)*b/dx; iy = f; f = f-iy; i = iy+1; g = 1.-f
						if(i >= nx)		cycle
	amp = (z / (t+dt)) * sqrt(nt*dt / (t+dt)) * (dt / db)

minx(1) = 1+i;  maxx(2) = nx-i
  do is= 1,2   {  iy = -iy; i = -i		# two branches of hyperbola
  do ix=	minx(is),	maxx(is) {
		if( adj)	
			modl(iz,ix)    = modl(iz,ix) + data(it,ix+iy)*amp*g +
						       data(it,ix+i )*amp*f
		else {
			data(it,ix+iy) = data(it,ix+iy) + modl(iz,ix)*amp*g
			data(it,ix+i ) = data(it,ix+i ) + modl(iz,ix)*amp*f
		}
	}}}
  do ib= 0, nx-1 {	   b = dx*ib*2./vrms(iz); iy = ib	# space shift
		   t = sqrt(z*z  + b*b); db = dx*b*2./(vrms(iz)*t)
						if(db > dt .or. sw == 2) exit

f = (t-t0)/dt; i = f; it = i+1; f = f-i ; i = it+1; g = 1.-f
						if( i > nt)		exit

amp = (z / (t+dt)) * sqrt(nt*dt / (t+dt)); if(ib == 0) amp = amp*0.5

minx(1) = 1+iy; maxx(2) = nx-iy
  do is= 1,2   {	iy = -iy		# two branches of hyperbola
  do ix=	minx(is),	maxx(is) {
		if( adj)	
			modl(iz,ix)    = modl(iz,ix) + data(it,ix+iy)*amp*g +
						       data(i ,ix+iy)*amp*f
		else {
			data(it,ix+iy) = data(it,ix+iy) + modl(iz,ix)*amp*g
			data(i ,ix+iy) = data(i ,ix+iy) + modl(iz,ix)*amp*f
		}
	}}}
  }
}