Also called Great-circle-distance --
    use long-ass formula from wikipedia (last in section):
    https://en.wikipedia.org/wiki/Great-circle_distance#Computational_formulas

    Space angle only makes sense in lon-lat, so convert zenith -> latitude.

def space_angle(phi1, theta1, phi2, theta2):
    """Also called Great-circle-distance --
    use long-ass formula from wikipedia (last in section):
    https://en.wikipedia.org/wiki/Great-circle_distance#Computational_formulas

    Space angle only makes sense in lon-lat, so convert zenith -> latitude.
    """
    from numpy import pi, sin, cos, arctan2, sqrt, square
    lamb1 = pi / 2 - theta1
    lamb2 = pi / 2 - theta2
    lambdelt = lamb2 - lamb1
    under = sin(phi1) * sin(phi2) + cos(phi1) * cos(phi2) * cos(lambdelt)
    over = sqrt(
        np.square((cos(phi2) * sin(lambdelt))) +
        square(cos(phi1) * sin(phi2) - sin(phi1) * cos(phi2) * cos(lambdelt))
    )
    angle = arctan2(over, under)
    return angle