Characteristic impedance & field pattern in arbitrary transmission lines

characteristicImpedanceTEM.py

# Transmission line characteristic impedance
# Amar, 12 July 2012
# version 3

# All lengths in cm, E-fields in V/cm, time in ns
import pylab as plt
import numpy as np
from numpy import pi
from scipy.weave import blitz
import cPickle, time
eta0 = 376.730313461    # Ohms

"""
# Blitz test:
a,b = 1,0
print blitz("b = a+a")
"""

################### Computation grid ##########################
size = 50.0             # mm
numberOfGridBits = 10 
numberOfGridPoints = 2**numberOfGridBits
dx = dy = 2*size/numberOfGridPoints

x = np.arange(-size/2.,size/2.,dx)
y = np.arange(-size/2.,size/2.,dy)

X,Y = np.meshgrid(x,y)
N = len(x)

################## Electrodes #################################

class CircularElectrode:
    def __init__(self,x0,y0,D,V0,innie=True):
        self.x0, self.y0, self.V0 = x0,y0,V0
        self.D = D
        if innie: self.mask = ( ((X-self.x0)**2 + (Y-self.y0)**2) <= self.D**2/4. )   
        else: self.mask = ( ((X-self.x0)**2 + (Y-self.y0)**2) >= self.D**2/4. )   

#electrodes = [ CircularElectrode(0,0,.24425,V0,innie=True), CircularElectrode(0,0,.5625,0,innie=False) ]   # GR900

class RectangularElectrode:
    def __init__(self,x0,y0,w,h,V0,innie=True):
        self.x0, self.y0, self.V0 = x0,y0,V0
        self.w,self.h = w,h
        if innie: self.mask = np.logical_and( np.abs(X-self.x0) <= self.w/2., np.abs(Y-self.y0) <= self.h/2. )   
        else: self.mask = np.logical_or( np.abs(X-self.x0) >= self.w/2., np.abs(Y-self.y0) >= self.h/2. )  
        # Note this subtleness between logical_and vs. logical_or for innies vs. outies


################# Calculation routine ########################

expr = "V[step:-step, step:-step] = ((V[0:-2*step, step:-step] + V[2*step:, step:-step]) + (V[step:-step,0:-2*step] + V[step:-step, 2*step:]))/4.0"
expr2 = "V = V_mask_inv*V + V_mask*V_electrode"

def calculateVoltage(electrodes, V0=1.0):
    """Calculation of voltage for an electrode configuration"""

    V_electrode = np.zeros(X.shape)
    for electrode in electrodes: V_electrode[electrode.mask] = electrode.V0
    V_mask = sum([electrode.mask for electrode in electrodes])
    V_mask_inv = np.logical_not(V_mask)

    V = V_mask_inv*np.average(V_electrode)   # initialized

    print "Preconditioning ...",
    for i in xrange(4,numberOfGridBits):
        step = 2**numberOfGridBits/2**i
        for j in xrange(int(5e1)):    
            blitz(expr, check_size=0)
            blitz(expr2,check_size=0)
    print "done"

    step = 1.0
    for i in xrange(int(4e3)):
        blitz(expr, check_size=0)
        blitz(expr2,check_size=0)
    V_unref = V.copy()
    for i in xrange(int(1e3)):
        blitz(expr, check_size=0)
        blitz(expr2,check_size=0)
    V_ref = V.copy()
    fractionalChange = np.sqrt( np.sum((V_ref - V_unref)**2)/np.sum(V_unref**2) )

    # Electric field
    ESquared = np.zeros(X.shape)
    ESquared[1:-1, 1:-1] = ((V[2:, 1:-1] - V[0:-2, 1:-1])/(2*dx))**2 + ((V[1:-1, 2:] - V[1:-1,0:-2])/(2*dy))**2
    
    # Characteristic impedance
    Z0 = eta0/( np.sum(ESquared)*dx*dy ) * V0**2
    return Z0, V, (fractionalChange, V_ref-V_unref)       # characteristic impedance, Ohms


##################### Test problem #############################
V0 = 1.0
outerConductor = CircularElectrode(0,0,46.06,0,innie=False)
innerConductor = CircularElectrode(0,0,20.00,V0,innie=True)
electrodes = [ outerConductor, innerConductor  ]
Z0,V,convergence = calculateVoltage(electrodes,V0)
print "Z0 = ", round(Z0,2), 
print "converged to ", convergence[0]

V0 = 2.0
outerConductor = RectangularElectrode(0,0,49.6,49.6,0,innie=False)
innerConductor = RectangularElectrode(0,0,19.8,19.8,V0,innie=True)
electrodes = [ outerConductor, innerConductor  ]
Z0,V,convergence = calculateVoltage(electrodes,V0)
print "Z0 = ", round(Z0,2), 
print "converged to ", convergence[0]


######### Plots ##############################
# Electric field meshgrid
Ex,Ey = np.zeros(X.shape),np.zeros(X.shape)
Ex[1:-1, 1:-1] = (V[2:, 1:-1] - V[0:-2, 1:-1])/(2*dx) 
Ey[1:-1, 1:-1] = (V[1:-1, 2:] - V[1:-1,0:-2])/(2*dy)
E = np.sqrt( Ex**2 + Ey**2 )


plt.figure()
plt.title(r"Voltage (V)")
plt.xlabel("x (mm)")
plt.ylabel("y (mm)")
C1 = plt.contourf(X,Y,V,30)
plt.colorbar(C1)

p.figure()
p.title(r"Electric field, $\mathcal{E}_x$ (V/cm)")
p.xlabel("x (mm)")
p.ylabel("y (mm)")
C2 = p.contourf(X,Y,E,30)
p.colorbar(C2)


##################### Batch problem setup ######################

outerConductor = RectangularElectrode(0,0,30,50,0,innie=False)

V0 = 1.0
w = np.arange(5,25,2)
h = np.arange(1,15,2)

W,H = np.meshgrid(w,h)
Z = np.zeros(W.shape)

for i in xrange(len(w)):
    for j in xrange(len(h)):
        innerConductor = RectangularElectrode(0,-15,w[i],h[j],V0,innie=True)
        electrodes = [ outerConductor, innerConductor  ] 
        Z0,V,junk = calculateVoltage(electrodes)        
        Z[j,i] = Z0
        print w[i],h[j],round(Z0)

p.figure()
p.title(r"Characteristic impedance ($\Omega$)")
p.xlabel("Center conductor width, w (mm)")
p.ylabel("Center conductor height, h (mm)")
C = p.contourf(W,H,Z,30)
p.colorbar(C)

p.figure()
p.title("Characteristic impedance ($\Omega$)\n"+"Inner conductor @ "+`(innerConductor.x0,innerConductor.y0)`)
p.xlabel("Center conductor width, w (mm)")
p.ylabel("Center conductor height, h (mm)")
C = p.contour(W,H,Z,np.arange(0,100,5), linewidth=10, colors='b')
p.clabel(C, inline=1, fontsize=12, linewidth=5)