It is straightforward to adapt the simple programs
`phasemig()` and
`phasemod()`
to depth variable velocity.
As you might suspect, the two processes are adjoint to each other,
and for reasons given at the end of chapter
it is desirable to code them to be so.
With differential equations and their boundary conditions,
the concept of adjoint is more subtle than previous examples.
Thus, I postponed till here the development of adjoint code
for phase-shift migration.
This coding is a little strenuous,
but it affords a review of many basic concepts,
so we do so here. (Feel free to skip this section.)
It is nice to have a high quality code for this fundamental operation.

Many situations in physics are expressed by the differential equation

(18) |

(19) |

(20) |

(21) |

A recursive solution begins at the bottom
with *u _{3}*=

The adjoint (complex conjugate) of the delay operator *Z*
is the time advance operator .The adjoint of equation (21)
is given by

(22) |

We observe the wavefield only on the surface *z*=0,
so the adjointness of
equations (21)
and (22)
is academic because it relates the wavefield at all depths
with the source at all depths.
We need to truncate to its first coefficient *u _{0}*
since the upcoming wave is known only at the surface.
This truncation changes the adjoint in a curious way.
We rewrite equation (21)
using a truncation operator that is
the row matrix getting
.Its adjoint is
or
which looks like

(23) |

beginning from and continuing downward with

(24) |

A final feature of the migration application
is that the image is formed from by summing over all frequencies.
Although I believe the mathematics above
and the code in subroutine `gazadj()` ,
I ran the dot product test to be sure!

# Phase shift modeling and migration. (Warning: destroys its input!) # subroutine gazadj( adj, dt,dx, v,nt,nx, modl, data ) integer adj, nt,nx, iw, ikx, iz,nz complex eiktau, cup, modl(nt,nx), data(nt,nx) real dt,dx, v(nt), pi, w,w0,dw, kx,kx0,dkx,qi call adjnull( adj, 0, modl,nt*nx*2, data,nt*nx*2 ) pi = 4.*atan(1.); w0 = -pi/dt; dw = 2.*pi/(nt*dt); qi=.5/(nt*dt) nz = nt; kx0 = -pi/dx; dkx= 2.*pi/(nx*dx) if( adj == 0) call ft2axis( 0, -1., nz, nx, modl) else { call ft2axis( 0, -1., nt, nx, data) call ft1axis( 0, 1., nt, nx, data) } do ikx = 2, nx { kx = kx0 + (ikx-1) * dkx do iw = 2, 1+nt/2 { w = w0 + (iw -1) * dw if( adj== 0) { data(iw,ikx) = modl(nz,ikx) do iz = nz-1, 1, -1 data(iw,ikx) = data(iw,ikx) * eiktau(dt,w,v(iz)*kx,qi) + modl(iz,ikx) } else { cup = data(iw,ikx) do iz = 1, nz { modl(iz,ikx) = modl(iz,ikx) + cup cup = cup * conjg( eiktau(dt,w,v(iz)*kx,qi)) } } }} if( adj == 0) { call ft1axis( 1, 1., nt, nx, data) call ft2axis( 1, -1., nt, nx, data) } else { call ft2axis( 1, -1., nz, nx, modl) } return; end

Finally, a few small details about the code.
The loop on spatial frequency `ikx`
begins at `ikx=2`.
The reason for the 2, instead of a 1,
is to omit the Nyquist frequency.
If the Nyquist frequency were to be included,
it should be divided into one half at positive Nyquist and one half
at negative Nyquist, which would clutter the code
without adding practical value.
Another small detail is that the loop
on temporal frequency `iw`
begins at `iw=1+nt/2`
which effectly omits negative frequencies.
This is purely an economy measure.
Including the negative frequencies would assure that the
final image be real, no imaginary part.
Omitting negative frequencies simply gives an imaginary part
that can be thrown away,
and gives the same real image, scaled by a half.
The factor of two speed up makes these tiny compromises well worthwhile.

12/26/2000