1"""Routines solving basic geometric problems in astronomy."""
2
3from numpy import nan, where
4from .functions import length_of
5
6def intersect_line_and_sphere(endpoint, center, radius):
7    """Compute distance to intersections of a line and a sphere.
8
9    Given a line through the origin (0,0,0) and an |xyz| ``endpoint``,
10    and a sphere with the |xyz| ``center`` and scalar ``radius``,
11    return the distance from the origin to their two intersections.
12
13    If the line is tangent to the sphere, the two intersections will be
14    at the same distance.  If the line does not intersect the sphere,
15    two ``nan`` values will be returned.
16
17    """
18    # See http://paulbourke.net/geometry/circlesphere/index.html#linesphere
19    # Names "b" and "c" designate the familiar values from the quadratic
20    # formula; happily, a = 1 because we use a unit vector for the line.
21
22    minus_b = 2.0 * (endpoint / length_of(endpoint) * center).sum(axis=0)
23    c = (center * center).sum(axis=0) - radius * radius
24    discriminant = minus_b * minus_b - 4 * c
25    dsqrt = discriminant ** where(discriminant < 0, nan, 0.5)  # avoid sqrt(<0)
26    return (minus_b - dsqrt) / 2.0, (minus_b + dsqrt) / 2.0
27