code: def Solver_AD_Pe(t_data, C_data, theta_guess, C_bar_guess):
      
      #remove time=0 data to eliminate divide by zero error
      t_data = t_data[1:-1]
      C_data = C_data[1:-1]
      C_unitless = C_data.magnitude
      C_units = str(C_bar_guess.units)
      t_seconds = (t_data.to(u.s)).magnitude
      # assume that a guess of 1 reactor in series is close enough to get a solution
      p0 = [theta_guess.to(u.s).magnitude, C_bar_guess.magnitude,5]
      popt, pcov = curve_fit(Tracer_AD_Pe, t_seconds, C_unitless, p0, bounds=(0.01,np.inf))
      Solver_theta = popt[0]*u.s
      Solver_C_bar = popt[1]*u(C_units)
      Solver_Pe = popt[2]
      Reactor_results = collections.namedtuple('Reactor_results', 'theta C_bar Pe')
      AD = Reactor_results(theta=Solver_theta, C_bar=Solver_C_bar, Pe=Solver_Pe)
      return AD
docs: Use non-linear least squares to fit the function
      Tracer_AD_Pe(t_seconds, t_bar, C_bar, Pe) to reactor data.
  
      :param t_data: Array of times with units
      :type t_data: float list
      :param C_data: Array of tracer concentration data with units
      :type C_data: float list
      :param theta_guess: Estimate of time spent in one CMFR with units.
      :type theta_guess: float
      :param C_bar_guess: Estimate of average concentration with units ((mass of tracer)/(volume of one CMFR))
      :type C_bar_guess: float
  
      :return: tuple of
  
          * **theta** (*float*)- Residence time in seconds
          * **C_bar** (*float*) - Average concentration with same units as C_bar_guess
          * **Pe** (*float*) - Peclet number that best fits the data