Pseudocode to working code

Next is the task of imaging. Conceptually, at each depth an inverse Fourier transform is followed by selection of its value at t = 0. (Reflectors explode at t=0). Since only the Fourier transform at one point, t = 0, is needed, other times need not be be computed. We know the $\omega =0$ Fourier component is found by the sum over all time, analogously, the t=0 component is found as the sum over all $\omega$.(This is easily shown by substituting t=0 into the inverse Fourier integral.) Finally, inverse Fourier transform k_x to x. The migration process, computing the image from the upcoming wave u, may be summarized in the following pseudo code:

$$\vbox{ \begin{tabbing} $U(\omega, k_x, \tau = 0) = FT[u(t, ... ... gt image$(x, \tau) = FT[$Image$(k_x, \tau)]$ \\ \end{tabbing}}$$

This pseudo code Fourier transforms a wavefield observed at the earth's surface $\tau =0$,and then it marches that wavefield down into the earth ($\tau\gt$)filling up a three-dimensional function, $U(\omega, k_x, \tau)$.Then it selects t=0, the time of the exploding reflectors by summing over all frequencies $\omega$.(Mathematically, this is like finding the signal at $\omega =0$ by summing over all t).

Turning from pseudocode to real code, an important practical reality is that computer memories are not big enough for the three-dimensional function $U(\omega, k_x, \tau)$.But it is easy to intertwine the downward continuation with the summation over $\omega$so a three-dimensional function need not be kept in memory. This is done in the real code in subroutine phasemig().  

subroutine phasemig( up, nt, nx, dt, dx, image, ntau, dtau, vel)
integer                  nt, nx,                ntau,         iw,nw,ikx,itau
real    dt,dx, w,w0,dw, kx,kx0,dkx,dtau, vel, sig1,sig2,pi, w2, vkx2
complex up(nt,nx), image(ntau,nx), cc
                        pi = 3.14159265;        sig1 = +1.;    sig2 = -1.
call ft1axis( 0, sig1, nt, nx, up)   
call ft2axis( 0, sig2, nt, nx, up)
        nw = nt;             w0  = -pi/dt;      dw = 2.*pi/(nt*dt)
                             kx0 = -pi/dx;      dkx= 2.*pi/(nx*dx)
call null( image, ntau*nx*2)
do iw   = 2, nw  {      w  = w0  + (iw -1) * dw
do ikx  = 2, nx  {      kx = kx0 + (ikx-1) * dkx
        w2   = w * w
        vkx2 = vel*vel * kx*kx / 4.
        if( w2 > vkx2 ) {
                cc = cexp( cmplx( 0., - w * dtau * sqrt( 1. - vkx2/w2 ) ) )
                do itau = 1, ntau {
                        up(iw,ikx) = up(iw,ikx) * cc
                        image(itau,ikx) = image(itau,ikx) + up(iw,ikx)
                        }
                }
        }}
call ft2axis( 1, sig2, ntau, nx, image)
return; end

Conjugate migration (modeling) proceeds in much the same way. Beginning from an upcoming wave that is zero at great depth, the wave is marched upward in steps by multiplication with $\exp ( i k_z \Delta z )$.As each level in the earth is passed, exploding reflectors from that level are added into the upcoming wave. Pseudo code for modeling the upcoming wave u is

$$\vbox{ \begin{tabbing} Image$(k_x, z) = FT[$image$(x, z)]$ \... ...gt \} \} \} \\ $u(t, x) = FT[U(\omega, k_x)]$ \\ \end{tabbing}}$$

Some real code for this job is in subroutine phasemod().  

subroutine phasemod( image, nz, nx, dz, dx, up, nt, dt, vel)
integer                     nz, nx,            nt,       iw,nw,ikx,iz
real    dt,dx,dz, w,w0,dw, kx,kx0,dkx, vel, sig1,sig2,pi, w2, vkx2
complex up(nt,nx), image(nz,nx), cc
                pi = 3.14159265;        sig1 = +1.;     sig2 = -1.
call ft2axis( 0, sig2, nz, nx, image)
        nw = nt;        w0  = -pi/dt;   dw = 2.*pi/(nt*dt)
                        kx0 = -pi/dx;   dkx= 2.*pi/(nx*dx)
call null( up, nw*nx*2)
do iw  = 2, nw  {  w  = w0  + (iw-1)  * dw     
do ikx = 2, nx  {  kx = kx0 + (ikx-1) * dkx
        w2   = w * w
        vkx2 = vel*vel * kx*kx / 4.
        if( w2 > vkx2 ) {
                cc = cexp( cmplx( 0., w * dz * sqrt(1. - vkx2/w2) ))
                do iz = nz, 1, -1 
                        up(iw,ikx) = up(iw,ikx) * cc  + image(iz,ikx)
                }
        }}
call ft1axis( 1, sig1, nt, nx, up)
call ft2axis( 1, sig2, nt, nx, up)
return; end

The positive sign in the complex exponential is a combination of two negatives, the up coming wave and the upward extrapolation. In principle, the three loops on $\omega$, k_x, and z are interchangeable, however, since this tutorial program uses a velocity v that is a constant function of depth, I speeded it by a large factor by putting the z-loop on the inside and pulling the complex exponential out of the inner loop. Figure 2 was made with subroutine phasemod() .