Caustics on droplets on glasses formed by streetlights Suppose you're out at night, and it's rainy, and your glasses are covered by water droplets, and you chance to look at a streetlight. 


*

*Are the ridged caustics of light seen at the edges of these droplets described by the Airy functions? 

*Is there a good way of simulating their appearance?

 A: Good question! There is an answer, but it's not as simple/satisfying as "yes they do!", and is pretty involved. There may be another way to do this, but this is how I approach similar problems.
Whether or not they are described well by the Airy functions will depend on the exact geometry of each of the droplets/splats and also on the amplitude and phase of the light passing through them, but having had a quick look I think that the Airy functions may well describe them nicely if you assume uniform illumination and plane wave incident light. You can try the following in python/matlab if you like, and it might give you some idea about whether or not you're right.
Near-field interference/propagation effects (i.e., Fresnel number > 1) are in my opinion best simulated not with the Fresnel integrals but by using the angular spectrum method (ASM). This involves splitting the wavefront into direction cosines and propagating each individually before recombining them at the screen. This sounds a lot more complicated than it is, but really it is very close to being the shift theorem in 2D.
\begin{equation}
U(x,y,z) = \frac{1}{\lambda^{2}}\iint\tilde{U}\left(\frac{\alpha}{\lambda},\frac{\beta}{\lambda}, 0\right)\exp\left[-ikz\sqrt{1-\alpha^{2}-\beta^{2}}\right]\exp[ik(\alpha x + \beta y)]d\alpha d\beta
\end{equation}
Where $U(x,y,z)$ is the field distribution at the screen, and $\tilde{U}(\frac{\alpha}{\lambda},\frac{\beta}{\lambda}, 0)$ is the field distribution at the object plane in terms of direction cosines. I mentioned before that this is effectively the shift theorem. This can be seen by removing $\exp\left[-ikz\sqrt{1-\alpha^{2}-\beta^{2}}\right]$, which just gives an expression for $U$ in terms of its Fourier transform. The the term we removed describes a propagation -- or shift -- along the $z$ axis and through a solid angle given by integration limits on $\alpha$ and $\beta$.
The ASM is implemented numerically as follows:
1) Define your complex field (i.e., amplitude distribution multiplied by the approximate amplitude and phase profile of one splat).
2) FFT the data from real space to spatial frequency space.
3) Calculate the propagation phase using the following (matlab):
k = 2*pi/wavelength
dx = x_range/(N-1) % N is # of grid points, x_range is screen size
kX = (2*pi/dx)*((1:N)/N - (N/2 +1)/N);
[kx, ky] = meshgrid(kx, kx);
mu = sqrt(k^2 - (kx.^2 + ky.^2));
phase = exp(-1i * z *mu);
4) Multiply the FFT of your field by this phase
5) Inverse FFT for the propagated field in real space.
This approach is valid into the far field (Fresnel number <1) but you need an increasingly large grid size if the propagation distance is large so it eventually becomes impractical.
Now that you have your field data, you can define N Airy functions in 2D (scaling these correctly will be tricky, but can be done with some approximation). If you reshape these and your field into column vectors, you can then take the inner product of your field data with each Airy function individually. The result will be a set of projection amplitudes which tells you the component of each Airy mode which is present in your field data. If you reshape these projection amplitudes into a matrix, calculating the sparsity of this matrix will tell you how efficiently (or how well) the Airy functions describe the field data: the more sparse, the better the fit.
