docs: def pack(R, L, pf=None, n=None, rng=None, periodic=False,
           beta_max=1e4, dL_max=0.02, dr_max=0.02):
      """Pack a number of non-intersecting spheres into a periodic system.
  
      Can specify packing by number of spheres or packing fraction.
  
      This implementation uses the Metropolis-Hastings algorithm for an
      NPT system.
  
      Parameters
      ----------
      R: float
          Sphere radius.
      L: float array, shape (d,)
          System lengths.
      pf: float or None
          Packing fraction
      n: integer or None
          Number of spheres.
      rng: RandomState or None
          Random number generator. If None, use inbuilt numpy state.
      periodic: bool
          Whether or not the system is periodic.
  
  
      Metropolis-Hastings parameters
      ------------------------------
      Playing with these parameters may improve packing speed.
  
      beta_max: float, greater than zero.
          Inverse temperature which controls how little noiseis in the system.
      dL_max: float, 0 < dL_max < 1
          Maximum fraction by which to perturb the system size.
      dr_max: float, 0 < dr_max < 1
          Maximum system fraction by which to perturb sphere positions.
  
      Returns
      -------
      r: float array, shape (n, d)
          Coordinates of the centres of the spheres for a valid configuration.
      R_actual: float
          Actual sphere radius used in the packing.
      """
      if pf is not None:
          if pf == 0.0:
              return np.array([]), R
          # If packing fraction is specified, find required number of spheres
          # and the actual packing fraction this will produce
          n, pf_actual = pf_to_n(L, pf, R)
      elif n is not None:
          if n == 0:
              return np.array([]), R
          # If n is specified, find packing fraction
          pf_actual = n_to_pf(L, n, R)
  
      # Calculate an initial packing fraction and system size
      # Start at at most 0.5%; lower if the desired packing fraction is very low
      pf_initial = min(0.005, pf_actual / 2.0)
      # Find system size that will create this packing fraction
      dim = L.shape[0]
      increase_initial_ratio = (pf_actual / pf_initial) ** (1.0 / dim)
      L_0 = L * increase_initial_ratio
  
      # Pack naively into this system
      r_0, R = pack_simple(R, L_0, n=n, rng=rng, periodic=periodic)
  
      mg = metro_rcp_factory(periodic, r_0, L_0, R, dr_max, dL_max, rng=rng)
  
      print('Initial packing done, Initial packing: {:g}'.format(mg.pf))
  
      t = 0
      while mg.pf < pf_actual:
          t += 1
          beta = beta_max * mg.pf
          mg.iterate(beta)
  
          if not t % every:
              print('Packing: {:g}%'.format(100.0 * mg.pf))
  
      print('Final packing: {:g}%'.format(100.0 * mg.pf))
  
      return mg.r, mg.R