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.oldnumeric.linear_algebra.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:\Documents and Settings\Chris &

Esther\Desktop\636_hw5_2\elen636_hw5_2.py", line 60, in <module>

matrix_inverse =

numpy.numarray.linear_algebra.generalized_inverse( matrix)

File

"C:\Python25\lib\site-packages\numpy\oldnumeric\linear_algebra.py", line

59, in generalized_inverse

return linalg.pinv(a, rcond)

File "C:\Python25\lib\site-packages\numpy\linalg\linalg.py", line

557, in pinv

u, s, vt = svd(a, 0)

File "C:\Python25\lib\site-packages\numpy\linalg\linalg.py", line

485, in svd

a = _fastCopyAndTranspose(t, a)

File "C:\Python25\lib\site-packages\numpy\linalg\linalg.py", line

107, in _fastCopyAndTranspose

cast_arrays = cast_arrays + (_fastCT(a.astype(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.py 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*****@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,y 1,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

# 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*****@tamu.edu

# 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_impedance, top_row_matrix, Toeplitz

from elen636_math import BlockToeplitz, h_array_over_gndplane

# 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_matrices = [] # a list of N Toeplitz matrices

for i in range(0, 2*M*N, N):

toeplitz_matrices.append( Toeplitz( z[i:i +N] ) )

# Now we fill our Block-Toeplitz matrix

matrix = BlockToeplitz(toeplitz_matrices, M, N)

# generate the voltage excitations of all the dipoles

h = h_array_over_gndplane(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_impedance_element25 = input_impedance_matrix[24]

# creates and opens files to output results

try:

file = open( 'mutual_impedance_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 - Tu***@python.org

http://mail.python.org/mailman/listinfo/tutor