# 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}^{2\pi} \int_{0}^{\pi} |F(\theta,\phi)|^2 sin(\theta) \, d\theta d\phi} }$$

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{2\pi }{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, 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])])])

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

pow_sum = 0
for i in range(vals_F.shape):
for j in range(vals_F.shape):
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*vals_F.shape)
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):
for j in range(vals_F.shape):
pow_sum = pow_sum + abs(vals_F[i, j]) * abs(vals_F[i, j])

integral = pow_sum/(vals_F.shape*vals_F.shape)  # 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.

• An approximate method using half power beam width. "The beam solid angle can be approximated for antennas with one narrow major lobe and very negligible minor lobes by simply multiplying the half-power beam widths (in radians) in two perpendicular planes." Half power beam width can be obtained approximately some measured data.
– AJN
Jul 25, 2020 at 3:00
• Thanks, but I'm planning to use it for omni-directional antennas so there's no main lobe.
– AMTK
Jul 25, 2020 at 4:31
• Can you give the exact steps /calculations which led you to "Naively dumping this data into the above formula unsurprisingly gives the wrong answer.". May be there was some error in the numerical integration.
– AJN
Jul 25, 2020 at 4:41
• Ok, I added the actual python code.
– AMTK
Jul 25, 2020 at 5:28
• Three great circles is too few points and the intersection points will be counted multiple times.
– AJN
Jul 25, 2020 at 5:31