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