old_contents: 
  import numpy as np
  from itertools import cycle
  import scipy.optimize as opt
  
  '''
  Krumholz model for the fraction of H2 as a function of surface density and
  metallicity.
  
  Following Schruba+11, the output is in terms of the H2-HI fraction rather than
  the the molecular fraction.
  '''
  
  
  def krumholz_ratio_model(Sigma, Z=0.5, c=1):
      '''
  
      Parameters
      ----------
      sigma : float or np.ndarray
          Surface Density in Msol/pc^2
      Z : float
          Metallicity in units of solar metallicity.
      c : float
          Clumping fraction. Expected to be near unity with a resolution of
          100 pc. c>=1.
      '''
  
      Sigma_comp = c * Sigma
  
      chi = 0.77 * (1 + 3.1 * np.power(Z, 0.365))
  
      s = np.log(1 + 0.6 * chi) / (0.04 * Sigma_comp * Z)
  
      delta = 0.0712 * np.power((0.1 / s) + 0.675, -2.8)
  
      frac = 1 - np.power(1 + np.power(0.75 * (s / (1 + delta)), -5.), -0.2)
  
      return frac / (1 - frac)
  
  
  def alternate_krumholz_ratio_model(Sigma, psi=1.0, c=1, Z=0.1):
      '''
  
      Parameters
      ----------
      sigma : float or np.ndarray
          Surface Density in Msol/pc^2
      psi : float
          Dust-adjusted radiation field (unitless). Related to metallicity
          (among other things). At Z=1, psi=1.6, and at Z=0.1, psi=1.0
      c : float
          Clumping fraction. Expected to be near unity with a resolution of
          100 pc. c>=1.
      '''
  
      Sigma_comp = c * Sigma
  
      s = Sigma_comp * Z / float(psi)
  
      term1 = (s / 11.) ** 3 * ((125 + s) / (96 + s)) ** 3
  
      return np.power(1 + term1, 1 / 3.) - 1
  
  
  def optimize_clump_factors(Sigma, R, Z=0.5):
      '''
      Solve for the clump factor needed to intersect each point.
  
      Parameters
      ----------
      Sigma : Quantity
          Surface densities.
      Z : float or array
          Metallicity values.
      '''
  
      if not isinstance(Z, np.ndarray):
          Z = cycle([Z])
  
      clump_values = []
  
      for sig, r, z in zip(Sigma, R, Z):
  
          def model(sigma, c):
              return krumholz_ratio_model(sigma, Z=z, c=c)
  
          popt, pcov = opt.curve_fit(model, sig, r)
  
          clump_values.append(popt)
  
      # return np.array(clump_values)
      return clump_values
  
new_contents: 
  import numpy as np
  from itertools import cycle
  import scipy.optimize as opt
  
  '''
  Krumholz model for the fraction of H2 as a function of surface density and
  metallicity.
  
  Following Schruba+11, the output is in terms of the H2-HI fraction rather than
  the the molecular fraction.
  '''
  
  
  def krumholz_ratio_model(Sigma, Z=0.5, c=1):
      '''
  
      Parameters
      ----------
      sigma : float or np.ndarray
          Surface Density in Msol/pc^2
      Z : float
          Metallicity in units of solar metallicity.
      c : float
          Clumping fraction. Expected to be near unity with a resolution of
          100 pc. c>=1.
      '''
  
      Sigma_comp = c * Sigma
  
      chi = 0.77 * (1 + 3.1 * np.power(Z, 0.365))
  
      s = np.log(1 + 0.6 * chi) / (0.04 * Sigma_comp * Z)
  
      delta = 0.0712 * np.power((0.1 / s) + 0.675, -2.8)
  
      frac = 1 - np.power(1 + np.power(0.75 * (s / (1 + delta)), -5.), -0.2)
  
      return frac / (1 - frac)
  
  
  def alternate_krumholz_ratio_model(Sigma, psi=1.0, c=1, Z=0.1):
      '''
  
      Parameters
      ----------
      sigma : float or np.ndarray
          Surface Density in Msol/pc^2
      psi : float
          Dust-adjusted radiation field (unitless). Related to metallicity
          (among other things). At Z=1, psi=1.6, and at Z=0.1, psi=1.0
      c : float
          Clumping fraction. Expected to be near unity with a resolution of
          100 pc. c>=1.
      '''
  
      Sigma_comp = c * Sigma
  
      s = Sigma_comp * Z / float(psi)
  
      term1 = (s / 11.) ** 3 * ((125 + s) / (96 + s)) ** 3
  
      return np.power(1 + term1, 1 / 3.) - 1
  
  
  def optimize_clump_factors(Sigma, R, Z=0.5):
      '''
      Solve for the clump factor needed to intersect each point.
  
      Parameters
      ----------
      Sigma : Quantity
          Surface densities.
      Z : float or array
          Metallicity values.
      '''
  
      if not isinstance(Z, np.ndarray):
          Z = cycle([Z])
  
      clump_values = []
  
      for sig, r, z in zip(Sigma, R, Z):
  
          def model(sigma, c):
              return krumholz_ratio_model(sigma, Z=z, c=c)
  
          popt, pcov = opt.curve_fit(model, sig, r)
  
          clump_values.append(popt[0])
  
      return np.array(clump_values)
  
subject: Fix output array shape