Solving the Alhazen–Ptolemy Problem: Determining Specular Points on Spherical Surfaces for Radiative Transfer of Titan's Seas

Perhaps the most notable use for finding the specular point for an arbitrary source–observer configuration is in Global Positioning Systems (GPS). Unfortunately, the tolerance of GPS requires Earth to be approximated as an ellipsoid instead of a sphere, rendering the analytical solution of the spherical-surface case less useful (Prakash et al. 1994; Southwell & Dempster 2018). However, our formulation can be extended to elliptical geometries more easily than previous approaches, arguing positively of its usefulness in this application.

Unlike Earth, Saturn's moon Titan—which harbors liquid hydrocarbons that allowed the Cassini mission to record numerous specular reflections (Stephan et al. 2010; Barnes et al. 2014)—can be well approximated as a sphere. The SRTC++ model (Barnes et al. 2018) does exactly that, and finding a faster method for computing the specular reflections due to scattering off of Titan's ubiquitous atmospheric haze aerosols in SRTC++ is the direct motivation for this work. There may also be utility for removing "Sun glint" from imaging data by employing a method similar to those in Kay et al. (2009), but with the higher precision enabled by our analytical solutions.

Another notable application is the rendering of computer graphics with ray tracing (Inakage 1986), a case where computation speed is critical, particularly when rendering a reflected image (instead of a simple point source) because the rendering fidelity is much more sensitive to the accuracy of the specular point calculation. The increased interest in ray-traced rendering, particularly with regard to the video game and movie industries, serves only to increase the usefulness of a faster method for handling even partially reflective surfaces.

For a planar surface, the specular point (Figure 1) can be found trivially. In the coordinate system centered on the specular point with the x-axis parallel to the surface, the positions of the source and observer must satisfy

$$\begin{eqnarray}&&{\theta }_{i}\mathop{=}\limits^{!}{\theta }_{e},\end{eqnarray} \tag{ 1 }$$

where θ_i is the incidence angle and θ_e is the emission angle. In the configuration shown in Figure 1, they are defined as

$$\begin{eqnarray*}&&{\theta }_{i}={\tan }^{-1}\left(\displaystyle \frac{{y}_{1}}{{x}_{1}}\right),\ \ \ \ \ \mathrm{for}\ \ \ \ \ 0\leqslant {\theta }_{i}\leqslant \displaystyle \frac{\pi }{2},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{\theta }_{e}={\tan }^{-1}\left(\displaystyle \frac{{y}_{2}}{{x}_{2}}\right),\ \ \ \ \ \mathrm{for}\ \ \ \ \ 0\leqslant {\theta }_{e}\leqslant \displaystyle \frac{\pi }{2},\ \ \ \ \ \mathrm{and}\end{eqnarray*}$$

$$\begin{eqnarray*}&&{x}_{s}={x}_{1}+{x}_{2},\end{eqnarray*}$$

where x_s is the separation between the source and observer in the x direction. Then, the position of the source and observer relative to the specular point can be found by combining the above and rearranging:

$$\begin{eqnarray*}&&{x}_{2}=\displaystyle \frac{{x}_{s}{y}_{2}}{{y}_{2}-{y}_{1}},\end{eqnarray*}$$

$$\begin{eqnarray*}&&{x}_{1}=\displaystyle \frac{{x}_{s}{y}_{1}}{{y}_{1}-{y}_{2}}.\end{eqnarray*}$$

Zoom In Zoom Out Reset image size

The case of a spherical surface (Figure 2) is far more complicated. First formulated in 150 CE by Ptolemy (Weisstein 2002), it is known as the Alhazen–Ptolemy problem because the first method of solving it (a geometric approach using conic sections) was provided—and proven—by Ibn al-Haythum in the 11th century (Katz 1995). Elkin (1965) provided the first algebraic solution and similar formulations by Riede (1989), Smith (1992), and Waldvogel (1992) confirmed Elkin's solution. Neumann (1998) proved that the ruler-and-compass construction of a general solution, a method that permits only circular arcs and straight lines, is impossible. That is, given only an unmarked and idealized compass and straight edge, there is no way to determine the specular point.

Zoom In Zoom Out Reset image size

1.3.1. Existing Methods

1.3.2. Our Approach

Because the incidence angle (θ_i ) must always be equal to emission angle (θ_e ), then in the limit where R_obs ≫ R_sph, we get

$$\begin{eqnarray}&&{\theta }_{2}=2{\theta }_{\mathrm{spec}}-{\theta }_{\mathrm{obs}}.\end{eqnarray} \tag{ 9 }$$

With the introduction of a constant,

$$\begin{eqnarray}&&c\equiv \displaystyle \frac{{R}_{\mathrm{sph}}}{{R}_{\mathrm{src}}},\end{eqnarray} \tag{ 10 }$$

and substituting in θ₂ from Equation (8b), Equation (9) becomes

$$\begin{eqnarray}&&\displaystyle \frac{\sin \left({\theta }_{\mathrm{src}}\right)-c\sin \left({\theta }_{\mathrm{spec}}\right)}{\cos \left({\theta }_{\mathrm{src}}\right)-c\cos \left({\theta }_{\mathrm{spec}}\right)}=\tan \left(2{\theta }_{\mathrm{spec}}-{\theta }_{\mathrm{obs}}\right).\end{eqnarray} \tag{ 11 }$$

Making use of the identities

$$\begin{eqnarray*}&&\begin{array}{l}\tan \left(x-y\right)=\displaystyle \frac{\sin \left(2x\right)-\sin \left(2y\right)}{\cos \left(2x\right)+\cos \left(2y\right)}\qquad \mathrm{and}\\ \sin x\cos 2y\pm \sin 2y\cos x=\sin \left(2y\pm x\right)\end{array}\end{eqnarray*}$$

and the fact that the Euler rotation in Equation (5a) results in θ_src = π/2, we get

$$\begin{eqnarray}&&\,\begin{array}{l}\cos \left(4{\theta }_{\mathrm{spec}}\right)+c\sin \left(3{\theta }_{\mathrm{spec}}\right)\\ \quad +\,\cos \left(2{\theta }_{\mathrm{obs}}\right)=c\sin \left({\theta }_{\mathrm{spec}}-2{\theta }_{\mathrm{obs}}\right).\end{array}\,\end{eqnarray} \tag{ 12 }$$

Equation (12) has 16 unique, nontrivial solutions, courtesy of Wolfram Mathematica (Wolfram Research 2020). We ignore eight of these because they represent special cases that are not physically meaningful for our purposes. The other eight are presented in Equation (14) using the coefficients we define in Equation (13a):

$$\begin{eqnarray}\,\begin{array}{rcl}{D}_{0} & \equiv & 1152\left({c}^{2}-4\right)({c}^{2}+{c}^{2}\cos \left(2{\theta }_{\mathrm{obs}}\right)\\ & & -\,846\cos \left(2{\theta }_{\mathrm{obs}}\right)-1)+{c}^{2}\sin {\left(2{\theta }_{\mathrm{obs}}\right)}^{2}\end{array}\,\end{eqnarray} \tag{ 13a }$$

$$\begin{eqnarray}&&\,\begin{array}{l}{D}_{1}\equiv \left(-4\left(160-128{c}^{2}+4{c}^{4}+96\cos \left(2{\theta }_{\mathrm{obs}}\right)\right.\right.\\ \ {\left.{\left.-96{c}^{2}\cos \left(2{\theta }_{\mathrm{obs}}\right)\right)}^{3}+{\left(-16{\left({c}^{2}-4\right)}^{3}-{D}_{0}\right)}^{2}\right)}^{1/2}\end{array}\,\end{eqnarray} \tag{ 13b }$$

$$\begin{eqnarray}&&{D}_{2}\equiv \sqrt[3]{16{\left({c}^{2}-4\right)}^{3}+{D}_{0}-{D}_{1}}\end{eqnarray} \tag{ 13c }$$

$$\begin{eqnarray}&&\,\begin{array}{l}{D}_{3}\equiv \displaystyle \frac{1}{3\sqrt[3]{4}{D}_{2}}\left(40-32{c}^{2}+{c}^{4}+24\cos \left(2{\theta }_{\mathrm{obs}}\right)\right.\\ \quad -\,\left.24{c}^{2}\cos \left(2{\theta }_{\mathrm{obs}}\right)\right)+\displaystyle \frac{{D}_{2}}{24\sqrt[3]{2}}\end{array}\,\end{eqnarray} \tag{ 13d }$$

$$\begin{eqnarray}&&{D}_{4}\equiv \sqrt{\displaystyle \frac{4-{c}^{2}}{6}+{D}_{3}}\end{eqnarray} \tag{ 13e }$$

$$\begin{eqnarray}&&{\theta }_{\mathrm{spec}}=\pm \left\{\begin{array}{c}{\cos }^{-1}\left(-\tfrac{{D}_{4}}{2}-\sqrt{-\tfrac{4-{c}^{2}}{3}-{D}_{3}-\tfrac{c\sin \left(2{\theta }_{\mathrm{obs}}\right)}{8{D}_{4}}}\right)\\ {\cos }^{-1}\left(-\displaystyle \frac{{D}_{4}}{2}+\sqrt{-\tfrac{4-{c}^{2}}{3}-{D}_{3}-\tfrac{c\sin \left(2{\theta }_{\mathrm{obs}}\right)}{8{D}_{4}}}\right)\\ {\cos }^{-1}\left(+\displaystyle \frac{{D}_{4}}{2}-\sqrt{-\tfrac{4-{c}^{2}}{3}-{D}_{3}+\tfrac{c\sin \left(2{\theta }_{\mathrm{obs}}\right)}{8{D}_{4}}}\right)\\ {\cos }^{-1}\left(+\displaystyle \frac{{D}_{4}}{2}+\sqrt{-\tfrac{4-{c}^{2}}{3}-{D}_{3}+\tfrac{c\sin \left(2{\theta }_{\mathrm{obs}}\right)}{8{D}_{4}}}\right)\end{array}\right.\end{eqnarray} \tag{ 14 }$$

The only differences between the solutions are the different signs associated with the four terms. The possible combinations are ± − − − , ± − + − , ± + − + , and ± +++. These solutions are plotted in Figure 3 for a source at 3000 km from the center of a sphere of radius 2575 km (i.e., Titan) and an observer at a much greater distance. After comparing these solutions against the numerical method, as in Figure 3, it is clear the branch that produces the physical specular point (on the exterior of the sphere) is

$$\begin{eqnarray}{\theta }_{\mathrm{spec}}=\left\{\begin{array}{ll}{\cos }^{-1}\left(\tfrac{{D}_{4}}{2}-\sqrt{-\tfrac{4-{c}^{2}}{3}-{D}_{3}+\tfrac{c\sin \left(2{\theta }_{\mathrm{obs}}\right)}{8{D}_{4}}}\right) & \ \ \ \ \ {\theta }_{\mathrm{spec}}\lt \tfrac{\pi ({\cot }^{-1}\left(c\right)-{\theta }_{\mathrm{obs}})}{\pi +2{\cot }^{-1}\left(c\right)}\\ {\cos }^{-1}\left(\tfrac{{D}_{4}}{2}+\sqrt{-\tfrac{4-{c}^{2}}{3}-{D}_{3}+\tfrac{c\sin \left(2{\theta }_{\mathrm{obs}}\right)}{8{D}_{4}}}\right) & \ \ \ \ \ {\theta }_{\mathrm{spec}}\geqslant \tfrac{\pi ({\cot }^{-1}\left(c\right)-{\theta }_{\mathrm{obs}})}{\pi +2{\cot }^{-1}\left(c\right)}\end{array}\right..\end{eqnarray} \tag{ 15 }$$

If ${R}_{\mathrm{obs}}\,\rlap{/}{\gg }\,{R}_{\mathrm{sph}}$ and ${R}_{\mathrm{src}}\,\rlap{/}{\gg }\,{R}_{\mathrm{sph}}$, then the solution becomes even more complicated. In this case, Equation (9) becomes

$$\begin{eqnarray}&&{\theta }_{2}+{\theta }_{1}=2{\theta }_{\mathrm{spec}},\end{eqnarray} \tag{ 16 }$$

and Equation (8a) becomes

$$\begin{eqnarray}&&{\theta }_{1}={\tan }^{-1}\left(\displaystyle \frac{{R}_{\mathrm{obs}}\sin \left({\theta }_{\mathrm{src}}\right)-{R}_{\mathrm{sph}}\sin \left({\theta }_{\mathrm{spec}}\right)}{{R}_{\mathrm{obs}}\cos \left({\theta }_{\mathrm{src}}\right)-{R}_{\mathrm{sph}}\cos \left({\theta }_{\mathrm{spec}}\right)}\right)\ \ \ \ \ \mathrm{and}\end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray}&&{\theta }_{2}={\tan }^{-1}\left(\displaystyle \frac{{R}_{\mathrm{src}}\sin \left({\theta }_{\mathrm{obs}}\right)-{R}_{\mathrm{sph}}\sin \left({\theta }_{\mathrm{spec}}\right)}{{R}_{\mathrm{src}}\cos \left({\theta }_{\mathrm{obs}}\right)-{R}_{\mathrm{obs}}\cos \left({\theta }_{\mathrm{spec}}\right)}\right),\ \ \ \ \ \end{eqnarray} \tag{ 17b }$$

which yields

$$\begin{eqnarray*}&&\begin{array}{l}{\tan }^{-1}\left(\displaystyle \frac{{R}_{\mathrm{obs}}\sin \left({\theta }_{\mathrm{src}}\right)-{R}_{\mathrm{sph}}\sin \left({\theta }_{\mathrm{spec}}\right)}{{R}_{\mathrm{obs}}\cos \left({\theta }_{\mathrm{src}}\right)-{R}_{\mathrm{sph}}\cos \left({\theta }_{\mathrm{spec}}\right)}\right)\\ \ \ +\,{\tan }^{-1}\left(\displaystyle \frac{{R}_{\mathrm{src}}\sin \left({\theta }_{\mathrm{obs}}\right)-{R}_{\mathrm{sph}}\sin \left({\theta }_{\mathrm{spec}}\right)}{{R}_{\mathrm{src}}\cos \left({\theta }_{\mathrm{obs}}\right)-{R}_{\mathrm{obs}}\cos \left({\theta }_{\mathrm{spec}}\right)}\right)=2{\theta }_{\mathrm{spec}}\,.\end{array}\end{eqnarray*}$$

Again making use of θ_src = π/2 and defining a new constant,

$$\begin{eqnarray}&&b\equiv \displaystyle \frac{{R}_{\mathrm{sph}}}{{R}_{\mathrm{obs}}},\end{eqnarray} \tag{ 18 }$$

the both-finite equivalent of Equation (12) is

$$\begin{eqnarray}&&\displaystyle \frac{c\sin \left({\theta }_{\mathrm{obs}}+{\theta }_{\mathrm{spec}}\right)-{cb}\sin \left(2{\theta }_{\mathrm{spec}}\right)+\cos {\theta }_{\mathrm{obs}}-b\cos {\theta }_{\mathrm{spec}}}{c\cos \left({\theta }_{\mathrm{obs}}+{\theta }_{\mathrm{spec}}\right)-{cb}\cos \left(2{\theta }_{\mathrm{spec}}\right)+\sin {\theta }_{\mathrm{obs}}-b\sin {\theta }_{\mathrm{spec}}}=\tan \left(2{\theta }_{\mathrm{spec}}\right)\,.\end{eqnarray} \tag{ 19 }$$

With b = 0 (and with some manipulation), we can recover Equation (12). Doing so is nontrivial and will be ignored for brevity, but it is clearly demonstrated by the agreement between the two methods when b ≈ 0 as shown in Section 5.2.

Similar to Equation (12), Equation (19) produces 16 unique, nontrivial solutions—only 8 of which are physically meaningful (the others represent special cases). These are presented in Equation (21) using the coefficients defined in Equation (20a):

$$\begin{eqnarray}&&{E}_{0}\equiv {b}^{2}+{c}^{2}+2{bc}\sin {\theta }_{\mathrm{obs}}-4\end{eqnarray} \tag{ 20a }$$

$$\begin{eqnarray}&&\,\begin{array}{l}{E}_{1}\equiv 24\cos {\left({\theta }_{\mathrm{obs}}\right)}^{2}\\ \quad \times \,\left({b}^{2}-{bc}\sin {\theta }_{\mathrm{obs}}-2{c}^{2}+2\right)+{E}_{0}^{2}\end{array}\,\end{eqnarray} \tag{ 20b }$$

$$\begin{eqnarray}&&\,\begin{array}{l}{E}_{2}\equiv 72{E}_{0}{\cos }^{2}{\theta }_{\mathrm{obs}}\left({b}^{2}-{bc}\sin {\theta }_{\mathrm{obs}}+4{c}^{2}-4\right)\\ \quad +\,432{\cos }^{2}{\theta }_{\mathrm{obs}}\left({(b-c\sin {\theta }_{\mathrm{obs}})}^{2}\right.\\ \quad \left.-\,{b}^{2}({c}^{2}-1){\cos }^{2}{\theta }_{\mathrm{obs}}\right)+2{E}_{0}^{3}\end{array}\,\end{eqnarray} \tag{ 20c }$$

$$\begin{eqnarray}&&{E}_{3}\equiv {\left(\sqrt{{E}_{2}^{2}-4{E}_{1}^{3}}+{E}_{2}\right)}^{1/3}\end{eqnarray} \tag{ 20d }$$

$$\begin{eqnarray}&&{E}_{4}\equiv \displaystyle \frac{\sqrt[3]{4}\,{E}_{3}^{2}+\sqrt[3]{16}\,{E}_{1}}{24{E}_{3}}\end{eqnarray} \tag{ 20e }$$

$$\begin{eqnarray}&&{E}_{5}\equiv \sqrt{\displaystyle \frac{{b}^{2}\cos {\left({\theta }_{\mathrm{obs}}\right)}^{2}}{4}+\displaystyle \frac{{E}_{3}}{12\sqrt[3]{2}}+\displaystyle \frac{{E}_{1}}{6\sqrt[3]{4}\ {E}_{3}}-\displaystyle \frac{{E}_{0}}{6}}\end{eqnarray} \tag{ 20f }$$

$$\begin{eqnarray}&&{E}_{6}\equiv \displaystyle \frac{{(b\cos {\theta }_{\mathrm{obs}})}^{3}-{{bE}}_{0}\cos {\theta }_{\mathrm{obs}}-4\cos {\theta }_{\mathrm{obs}}\left(b-c\sin {\theta }_{\mathrm{obs}}\right)}{16{E}_{5}}\end{eqnarray} \tag{ 20g }$$

$$\begin{eqnarray}&&\,{\theta }_{\mathrm{spec}}=\pm \left\{\begin{array}{l}{\cos }^{-1}\left(\tfrac{b\cos {\theta }_{\mathrm{obs}}-2{E}_{5}}{4}-\tfrac{1}{2}\sqrt{\tfrac{3{(b\cos {\theta }_{\mathrm{obs}})}^{2}-2{E}_{0}}{12}-{E}_{4}-{E}_{6}}\right)\\ {\cos }^{-1}\left(\tfrac{b\cos {\theta }_{\mathrm{obs}}-2{E}_{5}}{4}+\tfrac{1}{2}\sqrt{\tfrac{3{(b\cos {\theta }_{\mathrm{obs}})}^{2}-2{E}_{0}}{12}-{E}_{4}-{E}_{6}}\right)\\ {\cos }^{-1}\left(\tfrac{b\cos {\theta }_{\mathrm{obs}}+2{E}_{5}}{4}-\tfrac{1}{2}\sqrt{\tfrac{3{(b\cos {\theta }_{\mathrm{obs}})}^{2}-2{E}_{0}}{12}-{E}_{4}+{E}_{6}}\right)\\ {\cos }^{-1}\left(\tfrac{b\cos {\theta }_{\mathrm{obs}}+2{E}_{5}}{4}+\tfrac{1}{2}\sqrt{\tfrac{3{(b\cos {\theta }_{\mathrm{obs}})}^{2}-2{E}_{0}}{12}-{E}_{4}+{E}_{6}}\right)\end{array}\right..\end{eqnarray} \tag{ 21 }$$

As in Equation (14), the only differences between the solutions are the signs on four terms. The combinations of these signs are ± − − − , ± − + − , ± + − + , and ± +++ (the same as in Equation (14)). They are plotted in Figure 4. Comparison with the numerical method indicates the branch that produces the physical specular point (on the exterior of the sphere) is

$$\begin{eqnarray*}&&{E}_{7}\equiv \displaystyle \frac{3{(b\cos {\theta }_{\mathrm{obs}})}^{2}-2{E}_{0}}{12}-{E}_{4}+{E}_{6}\end{eqnarray*}$$

$$\begin{eqnarray}{\theta }_{\mathrm{spec}}=\left\{\begin{array}{ll}-{\cos }^{-1}\left(\tfrac{b\cos {\theta }_{\mathrm{obs}}+2{E}_{5}}{4}+\tfrac{1}{2}\sqrt{{E}_{7}}\right) & \ \ \ \ \ -\tfrac{\pi }{2}\leqslant {\theta }_{\mathrm{obs}}\lt {L}_{1}\\ {\cos }^{-1}\left(\tfrac{b\cos {\theta }_{\mathrm{obs}}+2{E}_{5}}{4}+\tfrac{1}{2}\sqrt{{E}_{7}}\right) & \ \ \ \ \ {L}_{1}\leqslant {\theta }_{\mathrm{obs}}\leqslant {L}_{2}\ \ \mathrm{and}\ \ \tfrac{\partial \,{E}_{7}}{\partial \,{\theta }_{\mathrm{obs}}}\leqslant 0\\ {\cos }^{-1}\left(\tfrac{b\cos {\theta }_{\mathrm{obs}}+2{E}_{5}}{4}-\tfrac{1}{2}\sqrt{{E}_{7}}\right) & \ \ \ \ \ {L}_{2}\lt {\theta }_{\mathrm{obs}}\leqslant \tfrac{\pi }{2}\ \ \mathrm{or}\ \ \tfrac{\partial \,{E}_{7}}{\partial \,{\theta }_{\mathrm{obs}}}\gt 0\end{array}\right.,\end{eqnarray} \tag{ 22 }$$

where the limits L₁ and L₂ are given by

$$\begin{eqnarray}&&{L}_{1}={\tan }^{-1}\left(\displaystyle \frac{c+b\sqrt{1-{b}^{2}+{c}^{2}}}{{b}^{2}-{c}^{2}}\right)\end{eqnarray} \tag{ 23a }$$

$$\begin{eqnarray}&&{L}_{2}={\tan }^{-1}\left(\displaystyle \frac{c-b\sqrt{1-{b}^{2}+{c}^{2}}}{{b}^{2}-{c}^{2}}\right).\end{eqnarray} \tag{ 23b }$$

It should be noted that L₁ may occur at less than −π/2, and in such cases, only the second two parts in Equation (22) are necessary. Additionally, while there is an exact solution for the partial derivative of E₇, in practical computation, it is robust to use the approximation

$$\begin{eqnarray*}&&\displaystyle \frac{\partial \,{E}_{7}}{\partial \,{\theta }_{\mathrm{obs}}}\approx {E}_{7}{\left|{}_{{\theta }_{\mathrm{obs}}}-{E}_{7}\right|}_{{\theta }_{\mathrm{obs}}-\zeta },\end{eqnarray*}$$

where ζ is directly related to the branch-deduction precision (assuming it is larger than the computation precision), e.g., ζ = 10⁻⁹ will result in the correct branch being chosen within 10⁻⁹. Our both-finite implementation in the following section makes use of this approach.

Excluding branch deductions (and after compiling with the gcc option -O3), our C++ implementation of the one-finite solution requires 76 floating-point operations, 6 trigonometric, and 5 sqrt calls with one expensive exponentiation operation to find a complex cube root. The corresponding both-finite code requires 131 floating-point operations, 5 trigonometric, and 5 sqrt calls, and one expensive cube root.

Benchmarked on a hexacore Intel i7-8750H, the one-finite solution averaged 23 ns per iteration and the both-finite solution averaged 121 ns per iteration for 10⁹ random combinations of b, c, and θ_obs. The code used for benchmarking is included in the linked GitHub repository. In cases where L₁ ≤ θ_obs ≤ L₂, the one-finite implementation has a more significant speed advantage because the both-finite implementation's branch deductions are most expensive in this region.

The agreement between the one-finite solution and the both-finite solution is shown in Figure 5. When using the 64 bit double type, the machine precision for either calculation is roughly 10⁻⁷ due to the trigonometric functions. The use of the long double type (which is implementation dependent but usually corresponds to between 80 and 128 bit storage) can be used to improve this precision to 10⁻¹⁵. It should be noted that any numerical methods will have the same or similar precision limitations. Figure 5 does not make use of the extra precision, and therefore, the point at which the average error exceeds 10⁻⁶ can be considered the point at which the disagreement between the two methods exceeds the machine precision. This occurs at around b/c = 2 × 10⁻⁵, i.e., when either R_obs or R_src is 50,000 times larger than the other. At a ratio of 500, the accuracy drops to ±1.72 × 10⁻² degrees, but only for the nearly worst-case-scenario of c = 0.95. As c decreases, the accuracy of θ_spec at the same b/c improves, e.g., for c = 0.1 the accuracy is ±1.72 × 10⁻⁶ degrees at b/c = 500.

Zoom In Zoom Out Reset image size

Our method is of particular use in simulating radiative transfer on Titan, i.e., in the SRTC++ model (Barnes et al. 2018). Specifically, calculating specular reflections due to atmospheric scattering allows us to compensate for the adjacency effect—when atmospheric scattering redirects light reflected off bright material adjacent to spectrally dark regions (e.g., lakes) and causes them to appear brighter than they should. This effect has been well documented on Earth (Odell & Weinman 1975; Tanre et al. 1981; Minomura et al. 2001; Sterckx et al. 2011; Kiselev et al. 2015), and Karkoschka & Schroder (2016) found it necessary to compensate for it in their reflectivity analysis of Huygens data.

Realistically simulating Titan's atmosphere and surface requires compensating for the adjacency effect, and many high-fidelity simulations. In such simulations, the specular point must be calculated for every scattering event (and there are usually billions). The effect of our method on the computation speed of these simulations is negligible (e.g., about 1 minute over 4 days of execution time). The same cannot be said of existing methods.

Another notable application is in rendering computer graphics. Computation speed is critical here as well, and often multiple light sources or reflections must be accounted for—each with a separate specular point. Sources with non-negligible spatial extent are also common, requiring the calculation of multiple specular points during each refresh. We have formulated this paper to emphasize spherical surfaces, but in fact, the analysis in the rotated reference frame can be applied to any reflective surfaces of revolution, provided the intersection between the surface and the plane formed by the source, observer, and specular point forms a circle.

Extending our approach to elliptical geometry would drastically broaden its usefulness by providing an analytical, branch-deducing method for calculating specular points on any closed conic section. As mentioned previously, this extension would enable our solution to be used in GPS communications. But it would also expand the usefulness for computer rendering because most planar slices of regular surfaces of revolution produce one or several noncircular conic sections. We leave the elliptical case to future work.

We have presented a formulation of the Alhazen–Ptolemy problem in which the physical specular point is predeterminable; we find this formulation to be superior to previous approaches because it does not rely on numerical methods for branch deduction. Furthermore, we have shown that the solution can be directly written as a piecewise function with only three branches in the both-finite case and only two branches in the one-finite case. Employing an Euler rotation, normalizing the radius of the sphere to unity, and translating the rotated coordinate system to center on the specular point are critical in this formulation and in enabling robust and efficient computational implementation. Our method is useful for studying the ethane composition of Titan's lakes with SRTC++, rendering 3D computer graphics, and lays the groundwork for extension to elliptical geometry.