How do I map my data to the Airy function in python?

My code takes an image of a hole in a hole and matches the data in Gaussian. Using a Gaussian fit, it computes the full width with half the maximum. This tells me the resolution of my imaging system.

Here is the fitting one I am getting with my code right now:

According to the theory, for diffraction images with holes, the data should correspond to the Airy disk function. For completeness, I want the data to be bound to a Bessel function or Airy disk pattern. I cannot find packages that will match these features.

Here is the image I'm using:

You can just make out the outer stripes around the central bright spot. These are the boundaries that I want to take into account in my physical form.

import numpy as np
import scipy.optimize as opt
import PIL
from PIL import ImageFilter
from pylab import *

#defining the Gaussian
def gauss(x, p): # p[0]==mean, p[1]==stdev
    return 1.0/(p[1]*np.sqrt(2*np.pi))*np.exp(-(x-p[0])**2/(2*p[1]**2))

im = PIL.Image.open('C:/Documents/User/3000.bmp').convert("L")  #convert to array
imArr = np.array(im, dtype=float)
bg = np.average(imArr)  #find the background, subtract it
imArr = imArr - bg

#get the approx coordinates of brightest spot by filtering
im2 = im.filter(ImageFilter.GaussianBlur(radius=2))
imArr2 = np.array(im2, dtype=float)
tuple = unravel_index(imArr2.argmax(), imArr2.shape)

#find and plot FWHM for the brightest spot
x = np.arange(tuple[1] - 100, tuple[1] + 100, dtype=np.float)
y = imArr[tuple[0], tuple[1] - 100:tuple[1] + 100]

y /= ((max(x) - min(x)) / len(x)) * np.sum(y) # renormalize to a proper Gaussian

p0 = [tuple[1], tuple[0]]
errfunc = lambda p, x, y: gauss(x, p) - y # distance to the target function
p1, success = opt.leastsq(errfunc, p0[:], args=(x, y))

fit_mu, fit_stdev = p1

FWHM = 2*np.sqrt(2*np.log(2))*fit_stdev
print "FWHM", FWHM

plt.plot(x,y)
plt.plot(x, gauss(x,p1), lw=3, alpha=.5, color='r')
plt.axvspan(fit_mu-FWHM/2, fit_mu+FWHM/2, facecolor='g', alpha=0.5)
plt.show()