The migration program is similar to the film loop program. But there are some differences. The film loop program has ``do loops'' nested four deep. It produces results for many values of t. Migration requires a value only at t=0. So one loop is saved, which means that for the same amount of computer time, the space volume can be increased. Unfortunately, loss of a loop seems also to mean loss of a movie. With -domain migration, it seems that the only interesting thing to view is the input and the output.
The input for this process will probably be field data, unlike for the film loop movie, so there will not be an analytic representation in the -domain. The input will be in the time domain and will have to be Fourier transformed. The beginning of the program defines some pulses to simulate field data. The pulses are broadened impulses and should migrate to approximate semicircles. Exact impulses were not used because the departure of difference operators from differential operators would make a noisy mess.
Next the program Fourier transforms the pseudodata from the time domain into the -frequency domain.
Then comes the downward continuation of each frequency. This is a loop on depth z and on frequency .Either of these loops may be on the inside. The choice can be made for machine-dependent efficiency.
For migration an equation for upcoming waves is required, unlike the downgoing wave equation required for the film loop program. Change the sign of the z-axis in equation (53). This affects the sign of aa and the sign of the phase of cshift.
Another difference with the film loop program is that the input now has a time axis whereas the output is still a depth axis. It is customary and convenient to reorganize the calculation to plot travel-time depth, instead of depth, making the vertical axes on both input and output the same. Using , equivalently , the chain rule gives
(61) |
(62) |
In the program, the time sample size dt and the travel-time depth sample dtau are taken to be unity, so the maximum frequency is the Nyquist. Notice that the frequency loop covers only the positive frequency axis. The negative frequencies serve only to keep the time function real, a task that is more easily done by simply taking the real part. A program accompanies.
#% Migration in the (omega,x,z)-domain subroutine kjartjac() real p(48,64), pi, alpha, dt, dtau, dw, w0, omega complex cp(48,64), cd(48), ce(48), cf(48), aa, a, b, c, cshift integer ix, nx, iz, nz, iw, nw, it, nt, esize nt= 64; nz= nt; nx= 48; pi= 3.141592 dt= 1.; dtau= 1.; w0=-pi/dt; dw= 2*pi/(dt*nt); nw= nt/2; alpha = .25 # alpha = v*v*dtau/(4*dx*dx) do iz= 1, nz { do ix=1,nx { p(ix,iz) = 0.; cp(ix,iz)=0. }} do it= nt/3, nt, nt/4{ do ix= 1, 4 # Broadened impulse source { cp(ix,it) = (5.-ix); cp(ix,it+1) = (5.-ix) }} call ft2axis( 0, 1., nx,nt, cp) do iz= 1, nz { do iw= 2, nw { omega = w0 + dw*(iw-1) aa = - alpha /( (0.,-1.)*omega ) a = -aa; b = 1.+2.*aa; c = -aa do ix= 2, nx-1 cd(ix) = aa*cp(ix+1,iw) + (1.-2.*aa)*cp(ix,iw) + aa*cp(ix-1,iw) cd(1) = 0.; cd(nx) = 0. call ctris( nx, -a, a, b, c, -c, cd, cp(1,iw)) cshift = cexp( cmplx( 0.,-omega*dtau)) do ix= 1, nx cp(ix,iw) = cp(ix,iw) * cshift do ix= 1, nx p(ix,iz) = p(ix,iz)+cp(ix,iw) # p(t=0) = Sum P(omega) }} esize=4 to history: integer n1:nx, n2:nz, esize call hclose() call srite( 'out', p, nx*nz*4 ) return; end
The output of the program is shown in Figure 12. Mainly, you see semicircle approximations. There are also some artifacts at late time that may be -domain wraparounds. The input pulses were apparently sufficiently broad-banded in dip that the figure provides a preview of the fact, to be proved later, that the actual semicircle approximation is an ellipse going through the origin.
kjartjac
Figure 12 Output of the program kjartjac: semicircle approximations. |
Notice that the waveform of the original pulses was a symmetric function of time, whereas the semicircles exhibit a waveform that is neither symmetric nor antisymmetric, but is a 45 phase-shifted pulse. Waves from a point in a three-dimensional world would have a phase shift of 90.Waves from a two-dimensional exploding reflector in a three-dimensional world have the 45 phase shift.