Migration for all constant-offset sections

Another way to greatly improve the speed of the algorithm is to process all the constant-offset sections at the same time. In this case the phase loop can be computed once for all the offsets. The need to keep the constant-offset sections separated is accomplished by observing that the DSR equation can be rewritten in an offset separable form  

$$\begin{array} {lcl} p(t=0,k_y,h=0,z) & = & \displaystyle{ {\... ...,h) \int dk_h \; e^{ik_z(\omega,k_y,k_h)z-ik_h h}}.\end{array}$$ (7)

The integral in h stacks the independently prestack migrated constant-offset sections. As we are interested in each constant-offset section, we rewrite the DSR in an offset separable form  

$$p(t=0,k_y,h,z)= {\int d\omega \; p(\omega,k_y,h) \int dk_h \; e^{ik_z(\omega,k_y,k_h)z-ik_h h}} \\$$ (8)

and observe that the integral in k_h can be replaced by a fast Fourier transform of the exponential

$$e^{ik_z(\omega,k_y,k_h)z}.$$

The algorithm becomes:

		FFT along t,y axes $p(t,y;h_0) \rightarrow P(\omega,k_y)$ 
		do z 
		 		do ky 
		 		do $\omega$ 
		 		 		do kh 
		 		 		$ph(\omega,k_y,k_h)=ph(\omega,k_y,k_h)e^{ik_z dz}$ 
		 		 		$phase(\omega,k_y,k_h)=ph(\omega,k_y,k_h)$ 
		 		FFT phase along kh axis
		 		$M(z,k_y,h)=M(z,k_y,h)+P(\omega,k_y,h)*phase(\omega,k_y,h)$

This algorithm can be executed in parallel for all the values of $\omega,k_y$.Note that for depth variable velocity the phase is carried from one depth level to another, corresponding to an integral in depth of the phase terms corresponding to each depth slice.