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

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.

source https://stackoverflow.com/questions/71995629/why-is-not-the-following-3d-polar-plot-of-array-factor-being-plotted