We apply the level set method to compute the three dimensional multivalued geometrical optics term
in a paraxial formulation.
The paraxial formulation is obtained from the 3-D stationary eikonal equation
by using one of the spatial directions as the arti.cial evolution direction.
The advection velocity field used to move level sets is obtained by the method of characteristics;
therefore the motion of level sets is defined in phase space.
The multivalued travel-time and amplitude-related quantity are obtained
from solving advection equations with source terms.
We derive an amplitude formula in a reduced phase space which is very convenient
to use in the level set framework.
By using a semi-Lagrangian method in the paraxial formulation,
the method has O(N2) rather than O(N4) memory storage requirement
for up to O(N2) multiple point sources in the five dimensional phase space,
where N is the number of mesh points along one direction.
Although the computational complexity is still O(MN4),
where M is the number of steps in the ODE solver for the semi-Lagrangian scheme,
this disadvantage is largely overcome by the fact that up to O(N2)
multiple point sources can be treated simultaneously.
Three dimensional numerical examples demonstrate the efficiency and accuracy of the method.