#!/usr/bin/env python

import matplotlib.pyplot as plt
import numpy as np
import os, sys

# This code is heavily drawn from the manual, "Passive Repeater Engineering,
# Manual 161A," 1984, Microflect Co, Inc., Salem, Oregon.  On page 49, an HP
# BASIC program is presented to calculate the antenna pattern for a passive
# repeater with rectangular aperture.  Output is an ascii plot of values and
# asterisks to generate a graphical plot.
#
# The code here references pages in the manual and the subroutine pattern()
# is take directly and changed only as needed to use python syntax and also
# store results in an array rather than print them.  The array values are
# then plotted using python's matplotlib library.
#
# Subroutines are in alphabetical order except that main() is at the
# very end.  Variable names are measurement_units.  That is, a descriptive
# name with an underscore separating it from a units suffix.  For example,
# height_m is some height in meters.
#
# Mike Markowski, mike.ab3ap@gmail.com
# Nov 2023

def cli(argv):
    '''Retrieve command line args provided by user.  Those not provided are
    prompted for.

    Inputs:
      argv[] (string): command as typed by user: progName -arg1 opt1...
    Outputs:
      f_GHz (float), frequency in GHz,
      v_ft (float), vertical length (height) in feet of repeater,
      h_ft (float), horizontal length (width) in feet of repeater,
      cut (string), 'Horizontal' or 'Vertical'.
      alpha_deg (float), half angle in degrees made by repeater, tx and rx.
    '''

    prog = argv[0]
    argv = argv[1:]
    cut = ''
    alpha_deg = f_GHz = h_ft = v_ft = -1

    # Retrieve command line args provided by user.
    i = 0
    while i < len(argv):
        match argv[i]:
            case '-a': # Included angle.
                i += 1
                alpha_deg = float(argv[i])/2
            case '-c': # Cut, v or h.
                i += 1
                cut = argv[i][0].upper()
            case '-f': # Freq in GHz.
                i += 1
                fStr = argv[i]
                f_GHz = float(fStr)
            case '-h': # Horizontal length (ft) of repeater.
                i += 1
                h_ft = float(argv[i])
            case '-v': # Vertical length (ft) of repeater.
                i += 1
                v_ft = float(argv[i])
            case _: # Error.
                usage(prog)
        i += 1

    # If required values not on cmd line, prompt user for them.
    if f_GHz == -1:
        fStr = input('Frequency in GHz: ')
        f_GHz = float(fStr)
    if v_ft == -1:
        v_ft = int(input('Passive vertical length in feet: '))
    if h_ft == -1:
        h_ft = int(input('Passive horizontal length in feet: '))
    if cut == '':
        cut = (input('\'H\'orizontal or \'V\'ertical cut: ')).upper()
    cut = ('Horizontal' if cut[0]=='H' else 'Vertical')
    if alpha_deg == -1:
        alpha_deg = float(input('%s included angle (deg) ' % cut))/2

    return f_GHz, fStr, v_ft, h_ft, cut, alpha_deg

def farField(aE_ft2, lambda_ft, dist_mi):
    '''Determine if repeater is in near or far field of station.
    Calculations are described on pp. 39 and 42.

    Inputs:
      aE_ft2 (float): effective area in square feet of repeater.
      lambda_ft (float): wavelength in feet if carrier.
      dist_mi (float): distance in miles between repeater and station.
    Output:
      (boolean): True if far-field, False if near-field.
    '''
    dist_ft = 5280*dist_mi
    return np.pi*lambda_ft*dist_ft/(4*aE_ft2) >= 2.5

def nearFieldRange_mi(aE_ft2, lambda_ft):
    # Formula from "Passive Repeater Engineering," on p. 39 graph.
    dist_mi = 10*aE_ft2/(np.pi*lambda_ft)/5280
    return dist_mi # Same as:  aE_ft2 * f_GHz / 1632

def nearFieldRangeCollins_mi(aE_ft2, f_GHz):
    '''This is the formula given on the Pickett Collins C18 microwave
    slide rule.
    '''
    return aE_ft2*f_GHz/2598

def passiveGain_dBi(area_ft2, lambda_ft):
    '''Calculate passive repeater gain using equation on p. 35.

    Inputs:
      area_ft2 (float): effective area in square feet of repeater.
      lambda_ft (float): wavelength in feet if carrier.
    Output:
      (float): passive gain in units of dBi.
    '''
    return 20*np.log10(4*np.pi*area_ft2/lambda_ft**2) # Passive repeater gain.

def pattern(Gp_dBi, a_ft, lambda_ft):
    '''Create far-field radiation pattern for a rectangular aperture passive
    repeater.  Converted from p.49 BASIC program.

    Inputs:
      Gp_dBi (float): max gain in dBi of repeater.
      a_ft (float): effective area in square ft of rectangular repeater.
      lambda_ft (float): wavelength in ft of carrier.
    Output:
      (float[[theta, gain], ...]): list of antenna pattern angle in degrees
        and corresponding gain in dBi.
    '''

    pattern = []
    for theta_deg in np.arange(0, 100, 0.02):
        if theta_deg == 0:
            pattern.append([0,0,0])
        else:
            u,P = Prect_dB(a_ft, lambda_ft, theta_deg)
            if u >= 14.1: # 9*pi/2 == 14.1, approx peak of 4th lobe.
                break
            pattern.append([theta_deg, u, P])

    hp_dB = Gp_dBi/2 # Half power drop.
    for i_deg in range(int(theta_deg)+1, 21, 2):
        uPrev = u
        u = np.pi*a_ft/lambda_ft*np.sin(np.radians(i_deg))
        u = 20*np.log10(u) # Pattern envelope.
        if u >= hp_dB:
            break
        pattern.append([i_deg, 0, -u])

    m = (uPrev - u)/2
    for j_deg in range(i_deg, 151, 5):
        # The curve now decreases linearly.
        u -= m*2.5
        if u > hp_dB:
            break
        pattern.append([j_deg, 0, -u])

    # The pattern now holds constant at the front-to-back ratio.
    for k_deg in range(j_deg, 180, 10):
        pattern.append([k_deg, 0, -hp_dB])
    pattern.append([180, 0, -hp_dB])
    return pattern

def plotCartesian(pattern, subtitle=''):
    '''
    Inputs:
      pattern (float[[theta, u, gain], ...]): long list of 3-item lists made
        up of angle in degrees, parameter u and gain in dBi.
      subtitle (string): additional information to put plot heading.
    Output:
      none (plot displayed on screen).
    '''

    theta_deg = pattern[:,0]
    u = pattern[:,1]
    gain = pattern[:,2]

    # Pattern only reliable out to 4 lobes.
    stop = np.where(u[1:]==0)[0][0] # Point beyond 4 lobes.
    theta_deg = theta_deg[:stop]
    u = u[:stop]
    gain = gain[:stop]

    fig, ax = plt.subplots(figsize=(10,6))
    ax.grid(True,which="minor",alpha=0.2)
    ax.grid(True,which="major",alpha=0.5)
    ax.minorticks_on()

    spanT = 2*theta_deg[-1] # -theta to +theta.
    spanU = 2*u[-1] # -mu to +mu.
    # Fwd and inverse tick fns, mu x-axis. Linear theta-mu mapping.
    def theta2mu(t):  return (t +t[-1])/spanT*spanU - u[-1]
    def mu2theta(mu): return (mu+u[-1])/spanU*spanT - theta_deg[-1]

    ax2 = ax.secondary_xaxis('top', functions=(theta2mu,mu2theta))
    ax2.set_xlim(-u[-1], u[-1])
    ax2.set_xlabel('$\mu$') # Top x-axis in unit of antenna parameter mu.
    ax2.set_xlim(theta_deg[0], theta_deg[-1])

    plt.xlabel('$\\theta$ (deg)') # Bottom x-axis in degrees.
    plt.ylabel('Gain (dBi)')
    plt.xlim(-theta_deg[-1], theta_deg[-1])
    plt.ylim(-35, 0)
    plt.plot(-theta_deg[::-1], gain[::-1], color='blue', lw=1)
    plt.plot(theta_deg, gain, color='blue', lw=1)

    fig.suptitle('Rectangular Passive Repeater Gain')
    plt.title(subtitle, fontsize=10)
    fig.tight_layout()

    plt.show()

def plotPolar(pattern, subtitle=''):
    '''
    Inputs:
      pattern (float[[theta, u, gain], ...]): long list of 3-item lists made
        up of angle in degrees, parameter u and gain in dBi.
      subtitle (string): additional information to put plot heading.
    Output:
      none (plot displayed on screen).
    '''

    theta = np.radians(pattern[:,0])
    u = pattern[:,1]
    gain = pattern[:,2]

    fig, ax = plt.subplots(figsize=(8,8), subplot_kw={'projection':'polar'})
    ax.grid(True, which="minor", alpha=0.2)
    ax.grid(True, which="major", alpha=0.5)
    ax.set_rlabel_position(330)
    ax.minorticks_on()
    ax.set_thetagrids(range(0,360,20))
    ax.set_theta_zero_location("N") # 0 deg pointed upward.

    plt.suptitle('Rectangular Passive Repeater Gain')
    plt.title(subtitle, fontsize=10)
    plt.xlabel('Gain (dBi)')
    plt.plot(theta, gain, color='blue', lw=1)
    plt.plot(-theta, gain, color='blue', lw=1) # Use symmetry.
    plt.show()

def Prect_dB(a_ft, lambda_ft, theta_deg):
    '''For rectangular apertures, calculate normalized power radiation
    pattern.  From p. 43 and 45, "Parameter u is related to the dimensions of
    the aperture and wavelength with consideration given to aperture shape."

    Inputs:
      a_ft (float): effective area in square ft of rectangular repeater.
      lambda_ft (float): wavelength in ft of carrier.
      theta_rad (float): angle in degrees off boresight.
    Output:
    '''
    theta_rad = np.radians(theta_deg)
    u = np.pi*a_ft/lambda_ft*np.sin(theta_rad)
    # Normalized power radiation pattern, rectangular apertures.
    P_dB = 10*np.log10((np.sin(u)/u)**2)
    return u, P_dB

def usage(prog):
    s = os.path.basename(prog) + '\n'
    s += '-a, included angle in degrees made by repeater, tx and rx,\n'
    s += '-c, v or h (vertical or horizontal),\n'
    s += '-f, frequency in GHz,\n'
    s += '-h, horizontal length (width) in feet of repeater,\n'
    s += '-v, vertical length (height) in feet of repeater.'
    print(s)
    sys.exit(1)

def wavelength_ft(freq_GHz):
    '''Use c=f x lambda to calculate wavelength, but use correction factor to
    accomdate f in GHz and wavelength in feet.  For c in m/s, f in GHz then
    wavelength in feet is:

      lambda_ft = (c * 39.37/12) / (f * 1e9) = 0.9836/f

    Input:
      freq_GHz (float): frequency whose wavelength is wanted.
    Output:
      (float): wavelength in feet of frequency.
    '''
    return 0.9836/freq_GHz

def main(argv):
    '''Test program to calculate a far-field antenna pattern for a passive
    repeater with rectangular aperture.
    '''

    f_GHz, fStr, v_ft, h_ft, cut, alpha_deg = cli(argv) # User inputs.

    alpha_rad = np.radians(alpha_deg)           # Half included angle.
    a_ft = np.cos(alpha_rad)                    # Resolve plane.
    a_ft *= (h_ft if cut[0] == 'H' else v_ft)   # Effective length.
    lambda_ft = wavelength_ft(f_GHz)            # Wavelength in ft.
    aE_ft2 = v_ft*h_ft*np.cos(alpha_rad)        # Effective area.
    Gp_dBi = passiveGain_dBi(aE_ft2, lambda_ft) # Passive repeater gain.

    d_mi = nearFieldRange_mi(aE_ft2, lambda_ft)
    print('Near/far field transition (Microflect): %.2f miles' % d_mi)
    d_mi = nearFieldRangeCollins_mi(aE_ft2, f_GHz)
    print('Near/far field transition (Collins):    %.2f miles' % d_mi)
    print('')
#    print('Frequency:                    %s GHz.' % fStr)
#    print('Passive repeater size:        %.0f\' x% .0f\'' % (v_ft, h_ft))
#    print('Plane of radiation pattern:   %s' % cut)
#    print('Included angle between paths: %.2f deg' % alpha_deg)
    print('Effective aperture: %.1f sq. ft.' % aE_ft2)
    print('Gain:               %.1f dBi' % Gp_dBi)

    # Create plot subtitle with repeater info.
    subtitle = '%.1f dBi, ' % Gp_dBi
    subtitle += '%.0f\'x%.0f\', ' % (v_ft, h_ft)
    subtitle += '%s plane, %s GHz, ' % (cut, fStr)
    subtitle += '$\\alpha$=%.1f deg, ' % alpha_deg
    subtitle += '$A_{eff}$=%.1f ft$^2$' % aE_ft2

    pat = pattern(Gp_dBi, a_ft, lambda_ft) # Create antenna pattern.
    plotPolar(np.array(pat), subtitle)     # Gain vs theta.
    plotCartesian(np.array(pat), subtitle) # Gain vs parameter u.

if __name__ == '__main__': # When run as program.
    main(sys.argv)