Message Boards Message Boards

2 Replies
1 Total Likes
View groups...
Share this post:

fitting the Airy diffraction pattern with FindFit

Posted 10 years ago

Hello! I have 2D data (images) which are results of Gaussian type beam diffraction. enter image description here According to the theory of diffraction the pattern can be discribed by function ~ J1(r)/r. I tried to fit my data with this function using the FIndFit routine, but had an error :

Blockquote Power::infy: "Infinite expression 1/0. encountered" Blockquote

data01 = Import["366.0000cm_660.0000nm.png"];
data02 = ImageData[data01];
min02 = Min[data02];
max02 = Max[data02];

(*Dimensionalities of images*)
dim01a = Import["366.0000cm_660.0000nm.png", "ImageSize"] [[1]];
dim01b = Import["366.0000cm_660.0000nm.png", "ImageSize"] [[2]];
dim01 = Min[dim01a, dim01b];

ReducPix = 5;
data03 = Table[ (data02[[i]][[j]][[1]] - min02)/max02, {i, 1, dim01b, 
    ReducPix}, {j, 1, dim01a, ReducPix}];

dim03a = Dimensions[data03][[1]];
dim03b = Dimensions[data03][[2]];
dim03 = Min[dim03a, dim03b]

ListPlot3D[data03, PlotRange -> All]
s[x_, y_] := Sqrt[a2 (x - a1)^2 + b2 (y - b1)^2]
Bessel3D = Piecewise[
  {{1/4, x == 0 && y == 0},
   {(a0 BesselJ[1, s[x, 0]]/s[x, 0]), y == 0},
   {(a0 BesselJ[1, s[0, y]]/s[0, y]), x == 0}},
  {(a0 BesselJ[1, s[x, y]]/s[x, y])}]

data3D = Flatten[MapIndexed[{#2[[1]], #2[[2]], #1} &, data03, {2}], 1];
fit = NonlinearModelFit[
   data3D, (Bessel3D)^2, {{a0, 1}, {a1, 50}, {a2, 10}, {b1, 50}, {b2, 
     10}}, {x, y}];

I tried to use different method in the function (NMinimize, Gradient, so on ) but it did not help. It is really puzzling for me. Will appreciate any suggestions!!

POSTED BY: Alex Mikh
2 Replies
Posted 10 years ago

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}]
POSTED BY: David Keith
Posted 10 years ago

Hi David! thanks for your interest to this problem! You are right, there is the singularity at zero point. I tried to eliminate it using a Picewise definiton of the model function in my code but for some reason it did no help and I still had the error. Looks like your definition of the function

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

airy[0.] = 1;`

is able to overcome that issue. Great! An addition specification of the range for parameters in NonlinearModelFit helps as well. However when I'm trying to fit my data the result is not accurate. A new error pops up:

NonlinearModelFit::nrlnum: The function value {2.3360210^119-3.4376510^103 I,3.2420410^118+0. I,4.6252810^117+6.1297110^101 I,6.7870110^116-1.18210^101 I,<<44>>,1.1806110^96+0. I,1.15086*10^96+0. I,<<12238>>} is not a list of real numbers with dimensions {12288} at {a0,a1,a2,b1,b2} = {0.206555,69.9985,-2.9634,50.0133,-2.99061}. >>

not sure what does it mean.. The 3D plots with parameters from the fitting are not very similiar to the original data at all and as you mentioned very sensitive to the initial values (guess) in the fitting. Thank you anyway!

Thanks for your suggestions!

POSTED BY: Alex Mikh
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
or Discard

Group Abstract Group Abstract