# Problems with numerical integration

I tried to plot beam intensity using a function that evaluates a numerical integral, but it didn't work.

Here is my code, which did not produce a result.

a1 = 0.6328*10^-6;
k = (2*π)/a1;
w = 0.02;
c = 0.02;
m = 1;
a = 800;
z = 1000;

M1[x_, r1_, r2_] =
NIntegrate[
Sum[k^2*z^(-2)*BesselJ[m, a *r1]*BesselJ[m, a* r2]*
BesselJ[b, k* r1* x/z]*BesselJ[b, k* r2* x/z]*
BesselI[m + b, 2*r1*r2*c^(-2)]*
Exp[-(c^(-2) + w^(-2))*(r1^2 + r2^2) + (I*k) /
(2*z)*(r2^2 - r1^2)]*r1*r2,
{r1, 0, ∞}, {r2, 0, ∞}],
{b, -∞, ∞}];
a1 = Table[{x, Abs[M1[x]]}, {x, -0.02, 0.02, 0.001}];
ListPlot[a1, PlotRange -> All] TimeUsed[]


Here is how appears in my notebook:

• Welcome on Mathematica.StackExchange. Please always provide copyable code in InputForm. This can be done by: (i) marking the code to copy (ii) right click (iii) selecting "Copy As" -> "Input Text". – Henrik Schumacher Jan 12 at 9:04
• Your first problem is that you define M1[x_,r1_,r2_]= but then you use M1[x]. I am guessing you want to change your function definition to M1[x_]:= so that your r1 and r2 don't disagree with your r1 and r2 variables of integration. – Bill Jan 12 at 9:15
• I'm sorry ,I have changed the error. – qb.Suo Jan 12 at 12:05

AbsoluteTiming[
Block[{b = 1000, s = 100},
Total@Flatten@
Table[k^2*z^(-2)*BesselJ[m, a*r1]*BesselJ[m, a*r2]*
BesselJ[b, k*r1*x/z]*BesselJ[b, k*r2*x/z]*
BesselI[m + b, 2*r1*r2*c^(-2)]*
Exp[-(c^(-2) + w^(-2))*(r1^2 + r2^2) + (I*k)/(2*z)*(r2^2 -
r1^2)]*r1*r2, {r1, 0, s}, {r2, 0, s}]]
]

(* {0.35228, 0. + 0. I} *)


I get zeroes for all 5-6 values of b I tried.

• Thank you.this expression was from a paper,and I think the author may do some approximation. – qb.Suo Jan 13 at 2:32