POLAR COORDINATES

The problems associated with overturned waves and initial conditions can be overcome by calculating the Green's functions in polar coordinates. Van Trier and Symes 1991 use polar coordinates for similar reasons in their finite difference solution to the eikonal equation.

I will follow exactly the same steps in deriving the algorithms in polar coordinates as I did in rectangular coordinates. In a 2-D polar coordinate frame, the wave equation is,

$$\frac{\partial^2 P }{ \partial r^2 } + \frac{1}{r}\frac{\pa... ...heta^2} = \frac {1}{v^2} \frac{\partial^2 P }{ \partial t^2 }$$

When this equation is Fourier transformed over time and $\theta$ we obtain,

$$\frac{\partial^2 P }{ \partial r^2 } + \frac{1}{r}\frac{\par... ...ial r} - \frac{1}{r^2} k_\theta^2 P = -\frac{\omega^2}{v^2} P$$

If we assume a solution of the form,

 

P(r) = P(r₀) e ^{i k_r ( r- r₀ )} , (2)

we find,

$$k_r = \frac{i}{2r} \pm \sqrt{\frac{\omega^2}{v^2}-\frac{k_\theta^2}{r^2}-\frac{1}{4r^2}}$$

Thus given a wavefield at some radius r we can propagate it to a new radius using,

$$P(r+\Delta r) = P(r) e^{-\frac{\Delta r}{2r}} e^{ i \sqrt{\... ...omega^2}{v^2}- \frac{k_\theta^2}{r^2}-\frac{1}{4r^2}} \Delta r}$$

The real part of this expression controls the amplitude decay due to increasing wavefront radius and the imaginary part gives the phase shift.

The initial conditions are not supplied at r=0, but rather at some finite initial radius r₀. The region between the origin and r₀ is assumed to be homogeneous. A suitable initial condition is,

$$P(r_0,\theta,\omega) = \frac{1}{\sqrt{r_0}} \frac{1}{\sqrt{\omega}} e^{ i \omega \frac{r_0}{v} }.$$

Figure  shows the amplitude and phase plots for a mono-frequency Green's function calculated for a constant velocity medium. The polar coordinate data has been re-interpolated onto the same rectangular grid as the rectangular coordinate data. The amplitude is now more uniform for all dips. Figure  shows one time slice through the full Green's function in the (x,z,t) domain. The result is a much more accurate representation of the wavefront.

 

pol-one
Figure 3 Amplitude (left) and phase (right) for a mono-frequency Green's function calculated in polar coordinates. The initial condition is a constant value on a semi-circle of radius 80m.

 

pol-all
Figure 4 A time slice through the (x,z,t) Green's function calculated by inverse transforming a complete range of mono-frequency Green's functions calculated in polar coordinates.