解决方案

有几个不同的事情，你需要改变，以修复你的代码。主要问题是，当你绘制你使用的数据不转换为carnival坐标。ax.scatter(...)函数采取的位置参数为x,y,z而不是r,θ,z。此外，你的θ的边界是错误的，需要移动几行代码来使它全部工作。

第一步

你应该做的第一件事是将f = np.sqrt(...)移动到def equation(...)函数中。这是因为f的值会随着theta的值而变化，而theta在函数中并不是恒定的。

## Incorrect
f = np.sqrt(np.math.factorial(p)/(np.math.factorial(m+p)))*((np.sqrt(2)*r/r0)**m)*(2*r**2/r0**2)*np.cos(m*(theta-theta0))
def equation(r, theta, z):
    return E0*f*np.exp(-r**2/r0**2)*np.cos(k0*z)

## Correct
def equation(r, theta, z):
    f = np.sqrt(np.math.factorial(p)/(np.math.factorial(m+p)))*((np.sqrt(2)*r/r0)**m)*(2*r**2/r0**2)*np.cos(m*(theta-theta0))
    return E0*f*np.exp(-r**2/r0**2)*np.cos(k0*z)

修改后的输出。x1c 0d1x

第二步

将您的圆柱坐标值转换为直角坐标！这使得ax.scatter(...)可以适当地绘制您的数据。

## Incorrect
X, Y, Z = np.meshgrid(r, theta, z)
E = equation(X, Y, Z)

# Create a 3D density plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Plot the density using a colormap
density_plot = ax.scatter(Z, Y, X, c=E, cmap='plasma', s=1)

## Correct
X, Y, Z = np.meshgrid(r, theta, z)
E = equation(X, Y, Z)

xs = X*np.cos(Y)
ys = X*np.sin(Y)
zs = Z

# Create a 3D density plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Plot the density using a colormap
density_plot = ax.scatter(xs, ys, zs, c=E, cmap='plasma', s=1)

修改后的输出。

第三步

更改θ的边界以创建完整的圆柱体而不是半个（2π而不是π）。

## Incorrect
theta = np.linspace(0, np.pi, 100)

## Correct
theta = np.linspace(0, 2*np.pi, 100)

修改后的输出。

最终产品

即使解决了阴谋问题，您的输出将不会像上面的图像。这是因为您用于计算E的公式是不正确的。这里是我尝试修复错误的一些代码。我还更改了E的公差级别，以便只绘制极值，更容易看到螺旋线。请参阅Laser Resonators and Beam Propagation的第222页Hodgson和Weber对所用方程的完整描述。

代码

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import time

start = time.perf_counter()

# Defined Parameters
ps = 1  # The order of the Laguerre polynomial
ms = 1 #The azimuthal order of the pulse
r0s = 1
z = np.linspace(-5, 5, 20)
x = np.linspace(-1.5,1.5,15)
y = np.linspace(-1.5,1.5,15)

a,b = np.meshgrid(x,y)

rs = np.sqrt(a**2 + b**2)
thetas = np.arctan2(a,b)
fs = np.sqrt( (2 * np.math.factorial(ps)) / (np.pi * np.math.factorial(ps+ms)))

# Define the equation E = f(r, theta, z)
def equation(r, theta, z, m, r0, f):
    multiplier = (1 / r0) * np.exp(-(r**2) / (r0**2)) * (((r * np.sqrt(2)) / r0)**m)
    laguerre1 = (m+1) - ((2 * (r**2)) / (r0**2))
    return f*multiplier*laguerre1* (np.cos(z+theta))

# Generate data for the plot
X, Y, Z = np.meshgrid(rs, thetas, z)
E = equation(X, Y, Z, ms, r0s, fs)

xs = (X*np.cos(Y))[abs(E)>0.25]
ys = (X*np.sin(Y))[abs(E)>0.25]
zs = Z[abs(E)>0.25]

# Create a 3D density plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Plot the density using a colormap
density_plot = ax.scatter(xs, ys, zs, c=E[abs(E)>0.25], cmap='plasma', s=1)

# Add a colorbar to the plot
cbar = plt.colorbar(density_plot)
cbar.set_label('E Value')

# Set axis labels
ax.set_xlabel('X Label')
ax.set_ylabel('Y Label')
ax.set_zlabel('Z Label')

# Show the plot
plt.show()

end = time.perf_counter()

print("\nTotal runtime:", (end-start), "s")