Pastebin

def _MassEnclosed(field,data):
    msink = data.get_field_parameter("sink")
    c = data.get_field_parameter("center")
    Radius = data['Radius']
    R = Radius.max()
    sp = pf.h.sphere(c,1.2*R)
    M = msink*data['Ones']

    nbin = int(ceil(R*pf['au']))
    Men = BinnedProfile1D(sp, nbin, 'Radius', 0., R, log_space=False)
    Men.add_fields('CellMass',weight=None,accumulation=True)
    Minterp = Men['CellMass'] # Rinterp is every AU so R[1] = 1 AU

    for i in range(Radius.size):
        M[i] += Minterp[int(np.round(Radius[i]*pf['au']))]

    return M
add_field("MassEnclosed", function=_MassEnclosed)