Mathematica/極座標曲面圖

沒有內建函式可以創建極座標3D曲面圖（即高度由半徑和角度控制）。但是，少量程式碼允許我們使用笛卡爾座標系和極座標系之間的恆等式來繪製一個。

$x=r\cos \left(\theta \right)\,$

$y=r\sin \left(\theta \right)\,$

我們的做法是構建一個值表，對應於一系列角度和半徑的高度。讓我們使用下面的函式作為示例

$f\left(r,\theta \right)=r\sin \left(\theta \right)\,$

構建表格的第一步是定義函式並說明我們希望的繪圖點之間的距離。

dtheta = Pi/20;                          (*Give a radial gridline spacing of Pi/20 radians*)
rmax = 1;                                (*Define the maximum radius*)
dr = rmax/10;                            (*Give 10 circumferential grid lines*)

f[r_, theta_] := r  Sin[theta];          (*This is the function definition*)

現在，我們根據以上資訊構建高度值的表格。

data = Table[
      f[r, theta],                       (*This is the function giving the value in the table*)
      {theta, 0, 2Pi, dtheta},           (*This is the increment for theta*)
      {r, 0, rmax, dr}];                 (*This is the increment for r*)

此時，我們可以使用ListPlot3D在笛卡爾座標系中繪製資料表，其中r在x軸上，θ在y軸上，f(r,θ)在z軸上。這將生成一個普通的SurfaceGraphics輸出。

gr1 = ListPlot3D[
    data,                                (*The array of height values*)
    DataRange -> {{0, rmax}, {0, 2*Pi}}];(*The range that this array covers*) 

下圖顯示了SurfaceGraphics輸出gr1。請注意，在此示例中，影像已進行抗鋸齒處理。有關抗鋸齒程式碼，請參閱影像頁面。

現在我們將SurfaceGraphics物件轉換為Graphics3D物件gr2，

gr2 = Graphics3D[gr1];

這是我們使用上述關係將笛卡爾圖轉換為極座標圖的點。我們首先定義一個規則，用於對由三個數字列表（形式為{x,y,z}）給定的點執行變換。我們將它稱為“替換”。

substitution = {r_, theta_, z_} -> {r Cos[theta], r Sin[theta], z};

現在，我們使用ReplaceAll函式遍歷Graphics3D物件gr2，並找到每個多邊形點。我們使用“：”運算子將多邊形點（三個數字的列表）設定為模式物件，併為其分配名稱“p”。現在我們使用RuleDelayed（：>）對每個點執行上述變換。我們不能使用Rule（->），因為這隻會在其輸入時進行評估，而不是在其使用時進行評估，並且不會變換點。

gr3 = ReplaceAll[gr2, p : Polygon[pts_] :> ReplaceAll[p,substitution] ]

最後一步是顯示生成的極座標圖gr3

Show[
  gr3,
  AxesLabel -> {"", "", z}]             (*Retitle the axes*)

由於我們不再使用笛卡爾座標系，因此我們無法將水平軸命名為x和y，但是，我們沒有更改z的值，因此我們可以將其保留為垂直軸的名稱。

對於Mathematica 6.0，請使用以下程式碼

MyListPolarPlot3D[data_, rRange_, thetaRange_, zRange_] :=
 Module[
  {},
  gr1 = ListPlot3D[
    data,
    DataRange -> {
      { rRange[[1]], rRange[[2]]},
      { thetaRange[[1]], thetaRange[[2]]}
      },
    DisplayFunction -> Identity,
    ColorFunction -> "SolarColors",
    ColorFunction -> Automatic,
    MeshFunctions -> {Function[{x, y, z}, x*Cos[y]], 
      Function[{x, y, z}, x*Sin[y]]},
    BoundaryStyle -> None,
    ColorFunctionScaling -> True,
    Mesh -> 30
    ];
  
  substitution = {r_, theta_, z_} -> {r Cos[theta], r Sin[theta], z};
  
  gr2 = gr1 /. 
    GraphicsComplex[p_List, rest__] :> 
     GraphicsComplex[ReplaceAll[p, substitution], rest] ;
  
  (*
  * Retitle the axes and show final graph
  *)
  Return[
   Show[
    gr2,
    AxesLabel -> {"X", "Y", "Z"},
    DisplayFunction -> $DisplayFunction ,
    BoxRatios -> {1, 1, 0.8},
    PlotRange -> {
      {-0.65*rRange[[2]], 0.65*rRange[[2]]},
      {-0.65 rRange[[2]], 0.65*rRange[[2]]},
      {zRange[[1]], zRange[[2]]}
      }
    ]
   ] ;  
  ]

這是一個示例

Fun[r_, t_] := 
  0.632 (0.710 \[ExponentialE]^(-1.166 (0.492+ r^2 - 
         1.403 r Cos[t])) + 
     0.710 \[ExponentialE]^(-1.166 (0.492+ r^2 + 1.403 r Cos[t]))) ;

dataPlot = 
  Table[ Fun[r, t], {t, 0.0, 2.0*Pi, 2*Pi/100}, {r, 0.0, 4.0, 
    0.08}] ;
MyListPolarPlot3D[ dataPlot, {0.0, 4.0}, {0.0, 2*Pi}, {0, 0.55}]