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

Given a light source, a spherical reflector, and an observer, where on the surface of the sphere will the light be directly reflected to the observer, i.e. where is the the specular point? This is known as the Alhazen-Ptolemy problem, and finding this specular point for spherical reflectors is useful in applications ranging from computer rendering to atmospheric modeling to GPS communications. Existing solutions rely upon finding the roots of a quartic equation and evaluating numerically which root provides the real specular point. We offer a formulation, and two solutions thereof, for which the correct root is predeterminable, thereby allowing the construction of the fully analytical solutions we present. Being faster to compute, our solutions should prove useful in cases which require repeated calculation of the specular point, such as Monte-Carlo radiative transfer, including reflections off of Titan's hydrocarbon seas.


INTRODUCTION
The specular point is the specific location on a surface where the outgoing vector of a mirror reflection coincides exactly with the vector to an observer. At such a point, the light from a point source is seen by the observer and, if the source has spatial extent, an image is produced in the reflection. Since superposition allows a source with spatial extent to be divided into point sources -each with a single specular point which in aggregate comprise the specular region -the same methods for finding a specular point can be repeatedly applied to determine the region of specular reflection corresponding to a source with spatial extent.  Perhaps the most notable use for finding the specular point for an arbitrary source-observer configuration is in Global Positioning Systems. Unfortunately, the tolerance of GPS requires the Earth be approximated as an ellipsoid instead of a sphere, rendering analytical solution of the spherical-surface case less useful (Southwell & Dempster 2018;Prakash et al. 1994). However, our formulation can be extended to elliptical geometries more easily than previous approaches, arguing positively of its usefulness in this application.
Unlike the 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 raytracing (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.

Planar Surface
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 where θ i is incidence angle and θ e is emission angle. In the configuration shown in Figure 1 they are defined as , and 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, x 2 = x s y 2 y 2 − y 1 ,

Alhazen-Ptolemy Problem
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 11 th 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 ruler-and-compass construction of a general solution, a method which 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.

Existing methods
The solution in Elkin (1965) produces a quartic equation for the distance from either the source or the observer to the specular point (though the he does not term it such), the roots of which can be used to retrieve a collection of distances from the specular point to the other point of interest. The resulting x-y pairs must then be analyzed to find the real specular point. There are several methods which are useful in retrieving the specular point from the generated x-y pairs: minimizing the source-specular-observer path length, employing a numerical root finding method with an initial 'guess' (e.g. Newton's method), or 'by probe'. Fujimura et al. (2019) provides a formulation which works well for the method of path-length minimization, but suffers numerical instability in cases where the radius of the sphere is very small relative to source or observer distance (e.g. in the one-finite case discussed below). Glaeser (1999) emonstrates the numerical root finding approach (using an algorithm from Schwarze 1990) and notes that it is subject to "numerical instabilities that may lead to a considerable loss of accuracy" under certain circumstances. Elkin (1965) shows the 'by probe' method, but this is almost impossible to capture algorithmically. These limitations beg the creation of a better, faster method.

Our approach
The salient problem, and in fact the only obstacle to a fully analytical solution, is that of deducing which root to use. We present a formulation that allows predetermination of the 'correct' root and enables the construction of a piecewise function which directly produces the location of the physical specular point. We find this formulation superior to previous work because of its analytical branch deduction -something lacking from every other formulation. The fully analytical approach avoids numerical instability and allows significantly faster computation.
We present two analytical solutions for this formulation: one that requires the distance to either the source or observer be approximateable as infinite (the "one-finite" case) and a second that works for any arbitrary configuration of source and observer (the "both-finite" case). Furthermore, we demonstrate that the complete solution can be expressed using no more than 3 roots in the both-finite case and no more than 2 roots in the one-finite case. We apply an Euler rotation and make use of the particulars of the resulting geometry to simplify the analysis, which also serves to minimize the associated computational expenses.

EULER ROTATION
In both cases, it is useful to place either the source or observer at polar angle π /2 (i.e. above the north pole) for two reasons: first, it reduces the subsequent analysis to two dimensions and second, it allows a simplification which will be used in Sections 3 and 4. The simplest way to accomplish this is to start in Cartesian coordinates with the origin at the center of the sphere and then apply two sequential rotations. The first rotation uses the angles to define new coordinates in a rotated frame (denoted by the prime, ), y k = x k cos (α) + y k sin (α) , and (5b) Here α is the azimuthal angle of the source (measured from the positive x-axis), and β is the inclination angle (measured from the positive z-axis). The second rotation is not strictly necessary; the critical portion is the placement of source or observer at π /2, accomplished by the first rotation. However, this second rotation further reduces the problem because it enforces − π /2 < θ obs < π /2 (in the rotated frame), allowing us to ignore some of the roots in Equations 13 and 20. This rotation takes γ = atan2 (y src , x src ) to establish coordinates in the twice-rotated, , frame to be With k ∈ {src, obs} (i.e. for x src , x obs , etc) and where atan2 (y, x) defines the angle from the positive x-axis to (x, y). At the end, once the specular point has been calculated in this new coordinate system, we will apply these two rotation in reverse for k ∈ {src, obs, spec} to recover the specular point in the original system. The rotated system (x k , y k , z k ) will be used for the remainder of this paper.
In the one-finite case, the rotation should be performed to place the coordinate of finite radius at π /2. This allows full use to be made of the approximation R obs R sph or R src R sph . In the both-finite case the rotation can place either coordinate at π /2, but we will place the source at π /2 for notational consistency.

ONE-FINITE CASE
If either the source or the observer can be considered infinitely distant, then after applying the Euler rotation in Equation 5 to place the finitely distant of the two at π /2 -as Figure 2 shows -θ 1 and θ 2 can be defined such that where θ e and θ i are then given by

Formulation
Since incidence angle (θ i ) must always be equal to emission angle (θ e ), then in the limit where R obs R sph we get With the introduction of a constant, and substituting in θ 2 from Equation 8b, Equation 9 becomes Making use of the identities tan (x − y) = sin (2x) − sin (2y) cos (2x) + cos (2y) and sin x cos 2y ± sin 2y cos x = sin (2y ± x) , and the fact that the Euler rotation in Equation 5 results in θ src = π /2, we get cos (4θ spec ) + c sin (3θ spec ) + cos (2θ obs ) = c sin (θ spec − 2θ obs ) .

Solution
Equation 12 has sixteen unique, non-trivial solutions, courtesy of Wolfram Mathematica (Wolfram Research 2020). We neglect 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 13. D 0 ≡ 1152 c 2 − 4 c 2 + c 2 cos (2θ obs ) − 846 cos (2θ obs ) − 1 + c 2 sin (2θ obs ) 2 (13a) The only differences between the solutions are the different signs associated with 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 observer at 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 4. BOTH-FINITE CASE
With b = 0 (and with some manipulation) we can recover Equation 12. Doing so is nontrivial and will be neglected for brevity, but it is clearly demonstrated by the agreement between the two methods when b ≈ 0 as shown in section 5.2. Solutions to cos(4θ spec ) + c sin(3θ spec ) + cos(2θ obs ) = c sin(θ spec − 2θ obs ) Figure 3. Plots of the solutions to Equation 12 where R sph /Rsrc (c) is 0.85. Note that the numerical method is only applied to [− π /2, π /2] because the Euler rotation in Section 2 enforces θ obs to be between − π /2 and π /2. θ obs is the angle of the observer, θspec is the angle of the specular point. Both are in the rotated reference frame and measured from the positive x-axis. An animated version is available which shows the solutions for various values of c.

Similar to Equation 12, Equation 19
produces sixteen unique, non-trivial solutions -only eight of which are physically meaningful (the others represent special cases). These are presented in Equation 21 using the coefficients defined in Equation 20 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 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

Solutions to
c sin(θ obs +θ spec )−cb sin(2θ spec )+cos(θ obs )−b cos(θ spec ) c cos(θ obs +θ spec )−cb cos(2θ spec )+sin(θ obs )−b sin(θ spec ) = tan(2θ spec ) Figure 4. Solutions to Equation 21 for R sph /Rsrc (c) of 0.85 and R sph /R obs (b) of 0.95. Note that the numerical method is only applied to [− π /2, π /2] because the Euler rotation in Section 2 enforces θ obs to be between − π /2 and π /2. θ obs is the angle of the observer, θspec is the angle of the specular point. Both are in the rotated reference frame and measured from the positive x-axis. An animated version is available which shows the solutions for various values of b and c.
where the limits L 1 and L 2 are given by It should be noted that L 1 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 7 , in practical computation it is robust to use the approximation where ζ is directly related to the branch-deduction precision (assuming it is larger than the computation precision), e.g. ζ = 10 −9 will result in the correct branch being chosen within 10 −9 . Our both-finite implementation in the following section makes use of this approach.

COMPUTATION
Implementations of these solutions in both Python and C++ can be found in this GitHub Repository 1 . Also in the repository are Wolfram Language Package files for easily reading the solutions into Mathematica. The C++ implementations were used for performance benchmarking. The solution for the one-finite case performs significantly better than that of the both-finite case, but in many applications the difference will be negligible because both routines are relatively minimal -especially compared to numerical methods.

Performance
Excluding branch deductions (and after compiling with 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 an 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 9 random combinations of b, c, and θ obs . The code used for benchmarking is included in the linked GitHub repository. In cases where L 1 ≤ θ obs ≤ L 2 the one-finite implementation has a more significant speed advantage because the both-finite implementation's branch deductions are most expensive in this region.

Precision
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 −7 due to the trigonometric functions. The use of the long double type (which is implementation dependent but usually corresponds to between 80-bit and 128-bit storage) can be used to improve this precision to 10 −15 . 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 −6 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 −5 , 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 −2 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 −6 degrees at b /c = 500. 6. CONCLUSION

SRTC++
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 Absolute error (radians) Absolute error of cos(4θ spec ) + c sin(3θ spec ) + cos(2θ obs ) = c sin(θ spec − 2θ obs ) θ obs 0 π 8 π 2 3π 4 π 10 −3 10 −4 10 −6 Maximum error Figure 5. The agreement between the one-finite and both-finite solutions as a function of the ratio between b and c with the maximum error across all θ obs for a given b/c. Note that even though θ obs is confined to [− π /2, π /2], we are showing the error for θ obs between 0 and π. We do so to demonstrate that for − π /2 > θ obs or θ obs < π /2 (in the reference frame before the second Euler rotation) the persistent error is symmetrically reversed.
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 compensation 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.

Computer rendering
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.

Extension to elliptical geometry
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 non-circular conic sections. We leave the elliptical case to future work.

In Summary
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 3 branches in the both-finite case and only 2 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.