为什么不绘制以下阵列因子的 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()
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()