Find all 4 possible normals to an ellipse

Question

Given a point p exterior to an axially aligned, origin centered ellipse E, find the (upto) four unique normals to E passing through p.

This is not a Mathematica question. Direct computation is too slow; I am willing to sacrifice precision and accuracy for speed.

I have searched the web, but all I found involved overly complex calculations which if implemented directly appear to lack the performance I need. Is there a more "programmatical" way to do this, like using matrices or scaling the ellipse into a circle?

Answer 1

Let's assume the ellipse E is in "standard position", center at the origin and axes parallel to the coordinate axes:

    (x/a)^2 + (y/b)^2 = 1   where a > b > 0

The boundary cases a=b are circles, where the normal lines are simply ones that pass through the center (origin) and are thus easy to find. So we omit discussion of these cases.

The slope of the tangent to the ellipse at any point (x,y) may be found by implicit differentiation:

    dy/dx = -(b^2 x)/(a^2 y)

For the line passing through (x,y) and a specified point p = (u,v) not on the ellipse, that is normal to ellipse E when its slope is the negative reciprocal of dy/dx:

    (y-v)/(x-u) * (-b^2 x)/(a^2 y) = -1       (N)

which simplifies to:

    (x - (1+g)u) * (y + gv) = -g(1+g)uv  where g = b^2/(a^2 - b^2)

In this form we recognize it is the equation for a right rectangular hyperbola. Depending on how many points of intersection there are between the ellipse and the hyperbola (2,3,4), we have that many normals to E passing through p.

By reflected symmetry, if p is assumed exterior to E, we may take p to be in the first quadrant:

    (u/a)^2 + (v/b)^2 > 1    (exterior to E)
          u,v > 0            (1'st quadrant)

We could have boundary cases where u=0 or v=0, i.e. point p lies on an axis of E, but these cases may be reduced to solving a quadratic, because two normals are the (coinciding) lines through the endpoints of that axis. We defer further discussion of these special cases for the moment.

Here's an illustration with a=u=5,b=v=3 in which only one branch of the hyperbola intersects E, and there will be only two normals:

If the system of two equations in two unknowns (x,y) is reduced to one equation in one unknown, the simplest root-finding method to code is a bisection method, but knowing something about the possible locations of roots/intersections will expedite our search. The intersection in the first quadrant is the nearest point of E to p, and likewise the intersection in the third quadrant is the farthest point of E from p. If the point p were a good bit closer to the upper endpoint of the minor axis, the branches of the hyperbola would shift together enough to create up to two more points of intersection in the fourth quadrant.

One approach would be to parameterize E by points of intersection with the x-axis. The lines from p normal to the ellipse must intersect the major axis which is a finite interval [-a,+a]. We can test both the upper and lower points of intersection q=(x,y) of a line passing through p=(u,v) and (z,0) as z sweeps from -a to +a, looking for places where the ellipse and hyperbola intersect.

In more detail:

1. Find the upper and lower points `q` of intersection of E with the
   line through `p` and `(z,0)` (amounts to solving a quadratic)

3. Check the sign of a^2 y(x-u) - b^2 x(y-v) at `q=(x,y)`, because it
   is zero if and only `q` is a point of normal intersection

Once a subinterval is detected (either for upper or lower portion) where the sign changes, it can be refined to get the desired accuracy. If only modest accuracy is needed, there may be no need to use faster root finding methods, but even if they are needed, having a short subinterval that isolates a root (or root pair in the fourth quadrant) will be useful.

** more to come comparing convergence of various methods **

Answer 2

I had to solve a problem similar to this, for GPS initialization. The question is: what is the latitude of a point interior to the Earth, especially near the center, and is it single-valued? There are lots of methods for converting ECEF cartesian coordinates to geodetic latitude, longitude and altitude (look up "ECEF to Geodetic"). We use a fast one with only one divide and sqrt per iteration, instead of several trig evaluations like most methods, but since I can't find it in the wild, I can't give it to you here. I would start with Lin and Wang's method, since it only uses divisions in its iterations. Here is a plot of the ellipsoid surface normals to points within 100 km of Earth's center (North is up in the diagram, which is really ECEF Z, not Y):

The star-shaped "caustic" in the figure center traces the center of curvature of the WGS-84 ellipsoid as latitude is varied from pole to equator. Note that the center of curvature at the poles is on the opposite side of the equator, due to polar flattening, and that the center of curvature at the equator is nearer to the surface than the axis of rotation.

Wherever lines cross, there is more than one latitude for that cartesian position. The green circle shows where our algorithm was struggling. If you consider that I cut off these normal vectors where they reach the axis, you would have even more normals for a given position for the problem considered in this SO thread. You would have 4 latitudes / normals inside the caustic, and 2 outside.

Answer 3

The problem can be expressed as the solution of a cubic equation which gives 1, 2, or 3 real roots. For the derivation and closed form solution see Appendix B of Geodesics on an ellipsoid of revolution. The boundary between 1 and 3 solutions is an astroid.