# Polar plotting Hankel Function with a lot of terms

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[],
Pr[.25, ap25m, \[Theta]]/normp25[],
Pr[2, a2m, \[Theta]]/norm2[],
Pr[10, a10m, \[Theta]]/norm10[],
Pr[20, a20m, \[Theta]]/norm20[]}, {\[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$$.

• The problem with your code is that you recompute the prefactor for each value of angle again and again. The idea would be to precompute the prefactors once and store them in the table. Or, if a minimal code modification is required you can use memoization option of MA reference.wolfram.com/language/tutorial/…. – yarchik Oct 11 '19 at 8:21
• @yarchik that's a great idea! I didn't know that! However, i still don't get the graphs for a = 10 and a = 20. – Rumman Oct 11 '19 at 21:50

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.

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}] • Thanks a lot for your response. I am going reiterate your explanation in my own words so that you can correct/check my understanding of your explanation. -Mathematica by default uses Machine Precision if the values are using decimal – Rumman Oct 15 '19 at 3:47
• - This was giving me precision issue. So you changed the list of radius values to "infinitely accurate precision" by using "/" and you controlled the precision later on in the calculation - I understand your notations. I played around with "@@@" and "##", "#1" Thanks for the example. That was helpful. - You used ParallelTable to evaluate the 5 sets of data – Rumman Oct 15 '19 at 4:02
• @Rumman Yes, exactly. There are different ways of controlling precision/accuracy (see 0.2520 or 0.2520 notation) but you can read about them in several posts here or in the documentation. Just look up Precision, Accuracy, N, SetAccuracy, ... and the standard options WorkingPrecision, AccuracyGoal`, and the links therein. – halirutan Oct 15 '19 at 6:45