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 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 at 9:15
  • $\begingroup$ I'm sorry ,I have changed the error. $\endgroup$ – qb.Suo Jan 12 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 at 2:32

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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