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

Question

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()

Although, 2D-polar plot of the gain/array factor with respect to theta and phi, individually, is working accurately. However, as I am trying to plot theta, phi, gain altogether it is not working. I found it to be bit tricky "3d-Polar-Plot". So, I humbly request if anyone could kindly suggest me or point out my mistakes. I would be obliged.

Answer 1

In your code, th[l], ph[o] are two arrays of length 100. They are passed to the function gain, which returns an array of length 100. G is a 2D array, so the element G[l, o] is expected to be a single number. Instead, with G[l][o] = g you were setting it to be an array, hence the error you would see on your terminal.

You need to change the for loop accordingly.

Also, to obtain a polar plot you need to apply the transformation in order to compute the X, Y coordinates, starting from theta, phi. I have used this spherical coordinate transformation: you might need to adapt it to your reference system.

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 is really terrible when it comes to 3D plots (as you can see from the shrink-ed scale), so you might want to switch to a different library. Here is 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()