是否可以在 Mathematica 中创建极坐标 CountourPlot/ListCountourPlot/DensityPlot？

Question

我正在寻找类似回音壁模式的东西——极坐标中的 2D 圆柱对称图。像这样的：

我在 Trott 的符号指南中找到了以下代码 sn-p。尝试在一个非常小的数据集上运行它；它吃掉了 4 GB 的内存并耗尽了我的内核：

(* add points to get smooth curves *)
addPoints[lp_][points_, \[Delta]\[CurlyEpsilon]_] := 
Module[{n, l}, Join @@ (Function[pair,
       If[(* additional points needed? *)
          (l = Sqrt[#. #]&[Subtract @@ pair]) < \[Delta]\[CurlyEpsilon], pair, 
          n = Floor[l/\[Delta]\[CurlyEpsilon]] + 1; 
          Table[# + i/n (#2 - #1), {i, 0, n - 1}]& @@ pair]] /@ 
          Partition[If[lp === Polygon, 
                       Append[#, First[#]], #]&[points], 2, 1])]

(* Make the plot circular *)
With[{\[Delta]\[CurlyEpsilon] = 0.1, R = 10}, 
 Show[{gr /. (lp : (Polygon | Line))[l_] :> 
     lp[{#2 Cos[#1], #2 Sin[#1]} & @@@(* add points *)
       addPoints[lp][l, \[Delta]\[CurlyEpsilon]]], 
   Graphics[{Thickness[0.01], GrayLevel[0], Circle[{0, 0}, R]}]}, 
  DisplayFunction -> $DisplayFunction, Frame -> False]]

这里，gr 是一个矩形 2D ListContourPlot，使用类似这样的东西生成（例如）：

data = With[{eth = 2, er = 2, wc = 1, m = 4}, 
   Table[Re[
     BesselJ[(Sqrt[eth] m)/Sqrt[er], Sqrt[eth] r wc] Exp[
       I m phi]], {r, 0, 10, .2}, {phi, 0, 2 Pi, 0.1}]];
gr = ListContourPlot[data, Contours -> 50, ContourLines -> False, 
  DataRange -> {{0, 2 Pi}, {0, 10}}, DisplayFunction -> Identity, 
  ContourStyle -> {Thickness[0.002]}, PlotRange -> All, 
  ColorFunctionScaling -> False]

有没有一种直接的方法来绘制这样的圆柱图？.. 我很难相信我不得不求助于 Matlab 来满足我的曲线坐标需求:)

Answer 1

以前的 sn-ps 已删除，因为这显然是我想出的最佳答案：

With[{eth = 2, er = 2, wc = 1, m = 4}, 
 ContourPlot[
  Re[BesselJ[(Sqrt[eth] m)/Sqrt[er], Sqrt[eth] r wc] Exp[I phi m]]/. 
                                         {r ->Norm[{x, y}], phi ->ArcTan[x, y]}, 
  {x, -10, 10}, {y, -10, 10}, 
  Contours -> 50, ContourLines -> False, 
  RegionFunction -> (#1^2 + #2^2 < 100 &), 
  ColorFunction -> "SunsetColors"
 ]
]

编辑

将ContourPlot 替换为Plot3D 并删除您获得的不受支持的选项：

Answer 2

这是一个相对简单的问题。关键是，如果你可以参数化它，你就可以绘制它。根据文档，ListContourPlot 和 ListDensityPlot 都接受两种形式的数据：高度值数组或坐标列表加上函数值 ({{x, y, f} ..})。第二种形式更容易处理，即使你的数据是第一种形式，我们也会将其转换为第二种形式。

简单地说，将{{r, t, f} ..} 形式的数据转换为{{x, y, f} ..} 形式的数据，当应用于从BesselJ[1, r/2] Cos[3 t] 获取的数据时，您会得到

如果你只有一个数据数组，比如this guy，那该怎么办？在这种情况下，您有一个二维数组，其中数组中的每个点都具有已知位置，为了绘制它，您必须将其转换为第二种形式。我偏爱MapIndexed，但还有其他方法可以做到这一点。假设您的数据存储在一个数组中，其中行对应于径向坐标，列是角坐标。然后要转换它，我会使用

R = 0.01;    (*radial increment*)
T = 0.05 Pi; (*angular increment*)
xformed = MapIndexed[ 
   With[{r = #2[[1]]*R, t = #2[[1]]*t, f = #1},
   {r Cos[t], r Sin[t], f}]&, data, {2}]//Flatten[#,1]&

给出相同的结果。

---

如果您有解析解，则需要将其转换为笛卡尔坐标，如上所示，但您使用替换规则。例如，

ContourPlot[ Evaluate[
    BesselJ[1, r/2]*Cos[3 t ] /. {r -> Sqrt[x^2 + y^2], t -> ArcTan[x, y]}], 
   {x, -5, 5}, {y, -5, 5}, PlotPoints -> 50, 
   ColorFunction -> ColorData["DarkRainbow"], Contours -> 25]

给予

有两点需要注意：1）Evaluate 用于确保正确执行替换，2）ArcTan[x, y] 考虑到点 {x,y} 所在的象限。