Polar plotting Hankel Function with a lot of terms

Question

I am trying to plot a normalized polar plot for the following function with different values of $a$

$$\left\lvert \sum_{n=1}^\infty i^n (2n+1) \frac {P_n^1(cos(\theta))}{\sqrt{\frac{\pi k a}{2}}[-H_{n+\frac{3}{2}}^2 (ka) + \frac{n+1}{ka}H_{n+\frac{1}{2}}^2(ka)]} \right\rvert^2 $$

where $P_n^1(cos(\theta))$ is the associated Legendre polynomial and $H_{n+\frac{3}{2}}^2 (ka)$ and $H_{n+\frac{1}{2}}^2(ka)$ are Hankel function of 2nd kind. Here $k=2\pi$ and the value of $a$ vary as $a = [.25, .05, 2, 10, 20]$

I can get plots upto $a=2$ with $n=100$ but for $a=10, 20$ I am having difficulty plotting. I am running into issues where machine precision is lost and the norm for Hankel function becomes too big. This is my attempt below

k = 2 \[Pi] ;
Pr[a_, m_, \[Theta]_] := 
Abs[Sum[I^n (2 n + 1) LegendreP[n, 1, Cos[\[Theta]]]/(
 Sqrt[(\[Pi] k a )/
   2] (-HankelH2[n + 3/2, k a ] + (n + 1)/(k a)
      HankelH2[n + 1/2, k a ])), {n, 1, m}]]^2;

ap05m = 2; normp05 = FindMaximum[Pr[.05, ap05m, \[Theta]], {\[Theta], 0, 2 \[Pi]}];
ap25m = 7 ; normp25 = FindMaximum[Pr[.25, ap25m, \[Theta]], {\[Theta], 0, 2 \[Pi]}];
a2m = 100; norm2 = FindMaximum[Pr[2.0, a2m, \[Theta]], {\[Theta], 0, 2 \[Pi]}];
a10m = 150; norm10 = FindMaximum[Pr[10.0, a10m, \[Theta]], {\[Theta], 0, 2 \[Pi]}];
a20m = 150; norm20 = FindMaximum[Pr[20.0, a20m, \[Theta]], {\[Theta], 0, 2 \[Pi]}];

PolarPlot[{Pr[.05, ap05m,\[Theta]]/normp05[[1]], 
           Pr[.25, ap25m, \[Theta]]/normp25[[1]],
           Pr[2, a2m, \[Theta]]/norm2[[1]], 
           Pr[10, a10m, \[Theta]]/norm10[[1]], 
           Pr[20, a20m, \[Theta]]/norm20[[1]]}, {\[Theta], 0, 2 \[Pi]},
PolarAxes -> True, PlotRange -> Automatic, 
PolarGridLines -> Automatic, PolarTicks -> {"Degrees", Automatic}, 
PolarAxesOrigin -> {0, 1}, PlotLegends -> "Expressions"]

The norm for $a=10$ and $a=20$ are

a10m = 150; norm10 = FindMaximum[Pr[10.0, a10m, \[Theta]], {\[Theta], 0, 2 \[Pi]}]

${7.24673*10^{27}, {\theta -> 0.121351}}$

a20m = 150; norm20 = FindMaximum[Pr[20.0, a20m, \[Theta]], {\[Theta], 0, 2 \[Pi]}]

${1.11063*10^{79}, {\theta -> 0.0809369}}$

These normalization factors are too big and as a result I don't see the graphs for $a=10$ and $a=20$

This is what the graph is supposed to look like given by the professor

This is what my current graph looks like

I have tried $Chop[]$ function, but it didn't work. I would appreciate any help in plotting the normalized graph for $a = 10, 20$ on the same graph with $a = .05, .25, 2$.

Answer 1

I am running into issues where machine precision is lost and the norm for Hankel function becomes too big.

Exactly. So you already know the answer: Don't use machine precision but ensure your computations are numerically sound by giving Mathematica enough digits to work with. I believe this answer about controlling precision is a good start.

In general, you have to understand that as soon as you say, e.g. 0.25, you told Mathematica that this is in machine precision. As long as you stick to infinitely accurate numbers like 1/4, you can tell Mathematica how precise it should evaluate numbers in certain algorithms later.

That being said, in the following I use ListPolarPlot and create the list of values myself. The adaptive plotting algorithms of Mathematica usually include many points to ensure your curve is smooth and calculating with higher precision takes a lot of time. We can speed this up by using ParallelTable but it will still require some minutes. I use the definitions of k and Pr you have already given.

as = {5/100, 1/4, 2, 10, 20};
ams = {2, 7, 100, 150, 150};
Block[{$MaxExtraPrecision = 100},
 norms = First@
     FindMaximum[Pr[##, \[Theta]], {\[Theta], 0, 2 \[Pi]}, 
      WorkingPrecision -> 30] & @@@ Transpose[{as, ams}]
 ]

First, note that I converted all your numbers to rationals and only in the final call to FindMaximum, I ask Mathematica to use a working precision of 30. Also, note that during the computation, Mathematica needs more precision and that's why we set $MaxExtraPrecision. If you don't set it, Mathematica will tell you that this setting is too low in a warning message.

Now, we can create the list of radii. For each curve, one list of radii, and we do this computation in parallel as it takes some minutes. With ParallelEvaluate, we distribute the setting of $MaxExtraPrecision to the parallel kernels.

ParallelEvaluate[$MaxExtraPrecision = 100];
data = ParallelTable[
     Pr[#1, #2, \[Theta]]/#3, {\[Theta], 0, 2 Pi, Pi/80}] & @@@ 
   Transpose[{as, ams, norms}];

As you can see, in above code blocks, I used @@@ and ## (or #1). Familiarize yourself how this works by looking at this example

f[##, {#2}] & @@@ {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}

and checking the documentation. Don't just copy the code.

After that, you can plot your data and recreate what your teacher had:

ListPolarPlot[data, Joined -> True, 
  PolarAxes -> True,
  PlotRange -> Automatic, 
  PolarGridLines -> Automatic, 
  PolarTicks -> {"Degrees", Automatic},
  PolarAxesOrigin -> {0, 1},
  PlotLegends -> (StringTemplate["a = ``\[Lambda]"] /@ N[as]), 
  PlotStyle -> {Blue, Red, Green, Magenta, Black}]