def max_ring_density(rho_pl_inf, f=1, verbose=True):
    """
    Use Roche formular to compute maximum density of ring material given inferred planet density in g/cm3.
    since rings cannot be sustained beyond roche limit, we assume the observed radius is the roche limit of possible ringed planet. 
    
    Parameters:
    ----------
    rho_pl_inf: ufloat, float, array-like;
        Inferred density of the planet given the observed radius and estimated mass. ufloat from uncertainties package allows the uncertainties in planet density to be propagated
        to the resulting ring density 
    
    f: float;
    	factor to account for 

    Returns
    -------
    
    rho_ring: Array-like;
        maximum ring density a planet can have. 
    """    
    rho_r = rho_pl_inf * (2.46* f**0.5)**3 
    if verbose: print(f"Any ring around this planet must have density below {rho_r} in order to be within the Roche limit. Check ```T_eq``` of this planet to ensure that ring of such density can exist.")
    return  rho_r