This pseudo code Fourier transforms a wavefield
observed at the earth's surface ,and then it marches that wavefield down into the earth ()filling up a three-dimensional function, .Then it selects *t*=0, the time of the exploding reflectors
by summing over all frequencies .(Mathematically,
this is like
finding the signal at 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 .But it is easy to intertwine the downward continuation
with the summation over 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 .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

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 *up*ward extrapolation.
In principle,
the three loops on , *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()` .

12/26/2000