// Find intersection(s) of a line segment and a sphere centered at the origin.
// There can be zero, one, or two results.

function line_sphere_intersection(p1, p2, center, r) =
    let (len = norm(p2 - p1))
        len == 0 ? (norm(p1 - center) == r ? [p1] : []) :
        let (dir = (p2 - p1) / len,
             t = (center - p1) * dir,  // Vec from 1st endpoint proj onto segment dir
             pc = p1 + t * dir,        // Point on line closest to center
             t_start = 0,              // Parameter boundaries (p1 = 0, p2 = len)
             t_end = len,
             dist = norm(pc - center), // Distance of point on line from center
             isq = r ^ 2 - dist ^ 2)
            isq < 0 ? [] :             // No intersection possible
                let (i = sqrt(isq),    // Complete the Pythagorean calculation
                     ti1 = t - i,      // Parameter at two intersections
                     ti2 = t + i,
                     pi1 = p1 + dir * ti1,   // Two points of intersection
                     pi2 = p1 + dir * ti2)   //   (if i == 0 then they're equal)
                    t_start <= ti1 ? (t_end < ti1 ? [] :
                                      t_end < ti2 || i == 0 ? [pi1] :
                                      [pi1, pi2]) :
                    t_start <= ti2 && t_end >= ti2 ? [pi2] :
                    [];