docs: def bodycomp(mass, tbw, method='reilly', simulate=False, n_rand=1000):
      '''Create dataframe with derived body composition values
  
      Args
      ----
      mass: ndarray
          Mass of the seal (kg)
      tbw: ndarray
          Total body water (kg)
      method: str
          name of method used to derive composition values
      simulate: bool
          switch for generating values with random noise
      n_rand: int
          number of density values to simulate
  
      Returns
      -------
      field: pandas.Dataframe
          dataframe containing columns for each body composition value
  
      References
      ----------
      Reilly, J.J., Fedak, M.A., 1990. Measurement of the body composition of
      living gray seals by hydrogen isotope dilution. Journal of Applied
      Physiology 69, 885–891.
  
      Gales, R., Renouf, D., Noseworthy, E., 1994. Body composition of harp
      seals. Canadian journal of zoology 72, 545–551.
      '''
      import numpy
      import pandas
  
      if len(mass) != len(tbw):
          raise SystemError('`mass` and `tbw` arrays must be the same length')
  
      bc = pandas.DataFrame(index=range(len(mass)))
  
      rnorm = lambda n, mu, sigma: numpy.random.normal(mu, sigma, n)
  
      if method == 'reilly':
          if simulate is True:
              bc['ptbw'] = 100 * (tbw / mass)
              bc['ptbf'] = 105.1 - (1.47 * bc['ptbw']) + rnorm(n_rand, 0, 1.1)
              bc['ptbp'] = (0.42 * bc['ptbw']) - 4.75 + rnorm(n_rand, 0, 0.8)
              bc['tbf'] = mass * (bc['ptbf'] / 100)
              bc['tbp'] = mass * (bc['ptbp'] / 100)
              bc['tba'] = 0.1 - (0.008 * mass) + \
                                 (0.05 * tbw) + rnorm(0, 0.3, n_rand)
              bc['tbge'] = (40.8 * mass) - (48.5 * tbw) - \
                            0.4 + rnorm(0, 17.2, n_rand)
          else:
              bc['ptbw'] = 100 * (tbw / mass)
              bc['ptbf'] = 105.1 - (1.47 * bc['ptbw'])
              bc['ptbp'] = (0.42 * bc['ptbw']) - 4.75
              bc['tbf'] = mass * (bc['ptbf'] / 100)
              bc['tbp'] = mass * (bc['ptbp'] / 100)
              bc['tba'] = 0.1 - (0.008 * mass) + (0.05 * tbw)
              bc['tbge'] = (40.8 * mass) - (48.5 * tbw) - 0.4
      elif method == 'gales':
          if simulate is True:
              raise ValueError('Random error simulation is currently only '
                  'implemented for `method` `reilly`. `simulate` must be passed '
                  'as `False` when using `method` `gales`.')
          else:
              bc['ptbw'] = 100 * (tbw / mass)
              bc['tbf'] = mass - (1.37 * tbw)
              bc['tbp'] = 0.27 * (mass - bc['tbf'])
              bc['tbge'] = (40.8 * mass) - (48.5 * tbw) - 0.4
              bc['ptbf'] = 100 * (bc['tbf'] / mass)
              bc['ptbp'] = 100 * (bc['tbp'] / mass)
      else:
          raise ValueError('`method` must be either `reilly` or `gales`, not '
              '`{}`'.format(method))
  
      return bc