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: [EMAIL PROTECTED] # 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
# 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: [EMAIL PROTECTED] # 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 - Tutor@python.org http://mail.python.org/mailman/listinfo/tutor
-- http://mail.python.org/mailman/listinfo/python-list