0
$\begingroup$

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:

enter image description here

$\endgroup$
  • $\begingroup$ 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". $\endgroup$ – Henrik Schumacher Jan 12 '19 at 9:04
  • 1
    $\begingroup$ 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. $\endgroup$ – Bill Jan 12 '19 at 9:15
  • $\begingroup$ I'm sorry ,I have changed the error. $\endgroup$ – qb.Suo Jan 12 '19 at 12:05
1
$\begingroup$

(Not an answer, extended comment.)

Please experiment with your sum, say, like this:

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.