python-3.x - Matplotlib 中的极坐标图，通过映射到笛卡尔坐标

Question

我有一个变量 (P)，它是角度 (theta) 的函数：

在这个等式中，K是一个常数，theta_p等于 0，I是第一类（0 阶）的修正贝塞尔函数，定义为：

现在，我想针对常数K的不同值绘制P与theta的关系。首先，我计算了参数I，然后将其代入第一个方程以计算不同角度 theta 的 P。我通过放置将其映射到笛卡尔坐标中：

x = P*cos(θ)

y = P*sin(θ)

这是我在常数 k=2.0 时使用 matplotlib 和 scipy 的 python 实现：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

def integrand(x, a, k):
   return a*np.exp(k*np.cos(x))

theta = (np.arange(0, 362, 2))
theta_p = 0.0

X = []
Y = []

for i in range(len(theta)):
    a = (1 / np.pi)
    k = 2.0
    Bessel = quad(integrand, 0, np.pi, args=(a, k))
    I = list(Bessel)[0]
    P = (1 / (np.pi * I)) * np.exp(k * np.cos(2 * (theta[i]*np.pi/180. - theta_p)))
    x = P*np.cos(theta[i]*np.pi/180.)
    y = P*np.sin(theta[i]*np.pi/180.)
    X.append(x)
    Y.append(y)

plt.plot(X,Y,  linestyle='-', linewidth=3, color='red')
axes = plt.gca()
plt.show()

对于不同的 K 值，我应该得到一组如下图所示的图表：

（请注意，为了便于可视化，将分布绘制在单元 1 的圆圈上）

但是，上面代码生成的图表似乎与上图不相似。知道上述实现有什么问题吗？在此先感谢您的帮助。

这是它的样子（对于 k=2）：

这些公式的参考是等式 5 和 6，您可以在此处找到

Answer 1

你的公式有误。

您的公式给出了单位圆上方函数的增量。所以在你的函数中得到你想要的情节，只需将 1 添加到它。

这是您想要的，带有一些整理好的python。...请注意，您可以将“P”值的整个计算作为一个 numpy 向量线，您不需要遍历索引。...您也可以直接在 matplotlib 中绘制极坐标图 - 您无需将其转换为笛卡尔坐标。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

theta = np.arange(0, 2*np.pi+0.1, 2*np.pi/100)

def integrand(x, a, k):
   return a*np.exp(k*np.cos(x))

for k in np.arange(0, 5, 0.5):
    a = (1 / np.pi)
    Bessel = quad(integrand, 0, np.pi, args=(a, k))
    I = Bessel[0]    
    P = 1 + (1/(np.pi * I)) * np.exp(k * np.cos(2 * theta))
    plt.polar(theta, P)

plt.show()