Phase-ray formulation

When a solution, $\Psi$, to the Helmholtz equation is known, obtaining a ray trajectory using equation (5) is relatively trivial. The expression for the wavefield gradient, $\nabla \Psi=\nabla(A {\rm e}^{{\rm i}\phi})$, divided by the wavefield, $\Psi$, is,  

$$\frac{\nabla \Psi}{\Psi} = \frac{\nabla A}{A} + {\rm i} {\bf K}.$$ (7)

An expression for the wavefield gradient vector, ${\bf K}$, is obtained by retaining the imaginary component of equation (7) and using the expression for ${\bf K}$ in equation (3),  

$$K \frac{ {\rm d}}{ {\rm d}s}{\bf r} = {\rm Im} \left( \frac{\nabla \Psi} {\Psi} \right)$$ (8)

Equation (8) may be rewritten explicitly as a system of two decoupled ordinary differential equations,  

$$\frac{ {\rm d}}{ {\rm d}s}\left[\begin{array}[pos] {c} x \\ ... ...\frac{\partial}{\partial z} \\ \end{array}\right]\Psi \right).$$ (9)

The solution for ray-path, ${\bf r}$, is computed through an initial evaluation of the right hand side of equations (9), and an iterative forward step by a constant interval using the precomputed quantity to determine the proper apportioning of the step along each coordinate. Note that because these differential equations are first-order, only one initial condition (position) is required for ray computation, and rays may be started from any location in the wavefield solution. Finally, although the two-dimensional formulation is presented here, the extension to three dimensions is trivial.