#!/usr/bin/env python

'''This program parses antenna pattern files as described in NSMA (nsma.org)
Recommendation WG 16.99.050, "Antenna Systems - Standard Format for
Digitized Antenna Patterns."

Recognizing that mistakes are made and people don't always strictly follow a
standard, the program attempts to be somewhat flexible.  The flexibility is
limited by my test cases, so please continue improve the code based on your
experiences.  Typical usage is to take output from the -s (summary)
argument to then make plots.  Alternatively, the -t (text) argument displays
the file with indentation to see what fields are used.  Example:

  $ nsma.py -s ant.adf 
  Line 402 warning: expected 180 points, got 179.
  ABC Antenna Company 800A-065-25-4N
   E-tilt: 4.0 deg.
   851 MHz         (-f 851):
    EL V/V         (-p V/V)
    AZ V/V         (-p V/V)

2D polar plots can be made using:

  $ nsma.py -f 851 -p v/v -2 ant.adf

Or interactive 3D with:

  $ nsma.py -f 851 -p v/v -3 ant.adf

If an input file specifies no frequencies, the -f option is not needed.
If there is only one frequency or polarization type in a file, each will be
assumed without need for command line args specifying them.  The following
will work because ant.adf has patterns for one frequency, and while there
are multiple cuts there is only one polarization type.

  $ nsma.py -3 ant.adf

It is possible to find the antenna gain given azimuth and elevation angles
in degrees.  For example, gain for an RF ray at azimuth of 0 and elevation
of 34 degrees:

  $ nsma.py -a 0 -e 36 ant.adf
  gain(az=0.0, el=36.0) = -31.2 dBi

The dBi units are drawn from the NSMA file.  When the RF ray is not directly
in a cut, gains from the nearest two planes are calculated and summed after
weighting based on ray closeness to each plane.  It is not exact, especially
in the typical case of only two cuts, but makes the most of sparse data.

Mike Markowski, mike.ab3ap@gmail.com
Mar 2023
'''

from matplotlib.collections import PatchCollection
from matplotlib.patches import Polygon
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import matplotlib.pyplot as plt
import numpy as np
import sys

apd = {} # Antenna pattern dictionary.
curSlice = None

def cli(argv):
	'''Retrieve command line arguments.  See usage() for available arguments.

	Input:
	  argv (string[]): command line arguments.
	Outputs:
	  freq (string): frequency in MHz as provided by user.
	  plot (string): '2', '3', '' for 2D, 3D, or no plot.
	  pol (string): polarizations of receiving and transmitting antennas.
	  text (string): 's', 't', '' for short, full or no text summary.
	'''

	az_deg = el_deg = None
	freq = plot = pol = text = ''
	fc = 'moccasin'

	prog = argv[0]
	argv = argv[1:]
	i = 0
	while i < len(argv):
		arg = argv[i]
		if arg == '-2': # 2D polar plots.
			plot = '2'
		elif arg == '-3': # 3D plot of cuts.
			plot = '3'
		elif arg == '-a': # Azimuth of RF ray.
			i += 1
			az_deg = float(argv[i])
		elif arg == '-c': # Color of 3d cuts.
			i += 1
			fc = argv[i]
		elif arg == '-e': # Elevation of RF ray.
			i += 1
			el_deg = float(argv[i])
		elif arg == '-f': # Frequency of cut to plot.
			i += 1
			freq = argv[i]
		elif arg == '-h': # Help message.
			usage()
		elif arg == '-p': # Polarization pair, rx/tx.
			i += 1
			pol = argv[i]
		elif arg == '-s': # Summary text of NSMA file.
			text = 's'
		elif arg == '-t': # Text detail of NSMA file.
			text = 't'
		elif arg[0] == '-':
			usage()
		i += 1
	if i > len(argv):
		usage()
	antFile = argv[i-1] # NSMA antenna data file to work with.
	return freq, plot, pol, text, antFile, az_deg, el_deg, fc

def findCuts(apd, freq, pol, RF, top=1):
	'''Find 'top' number of cuts in an NSMA antenna file at specified
	frequency and polarization pair that are closest to RF ray, a 3d vector.
	Since the spec does not impose an ordering scheme, none is assumed here.

	NOTE: untested for phi cuts.  Please send me a phi cut file
	if you have one and I will be glad to debug this code as needed.

	Input:
	  apd (dictionary): NSMA antenna dictionary describing antenna.
	  freq (float): frequency of cuts in NSMA file to plot.
	  RF (float[3]): unit vector of RF ray piercing pattern.
	  top (int): 1 or 2, number of nearest cuts to return.
	Output:
	  (string[]): list of closest cuts.
	  (float[:top], dictionary[:top], float[:top][3]): RF ray heights, cuts
	    and yz-plane unit vector edges of cuts nearest the RF vector.
	'''

	try:
		cuts = apd['patfre'][freq]['patcut']
	except KeyError:
		print('Error: no such frequency %s MHz.' % freq)
		sys.exit(1)

	h = np.array([]) # Heights above cut planes.
	c = np.array([]) # Cuts.
	e = [] # yz-plane unit vector edges of cuts.

	for cutName in cuts.keys():
		curC = cuts[cutName]
		if curC['polari'].lower() != pol.lower():
			continue # Only consider requested polarization pairs.

		# Convert cut names to pairs of angles.  A pair is used because the
		# RF ray can be closer to plane on either side of the x-axis.
		cName = cutName.lower()
		if cName in ['az', 'h']:
			angCut_rad = 0
		elif cName in ['el', 'v']:
			angCut_rad = np.pi/2
		else: # Phi cut.
			angCut_rad = np.radians(float(cName) % 180)

		# Create unit vector of cut angle on yz-plane.
		edge = np.array([0, np.cos(angCut_rad), np.sin(angCut_rad)])
		# Project RF onto pattern cut plane.
		rfProj = edge*np.dot(RF, edge) + np.array([RF[0],0,0])
		rfHgt = np.sqrt(1 - (rfProj**2).sum())
		# Save results.
		h = np.append(h, rfHgt)
		c = np.append(c, curC)
		e.append(edge) # Unit vector along yz-plane edge of cut.

	# Sort results.
	i = np.argsort(h) # Index order to sort h[] ascending.
	h = h[i]
	c = c[i]
	e = np.array(e)[i]
	return h[:top], c[:top], e[:top]

def findFreqs(apd):
	'''Find all frequencies mentioned in an NSMA antenna file.

	Input:
	  apd (dictionary): NSMA antenna dictionary describing antenna.
	Output:
	  (string[]): list of unique frequencies found in NSMA file.
	'''

	freqs = []
	for freq in apd['patfre'].keys(): # Each key is a frequency.
		if not freq in freqs:
			freqs.append(freq)
	return freqs

def findPols(apd):
	'''Find all polarization pairs mentioned in an NSMA antenna file.

	Input:
	  apd (dictionary): NSMA antenna dictionary describing antenna.
	Output:
	  (string[]): list of unique polarizations found in NSMA file.
	'''

	pols = []
	for freq in apd['patfre'].keys(): # Each key is a frequency.
		for cut in apd['patfre'][freq]['patcut'].keys(): # Each key is a cut.
			pol = apd['patfre'][freq]['patcut'][cut]['polari']
			if not pol in pols:
				pols.append(pol)
	return pols

def gain(apd, freq, pol, az_deg, el_deg):
	'''Given an antenna pattern dictionary, return antenna gain for
	specified frequency, polarization pair, and az/el.  It is an error
	if frequency or polarization pair do not exist in the NSMA file.
	If neccessary, gain for requested azimuth and elevation will be
	interpolated.  NSMA phi cuts provide a wealth of pattern data.
	It is harder to interpolate gain in the typical case when only
	vertical and horizontal cuts are given.  This subroutine does its
	best to do so with the minimal information.

	Input:
	  apd (dictionary): NSMA antenna dictionary describing antenna.
	  freq (float): frequency of cuts in NSMA file to plot.
	  pol (string): polarization pair, rx/tx, for file to plot.
	  az_deg (float): azimuth in degrees of RF ray piercing pattern.
	  el_deg (float): elevation in degrees of RF ray piercing pattern.
	Output:
	  (float): gain at az/el.
	'''

	xUnit = np.array([1,0,0])
	yUnit = np.array([0,1,0])
	zUnit = np.array([0,0,1])

	# Create RF unit vector, R, from az, el angles.
	az_rad = np.radians(az_deg)
	el_rad = np.radians(el_deg)
	Rxy = np.cos(el_rad)       # R projected onto xy-plane.
	Rx = Rxy*np.cos(az_rad)    # Rxy projected on x-axis.
	Ry = Rxy*np.sin(az_rad)    # Rxy projected on y-axis.
	Rz = np.sin(el_rad)        # R projected on z-axis.
	R = np.array([Rx, Ry, Rz]) # The RF unit vector!

	# Get nearest heights, cuts, and yz-plane unit vector.
	h, c, e = findCuts(apd, freq, pol, R, top=2)

	gain = np.zeros(h.size)
	for i in range(h.size):
		v = xUnit*np.dot(R, xUnit) + e[i]*np.dot(R, e[i]) # Map to cut.
		v /= np.linalg.norm(v) # Normalize.
		ang_deg = np.degrees(np.arccos(np.dot(xUnit, v)))
		sign = np.dot(e[i], v)
		if sign < 0:
			ang_deg = -ang_deg
		gain[i] = gainCut(c[i], ang_deg) # Gain on this cut.
	if h.size == 1: # NSMA file contained single cut.
		gainTot = gain[0]
	elif h[0] == h[1] == 0: # X-axis RF ray.
		gainTot = gain.mean()
	elif h[0] == 0: # Ray falls on one plane.
		gainTot = gain[0]
	elif h[1] == 0: # Ray falls on one plane.
		gainTot = gain[1]
	else: # Interpolate between two gain values weighted by planar distance.
		gainTot = h[1]/(h[0]+h[1])*gain[0] + h[0]/(h[0]+h[1])*gain[1]

	units    = 'dB' + apd['gunits'][0][2].lower() # dBi, dBd, or dBr.
	unitsPat = 'dB' + apd['gunits'][1][2].lower() # dBi, dBd, or dBr.
	if unitsPat == 'dBr':
		return gainTot + apd['mdgain'][0], units
	else:
		return gainTot, units

def gainCut(cut, ang_deg):
	'''Given a cut and an angle, return the interpolated gain.

	Inputs:
	  cut (dictionary): cut to be used to find gain.
	  ang_deg (float): angle in degrees to find in cut.
	Output:
	  (float): gain in cut at specified angle.
	'''

	theta = cut['angle']
	gain = cut['gain']
	i0 = np.where(theta<=ang_deg)[0][-1]
	i1 = np.where(theta>=ang_deg)[0][0]
	if i0 == i1: # Exact angle is in pattern.
		return gain[i0]
	d0 = abs(theta[i0] - ang_deg) # Angular difference.
	d1 = abs(theta[i1] - ang_deg) # Angular difference.
	gain = d1/(d0+d1)*gain[i0] + d0/(d0+d1)*gain[i1]
	return gain # Gain interpolated between nearest two angles.

def getFloat(s, err=None, die=True):
	'''Helper function to convert a string to a float.  If conversion fails,
	an error message can optonally be printed, and the program can optionally
	exit.

	Inputs:
	  s (string): string to convert to float.
	  err (string): if not None, err is printed upon conversion failure.
	  die (boolean): if True, exit program with status 1 upon conversion
	    failure.
	Output:
	  (float): the float converted from the input string.
	'''

	if s.strip() == '':
		return 0.
	try:
		n = float(s)
	except ValueError:
		if not err is None:
			print('%s' % err)
		if die:
			sys.exit(1)
		n = None
	return n

def getInt(s, err):
	'''Helper function to convert a string to a int.  If conversion fails,
	an error message can optonally be printed, and the program can optionally
	exit.

	Inputs:
	  s (string): string to convert to int.
	  err (string): if not None, err is printed upon conversion failure.
	  die (boolean): if True, exit program with status 1 upon conversion
	    failure.
	Output:
	  (int): the int converted from the input string.
	'''
	if s.strip() == '':
		return 0
	try:
		n = int(float(s))
	except ValueError:
		print('%s' % err)
		sys.exit(1)
	return n

def makeTitle(apd, freq, pol, cutName='', plot3d=True):
	'''Create a reasonable plot title from NSMA file data.

	Inputs:
	  apd (dictionary): NSMA antenna dictionary describing antenna.
	  freq (string): frequency in MHz of cuts to plot.
	  pol (string): polarization of cuts to plot.  E.g., H/V.
	  cutName (string): name of cut to plot.
	  plot3d (boolean): True/False for 3d/2d plot.
	Output:
	  (string): title that can be used for 2d or 3d plot.
	'''

	u0 = list(apd['gunits'][0].lower())
	u0[1] = u0[1].upper()
	u0 = "".join(u0) # dBi, dBd, or dBr.
	u1 = list(apd['gunits'][1].lower())
	u1[1] = u1[1].upper()
	u1 = "".join(u1) # dBi, dBd, or dBr.
	units = [u0, u1]
	s = ''
	if apd['modnum'] != '':
		s += '%s: ' % apd['modnum']       # Antenna model number.
	if freq != '':
		s += '%s Mhz' % freq              # Frequency of cuts plotted.
	if freq != '' and apd['mdgain'][0] != 0:
		s += ', '
	if apd['mdgain'][0] != 0:
		s += 'gain %s %s' % (apd['mdgain'][0], units[0]) # Gain at mid-band.
	if not plot3d and s != '':
		s += '\n%s' % cutName             # Name of cut.  'Azimuth,' etc.
	if pol != '':
		s += ', %s (rx/tx)' % pol.upper() # Polarization pair.  'V/V,' etc.
	return s, units

def miniNsma(apd):
	'''Print a summary of an antenna pattern.

	Input:
	  apd (dictionary): antenna pattern dictionary to be summarized.
	Output:
	  none, summary printed to stdout.
	'''

	for key in apd.keys():
		if key == 'antman':
			print('%s ' % apd[key], end='')
		elif key == 'eltilt' and apd[key][0] != 0:
			print(' E-tilt: %.1f deg.' % apd[key][0])
		elif key == 'modnum':
			print('%s' % apd[key])
		elif key == 'patfre':
			for f in apd[key].keys():
				print(' %s MHz         (-f %s):' % (f,f))
				for cut in apd['patfre'][f]['patcut'].keys():
					pol = apd['patfre'][f]['patcut'][cut]['polari']
					print('  %s %s         (-p %s)' % (cut, pol, pol))

def nsma(filename):
	'''Read an NSMA antenna file and return a dictionary of dictionaries
	contaoning the data.  The dictionary tree can be traversed for structured
	text output or to generate plots.

	Input:
	  filename (string): name of NSMA antenna file to be parsed.
	Output:
	  (dictionary): dictionary of NSMA data formatted as displayed using the 
	    command line -t option.
	'''

	apd = {} # Antenna pattern dictionary.
	cut = 0
	freq = 0
	pt = 0
	readingPts = False
	curC = curF = None

	n = 0
	f = open(filename)
	nsmaData = f.read().split('\n')
	f.close()

	for line in nsmaData:
		n += 1 # Line count for error messages.

		# Things to ignore.
		excl = line.find('!')
		if excl != -1: # Ignore comments.
			line = line[:excl].strip()
		if line == '': # Ignore blank lines.
			continue

		# Parse values.
		if line[6:8] == ':,': # Keyword.
			if readingPts: # Finished reading pattern data.
				if nPts != curC['angle'].size:
					print('Line %d warning: expected %d points, got %d.'
						% (n, nPts, curC['angle'].size))
				readingPts = False
			field = line[:6].lower()
			val = line[8:].strip()
			v = parse(field, val, n)
			if field == 'endfil':
				break
			elif field == 'fstlst':
				curC['fstlst'] = v
			elif field == 'nofreq':
				apd['patfre'] = {}
			elif field == 'numcut':
				if curF == None:
					curF = apd['patfre']['any'] = {}
				curF['patcut'] = {}
			elif field == 'nupoin':
				nPts = v
			elif field == 'patfre':
				curF = apd['patfre'][v] = {}
			elif field == 'patcut':
				if not v.lower() in ['az', 'el', 'h', 'v']:
					phi = getFloat(v, die=False)
					if phi is None:
						print('Line %d error: cut %s not in AZ, EL, H, V.'
							% (n, v))
						sys.exit(1)
				curC = curF[field][v] = {}
				curC['polari'] = ''
			elif field == 'polari':
				curC[field] = v
			else:
				apd[field] = v
		else: # Pattern data.
			if curC is None:
				print('Line %d parse error: %s' % (n, line))
				sys.exit(1)
			if not readingPts: # Create arrays for new cut.
				curC['angle'] = np.array([])
				curC['gain'] = np.array([])
				curC['phase'] = np.array([])
				readingPts = True
			v = parse(line, None, n)
			a = curC['angle']
			g = curC['gain']
			p = curC['phase']
			curC['angle'] = np.append(a, v[0])
			curC['gain'] = np.append(g, v[1])
			curC['phase'] = np.append(p, v[2])

	f.close()
	return apd

def parse(field, val, n):
	'''Most lines in an NSMA file are formatted as VAR:,VAL where ':,' is a
	separator.  This subroutine parses the VAL appropriately depending on
	which VAR keyword is read.

	Inputs:
	  field (string): variable whose value is to be parsed.
	  field (string): value to be parsed.
	  n (int): line number of NSMA file, used in error messages.
	'''

	global apd # Antenna pattern dictionary.

	if val == None: # Pattern data.
		# Parse angle, gain, phase where phase is optional.
		c1 = field.find(',')
		if c1 == -1:
			print('Line %d: pattern data does not contain a comma.' % n)
			sys.exit(1)
		err1 = "Line %d: can't convert angle %s to float." % (n, field[:c1])
		n1 = getFloat(field[:c1], err1)

		c2 = field.find(',', c1+1)
		c2p = len(field) if c2==-1 else c2
		err2 = "Line %d: can't convert gain %s to float." % (n, field[c1+1:c2p])
		n2 = getFloat(field[c1+1:c2p], err2)

		n3 = 0
		if c2 > -1:
			err3 = "Line %d: can't convert phase %s to float." \
				% (n, field[c2+1:])
			n3 = getFloat(field[c2+2:], err3)
		val = np.array([n1, n2, n3])
	elif field == 'revnum':
		if val != 'NSMA WG16.99.050':
			print('Warning: expected revnum "NSMA WG16.99.050".')
			print('         Line %d revnum is: "%s"' % (n, val))
	elif field == 'gunits': # dBi, dBd, dBr, or Lin.
		slash = val.find('/')
		gainBand = val[:slash] # Gain units for low-, mid-, high-band values.
		gainPat = val[slash+1:] # Gain units for pattern data.
		val = [gainBand, gainPat]
	elif field in ['hghfrq', 'lowfrq', 'lwgain', 'hggain', 'atvswr',
		'frtoba', 'radctr', 'potopo', 'maxpow', 'antlen', 'antwid', 'antdep',
		'antwgt']:
#		'antwgt', 'patfre']:
		err = "Line %d: can't convert %s to float." % (n, val)
		val = getFloat(val, err)
	elif field in ['nofreq', 'numcut', 'nupoin']:
		err = "Line %d: can't convert %s to int." % (n, val)
		val = getInt(val, err)
	elif field in ['azwidt', 'eltilt', 'elwidt', 'fstlst', 'mdgain']:
		# Parse num1,num2 where num2 is optional.
		comma = val.find(',')
		cp = len(field) if comma==-1 else comma
		err = "Line %d: can't convert %s to float." % (n, val[:cp])
		n1 = getFloat(val[:cp], err)
		n2 = 0
		if comma > -1: # Optional 2nd number exists.
			err = "Line %d: can't convert %s to float." % (n, val[comma+1:])
			n2 = getFloat(val[comma+1:], err)
		val = np.array([n1, n2])
	return val

def plotNsma2d(apd, freq, pol):
	'''Plot 2D antenna patterns for cuts for a given frequency and
	polarization pair.

	Inputs:
	  apd (dictionary): dictionary of NSMA antenna file values.
	  freq (string): frequency in MHz whose cuts are to be plotted.
	  pol (string): polarization pairs to be plotted.
	Output:
	  none, plots will appear on screen.
	'''

	try:
		freqs = apd['patfre']
	except KeyError:
		print('Error: no such frequency %s MHz.' % freq)
		sys.exit(1)
	try:
		cuts = apd['patfre'][freq]['patcut']
	except:
		print('Error: no pattern for %s MHz:\n' % freq)
		miniNsma(apd)
		sys.exit(1)

	for cutName in cuts.keys():
		curC = cuts[cutName]
		cName = cutName.lower()
		if cName == 'az':
			cutFullname = 'Azimuth'
		elif cName == 'el':
			cutFullname = 'Elevation'
		elif cName == 'h':
			cutFullname = 'Horizontal'
		elif cName == 'v':
			cutFullname = 'Vertical'
		else:
			cutFullname = 'Phi Cut %s deg' % cName

		if curC['polari'].lower() == pol.lower():
			theta = np.radians(curC['angle'])
			g = curC['gain']

			fig, ax = plt.subplots(figsize=(6,6),
				subplot_kw={'projection': 'polar'})
			ax.grid(True,which="minor",alpha=0.2)
			ax.grid(True,which="major",alpha=0.5)
			ax.minorticks_on()
			# ax.set_rmax(0)
			ax.set_thetagrids(range(0,360,30))
			ax.set_rlabel_position(210)
			# ax.set_rlabel_position(120)

			s, units = makeTitle(apd, freq, pol, cutFullname, plot3d=False)
			plt.title(s)
			plt.xlabel('Gain (%s)' % units[1])
			plt.plot(theta,g, color='blue', lw=1)
			plt.plot([theta[0],theta[-1]],[g[0],g[-1]], color='blue', lw=1)
			plt.show()

def plotNsma3d(apd, freq, pol, fc='moccasin'):
	'''Plot 2D antenna patterns for cuts for a given frequency and
	polarization pair.

	Inputs:
	  apd (dictionary): dictionary of NSMA antenna file values.
	  freq (string): frequency in MHz whose cuts are to be plotted.
	  pol (string): polarization pairs to be plotted.
	  fc (string): color of antenna pattern slices.
	Output:
	  none, plot will appear on screen.
	'''

	try:
		freqs = apd['patfre']
	except KeyError:
		print('Error: no such frequency %s MHz.' % freq)
		sys.exit(1)
	try:
		cuts = apd['patfre'][freq]['patcut']
	except:
		print('Error: no pattern for %s MHz:\n' % freq)
		miniNsma(apd)
		sys.exit(1)

	patterns = []
	for cutName in cuts.keys():
		curC = cuts[cutName]
		if curC['polari'].lower() == pol.lower():
			patterns.append(
				[cutName, curC['angle'], curC['gain'], curC['phase']])
	if len(patterns) == 0:
		print('Error: no patterns for %s MHz, polarity %s:\n' % (freq, pol))
		miniNsma(apd)
		sys.exit(1)

	# Plot all cuts in patterns[].
	x0 = y0 = z0 = np.Inf
	x1 = y1 = z1 = -np.Inf
	gainMin_dB = np.Inf
	for p in patterns:
		# Each p: [cutName angle gain phase].
		gainMin_dB = min(gainMin_dB, min(p[2]))
	cuts = []
	for p in patterns:
		# Each p: [cutName angle gain phase].
		theta = np.radians(p[1])
		g = p[2]
		g -= gainMin_dB      # Shift to remove negative values.
		x = g*np.cos(theta)  # Polar to Cartesian.
		y = g*np.sin(theta)  # Polar to Cartesian.
		z = np.zeros(x.size) # 2d to 3d on xy plane.
		pts = np.column_stack((x, y, z))
		cName = p[0].lower()
		if cName == 'az':
			ang_deg = 0
		elif cName == 'el':
			ang_deg = -90
		elif cName == 'h':
			ang_deg = 0
		elif cName == 'v':
			ang_deg = -90
		else: # Phi cut angle.
			ang_deg = p[0]
		Rx = rotX(ang_deg)
		rotPlane = np.matmul(pts, Rx) # Rotate plane about boresight.
		x = rotPlane.T[0]
		y = rotPlane.T[1]
		z = rotPlane.T[2]
		cuts.append(rotPlane)

		# Need these to scale plot.
		x0 = min(x0, min(x))
		x1 = max(x1, max(x))
		y0 = min(y0, min(y))
		y1 = max(y1, max(y))
		z0 = min(z0, min(z))
		z1 = max(z1, max(z))

	fig = plt.figure(figsize=(8,8))
	ax = fig.add_subplot(111, projection='3d')
	#ax.set_axis_off()

	for c in cuts:
		pattern = Poly3DCollection([c], alpha=0.4,
			facecolor=fc, edgecolor='black')
		ax.add_collection3d(pattern, zdir='z')

	s, units = makeTitle(apd, freq, pol)
	ax.set_title(s)
	ax.plot([x0,x1],[0,0],[0,0], lw=1, color='black') # Show boresight.
	ax.set_xlabel('Gain (%s)' % units[1])
	ax.set_ylabel(units[1])
	ax.set_zlabel(units[1])

	m0 = min(x0, y0, z0)
	m1 = max(x1, y1, z1)
	ax.set_xlim(m0,m1)
	ax.set_ylim(m0,m1)
	ax.set_zlim(m0,m1)

	xl0 = '%.1f' % gainMin_dB
	xl1 = '%.1f' % (round(x1 + gainMin_dB, 1) + 0)
	ax.set_xticks([0, x1], [xl0, xl1])
	ax.set_yticks([])
	ax.set_zticks([])

	if True:
		plt.show()
	else: # Make animated GIF frames moving smoothly along ellipse.
		a = 15 # Semi-major axis.
		b = 7.5 # Semi-minor axis.
		c = np.array([-45, 12.5]) # Center of ellipse.
		f = 0 # Frame number,
		for theta in np.linspace(0, 2*np.pi, 60, endpoint=False):
			ax.azim = c[0] + a*np.sin(theta)
			ax.elev = c[1] + b*np.cos(theta)
			plt.savefig('frame%03d.png' % f)
			f += 1

def printNsma(apd):
	'''Recursively print an antenna pattern data structure.

	Input:
	  apd (dictionary): antenna pattern dictionary to be summarized.
	Output:
	  none, summary printed to stdout.
	'''

	for key in apd.keys():
		if key == 'patfre':
			printNsmaFreq(apd[key])
		else:
			print('%s: %s' % (key, apd[key]))

def printNsmaCut(cut):
	'''Function supporting recursive printing of antenna pattern dictionary.

	Input:
	  cut (dictionary): cut dictionary to be printed.
	Output:
	  none, summary printed to stdout.
	'''

	print('    patcut:')
	for cutName in cut.keys():
		print('      %s: ' % cutName)
		for key in cut[cutName]:
			if key == 'angle':
				printNsmaPts(cut[cutName])
			elif key in ['gain', 'phase']:
				pass
			else:
				print('        %s: %s' % (key, cut[cutName][key]))

def printNsmaFreq(frq):
	'''Function supporting recursive printing of antenna pattern dictionary.

	Input:
	  frq (dictionary): frequency dictionary to be printed.
	Output:
	  none, summary printed to stdout.
	'''

	print('patfre:')
	for key in frq.keys():
		print('  %s MHz' % key)
		printNsmaCut(frq[key]['patcut'])

def printNsmaPts(cut):
	'''Function supporting recursive printing of antenna pattern dictionary.

	Input:
	  cut (dictionary): cut dictionary to be printed.
	Output:
	  none, summary printed to stdout.
	'''
	print('        pattern:')
	showedDots = False
	for i in range(cut['angle'].size):
		if 2 <= i < cut['angle'].size-2:
			if not showedDots:
				print('          ...')
				showedDots = True
		else:
			print('          %.3f, %.3f, %.3f'
				% (cut['angle'][i], cut['gain'][i], cut['phase'][i]))

def rotX(ang_deg):
	'''Generate a 3x3 matrix to rotate 3d points about the X axis.
	See https://en.wikipedia.org/wiki/Rotation_matrix#In_three_dimensions

	Input:
	  ang_deg (float): angle to revolve about x-axis.
	Output:
	  (float[][]): 3x3 rotation array.
	'''

	ang_rad = np.radians(ang_deg)
	Rx = np.zeros((3,3))
	Rx[0][0] = 1
	Rx[1][1] = Rx[2][2] = np.cos(ang_rad)
	Rx[1][2] = Rx[2][1] = np.sin(ang_rad)
	Rx[1][2] = -Rx[1][2]
	return Rx

def usage():
	print('nsma [-2] [-3] [-a az] [-e el] [-c color] [-f freq] [-p pol] [-s]'
		'[-t] file')
	print('     -2: plot all 2D cuts in directory \'cuts\'.')
	print('     -3: plot interactive 3D pattern.')
	print('     -a: azimuth of RF ray piercing pattern.')
	print('     -e: elevation of RF ray piercing pattern.')
	print('     -c: color of 3D antenna cuts.')
	print('     -f: frequency in MHz of 3D cuts to plot.')
	print('     -p: polarization.')
	print('     -s: summary text for -f and -p values.')
	print('     -t: full text antenna description.')
	print('     file: NSMA antenna data file.')
	sys.exit(1)

# m a i n
#
# A simple program to read and print or plot an NSMA antenna file.

def main(argv):

	freq, plot, pol, text, a, az_deg, el_deg, fc = cli(argv)
	apd = nsma(a)
	allFreqs = findFreqs(apd)
	allPols = findPols(apd)

	if freq == '':
		if len(allFreqs) == 0:
			freq = 'any'
		elif len(allFreqs) == 1:
			freq = allFreqs[0]
	if pol == '' and len(allPols) == 1:
		pol = allPols[0]

	if az_deg != None and el_deg != None:
		g,u = gain(apd, freq, pol, az_deg, el_deg)
		print('gain(az=%.1f, el=%.1f) = %.1f %s' % (az_deg, el_deg, g, u))
	if plot == '2':
		plotNsma2d(apd, freq, pol)
	if plot == '3':
		plotNsma3d(apd, freq, pol, fc)
	if text == 's':
		miniNsma(apd)
	if text == 't':
		printNsma(apd)

if __name__ == '__main__':
	main(sys.argv)