def predict_angles(p0_star, experiment, s0=None):
    """Predict Ewald sphere crossing angle for RLP x, returned in order
    (entering, exiting)."""
    from scitbx import matrix
    import math
    a = matrix.col(experiment.goniometer.get_rotation_axis())
    if s0 is None:
        b = matrix.col(experiment.beam.get_s0())
    else:
        b = s0
    m2 = a.normalize()
    m1 = m2.cross(b.normalize())
    m3 = m1.cross(m2)
    p0_sqr = p0_star.dot(p0_star)
    rho = math.sqrt(p0_sqr - p0_star.dot(m2) ** 2)
    p_star_m3 = (-0.5 * p0_sqr - p0_star.dot(m2) * b.dot(m2)) / b.dot(m3)
    p_star_m2 = p0_star.dot(m2)
    if rho ** 2 < p_star_m3 ** 2:
        return None
    p_star_m1 = math.sqrt(rho ** 2 - p_star_m3 ** 2)
    p0_star_m1 = p0_star.dot(m1)
    p0_star_m2 = p0_star.dot(m2)
    p0_star_m3 = p0_star.dot(m3)
    cp1 = +p_star_m1 * p0_star_m1 + p_star_m3 * p0_star_m3
    cp2 = -p_star_m1 * p0_star_m1 + p_star_m3 * p0_star_m3
    sp1 = +p_star_m1 * p0_star_m3 - p_star_m3 * p0_star_m1
    sp2 = -p_star_m1 * p0_star_m3 - p_star_m3 * p0_star_m1
    return math.atan2(sp1, cp1), math.atan2(sp2, cp2)
