numpy help

Howdy,

I'm a college student and for one of we are writing programs to
numerically compute the parameters of antenna arrays. I decided to use
Python to code up my programs. Up to now I haven't had a problem,
however we have a problem set where we are creating a large matrix and
finding it's inverse to solve the problem. To invert the matrix I've
tried using numpy.numarray. linear_algebra. inverse and
numpy.oldnumeri c.linear_algebr a.inverse which both give me the same
error ( I was hoping they called different routines but I think they
call the same one ).


This is the error message I receive:

Traceback (most recent call last):
File "C:\Documen ts and Settings\Chris &
Esther\Desktop\ 636_hw5_2\elen6 36_hw5_2.py", line 60, in <module>
matrix_inverse =
numpy.numarray. linear_algebra. generalized_inv erse(matrix)
File
"C:\Python25\li b\site-packages\numpy\ oldnumeric\line ar_algebra.py", line
59, in generalized_inv erse
return linalg.pinv(a, rcond)
File "C:\Python25\li b\site-packages\numpy\ linalg\linalg.p y", line
557, in pinv
u, s, vt = svd(a, 0)
File "C:\Python25\li b\site-packages\numpy\ linalg\linalg.p y", line
485, in svd
a = _fastCopyAndTra nspose(t, a)
File "C:\Python25\li b\site-packages\numpy\ linalg\linalg.p y", line
107, in _fastCopyAndTra nspose
cast_arrays = cast_arrays + (_fastCT(a.asty pe(type)),)
TypeError: can't convert complex to float; use abs(z)



I've tried inverting small complex matrices and it worked fine. Does
anyone know why it won't work for this larger matrix? Any ideas how I
can work around this problem and get the correct inverse matrix?

Chris

P.S. elen636_math.py is my personal library of functions I've create to
solve the problem while elen636_hw5_2.p y is the program that I'm
actually running


# 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.ed u
# 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)*factor ial(2*n+1))
sine_int_new = sine_int + ((-1.)**stop)*x**( 2.*stop +1.)/((2.*stop+1.)*f actorial(2.*sto p+1.))
converge_check = sine_int_new - sine_int
if abs(converge_ch eck) < 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)*factoria l(2*m))
cosine_int_new = cosine_int + ((-1)**stop)*x**(2 *stop)/((2*stop)*facto rial(2*stop))
converge_check = cosine_int_new - cosine_int
if abs(converge_ch eck) < convergence:
break
else:
start = stop
stop += 5.0
#print stop
return cosine_int_new


def mutual_impedanc e(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,v 0prime,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*(l ength1 -h))*(Ci(u0) +Ci(v0) -Ci(u1) -Ci(v1)) \
+sin(beta*(leng th1 -h))*(-Si(u0) +Si(v0) +Si(u1) -Si(v1)) \
+cos(beta*(leng th1 +h))*(Ci(u0prim e) +Ci(v0prime) -Ci(u2) -Ci(v2)) \
+sin(beta*(leng th1 +h))*(-Si(u0prime) +Si(v0prime) +Si(u2) -Si(v2)) \
+cos(beta*(leng th1 -2.0*length2 -h))*(-Ci(u1) -Ci(v1) +Ci(u3) +Ci(v3)) \
+sin(beta*(leng th1 -2.0*length2 -h))*(Si(u1) -Si(v1) -Si(u3) +Si(v3)) \
+cos(beta*(leng th1 +2.0*length2 +h))*(-Ci(u2) -Ci(v2) +Ci(u4) +Ci(v4)) \
+sin(beta*(leng th1 +2.0*length2 +h))*(Si(u2) -Si(v2) -Si(u4) +Si(v4)) \
+2.0*cos(beta*l ength1)*cos(bet a*h)*(-Ci(w1) -Ci(y1) +Ci(w2) +Ci(y2)) \
+2.0*cos(beta*l ength1)*sin(bet a*h)*(Si(w1) -Si(y1) -Si(w2) +Si(y2)) \
+2.0*cos(beta*l ength1)*cos(bet a*(2.0*length2 +h))*(Ci(w2) +Ci(y2) -Ci(w3) -Ci(y3)) \
+2.0*cos(beta*l ength1)*sin(bet a*(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*(l ength1 -h))*(-Si(u0) -Si(v0) +Si(u1) +Si(v1)) \
+sin(beta*(leng th1 -h))*(-Ci(u0) +Ci(v0) +Ci(u1) -Ci(v1)) \
+cos(beta*(leng th1 +h))*(-Si(u0prime) -Si(v0prime) +Si(u2) +Si(v2)) \
+sin(beta*(leng th1 +h))*(-Ci(u0prime) +Ci(v0prime) +Ci(u2) -Ci(v2)) \
+cos(beta*(leng th1 -2.0*length2 -h))*(Si(u1) +Si(v1) -Si(u3) -Si(v3)) \
+sin(beta*(leng th1 -2.0*length2 -h))*(Ci(u1) -Ci(v1) -Ci(u3) +Ci(v3)) \
+cos(beta*(leng th1 +2.0*length2 +h))*(Si(u2) +Si(v2) -Si(u4) -Si(v4)) \
+sin(beta*(leng th1 +2.0*length2 +h))*(Ci(u2) -Ci(v2) -Ci(u4) +Ci(v4)) \
+2.0*cos(beta*l ength1)*cos(bet a*h)*(Si(w1) +Si(y1) -Si(w2) -Si(y2)) \
+2.0*cos(beta*l ength1)*sin(bet a*h)*(Ci(w1) -Ci(y1) -Ci(w2) +Ci(y2)) \
+2.0*cos(beta*l ength1)*cos(bet a*(2.0*length2 +h))*(-Si(w2) -Si(y2) +Si(w3) +Si(y3)) \
+2.0*cos(beta*l ength1)*sin(bet a*(2.0*length2 +h))*(-Ci(w2) +Ci(y2) +Ci(w3) -Ci(y3)))

R12_in = R12/(sin(beta*lengt h1)*sin(beta*le ngth2))
X12_in = X12/(sin(beta*lengt h1)*sin(beta*le ngth2))
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_rang e = 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_rang e:
h = stagger*m
d = stagger*n
if d == 0: d = radius
trans = mutual_impedanc e(length1, length2, h, d)
z.append(comple x(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_rang e:
separation = sqrt((stagger_i mage)**2+(stagg er*n)**2)
trans = mutual_impedanc e(length1, length2, stagger*m, separation)
z.append(comple x(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_gn dplane(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


# Purpose:
# This program will calculate the impedance matrix
# for a planar antenna array over a ground plane.
# This is for ELEN 636 homework 5 problem 2
#
# Author: Christopher Smith
# E-mail: mc_anjo@tamu.ed u
# Date: 10/30/06

# library so arrays (matrices) can be used
import numpy.numarray, numpy
# library so we can perform linear algebra operations
# like inverse( )
import numpy.numarray. linear_algebra
# import special functions
from elen636_math import mutual_impedanc e, top_row_matrix, Toeplitz
from elen636_math import BlockToeplitz, h_array_over_gn dplane
# math library
from math import *


# length and radius of dipoles (wvlgths)
# horizontal (same in x and y) stagger of matrix (wvlgths)
# vertical separation between array and it's image (wvlgths)
# or twice the vertical separation between array and ground plane
# and the number of elements in the y and x directions m, n
length1 = 0.5
length2 = length1
radius = 0.0025
stagger = 0.55
stag_image = 0.5
M = 7
N = 7



# function to give us back the top row of the m x n matrix
z = top_row_matrix( length1, length2, stagger, stag_image, radius, M, N)

# We now have a list with 2*M*N elements. Using symmetry we can create
# N NxN arrays which correspond to N Toeplitz matrices for our N x N
# Block Toeplitz matrix of 2*M*N x 2*M*N elements. So from those N
# matrices we can completely fill the matrix.
toeplitz_matric es = [] # a list of N Toeplitz matrices
for i in range(0, 2*M*N, N):
toeplitz_matric es.append( Toeplitz( z[i:i +N] ) )

# Now we fill our Block-Toeplitz matrix
matrix = BlockToeplitz(t oeplitz_matrice s, M, N)

# generate the voltage excitations of all the dipoles
h = h_array_over_gn dplane(M, N)





# now we compute Cn by multiply the inverse mutual impedance matrix
# and multiplying it by the h matrix
matrix_inverse = numpy.numarray. linear_algebra. inverse(matrix)
Cn = numpy.numarray. matrixmultiply( matrix_inverse, h)

input_impedance _matrix = h/Cn

#input_impedanc e_element25 = input_impedance _matrix[24]

# creates and opens files to output results
try:
file = open( 'mutual_impedan ce_results.txt' , "w" )
except IOError, message:
print >sys.stderr, "File could not be opened:" , message
sys.exit(1)

print >file, input_impedance _matrix
file.close()









_______________ _______________ _______________ __
Tutor maillist - Tutor@python.or g