# The ATCA Sensitivity Calculator

import argparse
import os
import math
import sys
import json
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from refract import *

# Define some global parameters.
frequencyBands = {
    # Low and high frequencies for each of the ATCA receivers, and the file to
    # read to get Tsys information.
    '16cm': { 'low': 1730, 'high': 2999,
              'tsys': "systemps/ca02_21cm_x_polarisation.avg" },
    '4cm': { 'low': 4928, 'high': 10928,
             'tsys': "systemps/ca02_4cm_y_polarisation.avg" },
    '15mm': { 'low': 16000, 'high': 25000,
              'tsys': "systemps/12mm_recvtemps.avg" },
    '7mm': { 'low': 30000, 'high': 50000,
             'tsys': "systemps/ca02_7mm.avg" },
    '3mm': { 'low': 83857, 'high': 104785,
             'tsys': "systemps/nominal_3mm.avg" }
}
sideBands = {
    # Which sideband is used per frequency range.
    'USB': [ [ 4928, 10928 ], [ 41000, 50000 ], [ 97800, 104785 ] ],
    'LSB': [ [ 1730, 2999 ], [ 16000, 25000 ], [ 30000, 40999 ],
             [ 83857, 97799 ] ]
}
continuumBandwidth = 2049.0 # MHz bandwidth of the continuum bands.
nZoomChannels = 2048 # the number of channels each zoom has
antennaDiameter = 22.0 # metres
speedoflight = 299792458.0 # Speed of light, m/s
boltzmann = 1.3806488e-23 # J / K
latitudeATCA = -30.31288472 # degrees, latitude of the ATCA.
longitudeATCA = 149.5501388 # degrees, longitude of the ATCA.
eastAngle = longitudeATCA - 90.0
cangle = math.radians(eastAngle)
cosl = math.cos(math.radians(latitudeATCA))
sinl = math.sin(math.radians(latitudeATCA))
channelFlagging = {
    # The channels always flagged in the continuum band.
    'continuum': { 'CFB1M': [ 513, 1025, 1537 ],
                   'CFB64M': [ 9, 17, 25 ] },
    # The sampling clock birdies.
    'birdies': { 'CFB1M': [ 129, 157, 257, 641, 769, 1153, 1177, 1281, 1409, 1793, 1921 ] }
}
frequencyFlagging = {
    # Known RFI ranges.
    'rfi': [ [ 1059.0, 1075.0 ], [ 1103.0, 1117.0 ], [ 1145.0, 1159.0 ], [ 1165.0, 1191.0 ],
             [ 1217.0, 1239.0 ], [ 1240.0, 1252.0 ], [ 1380.0, 1382.0 ], [ 1428.0, 1432.0 ],
             [ 1436.0, 1450.0 ], [ 1456.0, 1460.0 ], [ 1493.0, 1495.0 ], [ 1499.0, 1511.0 ],
             [ 1525.0, 1628.0 ], [ 2489.0, 2496.0 ], [ 2879.0, 2881.0 ], [ 5622.0, 5628.0 ],
             [ 5930.0, 5960.0 ], [ 6440.0, 6480.0 ], [ 7747.0, 7777.0 ], [ 7866.0, 7896.0 ],
             [ 8058.0, 8088.0 ], [ 8177.0, 8207.0 ] ]
}
# Some conversion factors.
mhzToHz = 1.0e6 # Convert MHz to Hz
degreesToArcmin = 60.0 # Convert degrees to arcminutes
degreesToArcsec = 3600.0 # Convert degrees to arcseconds
jyToSI = 1e-26 # Convert Jy to W.m^-2.s^-1
mToKm = 1e-3 # Convert m to km

# Our error handling class.
class CalcError(Exception):
    def __init__(self, value):
        self.value = value
    def __str__(self):
        return repr(self.value)

def frequencyOverlap(inf1, inf2, ing1, ing2):
    # Given two frequency ranges f1 - f2, g1 - g2, returns the amount
    # of frequency overlap.
    overlap = 0.0
    
    # Check for correct frequency order.
    f1 = inf1
    f2 = inf2
    g1 = ing1
    g2 = ing2
    if (f2 < f1):
        f1 = inf2
        f2 = inf1
    if (g2 < g1):
        g1 = ing2
        g2 = ing1

    # Assess overlap.
    if ((g1 >= f1) and (g1 <= f2)):
        a1 = g2 - g1
        a2 = f2 - g1
        overlap += min(a1, a2)
    elif ((g2 >= f1) and (g2 <= f2)):
        a1 = g2 - g1
        a2 = g2 - f1
        overlap += min(a1, a2)
    elif ((f1 >= g1) and (f1 <= g2)):
        a1 = f2 - f1
        a2 = g2 - f1
        overlap += min(a1, a2)
    elif ((f2 >= g1) and (f2 <= g2)):
        a1 = f2 - f1
        a2 = f2 - g1
        overlap += min(a1, a2)

    return overlap

def frequencyBand(freq):
    # Which band does the frequency belong to?
    for b in frequencyBands:
        if (freq >= frequencyBands[b]['low'] and
            freq <= frequencyBands[b]['high']):
            return b
    return None

def channelResolution(corrMode):
    # Return the width of each type of channel for a given correlator mode.
    if (corrMode == 'CFB1M'):
        return { 'continuum': 1.0, 'zoom': (1.0 / 2048.0) }
    elif (corrMode == 'CFB64M'):
        return { 'continuum': 64.0, 'zoom': (64.0 / 2048.0) }
    else:
        return None

def checkArguments(args):
    # Check the frequency exists and is within a known band.
    if 'frequency' in args:
        b = frequencyBand(args.frequency)
        if b is None:
            raise CalcError("Frequency is not in an ATCA band.")

    # Check that any specified zoom band is within the frequency range.
    if ('zoomfreq' in args and args.zoomfreq is not None):
        if (args.zoomfreq < (args.frequency - (continuumBandwidth / 2)) or
            args.zoomfreq > (args.frequency + (continuumBandwidth / 2))):
            raise CalcError("Zoom frequency too far from continuum central frequency.")

    # Check that the elevation limit is reasonable.
    if (args.ellimit < 12 or args.ellimit >= 90):
        raise CalcError("Elevation limit out of range.")

    # Check that the hour angle limit is reasonable.
    if (args.halimit <= 0 or args.halimit > 12):
        raise CalcError("Hour-angle limit out of range.")
    # And the middle hour angle.
    if (abs(args.ha_middle) > 12):
        raise CalcError("Middle hour-angle out of range.")
    # And thus the derived maximum and minimum hour angles.
    hourAngle_min = args.ha_middle - abs(args.halimit)
    hourAngle_max = args.ha_middle + abs(args.halimit)
    # Or the actual specified minimum and maximum hour angles.
    if (('ha_min' in args and args.ha_min is not None) and
        ('ha_max' in args and args.ha_max is not None)):
        hourAngle_min = args.ha_min
        hourAngle_max = args.ha_max
    if (abs(hourAngle_min) > 12 or abs(hourAngle_max) > 12):
        raise CalcError("Interpreted hour-angle maximum/minimum out of range.")

    # Check that the declination is reasonable.
    if (args.dec < -90 or args.dec > (90 - 30.313 - args.ellimit)):
        raise CalcError("Declination not observable with specified elevation limit.")

    # Check we have a positive integration time.
    if (args.integration <= 0):
        raise CalcError("Integration time must be greater than 0 minutes.")

    # Check the number of integrations per hour angle is reasonable.
    if (args.per_ha <= 0):
        raise CalcError("Number of integrations per hour angle is invalid.")

    # Check for a specified rest frequency.
    if ('restfreq' in args and args.restfreq is not None):
        restfreq = args.restfreq
    elif ('zoomfreq' in args and args.zoomfreq is not None):
        # Use the centre frequency of the specified zoom band.
        restfreq = args.zoomfreq
    else:
        # Use the centre frequency of the continuum band.
        restfreq = args.frequency

    # Check that the number of zooms requested is between 1 and 16 inclusive.
    if (args.zoom_width < 1 or args.zoom_width > 16):
        raise CalcError("Specified zoom width is out of range.")

    # Check for consistency with time calculation mode.
    # First check how many target bands were specified.
    ntargetband = 0
    if (args.target_continuum):
        ntargetband += 1
    if (args.target_spectral):
        ntargetband += 1
    if (args.target_zoom):
        ntargetband += 1
    if (args.target_specific_zoom):
        ntargetband += 1
    # Check how many target modes were specified.
    ntargetmode = 0
    if (args.target_flux_density):
        ntargetmode += 1
    if (args.target_brightness_temperature):
        ntargetmode += 1
    # Check how many target conditions were specified.
    ntargetcondition = 0
    if (args.target_best):
        ntargetcondition += 1
    if (args.target_typical):
        ntargetcondition += 1
    if (args.target_worst):
        ntargetcondition += 1

    # Do we operation in sensitivity target mode?
    if (args.calculate_time):
        # Yes, so we check we have a clear picture of how to do the calculation.
        if ((ntargetband != 1) or (ntargetmode != 1) or (ntargetcondition != 1)):
            if (ntargetband != 1):
                raise CalcError("Time calculation undefined: no single target band selected.")
            elif (ntargetmode != 1):
                raise CalcError("Time calculation undefined: no single target mode selected.")
            elif (ntargetcondition != 1):
                raise CalcError("Time calculation undefined: no single target condition selected.")
        # Check that we actually have a target value.
        if (args.target <= 0.0):
            raise CalcError("Target sensitivity not specified.")

    # Return the interpreted values we made.
    return { 'hourAngle_min': hourAngle_min, 'hourAngle_max': hourAngle_max,
             'restfreq': restfreq }

def makeTemplate(centreFreq, bandwidth, channelWidth):
    # Make a blank spectrum that covers the specified frequency range with
    # the correct channel resolution.
    c = []
    v = []
    f = []
    n = []
    x = []
    p = []
    if (type(bandwidth) is float):
        lowFreq = centreFreq - (bandwidth - channelWidth) / 2 
        highFreq = centreFreq + (bandwidth - channelWidth) / 2
    elif (type(bandwidth) is list):
        lowFreq = centreFreq - bandwidth[0]
        highFreq = centreFreq + bandwidth[1]
    cFreq = lowFreq
    chanNum = 1
    while (cFreq <= highFreq):
        c.append(cFreq)
        v.append(0.0)
        n.append(0)
        f.append(False)
        p.append(0.0)
        x.append(chanNum)
        cFreq += channelWidth
        chanNum += 1
    nx = np.array(x)
    for i in xrange(0, len(sideBands['LSB'])):
        if ((centreFreq >= sideBands['LSB'][i][0]) and
            (centreFreq <= sideBands['LSB'][i][1])):
            # This band is LSB, so we flip the channel numbers.
            nx = nx[::-1]
            break

    return { 'centreFrequency': np.array(c), 'value': np.array(v),
             'count': np.array(n), 'flags': f,
             'flaggedBandwidth': p,
             'channelWidth': channelWidth, 'channelNumber': nx }

def averageTemplate(template):
    # Return the average unflagged value of a template.
    varr = np.array([])
    for i in xrange(0, len(template['value'])):
        if (template['flags'][i] == False):
            varr = np.append(varr, template['value'][i])
    return np.mean(varr)

def getFreq(item):
    return item[0]

def frequencyToWavelength(freq):
    # Convert a frequency in MHz to a wavelength in metres.
    wl = speedoflight / (freq * mhzToHz)
    return (wl)

def baselineToLambda(freq, baselineLength):
    # Convert a basline length in metre to lambda for a specified frequency.
    wl = frequencyToWavelength(freq)
    bl = baselineLength / wl
    return (bl)

def primaryBeamSize(freq, diam):
    # Calculate the size of the primary beam for an antenna with a specified diameter (m),
    # at the specified frequency (MHz), in arcmin.
    pbfwhm = math.degrees((speedoflight / diam) / (freq * mhzToHz)) * degreesToArcmin
    tpbfwhm = "%.1f" % pbfwhm
    return (float(tpbfwhm))

def synthesisedBeamSize(freq, baselineLength, dec, minHa, maxHa, weightFactor):
    # Calculate the synthesised beam size given a frequency, baseline length,
    # declination and the minimum and maximum hour angles observed.
    # Returns an array where minor axis is the first element, and major
    # FWHM is the second.
    blX = baselineToLambda(freq, baselineLength['dX'])
    blY = baselineToLambda(freq, baselineLength['dY'])
    blZ = baselineToLambda(freq, baselineLength['dZ'])

    # We need the hour angles in radians.
    minHaRad = math.radians(minHa * 15.0)
    maxHaRad = math.radians(maxHa * 15.0)

    # We use the full TMS equation 4.1 here.
    # Find the hour angle for maximum u.
    umaxHa = math.degrees(math.atan2(blX, blY)) / 15.0
    if (umaxHa < minHa):
        umaxHa = minHa
    if (umaxHa > maxHa):
        umaxHa = maxHa
    umaxHaRad = math.radians(umaxHa * 15.0)
    umax = abs(math.sin(umaxHaRad) * blX + math.cos(umaxHaRad) * blY)
    # Find the hour angle for maximum v.
    vmaxHa = math.degrees(math.atan2(-1 * blY, blX)) / 15.0
    if (vmaxHa < minHa):
        vmaxHa = minHa
    if (vmaxHa > maxHa):
        vmaxHa = maxHa
    vmaxHaRad = math.radians(vmaxHa * 15.0)
    vmax = abs(-1 * math.sin(math.radians(dec)) * math.cos(vmaxHaRad) * blX +
               math.sin(math.radians(dec)) * math.sin(vmaxHaRad) * blY +
               math.cos(math.radians(dec)) * blZ)

    # The resolution is simply the inverse of umax or vmax, but we also
    # convert to arcseconds and multiply by the image weighting factor.
    ures = math.degrees(1.0 / umax) * degreesToArcsec * weightFactor
    # Check for zero vmax, which occurs when observing a source on the
    # celestial equator withouth NS baselines.
    try:
        vres = math.degrees(1.0 / vmax) * degreesToArcsec * weightFactor
    except ZeroDivisionError:
        print "FATAL: Cannot observe a declination 0 source with an EW array"
        sys.exit(-1)

    # We now take care of significant figures.
    utmp = "%.2f" % ures
    vtmp = "%.2f" % vres
    ures = float(utmp)
    vres = float(vtmp)

    # Make the output with the minor axis first.
    res = [ ures, vres ]
    if (ures > vres):
        res = [ vres, ures ]
    return (res)

def ellipseArea(minor, major):
    # Returns the area of an ellipse with specified minor and major
    # axis lengths (in arcseconds), in sr.
    mn = (minor / 2.0) / degreesToArcsec
    mj = (major / 2.0) / degreesToArcsec
    a = math.pi * mn * mj # in square degrees
    asr = a / (180.0 / math.pi)**2 # in steradians
    return (asr)

def brightnessTemperatureSensitivity(rms, synthBeam, freq):
    # Calculate the brightness temperature sensitivity (K) of an observation
    # with an RMS noise level (mJy/beam), a synthesised beam size (arcsec),
    # and a frequency in MHz. This comes from TMS equation 1.2.
    A = ellipseArea(synthBeam[0], synthBeam[1])
    I = (rms / 1000.0) * jyToSI / A
    wl = frequencyToWavelength(freq)
    bt = wl * wl * I / (2.0 * boltzmann)
    return (bt)

def readTsys(filename, lf, hf):
    # Open the filename.
    if (os.path.isfile(filename)):
        # Read it.
        d = np.loadtxt(filename)
        # Sort it.
        ds = sorted(d, key=getFreq)
        # Split it.
        startIndex = 0
        endIndex = len(ds) - 1
        llf = lf / 1000.0
        hhf = hf / 1000.0
        while ((ds[startIndex][0] < llf) and
               (startIndex < (len(ds) - 1)) and
               (ds[startIndex + 1][0] < llf)):
            startIndex += 1
        while ((ds[endIndex][0] > hhf) and
               (endIndex > 0) and
               (ds[endIndex - 1][0] > hhf)):
            endIndex -= 1
        if (startIndex == endIndex):
            if (startIndex > 0):
                startIndex -= 1
            if (endIndex < (len(ds) - 1)):
                endIndex += 1
        dds = ds[startIndex:(endIndex + 1)]
        c = [np.around(row[0] * 1000.0) for row in dds]
        v = [(10 ** row[1]) for row in dds]
        n = []
        f = []
        for i in xrange(0, len(v)):
            n.append(1)
            f.append(False)
        return { 'centreFrequency': c, 'value': v, 'count': np.array(n),
                 'flags': f,
                 'channelWidth': 1.0 }
    else:
        raise CalcError("Can't find Tsys file %s." % filename)

def lowHigh(c, w):
    return { 'low': c - (w / 2.0), 'high': c + (w / 2.0) }

def overlaps(a, b):
    if (b['low'] >= a['high']):
        return False
    if (b['high'] <= a['low']):
        return False
    return True

def templateAverage(t):
    # Divide the values by the counts.
    for i in xrange(0, len(t['centreFrequency'])):
        if (t['count'][i] > 0):
            t['value'][i] /= float(t['count'][i])
            if (isinstance(t['flags'][i], list) == True):
                # Multiple input flags, we choose the most common.
                ft = 0
                ff = 0
                for j in xrange(0, len(t['flags'][i])):
                    if (t['flags'][i][j]):
                        ft += 1
                    else:
                        ff += 1
                if (ft > ff):
                    t['flags'][i] = True
                else:
                    t['flags'][i] = False

def linearInterpolate(p1, p2, pi):
    # Using information from p1 and p2, determine the value at pi.
    run = p2['frequency'] - p1['frequency']
    if (run != 0):
        slope = (p2['value'] - p1['value']) / run
    else:
        slope = 0
    nrun = pi['frequency'] - p1['frequency']
    return p1['value'] + nrun * slope

def templateInterpolate(t):
    # Interpolate values for channels with no counts.
    lg = 0
    for i in xrange(0, len(t['centreFrequency'])):
        if (t['count'][i] > 0):
            # The last good channel.
            lg = i
        else:
            # Find the next good channel.
            for j in xrange(i + 1, len(t['centreFrequency'])):
                if (t['count'][j] > 0):
                    hg = j
                    break
            t['value'][i] = linearInterpolate({ 'frequency': t['centreFrequency'][lg],
                                                'value': t['value'][lg] },
                                              { 'frequency': t['centreFrequency'][hg],
                                                'value': t['value'][hg] },
                                              { 'frequency': t['centreFrequency'][i] })
    
def templateFill(srcTemplate, destTemplate):
    # Fill in a template spectrum with values from another template, and do it with
    # a single pass of each array (no looping).
    i = 0 # The index of the destination template bin
    j = 0 # The index of the source template bin
    sfs = lowHigh(srcTemplate['centreFrequency'][j], srcTemplate['channelWidth'])
    dfs = lowHigh(destTemplate['centreFrequency'][i], destTemplate['channelWidth'])
    while (i < len(destTemplate['centreFrequency']) and
           j < len(srcTemplate['centreFrequency'])):
        if (overlaps(dfs, sfs)):
            destTemplate['value'][i] += srcTemplate['value'][j]
            destTemplate['count'][i] += 1
            if (destTemplate['count'][i] > 1):
                if (destTemplate['count'][i] == 2):
                    t = [ destTemplate['flags'][i], srcTemplate['flags'][j] ]
                    destTemplate['flags'][i] = t
                else:
                    destTemplate['flags'][i].append(srcTemplate['flags'][j])
            else:
                destTemplate['flags'][i] = srcTemplate['flags'][j]
            j += 1
            if (j < len(srcTemplate['centreFrequency'])):
                sfs = lowHigh(srcTemplate['centreFrequency'][j], srcTemplate['channelWidth'])
        elif (srcTemplate['centreFrequency'][j] < destTemplate['centreFrequency'][i]):
            j += 1
            if (j < len(srcTemplate['centreFrequency'])):
                sfs = lowHigh(srcTemplate['centreFrequency'][j], srcTemplate['channelWidth'])
        else:
            i += 1
            if (i < len(destTemplate['centreFrequency'])):
                dfs = lowHigh(destTemplate['centreFrequency'][i], destTemplate['channelWidth'])

    templateAverage(destTemplate)

    # Check that the edges aren't empty
    # Bottom edge.
    if (destTemplate['count'][0] == 0):
        # Have to interpolate from the source template.
        dfs = lowHigh(destTemplate['centreFrequency'][0], destTemplate['channelWidth'])
        j = 0
        sfs = lowHigh(srcTemplate['centreFrequency'][j], srcTemplate['channelWidth'])
        while (sfs['high'] < dfs['low']):
            j += 1
            sfs = lowHigh(srcTemplate['centreFrequency'][j], srcTemplate['channelWidth'])
        j -= 1 # Because the breaking point is when the frequency goes too high.
        i = 1
        while (destTemplate['count'][i] == 0):
            i += 1
        destTemplate['value'][0] = linearInterpolate({ 'value': srcTemplate['value'][j],
                                                       'frequency': srcTemplate['centreFrequency'][j] },
                                                     { 'value': destTemplate['value'][i],
                                                       'frequency': destTemplate['centreFrequency'][i] },
                                                     { 'frequency': destTemplate['centreFrequency'][0] } )
        destTemplate['count'][0] = 1
    # Top edge.
    if (destTemplate['count'][-1] == 0):
        dfs = lowHigh(destTemplate['centreFrequency'][-1], destTemplate['channelWidth'])
        j = 0
        sfs = lowHigh(srcTemplate['centreFrequency'][j], srcTemplate['channelWidth'])
        while (sfs['high'] <= dfs['low']):
            j += 1
            if (j < len(srcTemplate['centreFrequency'])):
                sfs = lowHigh(srcTemplate['centreFrequency'][j], srcTemplate['channelWidth'])
            else:
                j -= 1
                break
        # Break point is fine this time.
        i = -2
        while (destTemplate['count'][i] == 0 and (abs(i) < (len(destTemplate['centreFrequency'])))):
            i -= 1
        destTemplate['value'][-1] = linearInterpolate({ 'value': srcTemplate['value'][j],
                                                        'frequency': srcTemplate['centreFrequency'][j] },
                                                      { 'value': destTemplate['value'][i],
                                                        'frequency': destTemplate['centreFrequency'][i] },
                                                      { 'frequency': destTemplate['centreFrequency'][-1] } )
        destTemplate['count'][-1] = 1

    templateInterpolate(destTemplate)

def templateEfficiency():
    # The template returned by this routine contains all the efficiencies for
    # all the bands.
    c = [   900.0,   1200.0,   1500.0,   1800.0,  2100.0,  2300.0,  2500.0,  4400.0,  5900.0,
           7400.0,   8800.0,  10600.0,  16000.0, 16500.0, 17000.0, 17500.0, 18000.0, 18500.0,
          19000.0,  19500.0,  20000.0,  20500.0, 21000.0, 21500.0, 22000.0, 22500.0, 23000.0,
          23500.0,  24000.0,  24500.0,  25000.0, 25400.0, 30000.0, 31000.0, 32000.0, 33000.0,
          34000.0,  35000.0,  36000.0,  37000.0, 38000.0, 39000.0, 40000.0, 41000.0, 42000.0,
          43000.0,  44000.0,  45000.0,  46000.0, 47000.0, 48000.0, 49000.0, 50000.0, 83781.1,
          85556.2,  86834.3,  88680.5,  90526.6, 91946.7, 94005.9, 95852.1, 97272.2, 98976.3,
         100254.4, 102200.0, 102300.0, 106432.0 ]
    v = [  0.57,     0.57,     0.60,     0.53,    0.43,    0.42,    0.44,    0.65,    0.72,
           0.65,     0.64,     0.65,     0.58,    0.62,    0.63,    0.65,    0.67,    0.70,
           0.68,     0.64,     0.64,     0.60,    0.53,    0.55,    0.54,    0.51,    0.51,
           0.53,     0.49,     0.49,     0.46,    0.47,    0.60,    0.60,    0.60,    0.60,
           0.60,     0.60,     0.60,     0.60,    0.60,    0.60,    0.60,    0.59,    0.58,
           0.57,     0.56,     0.55,     0.54,    0.53,    0.52,    0.51,    0.50,    0.3297,
           0.3065,   0.3020,   0.2856,   0.2689,  0.2670,  0.2734,  0.2727,  0.2521,  0.2403,
           0.2336,   0.2322,   0.14,     0.14 ]
    f = []
    n = []
    for i in xrange(0, len(c)):
        n.append(1)
        f.append(False)
    return { 'centreFrequency': np.array(c), 'value': np.array(v),
             'count': np.array(n), 'flags': f,
             'channelWidth': 1.0 }

def fillAtmosphereTemplate(templateOpacity, templateTemperature, t, p, h):
    # Calculate the opacity and atmospheric temperature at the zenith for each frequency
    # in the template.
    atmos = calcOpacity(templateOpacity['centreFrequency'] * 1e6, math.radians(90.0), t, p, h)
    templateOpacity['value'] = np.array(atmos['tau'])
    templateTemperature['value'] = np.array(atmos['Tb'])

def plotTemplate(t, e, outname):
    plt.clf()
    plt.plot(t['centreFrequency'], t['value'])
    plt.plot(e['centreFrequency'], (e['value']), "green")
    plt.savefig(outname)

def plotSpectrum(template, conditions, outname):
    # Plot the template spectrum we are passed with frequency on the x-axis.
    # Initialise the plot.
    plt.clf()
    fig = plt.figure()
    ax = fig.add_subplot(111)

    # These are the colours we can use for the different lines.
    colours = [ "blue", "green", "red", "black", "yellow" ]

    # Go through the conditions (usually weather conditions) and plot a line for each.
    for i, c in enumerate(conditions):
        ax.plot(template[c]['centreFrequency'], template[c]['value'], colours[i], label=c)

    # Set the x-axis limits to be tight on the actual frequency range.
    plt.xlim(template[c]['centreFrequency'][0], template[c]['centreFrequency'][-1])
    plt.xlabel("Frequency [MHz]")
    plt.ylabel("RMS noise level [mJy/beam]")

    # Ensure that the x- and y-axes don't have an offset value.
    ax.get_yaxis().get_major_formatter().set_useOffset(False)
    ax.get_xaxis().get_major_formatter().set_useOffset(False)

    # Put the legend with the condition names at the top of the plot outside the border.
    plt.legend(bbox_to_anchor=(0., 1.02, 1., .102), loc=3,
               ncol=len(c), mode="expand", borderaxespad=0.)

    # Highlight the regions that are flagged in the template.
    for i in xrange(0, len(template[conditions[0]]['flags'])):
        if (template[conditions[0]]['flags'][i] == True):
            x1 = template[conditions[0]]['centreFrequency'][i] - template[c]['channelWidth'] / 2.0
            x2 = template[conditions[0]]['centreFrequency'][i] + template[c]['channelWidth'] / 2.0
            plt.axvspan(x1, x2, alpha=0.2, edgecolor='none', facecolor='red')
    plt.savefig(outname)

def flagTemplate(t, flagType, corrMode, edgeChan):
    # Set the flags in the template.
    if ((flagType == "continuum") or (flagType == "birdies")):
        # This is channel based flagging.
        flagSrc = channelFlagging[flagType]
        if (corrMode in flagSrc):
            flagSrc = flagSrc[corrMode]
            for c in xrange(0, len(flagSrc)):
                for i in xrange(0, len(t['flags'])):
                    if (t['channelNumber'][i] == flagSrc[c]):
                        t['flags'][i] = True
                        break
    elif (flagType == "rfi"):
        # This is frequency based flagging.
        oldi = 0
        flagSrc = frequencyFlagging[flagType]
        for r in xrange(0, len(flagSrc)):
            for i in xrange(oldi, len(t['flags'])):
                f1 = t['centreFrequency'][i] - t['channelWidth'] / 2.0
                f2 = t['centreFrequency'][i] + t['channelWidth'] / 2.0
                if (f1 > flagSrc[r][1]):
                    break
                t['flaggedBandwidth'][i] += frequencyOverlap(f1, f2, flagSrc[r][0], flagSrc[r][1])
                oldi = i
        # Check for which channels are above the percentage flagged cut.
        for i in xrange(0, len(t['flags'])):
            p = t['flaggedBandwidth'][i] / t['channelWidth']
            if (p > 0.5):
                t['flags'][i] = True
    elif (flagType == "edge"):
        # Flag edge channels.
        for i in xrange(0, len(t['channelNumber'])):
            cn = i + 1
            rcn = len(t['channelNumber']) - i
            if ((cn <= edgeChan) or
                (rcn <= edgeChan)):
                t['flags'][i] = True


def calculateSensitivity(rmsTemplate, nAnts):
    # Given a template filled with the RMS noise in each channel in the continuum
    # band, return the average sensitivity of the continuum band channels (the spectral
    # RMS), and the spectral RMS divided by the bandwidth (the continuum RMS). This
    # routine does this while respecting any flagging present in the template.
    rmsTotal = 0.0
    rmsZoomTotal = 0.0
    rmsN = 0
    rmsZoomN = 0
    for i in xrange(0, len(rmsTemplate['value'])):
        # We count flagged continuum channels in the zoom RMS calculation.
        rmsZoomTotal += rmsTemplate['value'][i]
        rmsZoomN += 1
        if (rmsTemplate['flags'][i] == False):
            rmsTotal += rmsTemplate['value'][i]
            rmsN += 1
    
    if (rmsN > 0):
        rmsSpectral = rmsTotal / rmsN
        totalBandwidth = rmsN * rmsTemplate['channelWidth']
        rmsContinuum = rmsSpectral / math.sqrt(float(rmsN))
        rmsZoom = rmsSpectral * math.sqrt(float(nZoomChannels))
    else:
        rmsSpectral = None
        totalBandwidth = None
        rmsContinuum = None
        rmsZoom = None

    # Calculate the SEFDs for each antenna and the array as a whole.
    # We use TMS equation 1.6 for this and convert to Jy.
    Aone = surfaceArea(antennaDiameter)
    Aall = nAnts * Aone
    sefdOne = 2.0 * boltzmann * rmsTemplate['systemp'] / (Aone * 1e-26)
    sefdAll = 2.0 * boltzmann * rmsTemplate['systemp'] / (Aall * 1e-26)

    return { 'rms': { 'spectral': rmsSpectral, 'continuum': rmsContinuum, 'zoom': rmsZoom },
             'bandwidth': { 'unflagged': totalBandwidth }, 'sefd': { 'antenna': sefdOne,
                                                                     'array': sefdAll } }

def calculateRms(tsys, efficiency, opacity, temperature, minHa, maxHa, perHa, nAntenna,
                 totalTime, weighting, sind, cosd):
    # Given the tsys and efficiency templates, the number of antennas involved in the
    # imaging, the total integration time and the image weighting scheme, this routine
    # will return another template with each channel being the RMS noise expected in
    # that channel. This comes from eqn 6.62 of TMS, where eta_Q is 1 (for CABB's
    # digitisation) but A is multiplied by our efficiency factor. That equation is
    # for only a single polarisation though, so for an unpolarised source, the noise
    # level is sqrt(2) lower, which is where the sqrt(2) factor in the numerator comes from
    # instead of the 2.

    # The variables for our output template.
    c = []
    v = []
    n = []
    f = []

    # Figure out how much time is spent in each integration period.
    haRange = maxHa - minHa
    nIntegrations = math.ceil(haRange * perHa)
    intTime = totalTime / nIntegrations

    Texcess = []
    systemTemperature = []
    for j in xrange(0, int(nIntegrations + 1)):
        # The hour angle at this integration.
        jHa = minHa + float(j) / perHa
        
        # The elevation at this hour angle.
        cosha = math.cos(math.radians(jHa * 15.0))
        sinel = sinl * sind + cosl * cosd * cosha
        
        # Calculate the excess temperature due to the atmosphere
        # and CMB.
        elFactor = np.exp(-1.0 * opacity['value'] / sinel)
        cbFactor = 2.7 * elFactor
        ivFactor = 1.0 - elFactor
        atFactor = temperature['value'] * ivFactor
        Texcess.append(atFactor + cbFactor)

    for i in xrange(0, len(tsys['centreFrequency'])):
        # Check that the frequencies are the same in both templates.
        if (tsys['centreFrequency'][i] == efficiency['centreFrequency'][i]):
            # All's good.
            c.append(tsys['centreFrequency'][i])

            # Get the average excess temperature now.
            exTemp = [ row[i] for row in Texcess ]
            excessTemp = math.fsum(exTemp) / float(len(exTemp))

            Tmeas = tsys['value'][i] + excessTemp
            TmeasEff = Tmeas / efficiency['value'][i]
            
            # The units of this is actually mJy since we keep the frequency
            # in MHz rather than converting to Hz (convenient isn't it!).
            v.append((math.sqrt(2.0) * boltzmann * Tmeas * weighting['avg']) /
                     (1e-26 * surfaceArea(antennaDiameter) *
                      efficiency['value'][i] *
                      math.sqrt(float(nAntenna) * float(nAntenna - 1) *
                                tsys['channelWidth'] * (totalTime * 60.0))))
            n.append(1)
            if (tsys['flags'][i] or efficiency['flags'][i]):
                f.append(True)
            else:
                f.append(False)
                systemTemperature.append(TmeasEff)

    return { 'centreFrequency': np.array(c), 'value': np.array(v),
             'count': np.array(n), 'flags': f,
             'channelWidth': tsys['channelWidth'], 'channelNumber': tsys['channelNumber'],
             'systemp': np.mean(systemTemperature),
             'systemTemperature': float("%.1f" % np.mean(systemTemperature))}

def surfaceArea(d):
    # Given the diameter of a dish (m), return its surface area (m^2).
    return (math.pi * ((d / 2.0) ** 2))

def weightingFactor(weighting, array, ant6):
    # Return the w_rms / w_mean weighting factors given the weighting scheme,
    # the array configuration and whether antenna 6 is included in the imaging.
    
    # The weighting factors we determined from simulations.
    factors = [
        { 'array': [ '6000', '6km', '3000', '3km' ],
          'factors': [
              { 'ca06': True,
                'R2': { 'avg': 1.039, 'min': 1.000, 'max': 1.079, 'beam': 1.32 },
                'R1': { 'avg': 1.040, 'min': 1.000, 'max': 1.080, 'beam': 1.32 },
                'R0': { 'avg': 1.871, 'min': 1.350, 'max': 2.781, 'beam': 0.84 },
                'R-1': { 'avg': 5.791, 'min': 3.685, 'max': 10.987, 'beam': 0.80 },
                'R-2': { 'avg': 5.847, 'min': 3.688, 'max': 11.240, 'beam': 0.80 } },
              { 'ca06': False,
                'R2': { 'avg': 1.000, 'min': 1.000, 'max': 1.000, 'beam': 0.97 },
                'R1': { 'avg': 1.002, 'min': 1.001, 'max': 1.004, 'beam': 0.89 },
                'R0': { 'avg': 1.882, 'min': 1.703, 'max': 1.943, 'beam': 0.66 },
                'R-1': { 'avg': 3.875, 'min': 2.543, 'max': 7.102, 'beam': 0.64 },
                'R-2': { 'avg': 3.908, 'min': 2.562, 'max': 7.222, 'beam': 0.64 } } ] },
        { 'array': [ '1500', '1.5km' ],
          'factors': [
              { 'ca06': True,
                'R2': { 'avg': 1.000, 'min': 1.000, 'max': 1.000, 'beam': 1.32 },
                'R1': { 'avg': 1.000, 'min': 1.000, 'max': 1.001, 'beam': 1.32 },
                'R0': { 'avg': 1.507, 'min': 1.181, 'max': 1.846, 'beam': 0.84 },
                'R-1': { 'avg': 7.925, 'min': 5.200, 'max': 16.732, 'beam': 0.80 },
                'R-2': { 'avg': 8.151, 'min': 5.163, 'max': 19.304, 'beam': 0.80 } },
              { 'ca06': False,
                'R2': { 'avg': 1.000, 'min': 1.000, 'max': 1.000, 'beam': 0.96 },
                'R1': { 'avg': 1.001, 'min': 1.000, 'max': 1.003, 'beam': 0.88 },
                'R0': { 'avg': 1.854, 'min': 1.576, 'max': 1.953, 'beam': 0.64 },
                'R-1': { 'avg': 3.900, 'min': 2.524, 'max': 8.218, 'beam': 0.62 },
                'R-2': { 'avg': 3.923, 'min': 2.506, 'max': 8.707, 'beam': 0.62 } } ] },
        { 'array': [ '750', '750m' ],
          'factors': [
              { 'ca06': True,
                'R2': { 'avg': 1.000, 'min': 1.000, 'max': 1.000, 'beam': 1.32 },
                'R1': { 'avg': 1.000, 'min': 1.000, 'max': 1.000, 'beam': 1.32 },
                'R0': { 'avg': 1.299, 'min': 1.143, 'max': 1.621, 'beam': 0.84 },
                'R-1': { 'avg': 12.893, 'min': 8.581, 'max': 17.674, 'beam': 0.80 },
                'R-2': { 'avg': 14.027, 'min': 8.882, 'max': 22.273, 'beam': 0.80 } },
              { 'ca06': False,
                'R2': { 'avg': 1.000, 'min': 1.000, 'max': 1.000, 'beam': 0.96 },
                'R1': { 'avg': 1.001, 'min': 1.000, 'max': 1.002, 'beam': 0.88 },
                'R0': { 'avg': 1.925, 'min': 1.850, 'max': 1.971, 'beam': 0.62 },
                'R-1': { 'avg': 3.557, 'min': 2.578, 'max': 5.255, 'beam': 0.59 },
                'R-2': { 'avg': 3.582, 'min': 2.583, 'max': 5.369, 'beam': 0.59 } } ] },
        { 'array': [ '367', 'EW352', 'EW367', 'EW352/367' ],
          'factors': [
              { 'ca06': True,
                'R2': { 'avg': 1.000, 'min': 1.000, 'max': 1.000, 'beam': 1.32 },
                'R1': { 'avg': 1.000, 'min': 1.000, 'max': 1.000, 'beam': 1.32 },
                'R0': { 'avg': 1.077, 'min': 1.029, 'max': 1.157, 'beam': 0.84 },
                'R-1': { 'avg': 18.304, 'min': 16.432, 'max': 17.498, 'beam': 0.80 },
                'R-2': { 'avg': 31.295, 'min': 20.574, 'max': 52.204, 'beam': 0.80 } },
              { 'ca06': False,
                'R2': { 'avg': 1.000, 'min': 1.000, 'max': 1.000, 'beam': 1.09 },
                'R1': { 'avg': 1.001, 'min': 1.000, 'max': 1.001, 'beam': 1.00 },
                'R0': { 'avg': 1.917, 'min': 1.838, 'max': 1.965, 'beam': 0.68 },
                'R-1': { 'avg': 3.271, 'min': 2.537, 'max': 4.639, 'beam': 0.64 },
                'R-2': { 'avg': 3.298, 'min': 2.550, 'max': 4.718, 'beam': 0.64 } } ] },
        { 'array': [ 'h214', 'H214', 'h168', 'H168', 'h75', 'H75' ],
          'factors': [
              { 'ca06': True,
                'R2': { 'avg': 1.000, 'min': 1.000, 'max': 1.000, 'beam': 1.32},
                'R1': { 'avg': 1.000, 'min': 1.000, 'max': 1.000, 'beam': 1.32},
                'R0': { 'avg': 1.106, 'min': 1.023, 'max': 1.186, 'beam': 0.84},
                'R-1': { 'avg': 16.865, 'min': 15.629, 'max': 18.294, 'beam': 0.80 },
                'R-2': { 'avg': 26.926, 'min': 18.717, 'max': 58.094, 'beam': 0.80 } },
              { 'ca06': False,
                'R2': { 'avg': 1.000, 'min': 1.000, 'max': 1.000, 'beam': 0.75 },
                'R1': { 'avg': 1.001, 'min': 1.001, 'max': 1.002, 'beam': 0.74 },
                'R0': { 'avg': 1.641, 'min': 1.529, 'max': 1.760, 'beam': 0.61 },
                'R-1': { 'avg': 1.984, 'min': 1.753, 'max': 2.281, 'beam': 0.59 },
                'R-2': { 'avg': 1.988, 'min': 1.755, 'max': 2.288, 'beam': 0.59 } } ] } ]

    # Search through the array.
    for f in xrange(0, len(factors)):
        arrayFound = False
        for i in xrange(0, len(factors[f]['array'])):
            if (array == factors[f]['array'][i]):
                arrayFound = True
                break
        if (arrayFound):
            g = 1
            if (ant6):
                g = 0
            if (weighting in factors[f]['factors'][g]):
                return (factors[f]['factors'][g][weighting])
    # We couldn't find the appropriate weighting factors.
    raise CalcError("Weighting factors are not available.")

def maximumBaseline(array):
    # Return the maximum baseline length for a named array, both on the track
    # only, and including CA06.

    stationLocations = {
        'W0': [ -4752438.459, 2790321.299, -3200483.747 ],
        'W2': [ -4752422.922, 2790347.675, -3200483.747 ],
        'W4': [ -4752407.385, 2790374.052, -3200483.747 ],
        'W6': [ -4752391.848, 2790400.428, -3200483.747 ],
        'W8': [ -4752376.311, 2790426.804, -3200483.747 ],
        'W10': [ -4752360.774, 2790453.181, -3200483.747 ],
        'W12': [ -4752345.237, 2790479.557, -3200483.747 ],
        'W14': [ -4752329.700, 2790505.934, -3200483.747 ],
        'W16': [ -4752314.163, 2790532.310, -3200483.747 ],
        'W32': [ -4752189.868, 2790743.321, -3200483.747 ],
        'W45': [ -4752088.877, 2790914.767, -3200483.747 ],
        'W64': [ -4751941.276, 2791165.342, -3200483.747 ],
        'W84': [ -4751785.907, 2791429.106, -3200483.747 ],
        'W98': [ -4751677.148, 2791613.741, -3200483.747 ],
        'W100': [ -4751661.611, 2791640.117, -3200483.747 ],
        'W102': [ -4751646.074, 2791666.493, -3200483.747 ],
        'W104': [ -4751630.537, 2791692.870, -3200483.747 ],
        'W106': [ -4751615.000, 2791719.246, -3200483.747 ],
        'W109': [ -4751591.695, 2791758.810, -3200483.747 ],
        'W110': [ -4751583.926, 2791771.999, -3200483.747 ],
        'W111': [ -4751576.158, 2791785.187, -3200483.747 ],
        'W112': [ -4751568.389, 2791798.375, -3200483.747 ],
        'W113': [ -4751560.621, 2791811.563, -3200483.747 ],
        'W124': [ -4751475.168, 2791956.633, -3200483.747 ],
        'W125': [ -4751467.399, 2791969.821, -3200483.747 ],
        'W128': [ -4751444.094, 2792009.386, -3200483.747 ],
        'W129': [ -4751436.325, 2792022.574, -3200483.747 ],
        'W140': [ -4751350.872, 2792167.644, -3200483.747 ],
        'W147': [ -4751296.492, 2792259.961, -3200483.747 ],
        'W148': [ -4751288.724, 2792273.149, -3200483.747 ],
        'W163': [ -4751172.197, 2792470.972, -3200483.747 ],
        'W168': [ -4751133.354, 2792536.913, -3200483.747 ],
        'W172': [ -4751102.281, 2792589.666, -3200483.747 ],
        'W173': [ -4751094.512, 2792602.854, -3200483.747 ],
        'W182': [ -4751024.596, 2792721.547, -3200483.747 ],
        'W189': [ -4750970.216, 2792813.865, -3200483.747 ],
        'W190': [ -4750962.448, 2792827.053, -3200483.747 ],
        'W195': [ -4750923.605, 2792892.994, -3200483.747 ],
        'W196': [ -4750915.837, 2792906.182, -3200483.747 ],
        'W392': [ -4749393.198, 2795491.050, -3200483.694 ],
        'N2': [ -4751628.291, 2791727.075, -3200457.305 ],
        'N5': [ -4751648.226, 2791738.818, -3200417.642 ],
        'N7': [ -4751661.517, 2791746.647, -3200391.200 ],
        'N11': [ -4751688.098, 2791762.304, -3200338.316 ],
        'N14': [ -4751708.034, 2791774.047, -3200298.653 ] }

    endStations = [
        { 'array': [ '6000', '6km', '3000', '3km' ],
          'stations': [ 'W2', 'W196' ] },
        { 'array': [ '1500', '1.5km' ],
          'stations': [ 'W98', 'W195' ] },
        { 'array': [ '750', '750m' ],
          'stations': [ 'W98', 'W148' ] },
        { 'array': [ '367', 'EW367', 'EW352/367' ],
          'stations': [ 'W104', 'W128' ] },
        { 'array': [ 'EW352' ],
          'stations': [ 'W102', 'W125' ] },
        { 'array': [ 'h214', 'H214' ],
          'stations': [ 'W98', 'W113', 'W104', 'N14' ] },
        { 'array': [ 'h168', 'H168' ],
          'stations': [ 'W100', 'W111', 'W104', 'N11' ] },
        { 'array': [ 'h75', 'H75' ],
          'stations': [ 'W104', 'W109', 'W104', 'N5' ] } ]

    # Search through the array.
    for s in xrange(0, len(endStations)):
        if (array in endStations[s]['array']):
            # Calculate the maximum dX, dY, dZ.
            mdX = 0
            mdY = 0
            mdZ = 0
            for i in xrange(0, len(endStations[s]['stations']) - 1):
                for j in xrange(i + 1, len(endStations[s]['stations'])):
                    dYp = abs(stationLocations[endStations[s]['stations'][i]][0] -
                              stationLocations[endStations[s]['stations'][j]][0])
                    dXp = abs(stationLocations[endStations[s]['stations'][i]][1] -
                              stationLocations[endStations[s]['stations'][j]][1])
                    dX = dXp * math.cos(cangle) - dYp * math.sin(cangle)
                    dY = dXp * math.sin(cangle) + dYp * math.cos(cangle)
                    if (dX > mdX):
                        mdX = dX
                    if (dY > mdY):
                        mdY = dY
                    dZ = abs(stationLocations[endStations[s]['stations'][i]][2] -
                             stationLocations[endStations[s]['stations'][j]][2])
                    if (dZ > mdZ):
                        mdZ = dZ
            mdX6 = 0
            mdY6 = 0
            mdZ6 = 0
            for i in xrange(0, len(endStations[s]['stations'])):
                dYp = abs(stationLocations[endStations[s]['stations'][i]][0] -
                          stationLocations['W392'][0])
                dXp = abs(stationLocations[endStations[s]['stations'][i]][1] -
                          stationLocations['W392'][1])
                dX = dXp * math.cos(cangle) - dYp * math.sin(cangle)
                dY = dXp * math.sin(cangle) + dYp * math.cos(cangle)
                if (dX > mdX6):
                    mdX6 = dX
                if (dY > mdY6):
                    mdY6 = dY
                dZ = abs(stationLocations[endStations[s]['stations'][i]][2] -
                         stationLocations['W392'][2])
                if (dZ > mdZ6):
                    mdZ6 = dZ
            return ( { 'track': { 'dX': mdX, 'dY': mdY, 'dZ': mdZ },
                       'ca06': { 'dX': mdX6, 'dY': mdY6, 'dZ': mdZ6 } } )

    # We couldn't find the appropriate maximum baseline lengths.
    raise CalcError("Baseline lengths are not available.")

def addToOutput(output, item, name, value, description, units):
    # Add some information to the JSON output object.
    if (isinstance(name, list) == False):
        output[item][name] = value
        output['description'][name] = description
        if (units is not None):
            output['units'][name] = units
    else:
        if (name[0] not in output[item]):
            output[item][name[0]] = {}
            output['description'][name[0]] = description
            if (units is not None):
                output['units'][name[0]] = units
        output[item][name[0]][name[1]] = value

def bandwidthToVelocity(lfreq, bw, restfreq):
    # Given a rest frequency in MHz, calculate the velocity span represented
    # by the specified bandwidth (MHz) starting at some low frequency (MHz).
    vr = bw * restfreq * speedoflight * mToKm / (lfreq * (lfreq + bw))
    vro = float("%.3f" % vr)
    return (vro)

def thingToString(t):
    if (type(t) is str):
        return (t)
    elif (type(t) is list):
        narr = []
        for i in xrange(0, len(t)):
            narr.append(thingToString(t[i]))
        return (', '.join(narr))
    elif (type(t) is int):
        return ("%d" % t)
    elif (type(t) is float or type(t) is np.float64):
        return ("%f" % t)

def humanOutputDict(d, description, units, l):
    for p in d:
        if (p != "description" and p != "units"):
            if (type(d[p]) is dict):
                print "%s%s:" % ((" " * l), p)
                humanOutputDict(d[p], description, units, (l + 1))
            else:
                od = p
                if (p in description):
                    od = description[p]
                ou = ""
                if (p in units):
                    ou = units[p]
                ov = thingToString(d[p])
                print "%s %s = %s %s" % ((" " * l), od, ov, ou)

def main(args):
    ####################################################################################################
    # Do some argument checking first.
    try:
        argsInterpreted = checkArguments(args)
    except CalcError as c:
        print "FATAL: ", c.value
        sys.exit(-1)

    # The number of antenna and baselines.
    nant = 5
    if (args.ca06):
        nant = 6
    nbaselines = nant * (nant - 1) / 2

    # Are we making calculations for a specific zoom band?
    specificZoomCalc = False
    if ('zoomfreq' in args and args.zoomfreq is not None):
        specificZoomCalc = True

    # Identify the constraints on the observed hour angles, by comparing them to the
    # specified elevation limit.
    sinel = math.sin(math.radians(args.ellimit))
    sind = math.sin(math.radians(args.dec))
    cosd = math.cos(math.radians(args.dec))

    hourAngle_min = argsInterpreted['hourAngle_min']
    hourAngle_max = argsInterpreted['hourAngle_max']

    coshaAtElLimit = sinel / (cosd * cosl) - (sind * sinl) / (cosd * cosl)
    if (abs(coshaAtElLimit) <= 1):
        haAtElLimit = math.degrees(math.acos(coshaAtElLimit)) # in degrees
        hahAtElLimit = haAtElLimit / 15.0 # in hours
        if (hahAtElLimit < abs(hourAngle_min)):
            hourAngle_min = -1 * hahAtElLimit
        if (hahAtElLimit < abs(hourAngle_max)):
            hourAngle_max = hahAtElLimit

    # Get the image weighting factors.
    try:
        imageWeights = weightingFactor(args.weighting, args.configuration, args.ca06)
    except CalcError as c:
        print "FATAL: ", c.value
        sys.exit(-1)

    # Get the lengths of the baselines and determine the maximum baseline length.
    try:
        baselineLengths = maximumBaseline(args.configuration)
    except CalcError as c:
        print "FATAL: ", c.value
        sys.exit(-1)
    maxBaselineLength = baselineLengths['track']
    if (args.ca06):
        maxBaselineLength = baselineLengths['ca06']

    # All parameters are valid.
    print "MESSAGE: Starting calculator."
    workArea = {}
    #\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\



    ####################################################################################################
    # Prepare the output, and begin by listing the parameters we were given.
    output = {
        'parameters': {},
        'continuum': {},
        'zoom': {},
        'specific_zoom': {},
        'sensitivities': {},
        'description': {},
        'source_imaging': {},
        'units': {} }
    # The array configuration.
    addToOutput(output, 'parameters', 'configuration', args.configuration, "Configuration", None)
    # The central frequency of the continuum band.
    addToOutput(output, 'parameters', 'central_frequency', args.frequency, "Central Frequency", "MHz")
    # The number of antenna in the array.
    addToOutput(output, 'parameters', 'n_antenna', nant, "Antenna Included", None)
    # The number of baselines formed by those antenna.
    addToOutput(output, 'parameters', 'n_baselines', nbaselines, "Baselines formed", None)
    # Whether the array is hybrid, or entirely East-West.
    addToOutput(output, 'parameters', 'hybrid', False, "Hybrid Array", None)
    if (args.configuration in [ 'h214', 'H214', 'h168', 'H168', 'h75', 'H75' ]):
        output['parameters']['hybrid'] = True
    # The rest frequency of the spectral line of interest.
    addToOutput(output, 'parameters', 'reference_rest_frequency', argsInterpreted['restfreq'],
                "Reference Rest Frequency", "MHz")

    # Source imaging parameters.
    # The declination of the source of interest.
    addToOutput(output, 'source_imaging', 'source_declination', args.dec, "Source Declination", "degrees")
    # The uv weighting scheme to use during imaging.
    addToOutput(output, 'source_imaging', 'weighting_scheme', args.weighting,
                "Weighting Scheme", None)
    # And the weighting factors associated with that weighting scheme.
    # The RMS noise level weighting factor, where the factor = 1 for natural weighting.
    addToOutput(output, 'source_imaging', 'weighting_factor', imageWeights['avg'],
                "RMS noise weighting factor", None)
    # The synthesised beam size weighting factor.
    addToOutput(output, 'source_imaging', 'beam_weighting_factor', imageWeights['beam'],
                "Synthesised beam weighting factor", None)
    #\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\



    ####################################################################################################
    # Calculate the frequency resolutions and put that in the output.
    workArea['resolutions'] = channelResolution(args.corrconfig)
    # These resolutions are the native resolutions made by the correlator, unsmoothed.
    # In the continuum band.
    addToOutput(output, 'continuum', 'correlator_channel_bandwidth', workArea['resolutions']['continuum'],
                "Correlator Channel Bandwidth", "MHz")
    # For a general zoom.
    addToOutput(output, 'zoom', 'correlator_channel_bandwidth', workArea['resolutions']['zoom'],
                "Correlator Channel Bandwidth", "MHz")
    if (specificZoomCalc):
        # For the specific zoom (which will be the same as the general zoom).
        addToOutput(output, 'specific_zoom', 'correlator_channel_bandwidth', workArea['resolutions']['zoom'],
                    "Correlator Channel Bandwidth", "MHz")
    #\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\



    ####################################################################################################
    # Make the continuum spectrum templates.
    print "MESSAGE: Generating template spectra..."
    workArea['continuum'] = makeTemplate(args.frequency, continuumBandwidth,
                                         workArea['resolutions']['continuum'])
    workArea['alternate'] = makeTemplate((args.frequency - workArea['resolutions']['continuum'] / 2),
                                         continuumBandwidth, workArea['resolutions']['continuum'])
    # Get the lowest and highest frequencies that we will need to read from the files.
    lowGlobalFreq = (workArea['alternate']['centreFrequency'][0] -
                     (workArea['resolutions']['continuum'] / 2.0))
    highGlobalFreq = (workArea['alternate']['centreFrequency'][-1] +
                      (workArea['resolutions']['continuum'] / 2.0))
    # Read the Tsys between those frequencies into memory if required.
    b = frequencyBand(args.frequency)
    try:
        cacheTsys = readTsys(frequencyBands[b]['tsys'], lowGlobalFreq, highGlobalFreq)
    except CalcError as c:
        print "FATAL: ", c.value
        sys.exit(-1)
    # Use the raw Tsys to fill in the required templates.
    templateFill(cacheTsys, workArea['continuum'])
    templateFill(cacheTsys, workArea['alternate'])
    # Make the specific zoom template here too.
    if (specificZoomCalc):
        # Determine which channel is closest in centre frequency to the nominated
        # specific zoom frequency.
        closestCentreFreq = 0.0
        closestCentreOffs = 3000.0
        for i in xrange(0, len(workArea['continuum']['centreFrequency'])):
            offs = abs(workArea['continuum']['centreFrequency'][i] - args.zoomfreq)
            if (offs < closestCentreOffs):
                closestCentreFreq = workArea['continuum']['centreFrequency'][i]
                closestCentreOffs = offs
        for i in xrange(0, len(workArea['alternate']['centreFrequency'])):
            offs = abs(workArea['alternate']['centreFrequency'][i] - args.zoomfreq)
            if (offs < closestCentreOffs):
                closestCentreFreq = workArea['alternate']['centreFrequency'][i]
                closestCentreOffs = offs
        # Compute how the bandwidth is distributed around the cente frequency.
        bandwidthAbove = (0.125 * (2.0 * float(args.zoom_width) + (-1.0) ** float(args.zoom_width) + 3.0) *
                          workArea['resolutions']['continuum'])
        bandwidthBelow = (0.125 * (2.0 * float(args.zoom_width) - (-1.0) ** float(args.zoom_width) + 1.0) *
                          workArea['resolutions']['continuum'])
        szBandwidths = [ bandwidthBelow, bandwidthAbove ]
        workArea['specificZoom'] = makeTemplate(closestCentreFreq, szBandwidths,
                                                workArea['resolutions']['zoom'])
        # Fill in the template from the continuum template.
        templateFill(workArea['continuum'], workArea['specificZoom'])
        
    # Put the frequency ranges of the templates in the output.
    # The continuum ranges.
    lowestFreq = (workArea['continuum']['centreFrequency'][args.edge] - 
                  (workArea['resolutions']['continuum'] / 2))
    highestFreq = (workArea['continuum']['centreFrequency'][-(args.edge + 1)] +
                   (workArea['resolutions']['continuum'] / 2))
    addToOutput(output, 'continuum', 'frequency_range', [ lowestFreq, highestFreq ],
                "Continuum Band Frequency Range", "MHz")
    if (specificZoomCalc):
        # Get the lowest and highest frequencies.
        szoomLowestFreq = (workArea['specificZoom']['centreFrequency'][args.zoom_edge] -
                           (workArea['resolutions']['zoom'] / 2))
        szoomHighestFreq = (workArea['specificZoom']['centreFrequency'][-(args.zoom_edge + 1)] +
                            (workArea['resolutions']['zoom'] / 2))
        addToOutput(output, 'parameters', 'zoom_frequency', args.zoomfreq, "Specific Zoom Frequency", "MHz")
        addToOutput(output, 'specific_zoom', 'frequency_range', [ szoomLowestFreq, szoomHighestFreq ],
                    "Specific Zoom Frequency Range", "MHz")
        # The number of zooms used.
        addToOutput(output, 'specific_zoom', 'n_zooms', args.zoom_width,
                    "Number of concatenated zoom channels", "channels")
    #\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\



    ####################################################################################################
    # Make antenna efficiency templates.
    # For the continuum bands.
    efficiencyMasterTemplate = templateEfficiency()
    workArea['continuum-efficiency'] = makeTemplate(args.frequency, continuumBandwidth,
                                                    workArea['resolutions']['continuum'])
    workArea['alternate-efficiency'] = makeTemplate((args.frequency - workArea['resolutions']['continuum'] / 2),
                                                    continuumBandwidth, workArea['resolutions']['continuum'])
    templateFill(efficiencyMasterTemplate, workArea['continuum-efficiency'])
    templateFill(efficiencyMasterTemplate, workArea['alternate-efficiency'])
    if (specificZoomCalc):
        # For the specific zoom band.
        workArea['specificZoom-efficiency'] = makeTemplate(closestCentreFreq, szBandwidths,
                                                           workArea['resolutions']['zoom'])
        templateFill(workArea['continuum-efficiency'], workArea['specificZoom-efficiency'])

    # Put the average computed efficiencies in the output.
    averageEfficiency = "%.1f" % (averageTemplate(workArea['continuum-efficiency']) * 100.0)
    addToOutput(output, 'parameters', 'antenna_efficiency', float(averageEfficiency),
                "Antenna Efficiency", "%")
    #\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\



    ####################################################################################################
    # Do template flagging.
    print "MESSAGE: Flagging..."
    flagTemplate(workArea['continuum'], 'continuum', args.corrconfig, 0)
    flagTemplate(workArea['continuum'], 'edge', args.corrconfig, args.edge)
    if (specificZoomCalc == True):
        flagTemplate(workArea['specificZoom'], 'edge', args.corrconfig, args.zoom_edge)
    if (args.birdies):
        flagTemplate(workArea['continuum'], 'birdies', args.corrconfig, 0)
    if (args.rfi):
        flagTemplate(workArea['continuum'], 'rfi', args.corrconfig, 0)
    #\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\



    ####################################################################################################
    # Do any bin smoothing now. We need to smooth the system temperature and efficiency templates.
    # The continuum frequency resolution after smoothing.
    contSmoothRes = workArea['resolutions']['continuum'] * float(args.smoothing)
    addToOutput(output, 'source_imaging', 'smoothing_window', args.smoothing,
                "Smoothing Window", "channels")
    # Check that we will end up with at least 2 channels, otherwise we die.
    if ((continuumBandwidth / contSmoothRes) < 2):
        print "FATAL: Smoothing factor too large."
        sys.exit(-1)
    # Make the smoothed templates and fill them from the unsmoothed templates.
    workArea['continuum-smooth'] = makeTemplate(args.frequency, continuumBandwidth, contSmoothRes)
    templateFill(workArea['continuum'], workArea['continuum-smooth'])
    workArea['continuum-efficiency-smooth'] = makeTemplate(args.frequency, continuumBandwidth, contSmoothRes)
    templateFill(workArea['continuum-efficiency'], workArea['continuum-efficiency-smooth'])
    addToOutput(output, 'continuum', 'channel_bandwidth', contSmoothRes,
                "Channel Bandwidth", "MHz")

    # The zooms are smoothed by a different argument.
    # The "general" zoom frequency resolution after smoothing.
    zoomSmoothRes = workArea['resolutions']['zoom'] * float(args.zoom_smoothing)
    addToOutput(output, 'source_imaging', 'zoom_smoothing_window', args.zoom_smoothing,
                "Zoom Smoothing Window", "channels")
    addToOutput(output, 'zoom', 'channel_bandwidth', zoomSmoothRes,
                "Channel Bandwidth", "MHz")
    if (specificZoomCalc):
        # The specific zoom will have the exact same frequency resolution.
        addToOutput(output, 'specific_zoom', 'channel_bandwidth', zoomSmoothRes,
                    "Channel Bandwidth", "MHz")
    # We don't want to smooth too much in the zooms either, so we check for that now.
    if ((workArea['resolutions']['continuum'] / zoomSmoothRes) < 2):
        print "FATAL: Zoom smoothing factor too large."
        sys.exit(-1)

    if (specificZoomCalc):
        # Make the specific zoom smoothed template.
        workArea['specificZoom-smooth'] = makeTemplate(closestCentreFreq, szBandwidths,
                                                       zoomSmoothRes)
        templateFill(workArea['specificZoom'], workArea['specificZoom-smooth'])
        workArea['specificZoom-efficiency-smooth'] = makeTemplate(closestCentreFreq, szBandwidths,
                                                                  zoomSmoothRes)
        templateFill(workArea['specificZoom-efficiency'], workArea['specificZoom-efficiency-smooth'])
    #\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
    


    ####################################################################################################
    # Calculate the properties of the primary beam.
    # Primary beam field of view, in arcminutes, for the continuum central frequency.
    pbfwhm = primaryBeamSize(args.frequency, antennaDiameter)
    addToOutput(output, 'source_imaging', 'field_of_view', pbfwhm,
                "Field of View (primary beam FWHM)", "arcmin")
    # And the low/high frequencies of the continuum band.
    pbfwhml = primaryBeamSize(workArea['continuum']['centreFrequency'][0], antennaDiameter)
    pbfwhmh = primaryBeamSize(workArea['continuum']['centreFrequency'][-1], antennaDiameter)
    addToOutput(output, 'source_imaging', 'field_of_view_range', [ pbfwhml, pbfwhmh ],
                "Field of View Range (primary beam FWHM)", "arcmin")
    if (specificZoomCalc):
        # And at the specific zoom frequency.
        zpbfwhm = primaryBeamSize(args.zoomfreq, antennaDiameter)
        addToOutput(output, 'source_imaging', 'field_of_view_zoom', zpbfwhm,
                    "Specific Zoom Band Field of View (primary beam FWHM)", "arcmin")
    #\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\



    ####################################################################################################
    # Calculate the properties of the synthesised beam.
    # Synthesised beam for the continuum central frequency.
    synthBeamContinuum = synthesisedBeamSize(args.frequency, maxBaselineLength, args.dec, hourAngle_min,
                                hourAngle_max, imageWeights['beam'])
    addToOutput(output, 'source_imaging', 'synthesised_beam_size', synthBeamContinuum,
                "Synthesised Beam Size (FWHM)", "arcsec")
    # The sizes at the high and low frequencies of the continuum band.
    synthBeamLowFreq = synthesisedBeamSize(lowestFreq, maxBaselineLength, args.dec, hourAngle_min,
                                           hourAngle_max, imageWeights['beam'])
    synthBeamHighFreq = synthesisedBeamSize(highestFreq, maxBaselineLength, args.dec, hourAngle_min,
                                            hourAngle_max, imageWeights['beam'])
    addToOutput(output, 'source_imaging', 'synthesised_beam_size_range', 
                [ synthBeamLowFreq, synthBeamHighFreq ],
                "Synthesised Beam Size Range (FWHM)", "arcsec")
    if (specificZoomCalc):
        # The synthesised beam for the frequency the user wanted in the specific zoom band.
        synthBeamZoom = synthesisedBeamSize(args.zoomfreq, maxBaselineLength, args.dec, hourAngle_min,
                                            hourAngle_max, imageWeights['beam'])
        addToOutput(output, 'source_imaging', 'synthesised_beam_size_zoom', synthBeamZoom,
                    "Specific Zoom Band Synthesised Beam Size (FWHM)", "arcsec")
    #\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\



    ####################################################################################################
    # Work out the velocity resolutions and ranges. Because the channels flagged at the
    # edge shouldn't be considered as usable bandwidth, we don't consider them as usable
    # for velocity range.
    # The usable frequency bandwidth over the continuum band.
    contVelBandwidth = continuumBandwidth - 2.0 * float(args.edge) * workArea['resolutions']['continuum']
    # We calculate the velocity span corresponding to that amount of bandwidth with the
    # actual lower frequency.
    addToOutput(output, 'continuum', 'spectral_bandwidth',
                bandwidthToVelocity(workArea['continuum']['centreFrequency'][args.edge],
                                    contVelBandwidth, argsInterpreted['restfreq']),
                "Spectral Bandwidth", "km/s")
    # The velocity width of a channel is dependant on its frequency, so we can't give
    # an exact figure for every channel here; we just divide the continuum range by the
    # number of channels here, excluding the edge channels.
    contRes = (output['continuum']['spectral_bandwidth'] / 
               float(len(workArea['continuum']['centreFrequency']) - 2 * args.edge))
    # We now compensate for continuum band smoothing.
    chanRes = float("%.3f" % (contRes / float(args.smoothing)))
    addToOutput(output, 'continuum', 'spectral_channel_resolution', chanRes,
                "Spectral Channel Resolution", "km/s")

    # The velocity width of a "general" zoom band needs to be calculated in the same way
    # since the user may not want the edge channels there.
    # The usable frequency bandwidth over a "general" zoom band.
    zoomVelBandwidth = (workArea['resolutions']['continuum'] - 
                        (2.0 * float(args.zoom_edge) * workArea['resolutions']['zoom']))
    addToOutput(output, 'zoom', 'effective_bandwidth', zoomVelBandwidth, 
                "Effective Bandwidth", "MHz")
    # But we can't actually calculate a "general" velocity width that way since it is again
    # frequency dependant. So we just divide up the continuum velocity resolution.
    zoomRes = contRes / float(nZoomChannels)
    # And compensate for zoom band smoothing.
    zoomChanRes = float("%.3f" % (zoomRes / float(args.zoom_smoothing)))
    addToOutput(output, 'zoom', 'spectral_channel_resolution', zoomChanRes, "Spectral Channel Resolution",
                "km/s")
    # The number of channels in the zoom, excluding the edge flagged channels.
    nZoomVelChans = nZoomChannels - 2 * args.zoom_edge
    # But we output only the number of smoothed channels.
    addToOutput(output, 'zoom', 'n_channels', (nZoomVelChans / args.zoom_smoothing),
                "# Channels", None)
    # The zoom velocity width is then just the velocity resolution of each channel multiplied
    # by the number of non-edge channels.
    zoomSpecBandwidth = zoomRes * float(nZoomVelChans)
    addToOutput(output, 'zoom', 'spectral_bandwidth', zoomSpecBandwidth, "Spectral Bandwidth", "km/s")

    if (specificZoomCalc == True):
        # We can actually calculate real velocity values for the specific zoom band, so we
        # do it basically the same as for the continuum band.
        # The frequency bandwidth of the specific zoom takes into account how many zoom
        # channels that the user will consolidate, and the number of channels flagged at the
        # edge.
        szoomVelBandwidth = (workArea['resolutions']['continuum'] * sum(szBandwidths) - 
                             (2.0 * float(args.zoom_edge) * workArea['resolutions']['zoom']))
        addToOutput(output, 'specific_zoom', 'effective_bandwidth', szoomVelBandwidth, 
                    "Effective Bandwidth", "MHz")
        # The velocity span corresponding to that amount of bandwidth is calculated with
        # respect to the lower frequency of the zoom.
        addToOutput(output, 'specific_zoom', 'spectral_bandwidth',
                    bandwidthToVelocity(workArea['specificZoom']['centreFrequency'][args.zoom_edge],
                                        szoomVelBandwidth, argsInterpreted['restfreq']),
                    "Spectral Bandwidth", "km/s")
        # We calculate the number of channels across the consolidated zoom.
        nszoomChans = int(sum(szBandwidths) * float(nZoomChannels) - float(2 * args.zoom_edge))
        # But we output only the number of smoothed channels.
        addToOutput(output, 'specific_zoom', 'n_channels', (nszoomChans / args.zoom_smoothing),
                    "# Channels", None)

        # The velocity resolution is just the velocity span divided by the number of channels,
        # excluding those flagged at the edge.
        szoomRes = (output['specific_zoom']['spectral_bandwidth'] / float(nszoomChans))
        # We compensate for zoom smoothing.
        szoomChanRes = float("%.3f" % (szoomRes / float(args.zoom_smoothing)))
        addToOutput(output, 'specific_zoom', 'spectral_channel_resolution', szoomChanRes,
                    "Spectral Channel Resolution", "km/s")
    #\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\



    ####################################################################################################
    # Get the atmospheric parameters for the zenith and store that in a template.
    print "MESSAGE: Calculating weather effects..."

    # The weather conditions that we will use for computing the atmosphere later.
    weatherConditions = {
        'JAN': { 'best': { 'temperature': 32.7, 'pressure': 986.8, 'humidity': 27.5 },
                 'typical': { 'temperature': 30.6, 'pressure': 989.8, 'humidity': 50.5 },
                 'worst': { 'temperature': 26.3, 'pressure': 1001.9, 'humidity': 91.0 } }, 
        'FEB': { 'best': { 'temperature': 29.6, 'pressure': 987.2, 'humidity': 36.0 },
                 'typical': { 'temperature': 24.2, 'pressure': 989.6, 'humidity': 65.0 },
                 'worst': { 'temperature': 23.7, 'pressure': 990.4, 'humidity': 86.0 } },
        'MAR': { 'best': { 'temperature': 27.4, 'pressure': 989.7, 'humidity': 33.0 },
                 'typical': { 'temperature': 18.9, 'pressure': 995.5, 'humidity': 82.0 },
                 'worst': { 'temperature': 28.2, 'pressure': 988.0, 'humidity': 69.7 } },
        'APR': { 'best': { 'temperature': 9.5, 'pressure': 1011.7, 'humidity': 76.0 },
                 'typical': { 'temperature': 16.8, 'pressure': 1013.4, 'humidity': 73.0 },
                 'worst': { 'temperature': 19.7, 'pressure': 1001.6, 'humidity': 85.0 } },
        'MAY': { 'best': { 'temperature': 19.6, 'pressure': 1008.0, 'humidity': 31.0 },
                 'typical': { 'temperature': 9.7, 'pressure': 1009.7, 'humidity': 86.0 },
                 'worst': { 'temperature': 18.7, 'pressure': 1012.3, 'humidity': 78.0 } },
        'JUN': { 'best': { 'temperature': -1.6, 'pressure': 1016.9, 'humidity': 95.0 },
                 'typical': { 'temperature': 8.1, 'pressure': 1002.7, 'humidity': 95.0 },
                 'worst': { 'temperature': 15.0, 'pressure': 997.7, 'humidity': 101.1 } },
        'JUL': { 'best': { 'temperature': 1.9, 'pressure': 1019.0, 'humidity': 91.0 },
                 'typical': { 'temperature': 18.8, 'pressure': 999.4, 'humidity': 50.5 },
                 'worst': { 'temperature': 15.6, 'pressure': 1004.2, 'humidity': 100.0 } },
        'AUG': { 'best': { 'temperature': 3.6, 'pressure': 1017.2, 'humidity': 73.0 },
                 'typical': { 'temperature': 8.2, 'pressure': 1010.5, 'humidity': 87.0 },
                 'worst': { 'temperature': 16.6, 'pressure': 1012.4, 'humidity': 93.0 } },
        'SEP': { 'best': { 'temperature': 19.9, 'pressure': 989.9, 'humidity': 27.0 },
                 'typical': { 'temperature': 15.4, 'pressure': 993.5, 'humidity': 61.0 },
                 'worst': { 'temperature': 20.3, 'pressure': 993.6, 'humidity': 66.0 } },
        'OCT': { 'best': { 'temperature': 26.6, 'pressure': 986.6, 'humidity': 22.0 },
                 'typical': { 'temperature': 28.6, 'pressure': 986.5, 'humidity': 33.0 },
                 'worst': { 'temperature': 25.8, 'pressure': 996.6, 'humidity': 57.0 } },
        'NOV': { 'best': { 'temperature': 32.5, 'pressure': 986.7, 'humidity': 19.7 },
                 'typical': { 'temperature': 22.1, 'pressure': 990.2, 'humidity': 58.0 },
                 'worst': { 'temperature': 24.7, 'pressure': 989.5, 'humidity': 71.9 } },
        'DEC': { 'best': { 'temperature': 29.5, 'pressure': 986.8, 'humidity': 30.0 },
                 'typical': { 'temperature': 21.8, 'pressure': 987.0, 'humidity': 70.0 },
                 'worst': { 'temperature': 27.9, 'pressure': 984.5, 'humidity': 71.0 } },
        'SUMMER': { 'best': { 'temperature': 29.5, 'pressure': 986.8, 'humidity': 30.0 },
                    'typical': { 'temperature': 29.9, 'pressure': 989.1, 'humidity': 49.0 },
                    'worst': { 'temperature': 26.3, 'pressure': 1001.9, 'humidity': 91.0 } },
        'AUTUMN': { 'best': { 'temperature': 19.5, 'pressure': 1008.0, 'humidity': 31.0 },
                    'typical': { 'temperature': 23.6, 'pressure': 1016.5, 'humidity': 51.0 },
                    'worst': { 'temperature': 28.2, 'pressure': 988.0, 'humidity': 69.7 } },
        'WINTER': { 'best': { 'temperature': -1.6, 'pressure': 1016.9, 'humidity': 95.0 },
                    'typical': { 'temperature': 9.7, 'pressure': 1004.6, 'humidity': 83.0 },
                    'worst': { 'temperature': 15.6, 'pressure': 1004.2, 'humidity': 100.0 } },
        'SPRING': { 'best': { 'temperature': 19.9, 'pressure': 989.9, 'humidity': 27.0 },
                    'typical': { 'temperature': 14.9, 'pressure': 997.5, 'humidity': 68.7 },
                    'worst': { 'temperature': 24.7, 'pressure': 989.5, 'humidity': 71.9 } },
        'APRS': { 'best': { 'temperature': 19.9, 'pressure': 989.9, 'humidity': 27.0 },
                  'typical': { 'temperature': 12.7, 'pressure': 1011.1, 'humidity': 72.0 },
                  'worst': { 'temperature': 19.7, 'pressure': 1001.6, 'humidity': 85.0 } },
        'OCTS': { 'best': { 'temperature': 26.6, 'pressure': 986.6, 'humidity': 22.0 },
                  'typical': { 'temperature': 22.1, 'pressure': 988.7, 'humidity': 66.7 },
                  'worst': { 'temperature': 26.3, 'pressure': 1001.9, 'humidity': 91.0 } },
        'ANNUAL': { 'best': { 'temperature': 19.9, 'pressure': 989.9, 'humidity': 27.0 },
                    'typical': { 'temperature': 16.6, 'pressure': 1010.3, 'humidity': 57.5 },
                    'worst': { 'temperature': 26.3, 'pressure': 1001.9, 'humidity': 91.0 } }
        }

    # This is the frequency resolution of the atmospheric corrections.
    atmosRes = max(workArea['resolutions']['continuum'], args.per_freq)
    addToOutput(output, 'parameters', 'atmosphere_frequency_resolution', atmosRes,
                "Frequency resolution of atmospheric parameters", "MHz")
    # Include the weather parameters we use in the output.
    addToOutput(output, 'parameters', 'atmospheric_conditions', weatherConditions[args.season],
                "Atmospheric conditions", None)
    
    # Form the opacity and temperature templates for each of the weather conditions
    # (best, typical, worst).
    workArea['opacity'] = {}
    workArea['temperature'] = {}
    workArea['sz-opacity'] = {}
    workArea['sz-temperature'] = {}
    for condition in weatherConditions[args.season]:
        # Make the continuum templates first.
        tempOpacity = makeTemplate(args.frequency, continuumBandwidth, atmosRes)
        tempTemperature = makeTemplate(args.frequency, continuumBandwidth, atmosRes)
        workArea['opacity'][condition] = makeTemplate(args.frequency, continuumBandwidth,
                                                      workArea['resolutions']['continuum'])
        workArea['temperature'][condition] = makeTemplate(args.frequency, continuumBandwidth,
                                                          workArea['resolutions']['continuum'])
        fillAtmosphereTemplate(tempOpacity, tempTemperature, 
                               (weatherConditions[args.season][condition]['temperature'] + 273.15),
                               (weatherConditions[args.season][condition]['pressure'] * 100.0),
                               (weatherConditions[args.season][condition]['humidity'] / 100.0))
        templateFill(tempOpacity, workArea['opacity'][condition])
        templateFill(tempTemperature, workArea['temperature'][condition])
        if (specificZoomCalc):
            # Make the specific zoom templates if we need to, and fill them from the already
            # computed templates for the continuum band.
            workArea['sz-opacity'][condition] = makeTemplate(closestCentreFreq, szBandwidths,
                                                             workArea['resolutions']['zoom'])
            workArea['sz-temperature'][condition] = makeTemplate(closestCentreFreq, szBandwidths,
                                                                 workArea['resolutions']['zoom'])
            templateFill(workArea['opacity'][condition], workArea['sz-opacity'][condition])
            templateFill(workArea['temperature'][condition], workArea['sz-temperature'][condition])
    #\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\



    ####################################################################################################
    # Compute the sensitivites with all the information we just collected.
    print "MESSAGE: Calculating sensitivities..."
    workArea['continuum-rms'] = {}
    workArea['continuum-smooth-rms'] = {}
    workArea['specificZoom-rms'] = {}
    computeBands = { 'continuum': 0, 'spectral': 0, 'zoom': 0 }
    workArea['rms'] = {
        'continuum': {},
        'spectral': {},
        'zoom': {},
        'specificZoom': {}
    }
    workArea['btrms'] = {
        'continuum': {},
        'spectral': {},
        'zoom': {},
        'specificZoom': {}
    }
    workArea['SEFD'] = {}

    # Sensitivities are computed for each of the different weather conditions we expect.
    sensitivityReached = False
    while (sensitivityReached == False):
        for condition in weatherConditions[args.season]:
            # The RMS noise in the smoothed continuum band, for each channel.
            workArea['continuum-smooth-rms'][condition] = calculateRms(workArea['continuum-smooth'],
                                                                       workArea['continuum-efficiency-smooth'],
                                                                       workArea['opacity'][condition],
                                                                       workArea['temperature'][condition],
                                                                       hourAngle_min, hourAngle_max, args.per_ha,
                                                                       nant, args.integration, imageWeights, sind, cosd)
            # Then derive the global average values in the continuum band.
            sensResSmooth = calculateSensitivity(workArea['continuum-smooth-rms'][condition], nant)
            # Check whether we have any unflagged continuum channels.
            if (sensResSmooth['bandwidth']['unflagged'] < 1.0):
                print "FATAL: No continuum bandwidth remains unflagged."
                sys.exit(-1)

            # We get the "general" zoom sensitivity from the unsmoothed continuum data, since smoothing
            # the continuum won't help improve the zoom sensitivity.
            workArea['continuum-rms'][condition] = calculateRms(workArea['continuum'],
                                                                workArea['continuum-efficiency'],
                                                                workArea['opacity'][condition],
                                                                workArea['temperature'][condition],
                                                                hourAngle_min, hourAngle_max, args.per_ha,
                                                                nant, args.integration, imageWeights, sind, cosd)
            # Then derive the global average values in the "general" zoom band.
            sensRes = calculateSensitivity(workArea['continuum-rms'][condition], nant)

            if (specificZoomCalc):
                # The RMS noise in the smoothed specific zoom band, for each channel.
                workArea['specificZoom-rms'][condition] = calculateRms(workArea['specificZoom-smooth'],
                                                                       workArea['specificZoom-efficiency-smooth'],
                                                                       workArea['sz-opacity'][condition],
                                                                       workArea['sz-temperature'][condition],
                                                                       hourAngle_min, hourAngle_max, args.per_ha,
                                                                       nant, args.integration, imageWeights, sind, cosd)
                # Then derive the global average values in the specific zoom band.
                szSensRes = calculateSensitivity(workArea['specificZoom-rms'][condition], nant)

            for t in computeBands:
                if (t == 'zoom'):
                    # Since we don't make a template for the "general" zoom, we have to manually
                    # adjust for the zoom smoothing factor now.
                    # The RMS noise goes down with the square-root of the smoothed bandwidth over the
                    # normal bandwidth.
                    sensRes['rms'][t] /= math.sqrt(float(args.zoom_smoothing))
                    workArea['rms'][t][condition] = float("%.3f" % sensRes['rms'][t])
                    # Calculate the brightness temperature sensitivity using the synthesised beam at
                    # the centre of the continuum band.
                    bts = brightnessTemperatureSensitivity(sensRes['rms'][t], synthBeamContinuum,
                                                           args.frequency)
                else:
                    # We calculate the continuum and continuum-spectral sensitivities in the
                    # same way.
                    workArea['rms'][t][condition] = float("%.3f" % sensResSmooth['rms'][t])
                    # Calculate the brightness temperature sensitivity using the synthesised beam at
                    # the centre of the continuum band.
                    bts = brightnessTemperatureSensitivity(sensResSmooth['rms'][t], synthBeamContinuum,
                                                           args.frequency)
                # We put the brightness sensitivity in mK for readability.
                workArea['btrms'][t][condition] = float("%.2f" % (bts * 1000.0))
                # We keep the SEFDs in Jy.
                workArea['SEFD'][condition] = {
                    'antenna': float( "%.1f" % (sensRes['sefd']['antenna']) ),
                    'array': float( "%.1f" % (sensRes['sefd']['array'])) }

            if (specificZoomCalc):
                # The spectral RMS of the specific zoom band.
                workArea['rms']['specificZoom'][condition] = float("%.3f" % szSensRes['rms']['spectral'])
                # Calculate the brightness temperature sensitivity using the synthesised beam at the
                # nominated specific zoom band frequency.
                bts = brightnessTemperatureSensitivity(szSensRes['rms']['spectral'], synthBeamZoom,
                                                       args.zoomfreq)
                # The brightness sensitivity is again in mK.
                workArea['btrms']['specificZoom'][condition] = float("%.2f" % (bts * 1000.0))

        # Check if we need to adjust the integration time.
        if (args.calculate_time == False):
            # We operate from time to sensitivity, so we exit now.
            sensitivityReached = True
        else:
            # We compare the sensitivity we obtained with the target.
            if (args.target_best):
                conditionTarget = 'best'
            elif (args.target_typical):
                conditionTarget = 'typical'
            elif (args.target_worst):
                conditionTarget = 'worst'

            if (args.target_continuum):
                bandTarget = 'continuum'
            elif (args.target_spectral):
                bandTarget = 'spectral'
            elif (args.target_zoom):
                bandTarget = 'zoom'
            elif (args.target_specific_zoom):
                bandTarget = 'specificZoom'

            if (args.target_flux_density):
                modeTarget = 'rms'
            elif (args.target_brightness_temperature):
                modeTarget = 'btrms'

            compareSensitivity = workArea[modeTarget][bandTarget][conditionTarget]
            # Calculate the ratio of the obtained sensitivity to that desired.
            sensRatio = compareSensitivity / args.target
            # We stop if we're within 1% of the target sensitivity.
            if (abs(sensRatio - 1.0) < 0.01):
                sensitivityReached = True
            else:
                # Change the integration time appropriately.
                args.integration *= sensRatio * sensRatio
            

    # We now stick all this information into the output.
    # The effective bandwidth of the continuum band depends on what the user chose to flag.
    addToOutput(output, 'continuum', 'effective_bandwidth', sensResSmooth['bandwidth']['unflagged'],
                "Effective Bandwidth", "MHz")
    # The effective number of channels in the continuum band depends on the flagging and the
    # smoothing factor.
    addToOutput(output, 'continuum', 'n_channels',
                (int(sensResSmooth['bandwidth']['unflagged'] / contSmoothRes)),
                "# Channels", None)
    # The computed system temperatures over the continuum band, for each weather condition, in K.
    addToOutput(output, 'sensitivities', 'system_temperature',
                [ workArea['continuum-rms']['best']['systemTemperature'],
                  workArea['continuum-rms']['typical']['systemTemperature'],
                  workArea['continuum-rms']['worst']['systemTemperature'] ],
                "System Temperature", "K")

    # The continuum sensitivities, for each weather condition, in mJy/beam.
    addToOutput(output, 'sensitivities', [ 'rms_noise_level', 'continuum' ],
                [ workArea['rms']['continuum']['best'], workArea['rms']['continuum']['typical'],
                  workArea['rms']['continuum']['worst'] ],
                "RMS noise level", "mJy/beam")
    # The continuum brightness temperature sensitivity, for each weather condition, in mK.
    addToOutput(output, 'sensitivities', [ 'brightness_temperature_sensitivity', 'continuum' ],
                [ workArea['btrms']['continuum']['best'], workArea['btrms']['continuum']['typical'],
                  workArea['btrms']['continuum']['worst'] ],
                "Brightness Temperature Sensitivity", "mK")

    # The RMS spectral noise in the continuum band, for each weather condition, in mJy/beam.
    addToOutput(output, 'sensitivities', [ 'rms_noise_level', 'spectral' ],
                [ workArea['rms']['spectral']['best'], workArea['rms']['spectral']['typical'],
                  workArea['rms']['spectral']['worst'] ],
                "RMS noise level", "mJy/beam")
    # The RMS spectral brightness noise in the continuum band, for each weather condition, in mK.
    addToOutput(output, 'sensitivities', [ 'brightness_temperature_sensitivity', 'spectral' ],
                [ workArea['btrms']['spectral']['best'], workArea['btrms']['spectral']['typical'],
                  workArea['btrms']['spectral']['worst'] ],
                "Brightness Temperature Sensitivity", "mK")

    # The RMS spectral noise in a "general" zoom band, for each weather condition, in mJy/beam.
    addToOutput(output, 'sensitivities', [ 'rms_noise_level', 'zoom' ],
                [ workArea['rms']['zoom']['best'], workArea['rms']['zoom']['typical'],
                  workArea['rms']['zoom']['worst'] ],
                "RMS noise level", "mJy/beam")
    # The RMS spectral brightness noise in a "general" zoom band, for each weather condition, in mK.
    addToOutput(output, 'sensitivities', [ 'brightness_temperature_sensitivity', 'zoom' ],
                [ workArea['btrms']['zoom']['best'], workArea['btrms']['zoom']['typical'],
                  workArea['btrms']['zoom']['worst'] ],
                "Brightness Temperature Sensitivity", "mK")

    # The SEFD of a single antenna, for each weather condition, in Jy.
    addToOutput(output, 'sensitivities', 'antenna_sensitivity',
                [ workArea['SEFD']['best']['antenna'], workArea['SEFD']['typical']['antenna'],
                  workArea['SEFD']['worst']['antenna']],
                "Antenna SEFD", "Jy")
    # The SEFD of the entire array combined, for each weather condition, in Jy.
    addToOutput(output, 'sensitivities', 'array_sensitivity',
                [ workArea['SEFD']['best']['array'], workArea['SEFD']['typical']['array'],
                  workArea['SEFD']['worst']['array']],
                "Array SEFD", "Jy")

    if (specificZoomCalc):
        # The computed system temperatures over the specific zoom band, for each weather condition, in K.
        addToOutput(output, 'sensitivities', 'specific_zoom_system_temperature',
                    [ workArea['specificZoom-rms']['best']['systemTemperature'],
                      workArea['specificZoom-rms']['typical']['systemTemperature'],
                      workArea['specificZoom-rms']['worst']['systemTemperature'] ],
                    "System Temperature", "K")
        
        # The RMS spectral noise in the specific zoom band, for each weather condition, in mJy/beam.
        addToOutput(output, 'sensitivities', [ 'rms_noise_level', 'specific_zoom' ],
                    [ workArea['rms']['specificZoom']['best'], workArea['rms']['specificZoom']['typical'],
                      workArea['rms']['specificZoom']['worst'] ],
                    "RMS noise level", "mJy/beam")
        # The RMS spectral brightness noise in the specific zoom band, for each weather condition, in mK.
        addToOutput(output, 'sensitivities', [ 'brightness_temperature_sensitivity', 'specific_zoom' ],
                    [ workArea['btrms']['specificZoom']['best'], workArea['btrms']['specificZoom']['typical'],
                      workArea['btrms']['specificZoom']['worst'] ],
                    "Brightness Temperature Sensitivity", "mK")
    # The integration time is added now, since it may have changed if the user asked for a
    # particular sensitivity target.
    addToOutput(output, 'source_imaging', 'integration_time', args.integration,
                "Time on Source", "minutes")

    #\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\



    ####################################################################################################
    # Plot the RMS spectral noise templates into output files.
    # The name of the output file; we make sure it doesn't end with ".png.png" first.
    outfile = args.output.replace('.png', '') + '.png'
    plotSpectrum(workArea['continuum-smooth-rms'], [ 'typical', 'best', 'worst' ], outfile)
    # Put the name of the output file in the output dictionary.
    output['output_plot'] = outfile

    if (specificZoomCalc):
        # We make an output plot of the RMS spectral noise of the specific zoom as well.
        szoutfile = args.output.replace('.png', '') + '.sz.png'
        plotSpectrum(workArea['specificZoom-rms'], [ 'typical', 'best', 'worst'], szoutfile)
        # Put the name of this output file in the output dictionary.
        output['output_zoom_plot'] = szoutfile
    #\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\



    ####################################################################################################
    # Output all the quantities we have computed and stored.
    if (args.human_readable):
        # Output the object in a human readable format.
        humanOutputDict(output, output['description'], output['units'], 0)
    else:
        # Output the JSON.
        print json.dumps(output)
    #\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
    # SENSITIVITY CALCULATOR ENDS

if __name__ == "__main__":
    # Set up the arguments.
    for i, arg in enumerate(sys.argv):
        # Deal with the possibly negative number that is declination.
        if (arg[0] == '-') and arg[1].isdigit():
            sys.argv[i] = ' ' + arg
            
    parser = argparse.ArgumentParser() 
    parser.add_argument("--ca06", action="store_true",
                        help="include CA06 for the calculations")
    parser.add_argument("-a", "--halimit", type=float, default=6,
                        help="the hour angle limit to use (decimal hours)")
    parser.add_argument("-b", "--corrconfig", default="CFB1M",
                        help="the correlator configuration",
                        choices=[ "CFB1M", "CFB64M" ])
    parser.add_argument("-B", "--birdies", action="store_true",
                        help="flag the birdies in the continuum band")
    parser.add_argument("-c", "--configuration",
                        help="the array configuration", default="6000",
                        choices=[ "6000", "6km", "3000", "3km", "1500", "1.5km",
                                  "750", "750m", "367", "EW352", "EW367", "EW352/367",
                                  "h214", "H214", "h168", "H168", "h75", "H75" ]);
    parser.add_argument("-C", "--calculate-time", action="store_true",
                        help="the calculator will determine the time required to reach the target sensitivity")
    parser.add_argument("-d", "--dec", type=float, default=-30,
                        help="the declination of the source (decimal degrees)")
    parser.add_argument("-e", "--ellimit", type=float, default=12,
                        help="the elevation limit to use (decimal degrees)")
    parser.add_argument("-E", "--edge", type=int, default=0,
                        help="the number of edge channels to flag")
    parser.add_argument("-f", "--frequency", type=int,
                        help="the central frequency of the observations (MHz)")
    parser.add_argument("-F", "--per-freq", type=float, default=50.0,
                        help="the minimum frequency spacing between atmospheric corrections (MHz)")
    parser.add_argument("-H", "--ha-min", type=float,
                        help="the lowest hour angle observed (decimal hours)")
    parser.add_argument("-K", "--ha-max", type=float,
                        help="the highest hour angle observed (decimal hours)")
    parser.add_argument("-M", "--ha-middle", type=float, default=0.0,
                        help="the middle hour angle observed (decimal hours)")
    parser.add_argument("-o", "--output", default="spectrum",
                        help="the name of the PNG file to output the RMS v frequency spectrum to (no extension)")
    parser.add_argument("-p", "--human-readable", action="store_true",
                        help="make human readable output instead of the default JSON")
    parser.add_argument("-P", "--per-ha", type=float, default=5.0,
                        help="the number of calculations of atmosphere made per hour")
    parser.add_argument("-r", "--restfreq", type=float,
                        help="use this rest frequency to calculate the velocity parameters (GHz)")
    parser.add_argument("-R", "--rfi", action="store_true",
                        help="flag known RFI-affected regions of the continuum spectrum")
    parser.add_argument("-s", "--smoothing", type=int, default=1,
                        help="the number of continuum spectral channels to bin together in the output")
    parser.add_argument("-S", "--season", default="ANNUAL",
                        choices=[ "JAN", "FEB", "MAR", "APR", "MAY", "JUN",
                                  "JUL", "AUG", "SEP", "OCT", "NOV", "DEC",
                                  "SUMMER", "AUTUMN", "WINTER", "SPRING",
                                  "APRS", "OCTS", "ANNUAL" ],
                        help="the conditions to assume for weather dependence")
    parser.add_argument("-t", "--integration", type=float, default=720,
                        help="the amount of on-source integration time (min)")
    parser.add_argument("-T", "--target", type=float, default=0.0,
                        help="the target sensitivity to reach")
    parser.add_argument("--target-continuum", action="store_true",
                        help="the target sensitivity is for the continuum image")
    parser.add_argument("--target-spectral", action="store_true",
                        help="the target sensitivity is for the spectral RMS in the continuum band")
    parser.add_argument("--target-zoom", action="store_true",
                        help="the target sensitivity is for the spectral RMS in a typical zoom band")
    parser.add_argument("--target-specific-zoom", action="store_true",
                        help="the target sensitivity is for the spectral RMS in the specific zoom band")
    parser.add_argument("--target-flux-density", action="store_true",
                        help="the target sensitivity is specified as a flux density in mJy")
    parser.add_argument("--target-brightness-temperature", action="store_true",
                        help="the target sensitivity is specified as a brightness temperature in mK")
    parser.add_argument("--target-best", action="store_true",
                        help="the target sensitivity is to be obtained in the best weather conditions")
    parser.add_argument("--target-typical", action="store_true",
                        help="the target sensitivity is to be obtained in typical weather conditions")
    parser.add_argument("--target-worst", action="store_true",
                        help="the target sensitivity is to be obtained in the worst weather conditions")
    parser.add_argument("-w", "--weighting", default="R2",
                        help="the image weighting scheme",
                        choices=[ "R2", "R1", "R0", "R-1", "R-2" ])
    parser.add_argument("-W", "--zoom-width", default=1, type=int,
                        help="the number of concatenated zoom bands to use for the specific zoom")
    parser.add_argument("-y", "--zoom-smoothing", type=int, default=1,
                        help="the number of zoom spectral channels to bin together in the output")
    parser.add_argument("-z", "--zoomfreq", type=int,
                        help="a specific zoom frequency to calculate parameters for")
    parser.add_argument("-Z", "--zoom-edge", type=int, default=0,
                        help="the number of edge channels to flag in the zoom band")
    args = parser.parse_args()
    main(args)
