def Hernquist(r,a,rc=0,rho0=1.,beta=0):
    """return unscaled Herquist density

    r      : the radius to evaluate
    a      : the scale radius
    rc     : the core radius
    rho0   : the central density


    """
    ra = r/a
    return rho0/((ra+rc)*((1+ra)**3.))