Inverse correction of refractive index and radius for Rytov

    This method returns the inverse of :func:`correct_rytov_output`.

    Parameters
    ----------
    radius_sc: float
        Systematically corrected radius of the sphere [m]
    sphere_index_sc: float
        Systematically corrected refractive index of the sphere
    medium_index: float
        Refractive index of the surrounding medium
    radius_sampling: int
        Number of pixels used to sample the sphere radius when
        computing the Rytov field.

    Returns
    -------
    radius: float
        Fitted radius of the sphere [m]
    sphere_index: float
        Fitted refractive index of the sphere

    See Also
    --------
    correct_rytov_output: the inverse of this method

def correct_rytov_sc_input(radius_sc, sphere_index_sc, medium_index,
                           radius_sampling):
    """Inverse correction of refractive index and radius for Rytov

    This method returns the inverse of :func:`correct_rytov_output`.

    Parameters
    ----------
    radius_sc: float
        Systematically corrected radius of the sphere [m]
    sphere_index_sc: float
        Systematically corrected refractive index of the sphere
    medium_index: float
        Refractive index of the surrounding medium
    radius_sampling: int
        Number of pixels used to sample the sphere radius when
        computing the Rytov field.

    Returns
    -------
    radius: float
        Fitted radius of the sphere [m]
    sphere_index: float
        Fitted refractive index of the sphere

    See Also
    --------
    correct_rytov_output: the inverse of this method
    """
    params = get_params(radius_sampling)

    # sage script:
    # var('sphere_index, sphere_index_sc, na, nb, medium_index')
    # x = sphere_index / medium_index - 1
    # eq = sphere_index_sc == sphere_index + ( na*x^2 + nb*x) * medium_index
    # solve([eq], [sphere_index])
    # (take the positive sign solution)
    na = params["na"]
    nb = params["nb"]

    prefac = medium_index / (2 * na)
    sm = 2 * na - nb - 1
    rt = nb**2 - 4 * na + 2 * nb + 1 + 4 / medium_index * na * sphere_index_sc
    sphere_index = prefac * (sm + np.sqrt(rt))

    x = sphere_index / medium_index - 1

    radius = radius_sc / (params["ra"] * x**2
                          + params["rb"] * x
                          + params["rc"])

    return radius, sphere_index