function y = airy_func(lambda, diam, theta)
//
// Compute the Airy Function for light passing through a circular aperture.
//
//                                               2
//                  pi * diam^2   [  2 * J1(x)  ]
//     I(theta)  =  ------------  [  ---------  ]
//                  4 * lambda^2  [     x       ]
//
// where
//                     x    =  pi * diam * theta / lambda
//
//                   J1(x)  =  Bessell function of the first kind, order 1 
//
// Arguments:
//                 lambda    (input)     wavelength of light (meters)
//
//                 diam      (input)     diameter of aperture (meters)
//
//                 theta     (input)     angle from center of diffraction
//                                          pattern at which to evaluate
//                                          intensity (radians)
//
//                 y         (output)    normalized intensity of light
//                                          at that angle (= 1.0 at center)
//
// MWR 4/9/2007

  verbose = 1;

  // sanity checks
  if (lambda <= 0) 
    mprintf("airy_func: wavelength %9.4e must be positive \n", lambda);
    error("airy_func: quitting now \n");
  end
  if (diam <= 0) 
    mprintf("airy_func: aperture diameter %9.4e must be positive \n", diam);
    error("airy_func: quitting now \n");
  end
  if (theta < 0) 
    mprintf("airy_func: angle %9.4e must be positive \n", theta);
    error("airy_func: quitting now \n");
  end

  // in theory, the Airy function has this amplitude at center
  //        theoretical amp = (pi*diam*diam)/(4*lambda*lambda);
  // but we are computing a normalized value, so we'll set it equal to 1
  amp = 1.0;

  // special case for very center of Airy pattern, theta = 0 --
  //   we return the amplitude (which we are normalizing)
  if (theta == 0)
    y = amp;
  else
    // this auxiliary quantity simplifies the calculation
    pi = 3.14159265358979;
    x = (pi*diam*theta) / (lambda);
    xx = (2.0*besselj(1.0, x)) / (x);
    y = amp*xx*xx;
  end

endfunction
