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 **