def powerhalo(r,rs=1.,rc=0.,alpha=1.,beta=1.e-7):
    """return generic twopower law distribution

    inputs
    ----------
    r      : (float) radius values
    rs     : (float, default=1.) scale radius
    rc     : (float, default=0. i.e. no core) core radius
    alpha  : (float, default=1.) inner halo slope
    beta   : (float, default=1.e-7) outer halo slope

    returns
    ----------
    densities evaluated at r

    notes
    ----------
    different combinations are known distributions.
    alpha=1,beta=2 is NFW
    alpha=1,beta=3 is Hernquist
    alpha=2.5,beta=0 is a typical single power law halo


    """
    ra = r/rs
    return 1./(((ra+rc)**alpha)*((1+ra)**beta))