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

Question

我想绘制类似于耳语画廊模式的东西 - 在极坐标下的二维圆柱对称图。像这样的东西:

我在Trott的符号指南中找到了以下代码片段。尝试在一个非常小的数据集上运行它；它吃掉了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来满足我的曲线坐标需求 :)

- Leo Alekseyev

---

既然您更喜欢解析解而不是数值解，即ContourPlot与ListContourPlot，那么您应该更改标题以反映这一事实。 - rcollyer

Answer 1

之前的代码片段已删除，因为这明显是我想到的最佳答案：

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并删除不支持的选项，你会得到：

- Dr. belisarius

---

@belisarius 我把 DensityPlot 改成了 ContourPlot，因为它对我来说运行速度快得多，并且产生类似的结果。两者都可以正常工作。 - Leo Alekseyev

@belisarius，我会用Evaluate包装Replace以避免在没有时出现Plot问题。 - rcollyer

@rcollyer 我已经阅读了你在答案中的评论，我认为它是正确的。但在这种情况下似乎并不必要。我真的没有仔细考虑过，但如果你认为那里存在潜在的问题，请随意编辑我的代码。关于何时需要以及何时不需要的解释也将是很好的。 - Dr. belisarius

@belisarius，好问题。我自己也不确定，尽管Plot本身似乎是最糟糕的罪犯。 - rcollyer

Answer 2

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

简单地说，要将形式为{{r，t，f} ..}的数据转换为形式为{{x，y，f} ..}的数据，您可以执行N[{#[[1]] Cos[ #[[2]] ]，#[[1]] Sin[ #[[2]] ]，#[[3]]}]& /@ data，当应用于从BesselJ[1，r/2] Cos[3 t]获取的数据时，你会得到

这个图像： 那么如果您只有一个数据数组呢？例如这个人所说的那样？在这种情况下，您有一个二维数组，其中数组中的每个点都有已知位置，为了绘制它，您必须将其转换为第二种形式。我倾向于使用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}所在象限。

- rcollyer

---

数据可以通过以下方式生成：data = With[{eth = 2, er = 2, wc = 1, m = 4}, Flatten[#, 1] &@ Table[{r, phi, Re[BesselJ[(Sqrt[eth] m)/Sqrt[er], Sqrt[eth] r wc] Exp[ I m phi]]}, {r, 0, 10, .15}, {phi, 0, 2 Pi, 0.035}]] 这是一个合理的方法，可以给您很多控制权，但不像belisarius的最终答案那样简洁。 - Leo Alekseyev