title: Fitting Voigt function to data in Python
body: <p>I recently got a script running to fit a gaussian to my absorption profile with <a href="https://stackoverflow.com/questions/59047395/fitting-gaussian-to-absorbtion-line-in-python">help of SO</a>. My hope was that things would work fine if I simply replace the Gauss function by a Voigt one, but this seems not to be the case. I think mainly due to the fact that it is a shifted voigt. </p>
  
  <p>Edit: The profiles are absorption lines that vary in optical thickness. In practice they will be a mix between optically thick and thin features. Like the bottom part in this diagram. The current data will be more like the top image, but maybe the bottom is already flattened a bit. (And we only see the left side of the profile, a bit beyond the center)
  <a href="https://i.stack.imgur.com/mL0DU.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/mL0DU.png" alt="enter image description here"></a></p>
  
  <p>For a Gauss it looks like this and as predicted the bottom seems to be less deep than the fit wants it to be, but still quite close. The profile itself should still be a voigt though. But now I realize that the central points might throw off the fit. So maybe a weight should be added based on wing position? 
  <a href="https://i.stack.imgur.com/8wwwF.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/8wwwF.png" alt="enter image description here"></a></p>
  
  <p>I'm mostly wondering if the shifted function could be mis-defined or if its my starting values.</p>
  
  <pre><code>import matplotlib.pyplot as plt
  from scipy.optimize import curve_fit
  import numpy as np
  from scipy.special import wofz
  
  x = np.arange(13)
  xx = xx = np.linspace(0, 13, 100)
  y = np.array([19699.959 , 21679.445 , 21143.195 , 20602.875 , 16246.769 ,
          11635.25  ,  8602.465 ,  7035.493 ,  6697.0337,  6510.092 ,
                7717.772 , 12270.446 , 16807.81  ])
  # weighted arithmetic mean (corrected - check the section below)
  #mean = 2.4
  sigma = 2.4
  gamma = 2.4
  
  def Gauss(x, y0, a, x0, sigma):
      return y0 + a * np.exp(-(x - x0)**2 / (2 * sigma**2))
  
  
  def Voigt(x, x0, y0, a, sigma, gamma):
      #sigma = alpha / np.sqrt(2 * np.log(2))
  
      return y0 + a * np.real(wofz((x - x0 + 1j*gamma)/sigma/np.sqrt(2))) / sigma /np.sqrt(2*np.pi)
  
  
  popt, pcov = curve_fit(Voigt, x, y, p0=[8, np.max(y), -(np.max(y)-np.min(y)), sigma, gamma])
  #p0=[8, np.max(y), -(np.max(y)-np.min(y)), mean, sigma])
  
  plt.plot(x, y, 'b+:', label='data')
  plt.plot(xx, Voigt(xx, *popt), 'r-', label='fit')
  
  plt.legend()
  plt.show()
  </code></pre>