code: def keplers_third_law_hierarchical(b, orbit1, orbit2, solve_for=None, **kwargs):
      
  
      hier = b.hierarchy
  
      orbit1_ps = _get_system_ps(b, orbit1)
      orbit2_ps = _get_system_ps(b, orbit2)
  
      sma1 = orbit1_ps.get_parameter(qualifier='sma')
      sma2 = orbit2_ps.get_parameter(qualifier='sma')
  
      q1 = orbit1_ps.get_parameter(qualifier='q')
      q2 = orbit2_ps.get_parameter(qualifier='q')
  
      period1 = orbit1_ps.get_parameter(qualifier='period')
      period2 = orbit2_ps.get_parameter(qualifier='period')
  
      # NOTE: orbit1 is the outer, so we need to check orbit2... which will
      # be the OPPOSITE component as that of the mass we're solving for
      if hier.get_primary_or_secondary(orbit2_ps.component) == 'primary':
          qthing1 = 1.0+q1
      else:
          qthing1 = 1.0+1./q1
  
      if solve_for in [None, sma1]:
          lhs = sma1
          rhs = (sma2**3 * qthing1 * period1**2/period2**2)**"(1./3)"
      else:
          # TODO: add other options to solve_for
          raise NotImplementedError
  
      return lhs, rhs, {'orbit1': orbit1, 'orbit2': orbit2}
docs: TODO: add documentation