docs: def _do_two_gaussian_fit(freqs, signal, bounds=None):
     """
     Helper function for the two gaussian fit
     """
     initial = _two_func_initializer(freqs, signal)
     # Edit out the ones we want in the order we want them: 
     initial = (initial[0], initial[1],
                initial[6], initial[7],
                initial[2], initial[3],
                initial[10], initial[11])
  
     # We want to preferntially weight the error on estimating the height of the
     # individual peaks, so we formulate an error-weighting function based on
     # these peaks, which is simply a two-gaussian bumpety-bump:
     w = (ut.gaussian(freqs, initial[0], 0.075, 1, 0, 0) +
          ut.gaussian(freqs, initial[1], 0.075, 1, 0, 0))
  
     # Further, we want to also optimize on the individual gaussians error, to
     # restrict the fit space a bit more. For this purpose, we will pass a list
     # of gaussians with indices into the parameter list, so that we can do
     # that (see mopt.err_func for the mechanics).
     func_list = [[ut.gaussian, [0,2,4,6,7],
                   ut.gaussian(freqs, initial[0], 0.075, 1, 0, 0)],
                  [ut.gaussian, [1,3,5,6,7],
                   ut.gaussian(freqs, initial[1], 0.075, 1, 0, 0)]]
  
     params, _ = lsq.leastsqbound(mopt.err_func, initial,
                                  args=(freqs, np.real(signal),
                                  ut.two_gaussian, w, func_list),
                                  bounds=bounds)
  
     return params