Important: For homework #3 it’s not mandatory that you fit the depth of your lines. This quick tutorial is meant to be informative for those who do decide to fit the depth of the Na I D lines.

Hey class 👋

I wanted to make a very quick tutorial on how to generally fit a curve with scipy

Let’s begin with three function models you’re given in this tutorial of a Gaussian, Lorentzian, and a Vogt profile:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def G(x, alpha):
    """ Return Gaussian line shape at x with HWHM alpha """
    return np.sqrt(np.log(2) / np.pi) / alpha\\
                             * np.exp(-(x / alpha)**2 * np.log(2))

def L(x, gamma):
    """ Return Lorentzian line shape at x with HWHM gamma """
    return gamma / np.pi / (x**2 + gamma**2)

def V(x, alpha, gamma):
    """
    Return the Voigt line shape at x with Lorentzian component HWHM gamma
    and Gaussian component HWHM alpha.

    """
    sigma = alpha / np.sqrt(2 * np.log(2))

    return np.real(wofz((x + 1j*gamma)/sigma/np.sqrt(2))) / sigma\\
                                                           /np.sqrt(2*np.pi

Given these functions now, one thing we can do is fit these models to real data. This can be done using the curve_fit function from scipy. Here’s how it’s generally done:

Let’s say I have some random data (x,y) and I want do fit a Gaussian profile (i.e the G function from above). This is what my data looks like:

To fit the data, I will use the scipy function curve_fit

from scipy.optimize import curve_fit

# NOTES: 
# 1) Remember that I decided to call my data variables: x, y 
# 2) curve_fit also takes uncertainties for the y-variable (i.e use sigma argument in the curve_fit function)
# 3) In curve_fit, I will pass in the following arguments: 
# Model (here I'm fitting my Gaussian function)
# x,y data
# bounds: i.e the boundary at which you want your parameter to vary
# bounds (for 2 parameters): ((min_param_1, min_param_2), (max_param_1, max_param_2))

# Optimize and fit data
best_fit, covariance = curve_fit(G, x, y, bounds=((-5, 0.1), (5, 3))

Cool, now we can plot the best_fit parameters to see what it looks like

plt.figure(figsize=(3,3))
plt.scatter(x, y, color='tomato')
plt.plot(x, G(x, *best_fit), color='k', lw=2)