为什么不绘制以下阵列因子的 3D 极坐标图?

Why is not the following 3D polar plot of Array Factor being plotted?

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

####################################################################################
fig = plt.figure()
ax = fig.gca(projection='3d')   #fig.add_subplot(1,1,1, projection='3d')
####################################################################################

freq = 30e+6        # in Hz     # 30 MHz
lam = 3e+8/freq     # in m
k_o = 2*np.pi/lam

theta = np.linspace(0, np.pi, 100)  # in radians
phi = np.linspace(0, 2*np.pi, 100)  # in radians


###############################################################################################################################
N_E = 6     # number of elements # positions on the vertices of a pentagon
###############################################################################################################################
D_XY = 3.5  # in m
R_XY = 1.7*D_XY
C1 = 0.309
C2 = 0.809
S1 = 0.951
S2 = 0.5878
XY = np.array([[0, 0], [0, R_XY], [-R_XY*S1, R_XY*C1], [-R_XY*S2, -R_XY*C2], [R_XY*S2, -R_XY*C2], [R_XY*S1, R_XY*C1]])
#print(len(XY))
w = np.array([3, 1, 1, 1, 1, 1])
###############################################################################################################################

def gain(XY, k_o, w, theta, phi):
    """Return the power as a function of azimuthal angle, phi."""

    AF = 0.0
    for n in range(len(XY)):

        relative_phase = k_o*( (XY[n][0])*np.sin(theta)*np.cos(phi) + (XY[n][1])*np.sin(theta)*np.sin(phi) )    #Relative phase calculation
        beta_n = -k_o*( (XY[n][0])*np.sin(mt.radians(30))*np.cos(mt.radians(60)) + (XY[n][1])*np.sin(mt.radians(30))*np.sin(mt.radians(60)) )   #beta

        psi = relative_phase + beta_n   # Progressive phase-shift

        AF = AF + w[n]*np.exp(1j*psi)

    g = np.abs(AF)**2
    return g

###############################################################################################################################

def get_directive_gain(g, minDdBi=-20):
    """Return the "directive gain" of the antenna array producing gain g."""
    DdBi = 10 * np.log10(g / np.max(g))
    return np.clip(DdBi, minDdBi, None)

###############################################################################################################################

th, ph = np.meshgrid(theta, phi)
G = np.zeros((100,100))
for l in range(100):
    for o in range(100):
        g = gain(XY, k_o, w, th[l], ph[o])
        G[l][o] = g # get_directive_gain(g)

plot = ax.plot_surface(th, ph, G, cmap='viridis', edgecolor='none')
plt.show()

不过,gain/array 因子相对于 theta 和 phi 的二维极坐标图,单独地,工作准确。然而,当我试图绘制 theta、phi、gain 时,它是行不通的。我发现“3d-Polar-Plot”有点棘手。所以,虚心请教各位高人指点或指出错误。我将有义务。

在您的代码中,th[l], ph[o] 是两个长度为 100 的数组。它们被传递给函数 gain,其中 returns 是一个长度为 100 的数组。G是一个二维数组,所以元素 G[l, o] 应该是一个数字。相反,使用 G[l][o] = g 您将其设置为一个数组,因此您会在终端上看到错误。

您需要相应地更改 for 循环。

此外,要获得极坐标图,您需要应用变换以计算 X、Y 坐标,从 theta, phi 开始。我用过这个spherical coordinate transformation:你可能需要根据你的参考系统调整它。

th, ph = np.meshgrid(theta, phi)
G = np.zeros_like(th)
for l in range(len(theta)):
    g = gain(XY, k_o, w, theta[l], phi)
    # set all the rows at column `l` with g
    G[:, l] = g #get_directive_gain(g)

from matplotlib.colors import Normalize
import matplotlib.cm as cm
norm = Normalize(vmin=G.min(), vmax=G.max())
colors = norm(G)
cmap = cm.get_cmap("viridis")

# spherical projection
X = G * np.sin(th) * np.cos(ph)
Y = G * np.sin(th) * np.sin(ph)
Z = G * np.cos(th) 
surf = ax.plot_surface(X, Y, Z, cmap=cmap, facecolors=cmap(colors))
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("G")

mappable = cm.ScalarMappable(norm=norm, cmap=cmap)
mappable.set_array(colors)

cb = fig.colorbar(mappable)
cb.set_label("Gain")
plt.show()

Matplotlib 在 3D 绘图方面真的很糟糕(从 shrink-ed 比例可以看出),因此您可能想切换到其他库。这里是 K3D-Jupyter:

### K3D-related stuff

def ij2k(cols, i, j):
    """Create the connectivity for the mesh.
    https://github.com/K3D-tools/K3D-jupyter/issues/273
    """
    return cols * i + j

def get_vertices_indices(x, y, z):
    """Compute the vertices matrix (Nx3) and the connectivity list for
    triangular faces.
    Parameters
    ==========
        x, y, z : np.array
            2D arrays
    """

    rows, cols = x.shape
    x = x.flatten()
    y = y.flatten()
    z = z.flatten()
    vertices = np.vstack([x, y, z]).T
    indices = []
    for i in range(1, rows):
        for j in range(1, cols):
            indices.append(
                [ij2k(cols, i, j), ij2k(cols, i - 1, j), ij2k(cols, i, j - 1)]
            )
            indices.append(
                [ij2k(cols, i - 1, j - 1), ij2k(cols, i, j - 1), ij2k(cols, i - 1, j)]
            )
    return vertices, indices

vertices, indices = get_vertices_indices(X, Y, Z)

import k3d
fig = k3d.plot(grid_visible=False)
mesh = k3d.mesh(
    vertices, indices,
    side="double",
    flat_shading=False,
    attribute=Z.astype(np.float32),
    color_map=k3d.colormaps.matplotlib_color_maps.Viridis
)
fig += mesh
fig.display()