title: Multiple gaussian auto-fitting in Python
body: <p>The problem is to fitter on all my wavelength peaks a Gaussian in order to make a medium adjustment as accurate as possible</p>
  
  <p>My question is how to make the Gaussian adjustment on all my peaks automatically without having to manually specify the coordinates of the peaks</p>
  
  <p>For that, I realized the Gaussian adjustment of the brightest peaks, but I would like to generalize it to the following peaks. Subsequently, the Gaussian adjustment will allow me to obtain a polynomial adjustment fine enough to stagger pixels in wavelength</p>
  
  <pre class="lang-py prettyprint-override"><code>
  import numpy as np
  from astropy.io import fits
  import matplotlib.pyplot as plt
  from scipy import interpolate
  from tqdm import tqdm
  from scipy import ndimage
  import peakutils
  from scipy.optimize import curve_fit
  
  def gauss(x, x0, amp, wid):
      return amp * np.exp( -((x - x0)/wid)**2)
  
  def multi_gauss(x, *params):
      y = np.zeros_like(x)
      for i in range(0, len(params), 3):
          x0, amp, wid = params[i:i+3]
          y = y + gauss(x, x0, amp, wid)
      return y
  
  
  neon = fits.getdata(data_directory + wave_filename + '.fits')
  
  neon_sp = np.mean(neon, axis= 0)
  n_pix = len(neon_sp)
  
  peaks_index = peakutils.peak.indexes(neon_sp, thres=0.05, min_dist=2)
  
  ### peals around the brightest peak
  bright_index = peaks_index[np.argmax(neon_sp[peaks_index])]
  
  delta_pix = 20
  ind_min = bright_index - delta_pix
  ind_max = bright_index + delta_pix
  
  peak_select = peaks_index[np.where((peaks_index &gt; ind_min) &amp; (peaks_index &lt; ind_max))]
  peak_select_sort = peak_select[np.argsort(-neon_sp[peak_select])]
  
  if peak_select_sort[1] &gt; peak_select_sort[0] :
      ind_max = bright_index + 40
  else : 
      ind_min = bright_index - 40
  
  peak_select = peaks_index[np.where((peaks_index &gt; ind_min) &amp; (peaks_index &lt; ind_max))]
  peak_select_sort = peak_select[np.argsort(-neon_sp[peak_select])]
  
  plt.figure(num=0)
  plt.clf()
  plt.plot(neon_sp)
  plt.plot(peaks_index,neon_sp[peaks_index], 'r+')
  plt.plot(peak_select,neon_sp[peak_select], 'ro')
  
  ### Gaussian fit
  x = np.arange(n_pix)
  xx = np.arange(0, n_pix, .1)
  
  
  n_peak = 4
  bright_index_fit = np.zeros(n_peak)
  
  for i in range(n_peak):
      p = peak_select_sort[i]
      guess = [p, neon_sp[p], .5]
      popt, pcov = curve_fit(gauss, x, neon_sp, p0=guess)
      fit = gauss(xx, *popt)
  
      bright_index_fit[i] = popt[0]
  
      plt.plot(xx,fit, '--')
  
  bright_wave = [703.2, 724.5, 693.0, 743.9]
  
  </code></pre>