# Purpose:
#         This is a library of functions for ELEN 636 that
#         so far has the ability to calculate the Sine and
#         Cosine integrals as well as the mutual impedance
#         between two parallel antennas.
#
# Author: Christopher Smith
# E-mail: mc_anjo@tamu.edu
# Date:   10/30/2006

###############
### NOTE: The functions below for the sine and cosine integrals are similar
###       to the functions I turned in for homework assignment 4 problem 6
###       except that I added the ability to check for convergence.
###       I also added the factor into the mutual impedance formula so that the
###       answer is given in terms of the terminal input impedance instead of
###       the loop impedance as it was formerly giving.
###############

# depends on the math library
from math import *
import numpy.numarray, numpy

def factorial(n):
    """
    This function calculates the factorial of a number.
    """
    sum = 1.0
    for m in range(1, int(n)+1):
        sum = float(m)*sum
    return sum

def Si(x):
    """
    This function computes the sine integral. It uses a power series
    expansion that can be found in Abramowitz and Stegun's math
    functions reference book.
    """
    start = 0.0
    stop = 10.0
    sine_int = 0.0
    convergence = 1.0*10**(-6)   # want to have the difference between
                                 # the last run and this run below
                                 # this value
    while 1:
        for n in range(int(start), int(stop)):
            n = float(n)
            sine_int += ((-1)**n)*x**(2*n +1)/((2*n+1)*factorial(2*n+1))
        sine_int_new = sine_int + ((-1.)**stop)*x**(2.*stop +1.)/((2.*stop+1.)*factorial(2.*stop+1.))
        converge_check = sine_int_new - sine_int
        if abs(converge_check) < convergence:
            break
        else:
            start = stop
            stop += 5.0
    return sine_int_new

def Ci(x):
    """
    This function computes the cosine integral. It uses a power series
    expansion that can be found in Abramowitz and Stegun's math
    functions reference book.
    """
    start = 1.0
    stop = 10.0
    convergence = 1.0*10.**(-6)   # want to have the difference between
                                 # the last run and this run below
                                 # this value
    # The first number in the sum is Euler's constant to 10 digits
    cosine_int = 0.5772156649 + log(x)
    while 1:
        for n in range(int(start), int(stop)):
            m = float(n)
            cosine_int = cosine_int +((-1)**m)*x**(2*m)/((2*m)*factorial(2*m))
        cosine_int_new = cosine_int + ((-1)**stop)*x**(2*stop)/((2*stop)*factorial(2*stop))
        converge_check = cosine_int_new - cosine_int
        if abs(converge_check) < convergence:
            break
        else:
            start = stop
            stop += 5.0
        #print stop
    return cosine_int_new


def mutual_impedance(length1_tot, length2_tot, stagger, d):
    """
    This function computes the mutual impedance between two antennas
    for the Parallel in Echelon Configuration. The formulas are taken
    from a paper by Howard King, "Mutual Impedance of Unequal Length
    Antennas in Echelon"
    
    NOTE: all measurements should be entered in wavelengths
    """

    # stagger (this is the vertical separation between antenna centers)
    # d (this is the horizontal separation between the antennas)
    
    # length1 and length2 are the half length of the antennas, this is
    # to conform to King's formulas
    length1 = length1_tot/2.0
    length2 = length2_tot/2.0

    # vertical separation between center of antenna 1 and bottom of antenna 2
    h = stagger - length2

    # wave propagation constant
    beta = 2.0*pi

    # formulas to put into mutual impedance equation
    u0 = beta*(sqrt(d**2 +(h -length1)**2) +(h -length1))
    v0 = beta*(sqrt(d**2 +(h -length1)**2) -(h -length1))
    u0prime = beta*(sqrt(d**2 +(h +length1)**2) -(h +length1))
    v0prime = beta*(sqrt(d**2 +(h +length1)**2) +(h +length1))
    u1 = beta*(sqrt(d**2 +(h -length1 +length2)**2) +(h -length1 +length2))
    v1 = beta*(sqrt(d**2 +(h -length1 +length2)**2) -(h -length1 +length2))
    u2 = beta*(sqrt(d**2 +(h +length1 +length2)**2) -(h +length1 +length2))
    v2 = beta*(sqrt(d**2 +(h +length1 +length2)**2) +(h +length1 +length2))
    u3 = beta*(sqrt(d**2 +(h -length1 +2.0*length2)**2) +(h -length1 +2.0*length2))
    v3 = beta*(sqrt(d**2 +(h -length1 +2.0*length2)**2) -(h -length1 +2.0*length2))
    u4 = beta*(sqrt(d**2 +(h +length1 +2.0*length2)**2) -(h +length1 +2.0*length2))
    v4 = beta*(sqrt(d**2 +(h +length1 +2.0*length2)**2) +(h +length1 +2.0*length2))
    w1 = beta*(sqrt(d**2 +h**2) -h)
    y1 = beta*(sqrt(d**2 +h**2) +h)
    w2 = beta*(sqrt(d**2 +(h +length2)**2) -(h +length2))
    y2 = beta*(sqrt(d**2 +(h +length2)**2) +(h +length2))
    w3 = beta*(sqrt(d**2 +(h +2.0*length2)**2) -(h +2.0*length2))
    y3 = beta*(sqrt(d**2 +(h +2.0*length2)**2) +(h +2.0*length2))
    #print u0,v0,u0prime,v0prime,u1,v1,u2,v2,u3,v3,u4,v4,w1,y1,w2,y2,w3,y3

    # real part of the mutual impedance between two antennas
    R12 = 15*(cos(beta*(length1 -h))*(Ci(u0) +Ci(v0) -Ci(u1) -Ci(v1)) \
              +sin(beta*(length1 -h))*(-Si(u0) +Si(v0) +Si(u1) -Si(v1)) \
              +cos(beta*(length1 +h))*(Ci(u0prime) +Ci(v0prime) -Ci(u2) -Ci(v2)) \
              +sin(beta*(length1 +h))*(-Si(u0prime) +Si(v0prime) +Si(u2) -Si(v2)) \
              +cos(beta*(length1 -2.0*length2 -h))*(-Ci(u1) -Ci(v1) +Ci(u3) +Ci(v3)) \
              +sin(beta*(length1 -2.0*length2 -h))*(Si(u1) -Si(v1) -Si(u3) +Si(v3)) \
              +cos(beta*(length1 +2.0*length2 +h))*(-Ci(u2) -Ci(v2) +Ci(u4) +Ci(v4)) \
              +sin(beta*(length1 +2.0*length2 +h))*(Si(u2) -Si(v2) -Si(u4) +Si(v4)) \
              +2.0*cos(beta*length1)*cos(beta*h)*(-Ci(w1) -Ci(y1) +Ci(w2) +Ci(y2)) \
              +2.0*cos(beta*length1)*sin(beta*h)*(Si(w1) -Si(y1) -Si(w2) +Si(y2)) \
              +2.0*cos(beta*length1)*cos(beta*(2.0*length2 +h))*(Ci(w2) +Ci(y2) -Ci(w3) -Ci(y3)) \
              +2.0*cos(beta*length1)*sin(beta*(2.0*length2 +h))*(-Si(w2) +Si(y2) +Si(w3) -Si(y3)))

    # imaginary part of the mutual impedance between two antennas
    X12 = 15*(cos(beta*(length1 -h))*(-Si(u0) -Si(v0) +Si(u1) +Si(v1)) \
              +sin(beta*(length1 -h))*(-Ci(u0) +Ci(v0) +Ci(u1) -Ci(v1)) \
              +cos(beta*(length1 +h))*(-Si(u0prime) -Si(v0prime) +Si(u2) +Si(v2)) \
              +sin(beta*(length1 +h))*(-Ci(u0prime) +Ci(v0prime) +Ci(u2) -Ci(v2)) \
              +cos(beta*(length1 -2.0*length2 -h))*(Si(u1) +Si(v1) -Si(u3) -Si(v3)) \
              +sin(beta*(length1 -2.0*length2 -h))*(Ci(u1) -Ci(v1) -Ci(u3) +Ci(v3)) \
              +cos(beta*(length1 +2.0*length2 +h))*(Si(u2) +Si(v2) -Si(u4) -Si(v4)) \
              +sin(beta*(length1 +2.0*length2 +h))*(Ci(u2) -Ci(v2) -Ci(u4) +Ci(v4)) \
              +2.0*cos(beta*length1)*cos(beta*h)*(Si(w1) +Si(y1) -Si(w2) -Si(y2)) \
              +2.0*cos(beta*length1)*sin(beta*h)*(Ci(w1) -Ci(y1) -Ci(w2) +Ci(y2)) \
              +2.0*cos(beta*length1)*cos(beta*(2.0*length2 +h))*(-Si(w2) -Si(y2) +Si(w3) +Si(y3)) \
              +2.0*cos(beta*length1)*sin(beta*(2.0*length2 +h))*(-Ci(w2) +Ci(y2) +Ci(w3) -Ci(y3)))

    R12_in = R12/(sin(beta*length1)*sin(beta*length2))
    X12_in = X12/(sin(beta*length1)*sin(beta*length2))
    impedance = (R12_in, X12_in)
    return impedance


def top_row_matrix(length1, length2, stagger, stagger_image, radius, m, n):
    """
    This function will find the top row of a mutual impdedance matrix
    over a ground plane. From the top row we can find the overall matrix
    since it is a block Toeplitz matrix
    """

    z = []  # list to store our impedance values

    # index to step over for the staggering and separation between antennas
    stagger_range = range(0, m)
    separation_range = range(0, n)

    # calculate the mutual impedance values for the real planar array
    # the first loop gives us the stagger between rows while the second
    # loop gives us the separation between dipoles on the same row
    for m in stagger_range:
       for n in separation_range:
           h = stagger*m
           d = stagger*n 
           if d == 0: d = radius
           trans = mutual_impedance(length1, length2, h, d)
           z.append(complex(trans[0],trans[1]))

    # Suppose the real antenna array is a plane in x-y at z = 0 and the
    # imaginary antenna array is a plane in x-y at z = some spacing
    # Since the mutual impedance would include spacing in the
    # x,y and z directions we need to condense it down to two spacings
    # so that we can input the spacings into our formulas. For the
    # spacing between dipoles we will use the hypotenuse formed by the
    # x and z coordinates. The y coordinate will be the staggering between
    # dipoles.
    for m in stagger_range:
        for n in separation_range:
            separation = sqrt((stagger_image)**2+(stagger*n)**2)
            trans = mutual_impedance(length1, length2, stagger*m, separation)
            z.append(complex(trans[0],trans[1]))

    return z


def Toeplitz(a):
    """
    This function takes a list in. The list represents
    the top row in a Toeplitz matrix. From this row we
    can fill the rest of the matrix and then output it.

    Note: This function has a dependency on the numpy
          library.
    """

    # need the length of the list being input
    n = len(a)

    # fill a matrix the size we need with zeros
    matrix = numpy.zeros((n,n), dtype = type(a))

    a_new = []
    for i in range(-n+1,n):
        a_new.append( a[ abs(i) ] )

    for i in range(0,n):
        matrix[i,:] = a_new[ (n-1) -i : (2*n -1) -i ]

    return matrix


def BlockToeplitz(a, M, N):
    """
    This function takes a list of N Toeplitz matrices in
    and creates a Block Toeplitz matrix from them.

    Note: This function has a dependency on the numpy
          library.
          This function is also specifically for an array
          over a groundplane problem. Which means that
          the matrix is twice the length and width of a
          normal array matrix. The code would need to be
          modified to apply to a normal Block-Toeplitz.
    """


    # fill a matrix the size we need with zeros
    matrix = numpy.zeros((2*M*N,2*M*N), dtype = type(a))

    a_new = []
    for i in range(-2*N+1,2*N):
        a_new.append( a[ abs(i) ] )

    for i in range(0, 2*N):
        i_new = i*M
        for j in range(0, 2*N):
            j_new = j*M
            matrix[ i_new:i_new+N, j_new:j_new+N ] = a_new[2*N-1-j+i]

    return matrix


def h_array_over_gndplane(M, N):
    """
    This function generates the voltage excitation matrix for an
    array over a ground plane. The antennas are parallel to the
    ground plane which means that according to image theory the
    current along the antenna images are in the opposite direction
    from those of the real array.

    Note: This function has a dependency on the numpy library.
    """
    
    h = []
    for n in range(0, 2*M*N):
        if n < M*N:
            h.append(1.0)
        else:
            h.append(-1.0)
    h = numpy.numarray.array((h), shape = (2*M*N,1))
    return h

