Calculating Antenna Directivity Numerically

The equation for peak directivity of an antenna is this (example source)

$${D = frac{1 }{ displaystyle frac{1 }{4 pi} int_{0}^{2pi} int_{0}^{pi} |F(theta,phi)|^2 sin(theta) , dtheta dphi} }$$

Where F is the amplitude at the specified spherical coordinates ${(theta,phi) }$, normalized to a peak value of 1.

This can be approximated numerically as follows:

$${D = frac{1 }{ displaystyle frac{1 }{4 pi} frac{pi }{N} frac{2pi }{M} sum_{j=1}^{M} Big[ sum_{i=1}^{N} |F(theta _{i},phi _{j})|^2 sin(theta_{i}) Big] } }$$

I’ve satisfied myself that this is valid by setting ${ F = sqrt{sin(theta)} }$, and checking the integral against the numerical calculation with $theta$ and $phi$ in 1 degree increments.

My question is, what can I do with data that is not organized by $theta$ and $phi$? For instance say I obtain measurements around the equator ($theta = 90^o$) and vertical meridians at $phi = 0^o / 180^o$ and $phi = 90^o / 270^o$ (I’m visualizing the device-under-test as a little Earth if it’s not obvious!) I tried manually creating the data arrays $F,theta $ and $phi$ each with 3 rows of 360 points, each row representing a trip around the specified circumference. Naively dumping this data into the above formula unsurprisingly gives the wrong answer.

In my example I can actually just throw out the row of equator data, since I’ve defined ${ F = sqrt{sin(theta)} }$ and the equator has constant ${theta = 90^o}$. With the remaining two rows I actually do get the right answer. But I need a method that will work for an arbitrary pattern (ie. a real antenna measurement).

---

By request here is the numerical integration python code, it’s a bit ugly. You can see the commented out sections that do work, but the code as written gives the wrong answer:

import numpy as np

# # Standard phi and theta distribution **working**
# THETA = np.deg2rad(np.arange(0,181,1))
# PHI = np.deg2rad(np.arange(-180,181,1))
# THETA, PHI = np.meshgrid(THETA, PHI)

# 2 meridians (front and side) **working**
# vals_theta = np.hstack([np.vstack([np.arange(0,181,1).reshape(-1,1),np.arange(179,-0.1,-1).reshape(-1,1)]),np.vstack([np.arange(0,181,1).reshape(-1,1),np.arange(179,-0.1,-1).reshape(-1,1)]),np.vstack([np.arange(0,181,1).reshape(-1,1),np.arange(179,-0.1,-1).reshape(-1,1)]),np.vstack([np.arange(0,181,1).reshape(-1,1),np.arange(179,-0.1,-1).reshape(-1,1)])])
# vals_phi = np.hstack([np.vstack([0*np.ones([181,1]),180*np.ones([180,1])]),np.vstack([90*np.ones([181,1]),270*np.ones([180,1])]),np.vstack([180*np.ones([181,1]),0*np.ones([180,1])]),np.vstack([270*np.ones([181,1]),90*np.ones([180,1])])])

# 3 main circumferences (equator, front meridian, side meridian) - need to remove overlaps to increase accuracy
vals_theta = np.hstack([90*np.ones([361,1]),np.vstack([np.arange(0,181,1).reshape(-1,1),np.arange(179,-0.1,-1).reshape(-1,1)]),np.vstack([np.arange(0,181,1).reshape(-1,1),np.arange(179,-0.1,-1).reshape(-1,1)])])
vals_phi = np.hstack([np.arange(0,360.1,1).reshape(-1,1),np.vstack([0*np.ones([181,1]),180*np.ones([180,1])]),np.vstack([90*np.ones([181,1]),270*np.ones([180,1])])])
THETA = np.deg2rad(vals_theta)
PHI = np.deg2rad(vals_phi)

vals_F = np.sqrt(np.sin(THETA)) # should have a directivity of 1.27

pow_sum = 0
for i in range(vals_F.shape[0]):
    for j in range(vals_F.shape[1]):
        pow_sum = pow_sum + abs(vals_F[i,j])*abs(vals_F[i,j])*np.sin(THETA[i,j])

integral = pow_sum*np.pi/(vals_F.shape[0]*vals_F.shape[1])
directivity = 1/(1/2*integral)

print(directivity)

Despite the messy expressions for vals_theta and vals_phi they are quite simple. But I didn’t expect this to work because the values being summed no longer align with the limits of the integral/summation.

vals_theta
array([[90.,  0.,  0.],
       [90.,  1.,  1.],
       [90.,  2.,  2.],
       ...,
       [90.,  2.,  2.],
       [90.,  1.,  1.],
       [90.,  0.,  0.]])

vals_phi
array([[  0.,   0.,  90.],
       [  1.,   0.,  90.],
       [  2.,   0.,  90.],
       ...,
       [358., 180., 270.],
       [359., 180., 270.],
       [360., 180., 270.]])

---

One other idea I had is to simply average the power at all points, as per the definition of directivity. I guess this assumes that the points being sampled are evenly distributed around the sphere in some way which seems possibly true for the three great circles.

Here is the modified calculation in python:

for i in range(vals_F.shape[0]):
    for j in range(vals_F.shape[1]):
        pow_sum = pow_sum + abs(vals_F[i, j]) * abs(vals_F[i, j])

integral = pow_sum/(vals_F.shape[0]*vals_F.shape[1])  # average power
directivity = 1/integral

This actually gives a pretty close answer when ${ F = sqrt{sin(theta)} }$ (it gives 1.32 compared to the true value of 1.27). But for ${ F = sin^5(theta) }$ it’s a ways off, 2.01 compared to the true value of 2.71, so maybe I just got lucky the first time.