Matplotlib - 投影された極座標での等高線と震度プロット

Question

（r、theta）座標の不均一なグリッド上に定義されたスカラーフィールドとベクトルフィールドの輪郭と震度プロットをプロットする必要があります。

私が持っている問題の最小の例として、磁気双極子ためStream functionの等高線プロットを考える、そのような機能の輪郭は、（この場合、磁界）corresponedingベクトル場の流線です。

以下のコードは、（r、theta）座標で不均一なグリッドを取り、それをデカルト平面にマッピングし、ストリーム関数の等高線図をプロットします。何らかの理由で

import numpy as np 
import matplotlib.pyplot as plt 

r = np.logspace(0,1,200) 
theta = np.linspace(0,np.pi/2,100) 

N_r = len(r) 
N_theta = len(theta) 

# Polar to cartesian coordinates 
theta_matrix, r_matrix = np.meshgrid(theta, r) 
x = r_matrix * np.cos(theta_matrix) 
y = r_matrix * np.sin(theta_matrix) 

m = 5 
psi = np.zeros((N_r, N_theta)) 

# Stream function for a magnetic dipole 
psi = m * np.sin(theta_matrix)**2/r_matrix 

contour_levels = m * np.sin(np.linspace(0, np.pi/2,40))**2. 

fig, ax = plt.subplots() 
# ax.plot(x,y,'b.') # plot grid points 
ax.set_aspect('equal') 
ax.contour(x, y, psi, 100, colors='black',levels=contour_levels) 
plt.show() 

しかし、私が手プロットは、右を見ていない：

私は輪郭関数呼び出しでxとyを入れ替える場合は、私は望ましい結果を得る：

を

ベクトルフィールドの震度プロットを同じグリッドに定義し、xy平面にマッピングしようとすると、同じことが起こりますが、関数呼び出しのxとyを交換することはもう機能しません。

私はどこかで愚かな間違いをしたようですが、私はそれが何であるか理解できません。

Answer 1

psi = m * np.sin(theta_matrix)**2/r_matrix の場合、θが0からπ/ 2になるとpsiが増加し、rが増加するとpsiは減少します。

したがって、シータが増加するにつれて、psiの等高線はrで増加するはずです。その結果、 は、中心から放射するときに反時計回りに進む曲線になります。これは、あなたが投稿最初のプロットと一致 あり、その結果は、結果の妥当性を確認するための別の方法は、表面プロットを見ることである

ax.contour(x, y, psi, 100, colors='black',levels=contour_levels) 

---

であなたのコードの最初のバージョンで返されますpsiの：

import numpy as np 
import matplotlib.pyplot as plt 
import mpl_toolkits.mplot3d.axes3d as axes3d 

r = np.logspace(0,1,200) 
theta = np.linspace(0,np.pi/2,100) 

N_r = len(r) 
N_theta = len(theta) 

# Polar to cartesian coordinates 
theta_matrix, r_matrix = np.meshgrid(theta, r) 
x = r_matrix * np.cos(theta_matrix) 
y = r_matrix * np.sin(theta_matrix) 

m = 5 

# Stream function for a magnetic dipole 
psi = m * np.sin(theta_matrix)**2/r_matrix 

contour_levels = m * np.sin(np.linspace(0, np.pi/2,40))**2. 

fig = plt.figure() 
ax = fig.add_subplot(1, 1, 1, projection='3d') 
ax.set_aspect('equal') 

ax.plot_surface(x, y, psi, rstride=8, cstride=8, alpha=0.3) 
ax.contour(x, y, psi, colors='black',levels=contour_levels) 
plt.show()