def king(r, rc, rt, sigma_0, alpha=2):
    """
    See http://iopscience.iop.org/1538-3881/139/6/2097/fulltext/

    Parameters
    ----------
    r: float
        radius
    rc: float
        core radius
    rt: float
        truncation radius
    sigma_0: float
        central density
    """

    def z(x):
        return 1/(1+(x/rc)**2)**(1./alpha)

    term1 = (1 - z(rt))**-alpha
    term2 = (z(r) - z(rt))**alpha
    sigma = sigma_0 * term1 * term2

    return sigma