How do I draw a DensityPlot in polar coordinates?

Question

I am trying to plot $J_m(\chi_{mn} r) \cos(m\phi),$ where $J_m(\chi_{mn}) =0,$ as a DensityPlot in polar coordinates (solution to 2D wave equation in polar coordinates). First, some definitions

besselzero = N@Table[BesselJZero[l, n], {l, 0, 25}, {n, 1, 25}];
bzeros[l_, n_] := besselzero[[l + 1, n]];

for speed purposes.

For the purposes of this question consider $m=3$ (but there's essentially the same problem for any $m$):

DensityPlot[
 Cos[3 ArcTan[y/x] ] BesselJ[3, bzeros[3,1]Sqrt[x^2+y^2]],{x,-1,1},{y,-1,1},
 Epilog->{Thick,Circle[]},PlotRange->{-1,1},
 RegionFunction->Function[{x,y,z},Norm[{x,y}]<1.]]

This is nearly right, except that the sign on the left-hand side ($y<0$) of the plot is wrong (it should not be mirror symmetric around $y=0$, rather continue alternating positive and negative after each $\Delta\phi=\pi/3$).

Looking at the Plot3D of ArcTan[y/x] I realized the cause of the problem was because of how ArcTan is defined, so I tried a workaround:

DensityPlot[
 Cos[3 Piecewise[{{0, x > 0 && y == 0}, {π/2, 
       x == 0 && y > 0}, {(3 π)/2, x == 0 && y < 0}, {(2 π)/2,
        x < 0 && y == 0}, {ArcTan[y/x], 
       x > 0 && y > 0}, {ArcTan[y/x] + π, 
       x < 0 && y > 0}, {ArcTan[y/x] + π, 
       x < 0 && y < 0}, {ArcTan[y/x] + 2 π, 
       x > 0 && y < 0}}] ] BesselJ[3, 
   bzeros[3, 1] Sqrt[x^2 + y^2]], {x, -1, 1}, {y, -1, 1},
 Epilog -> {Thick, Circle[]}, PlotRange -> {-1, 1},
 RegionFunction -> Function[{x, y, z}, Norm[{x, y}] < 1.]]

This is nearly right, except for those white regions. Is there a better way to plot functions of this sort?

Answer 1

Use the two-argument form ArcTan[x, y] instead of ArcTan[y/x]:

ArcTan[x,y] gives the arc tangent of y/x, taking into account which quadrant the point $(x,y)$ is in.

It just works:

Answer 2

Quickly, now, which values are positive and which are negative?

Even when a legend may be present, it can help to have additional visual cues concerning the values in a density plot. ReliefPlot will do this with a "hillshading" algorithm. It also reverses the default colors, because it is more natural to think of light colors as being high or closer to the viewer.

To overcome the technical computing problems, isolate the calculations of the arctangent and the Bessel function:

ClearAll[arctan, besselj];
arctan[0., 0.] = 0;
arctan[x_, y_] := ArcTan[x, y];
besselj[m_, \[Chi]_, x_] /; x > 1 := 0;
besselj[m_, \[Chi]_, x_] := BesselJ[m, \[Chi] x];

This takes care of the special case at the origin and expedites the calculations beyond the unit circle (which are not wanted).

Because ReliefPlot does not appear to have the analog of RegionFunction, a reasonable workaround is to mask the undesired areas:

mask = RegionPlot[x^2 + y^2 > 1, {x, -1, 1}, {y, -1, 1}, 
   PlotStyle -> White, BoundaryStyle -> None];

Finally, we need to remember that ReliefPlot builds the image from the bottom to the top. To generate its array, then, we need to go in reverse order of the first coordinate. Whence the work is done as

Block[{m = 3, n = 1, \[Chi], data},
 \[Chi] = BesselJZero[m, n];
 data = ParallelTable[ 
   besselj[m, \[Chi], Sqrt[x^2 + y^2]] Cos[m arctan[x, y]], 
   {x, 1, -1, -.01}, {y, -1, 1, .01}];
 Show[
  ReliefPlot[data, DataRange -> {{-1, 1}, {-1, 1}}, 
   Axes -> {True, True}, AxesOrigin -> {-1.05, -1.05}, Frame -> False],
  mask
  ]
 ]

The improvements are subtle but helpful: it is now visually clear where the peaks and valleys are located. If the effect is not strong enough, exaggerate it by lowering the second LightingAngle value (the help page shows how).

Answer 3

As to the white lines:

Those are caused by small gaps in the plot introduced by the Exclusions option, which is trying to avoid jumps caused by discontinuities in your function. (In this case it's because your function uses Piecewise)

You can turn them off if necessary by using Exclusions -> None:

DensityPlot[
 Cos[3 Piecewise[{
          {0, x > 0 && y == 0}, 
          {Pi/2, x == 0 && y > 0}, 
          {(3 Pi)/2, x == 0 && y < 0}, 
          {(2 Pi)/2, x < 0 && y == 0},        
          {ArcTan[y/x], x > 0 && y > 0}, 
          {ArcTan[y/x] + Pi, x < 0 && y > 0}, 
          {ArcTan[y/x] + Pi, x < 0 && y < 0}, 
          {ArcTan[y/x] + 2 Pi, x > 0 && y < 0}
         }]] BesselJ[3, bzeros[3, 1] Sqrt[x^2 + y^2]], 
 {x, -1, 1}, {y, -1, 1}, 
 Epilog -> {Thick, Circle[]}, PlotRange -> {-1, 1}, 
 RegionFunction -> Function[{x, y, z}, Norm[{x, y}] < 1.], 
 Exclusions -> None]

Answer 4

The difficulties with either single-argument Arctan[y/x] or two-argument ArcTan[x,y] (the latter being undefined at the origin, as was addressed by @whuber) were also discussed on MathGroup. Heike's suggestion there was to use Arg[x + I y] to get the polar angle. But in fact it turns out to be slightly more efficient (and perhaps more readable) to use the function CoordinatesFromCartesian in the VectorAnalysis package:

Needs["VectorAnalysis`"]

SetCoordinates[Cylindrical];

DensityPlot[
 Apply[Cos[3 #2] BesselJ[3, bzeros[3, 1] #1] &, 
  CoordinatesFromCartesian[{x, y, 0}]], {x, -1, 1}, {y, -1, 1}, 
 Epilog -> {Thick, Circle[]}, PlotRange -> {-1, 1}, 
 RegionFunction -> Function[{x, y}, Norm[{x, y}] < 1.]]

(there was an unnecessary z argument in the RegionFunction of your example)

One final remark, related to your setup with BesselJZeros: instead of preparing them in a table, you may want to use the fact that BesselJZeros threads over lists! That makes the initial calculation faster than with a table. Of course you only do that loop once in this application, but I thought it's worth pointing out for the benefit of more complex applications.