A system of ray tracing equations

The ray tracing method I have implemented is based on the ray tracing system of first-order partial-differential equations, derived from the Eikonal equation by the method of characteristics (see Cerveny1987 for more details):

$$\begin{eqnarray} {dx_i \over ds } & = & {p_i \over w} \nonumber \\ {dp_i \over ds } & = & {\partial w \over \partial x_i }\end{eqnarray}$$ (37)

Here w is the slowness, x_i=x_i(s) (i=1,2,3) are the coordinates of the ray as a function of s, the arclength along the ray, and p_i(s) are the components of the slowness vector ${\bf \bar p}$.The traveltime t(s) along the ray is given by  

$${dt \over ds } = w.$$ (38)

The slowness vector is perpendicular to the wavefront (${\bf \bar p} = \nabla t$), and must satisfy

$$\vert{\bf \bar p}\vert^2 = \sum^3_{i=3}p_i^2 = w^2.$$ (39)

The system of ray equations together with equation  () can be solved by a standard numerical integration method. I use a fourth-order Rungge-Kutta method. It propagates the properties of the ray (x_i,p_i, and t) over an increment in arclength by combining the information from several Euler-style steps (each involving one evaluation of the right-hand side of equations () and ()), and then uses the information obtained to match a fourth-order Taylor series expansion of the ray variables at the current position.

Figure  shows some of rays traced through the reference model shown in Figure .

 

ray0
Figure 15 Some of the ray trajectories traced through the reference medium to calculate the operator L.