docs: def eval(e, amplitude, e_0, e_break, alpha_1, alpha_2):
          """One dimensional broken power law model function"""
          K = np.where(e < e_break, 1, (e_break / e_0) ** (alpha_2 - alpha_1))
          alpha = np.where(e < e_break, alpha_1, alpha_2)
          return amplitude * K * (e / e_0) ** -alpha