airy_abberations

from scipy.special import j1
import numpy as np

def airy_pattern_theta(theta, I0=1, wavelength=5e-7, aperture_radius=0.3e-2):
  """
  Airy function, with maximum intensity normalised to 1.
  default values for the wavelength and aperture radius are the wavelength of green light,
  and an approximate radius of a pupil.
  """
  return I0 * np.power(wavelength * j1(aperture_radius * np.sin(theta) * 2*np.pi/wavelength) / (np.pi * aperture_radius * np.sin(theta)), 2)

def airy_pattern(x, d, I0=1, wavelength=5e-7, aperture_radius=0.2e-2):
  theta = np.arctan2(x,d)
  return airy_pattern_theta(theta, I0=I0, wavelength=wavelength, aperture_radius=aperture_radius)


def gaussian(x, mu=0, sigma=1):
  return np.exp(-np.power(x-mu,2) / (2*np.power(sigma,2)))


from matplotlib import pyplot
import colorpy.ciexyz
import colorpy.colormodels

fig = pyplot.figure()
ax = fig.add_subplot(1,1,1)

aperture_radius = 3e-3
focal_distance = 3e-2
n_mins = 5
airy_minimum = np.arcsin(0.61 * 545e-9 / aperture_radius)
r_min_max = n_mins * airy_minimum

n_rs = 10000

rs = np.linspace(-r_min_max, r_min_max, n_rs)



wavelengths = range(380,700)
spectrums = np.zeros((len(wavelengths), n_rs))

for index, wavelength in enumerate(wavelengths):
  airy = airy_pattern(rs, 1, wavelength=wavelength * 1e-9)

  spectrums[index] = airy

  ax.plot(rs, airy, color=[i/255.0 for i in colorpy.colormodels.irgb_from_xyz(colorpy.ciexyz.xyz_from_wavelength(wavelength))])


for index, r in enumerate(rs):

  spectrum = np.array(zip(wavelengths, spectrums[:,index]))

  colour = [i/255. for i in colorpy.colormodels.irgb_from_xyz(colorpy.ciexyz.xyz_from_spectrum(spectrum))]

  ax.scatter(r, 1, color=colour)


ax.plot(rs, gaussian(rs, sigma=airy_minimum/2.0)+0.0015, color=(1,1,1), lw=4)
ax.plot(rs, gaussian(rs, sigma=airy_minimum/2.0)+0.0015, color=(0,0,0), lw=2)

ax.set_xlim(-r_min_max, r_min_max)
ax.set_ylim(0, 1.01)

pyplot.show()