def late_filling(target, pressure='pore.pressure',
                 Pc_star='pore.pc_star',
                 Swp_star=0.2, eta=3):
    r"""
    Calculates the fraction of a pore or throat filled with invading fluid
    based on the capillary pressure in the invading phase.  The invading phase
    volume is calculated from:

        .. math::
            S_{nwp} = 1 - S_{wp}^{*} (P^{*}/P_{c})^{\eta}

    Parameters
    ----------
    pressure : string
        The capillary pressure in the non-wetting phase (Pc > 0).

    Pc_star : string
        The minimum pressure required to create an interface within the pore
        body or throat.  Typically this would be calculated using the Washburn
        equation.

    Swp_star : float
        The residual wetting phase in an invaded pore or throat at a pressure
        of ``pc_star``.

    eta : float
        Exponent controlling the rate at which wetting phase is displaced with
        increasing pressure.

    Returns
    -------
    An array containing the fraction of each pore or throat that would be
    filled with non-wetting phase at the given phase pressure.  This does not
    account for whether or not the element is actually invaded, which requires
    a percolation algorithm of some sort.

    """
    element = pressure.split('.')[0]
    network = target.project.network
    phase = target.project.find_phase(target)
    pc_star = phase[Pc_star]
    Pc = phase[pressure]
    # Remove any 0's from the Pc array to prevent numpy div by 0 warning
    Pc = sp.maximum(Pc, 1e-9)
    Swp = Swp_star*((pc_star/Pc)**eta)
    values = sp.clip(1 - Swp, 0.0, 1.0)
    # Now map element onto target object
    if element == 'throat':
        Ts = network.map_throats(throats=target.Ts, origin=target)
        values = values[Ts]
    else:
        Ps = network.map_pores(pores=target.Ps, origin=target)
        values = values[Ps]
    return values