def next_point(
    point_x: float, point_y: float, incoming_gradient: float
) -> tuple[float, float, float]:
    # normal_gradient = gradient of line through which the beam is reflected
    # outgoing_gradient = gradient of reflected line
    normal_gradient = point_y / 4 / point_x
    s2 = 2 * normal_gradient / (1 + normal_gradient * normal_gradient)
    c2 = (1 - normal_gradient * normal_gradient) / (
        1 + normal_gradient * normal_gradient
    )
    outgoing_gradient = (s2 - c2 * incoming_gradient) / (c2 + s2 * incoming_gradient)

    # to find the next point, solve the simultaeneous equations:
    # y^2 + 4x^2 = 100
    # y - b = m * (x - a)
    # ==> A x^2 + B x + C = 0
    quadratic_term = outgoing_gradient**2 + 4
    linear_term = 2 * outgoing_gradient * (point_y - outgoing_gradient * point_x)
    constant_term = (point_y - outgoing_gradient * point_x) ** 2 - 100

    x_minus = (
        -linear_term - sqrt(linear_term**2 - 4 * quadratic_term * constant_term)
    ) / (2 * quadratic_term)
    x_plus = (
        -linear_term + sqrt(linear_term**2 - 4 * quadratic_term * constant_term)
    ) / (2 * quadratic_term)

    # two solutions, one of which is our input point
    next_x = x_minus if isclose(x_plus, point_x) else x_plus
    next_y = point_y + outgoing_gradient * (next_x - point_x)

    return next_x, next_y, outgoing_gradient