Hi Alex,

This is an interesting problem. The airy function has a singularity at 0 which can be removed by defining the value there as 1, which is the limit. I do that in the code below, and also make use of the radial symmetry to fit for only the central amplitude and the center. (Which is not all that is needed.)

I have found that the code below works intermittently. It seems still dependent on choices of values made during the fit function, even though I think I have removed the singularity. I Hopefully someone with more knowledge will speak up.

In[1]:= (2 BesselJ[1, r]/r)^2 /. r -> 0

During evaluation of In[1]:= Power::infy: Infinite expression 1/0^2 encountered. >>

During evaluation of In[1]:= Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered. >>

Out[1]= Indeterminate

In[2]:= Limit[(2 BesselJ[1, r]/r)^2, r -> 0]

Out[2]= 1

In[3]:= airy[r_] := (2 BesselJ[1, r]/r)^2 /; r != 0.

In[4]:= airy[0.] = 1

Out[4]= 1

In[5]:= airyData = 
  Flatten[Table[{x, y, 5 airy[Sqrt[(x - 1)^2 + (y - 2)^2]]}, {x, -10, 
     10, .1}, {y, -10, 10, .1}], 1];

In[6]:= ListPlot3D[airyData];

In[7]:= NonlinearModelFit[airyData, 
 a airy[Sqrt[(x - h)^2 + (y - k)^2]], {a, h, k}, {x, y}]