_validsolvefor['mass'] = ['mass', 'sma', 'period', 'q']
def mass(b, component, solve_for=None, **kwargs):
    """
    Create a constraint for the mass of a star based on Kepler's third
    law from its parent orbit.

    This constraint is automatically created and attached for all stars
    in binary orbits via <phoebe.frontend.bundle.Bundle.set_hierarchy>.

    This is usually passed as an argument to
     <phoebe.frontend.bundle.Bundle.add_constraint> as
     `b.add_constraint('mass', component='primary')`, where `component` is
     one of <phoebe.parameters.HierarchyParameter.get_stars>.

    Arguments
    -----------
    * `b` (<phoebe.frontend.bundle.Bundle>): the Bundle
    * `component` (string): the label of the orbit or component in which this
        constraint should be built.
    * `solve_for` (<phoebe.parameters.Parameter>, optional, default=None): if
        'mass' should not be the derived/constrained parameter, provide which
        other parameter should be derived (ie 'period', 'sma', 'q').  Note:
        you cannot solve_for 'period' and 'sma' in the same orbit as the solution
        will not be unique.

    Returns
    ----------
    * (<phoebe.parameters.Parameter>, <phoebe.parameters.ConstraintParameter>, list):
        lhs (Parameter), rhs (ConstraintParameter), addl_params (list of additional
        parameters that may be included in the constraint), kwargs (dict of
        keyword arguments that were passed to this function).

    Raises
    --------
    * NotImplementedError: if the value of `solve_for` is not implemented.
    """
    # TODO: optimize this - this is currently by far the most expensive constraint (due mostly to the parameter multiplication)

    hier = b.get_hierarchy()
    if not len(hier.get_value()):
        # TODO: change to custom error type to catch in bundle.add_component
        # TODO: check whether the problem is 0 hierarchies or more than 1
        raise NotImplementedError("constraint for mass requires hierarchy")

    component_ps = _get_system_ps(b, component)

    parentorbit = hier.get_parent_of(component)
    parentorbit_ps = _get_system_ps(b, parentorbit)

    m1 = component_ps.get_parameter(qualifier='mass', **_skip_filter_checks)
    masses = []
    masses.append(m1)

    siblings = hier.get_stars_of_sibling_of(component)
    siblings = siblings if isinstance(siblings, list) else [siblings]

    for i, sibling in enumerate(siblings):
        sibling_ps = _get_system_ps(b, sibling)
        m2 = sibling_ps.get_parameter(qualifier='mass', **_skip_filter_checks)
        if i==0:
            msum = m2
        else:
            msum += m2
        masses.append(m2)

    mass = m1
    mass_sibling = msum

    # we need to find the constraint attached to the other component... but we
    # don't know who is constrained, or whether it belongs to the sibling or parent
    # orbit, so we'll have to do a bit of digging.
    mass_constraint_sibling = None
    for p in b.filter(constraint_func='mass', component=[parentorbit, sibling], context='constraint', **_skip_filter_checks).to_list():
        if p.constraint_kwargs['component'] == sibling:
            mass_constraint_sibling = p
            break
    if mass_constraint_sibling is not None:
        sibling_solve_for = mass_constraint_sibling.qualifier
        logger.debug("constraint.mass for component='{}': sibling ('{}') is solved for '{}'".format(component, sibling, sibling_solve_for))
    else:
        # this could happen when we build the first constraint, before the second has been built
        sibling_solve_for = None

    sma = parentorbit_ps.get_parameter(qualifier='sma', **_skip_filter_checks)
    # NOTE: sidereal period
    period = parentorbit_ps.get_parameter(qualifier='period', **_skip_filter_checks)
    q = parentorbit_ps.get_parameter(qualifier='q', **_skip_filter_checks)

    G = c.G.to('solRad3 / (solMass d2)')
    G.keep_in_solar_units = True

    if hier.get_primary_or_secondary(component) == 'primary':
        # NOTE: 1.0 must be on the right for distribution math to behave
        qthing = q+1.0
    else:
        # NOTE: 1.0 must be on the right for distribution math to behave
        qthing = 1./q+1.0

    if solve_for in [None, mass]:
        lhs = mass
        rhs = (4*np.pi**2 * sma**3 ) / (period**2 * qthing * G)

    elif solve_for==sma:
        if sibling_solve_for in ['period', 'sma']:
            raise ValueError("cannot solve for '{}' when sibling ('{}') is solved for '{}'".format(solve_for.twig, sibling, sibling_solve_for))
        lhs = sma
        rhs = ((mass * period**2 * qthing * G)/(4 * np.pi**2))**"(1./3)"

    elif solve_for==period:
        if sibling_solve_for in ['period', 'sma']:
            raise ValueError("cannot solve for '{}' when sibling ('{}') is solved for '{}'".format(solve_for.twig, sibling, sibling_solve_for))
        lhs = period
        rhs = ((4 * np.pi**2 * sma**3)/(mass * qthing * G))**"(1./2)"

    elif solve_for==q:
        lhs = q

        if hier.get_primary_or_secondary(component) == 'primary':
            rhs = mass_sibling / mass
        else:
            rhs = mass / mass_sibling

    else:
        raise NotImplementedError

    return lhs, rhs, [period, sma, q] + masses, {'component': component}