Pythagon: April 2013

I couldn't find any good Smith chart plotters for Python, so I coded one up using matplotlib. The new matplotlib widgets are great!

Here's a screenshot of a coax-waveguide match tuner.

Here's the code for the Smith chart plot:

	# Interactive tuner for broadband coax-waveguide launcher
	# Analytic formulae from Slater, / Microwave Transmission /
	# Amar, version 7, 2013-04-09
	# Changes:
	# v2: S-parameters and Smith plot
	# v3: polar smith plot
	# v4: S and T matrices to get at S12 phase, objectized circuit description. Reverted to cartesian smith plot.
	# v5: double launcher configuration, freed up 'a' as parameter
	# v6: cleaned up code, parameters are general and swappable now, plotting and tuning is modular
	# v7: phase and amplitude sensitivity plots

	# All lengths in cm
	# All frequencies in MHz
	# L is in uH, C is in uF

	import numpy as np
	from numpy import pi, sin, cos, tan
	from scipy import optimize
	import scipy.linalg as lg
	from functools import partial
	import matplotlib.pyplot as plt
	from matplotlib.widgets import Slider, Button, RadioButtons

	########## Microwave formulae ####################

	Z0 = 50. # Ohms, characteristic impedance of feedline
	eta0 = 376.730313461 # Ohms, impedance of free space
	speedOfLight = 29.97e3 # MHz cm, speed of light
	I = np.mat(np.eye(2))

	# Convention used for T-matrix definition
	# dicke: (port 2) = T (port 1)
	# agilent: (port 1) = T (port 2)
	convention = "agilent"

	def SfromZ(Z): return (Z-I) * lg.inv(Z+I)
	def ZfromS(S): return (I+S) * lg.inv(I-S)

	if convention == "dicke":
	def TfromS(S): return (1/S[1,0]) * np.mat( [ [-lg.det(S), S[1,1]],[-S[0,0],1] ] ) # definition of T to match Dicke convention: (b2 a2)^T = T (a1 b1)^T
	def SfromT(T): return (1/T[1,1]) * np.mat( [ [T[1,0], 1],[1,T[0,1]] ] ) # definition of T to match Dicke convention: (b2 a2)^T = T (a1 b1)^T
	elif convention == "agilent":
	def TfromS(S): return (1/S[1,0]) * np.mat( [ [-lg.det(S), S[0,0]],[-S[1,1],1] ] ) # definition of T to match Agilent convention: (b1 a1)^T = T (a2 b2)^T
	def SfromT(T): return (1/T[1,1]) * np.mat( [ [T[0,1], lg.det(T)],[1,-T[1,0]] ] ) # definition of T to match Agilent convention: (b1 a1)^T = T (a2 b2)^T


	class TwoPort:
	def __init__(self,S,label=''):
	self.S = S
	self.label = label

	class Line(TwoPort):
	def __init__(self,beta,l,label=''):
	"""propagation constant = beta, length = l"""
	self.label = label
	self.S = np.mat([ [0,np.exp(-1jbetal)],[np.exp(-1jbetal),0] ])

	class TSection(TwoPort):
	def __init__(self,Zs1,Z_12,Zs2,label=''):
	""" T-section impedance network
	-------[Zs1]------------[Zs2]-------
	\|
	\|
	[Z_12]
	\|
	\|
	------------------------------------- """
	self.label = label
	Z_11 = Zs1 + Z_12
	Z_22 = Zs2 + Z_12
	self.Z = np.mat( [[Z_11,Z_12],[Z_12,Z_22]] )/Z0 # normalizes to Z0
	self.S = SfromZ(self.Z)

	class Joint(TwoPort):
	def __init__(self,Z1,Z2,label=''):
	""" (zero length) joint between different char. impedance transmission lines
	--------------*=================
	Z1 Z2
	--------------*================= """
	self.label = label
	Gamma = (Z2-Z1)/(Z2+Z1)
	self.S = np.mat([ [Gamma,np.sqrt(1-Gamma2)],[np.sqrt(1-Gamma2),-Gamma] ])


	class Cascade(TwoPort):
	def __init__(self,cascadeList,label=''):
	""" 2 port elements in order from left to right, inputs at left """
	self.network = cascadeList
	self.label = label
	T = I
	if convention == "dicke":
	for dev_i in cascadeList: T = TfromS(dev_i.S) * T # Dicke convention
	elif convention == "agilent":
	for dev_i in cascadeList: T = T * TfromS(dev_i.S) # Agilent convention
	self.S = SfromT(T)


	class SymmetryPlane(TwoPort):
	def __init__(self,symmetry,label='symmetry plane'):
	"""Symmetry plane: symmetry = \pm 1"""
	self.S = np.mat([ [0,symmetry],[symmetry,0] ])


	def makeLauncher(y,x,f=910.,verbose=False):
	""" Description of the coax-waveguide launcher """

	c,l_seriesStub,z_shuntStub,l_shuntStub,l_WG = x # tunable parameters
	a,b,probeDiameter = y # fixed parameters
	d = a/2.
	l_probe = b
	Z_seriesStub = 50.00
	l_coax = 5.

	fc = speedOfLight/(2*a)

	alpha = 2(l_probe/b)2 np.sin(pid/a)*2 # coupling strength (Slater)
	#Z_seriesStub = eta0 * np.log(D_seriesStub/d_seriesStub)/(2*pi)

	gamma = 1/np.sqrt( 1- (fc/f)**2 )
	beta0 = 2pif/speedOfLight
	beta_WG = beta0/gamma
	Z_eq = eta0 * gamma * (b/a)

	# Circuit model for probe antenna, from Marcuvitz
	X_b = Z_eq * (beta_WG * a/(2pi)) (piprobeDiameter/a)2 /( 1+ 11./24 (pid/a)*2 ) # 5.11 (4b)
	s = 0
	for n in range(3,100,2): s+= (1/np.sqrt(n*2 - (beta0a/pi)**2) - 1./n)
	S0 = np.log(4a/(piprobeDiameter)) - 2. + 2*s
	X_a = 0.5X_b + Z_eq beta_WGa/(4pi) * ( S0 -(beta0probeDiameter/4)*2 ) # 5.11 (3)

	#X_a = 0.55 * Z_eq * (abeta_WG/pi) /np.sin(pid/a)**2 # from Marcuvitz
	#C = 88.54e-9 * 10. * (pi/4.)0.62*2 / probeGap # epsilon_0 = 8.854 pF/m = 88.54 fF/cm

	X_series = Z_seriesStub * np.tan(beta0 * l_seriesStub)
	X_backshort = alphaZ_eq np.tan(beta_WG * c) # prescription from Slater

	# launcher
	Z_12 = -1jX_b + 1jX_backshort
	Z_s1 = 1jX_a + 1jX_series
	Z_s2 = -1j*X_b
	launcher = TSection(Z_s1,Z_12,Z_s2,label="launch antenna")
	reverseLauncher = TSection(Z_s2,Z_12,Z_s1,label="pickup antenna")

	# waveguide joint
	Gamma = (alphaZ_eq - Z0)/(alphaZ_eq + Z0)
	TLtoWGJoint = Joint(Z0,alpha*Z_eq,label='coax-waveguide joint')
	WGtoTLJoint = Joint(alpha*Z_eq,Z0,label='waveguide-coax joint')

	# lines
	waveguideLine = Line(beta_WG,l_WG,label='waveguide line') # length of waveguide
	coaxLine = Line(beta0,l_coax,label='coax line') # length of coax between launcher and shunt stub

	# line + shunt stub
	X_shuntStub = Z0np.tan(beta0 l_shuntStub)
	shuntStub = TSection(0,1j*X_shuntStub,0,label='shunt stub')
	stubDistance = Line(beta0,z_shuntStub,label='shunt stub location')

	# shorting plane
	shortingPlane = SymmetryPlane(-1)

	# Full network
	network = Cascade( [coaxLine,shuntStub,stubDistance,launcher,TLtoWGJoint,waveguideLine,shortingPlane,waveguideLine,WGtoTLJoint,reverseLauncher,stubDistance,shuntStub,coaxLine], label='full launcher' )

	if verbose:
	print "fc = ", fc, "MHz"
	print "gamma = ",gamma
	print "alpha = ",alpha
	print "Z_eq = ",Z_eq,"Ohms"
	print "lambda_g = ",2*pi/beta_WG,"cm","\n"
	print "X_bs = ",X_backshort,"Ohms"
	print "L_probe = ",X_a/(2pif1e6) 1e9,"nH"
	print "L_seriesStub = ",X_series/(2pif1e6) 1e9,"nH"
	#print "L_shuntStub = ",X_shuntStub/(2pif1e6) 1e9,"nH"

	return network


	def S(y,x,f=910.): return makeLauncher(y,x,f).S # return S-matrix of cascaded elements in full network


	def targetFunction(y,x,i=0):
	""" Weighted average of S11 over a band of frequencies, to get wideband optimum """
	frequenciesWeights = [ (870,0.7), (890,0.9), (900,0.95), (905,1.1), (910,1.1), (915,1.1), (920,0.95), (930,0.9), (950,0.7) ]
	weightedSiisquared = [ np.abs(S(y,x,f=fi))[i,i]*2 wi for (fi,wi) in frequenciesWeights ]
	return np.average(weightedSiisquared)



	################ Parameter data structure #########
	class Parameter:
	def __init__(self,value,label='',symbol='',limit=()):
	self.value = value
	self.limit = limit
	self.label = label
	self.symbol = symbol



	################ Plot and tune ####################

	def tuner(S,fixedParameters,tunableParameters,freqs):

	""" Interactive plot and tuner """

	y0 = np.array( [p.value for p in fixedParameters] )
	x_opt = np.array( [p.value for p in tunableParameters] )

	figure = plt.figure(figsize=(16,9))

	ax0 = plt.subplot(221)
	ax1 = plt.subplot(223)
	ax2 = plt.subplot(122)
	plt.tight_layout()

	plt.subplots_adjust(bottom=0.08*len(tunableParameters))
	plt.subplots_adjust(left=0.10)
	plt.subplots_adjust(right=1.00)


	Gamma = [S(y0,x_opt,fi) for fi in freqs]
	S11_dB = np.array([10 * np.log10(np.abs(Gamma_i[0,0])**2) for Gamma_i in Gamma])
	S12_dB = np.array([10 * np.log10(np.abs(Gamma_i[0,1])**2) for Gamma_i in Gamma])
	S22_dB = np.array([10 * np.log10(np.abs(Gamma_i[1,1])**2) for Gamma_i in Gamma])
	phase11 = np.array( [ np.angle(Gamma_i[0,0]) for Gamma_i in Gamma] )
	phase12 = np.array( [ np.angle(Gamma_i[0,1]) for Gamma_i in Gamma] )

	##### S-parameters plot ######
	ax0.axis([800, 1000, -60, 3])
	ax1.axis([800, 1000, -180, 180])

	l1, = ax0.plot(freqs,S11_dB, lw=2, color='red', label=" $\|S_{11}\|$ ")
	l2, = ax0.plot(freqs,S12_dB, lw=2, color='blue', label=" $\|S_{12}\|$ ")
	#l3, = ax0.plot(freqs,S22_dB, lw=2, color='green', label=" $\|S_{22}\|$ ")
	l4, = ax1.plot(freqs,phase11 * 180/np.pi, lw=2, color='red', label=" $\phi(S_{11})$ ")
	l5, = ax1.plot(freqs,phase12 * 180/np.pi, lw=2, color='blue', label=" $\phi(S_{12})$ ")

	ax0.set_ylabel("S-parameters [dB]")
	ax0.legend(loc="lower right")
	ax0.grid()

	ax1.set_xlabel("Frequency, $f$ [MHz]")
	ax1.set_ylabel("Phase [deg]")
	ax1.legend(loc="lower right")
	ax1.grid()


	########## Smith plot ##########

	numPoints = 20

	ax2.axis([-1,1,-1,1])

	ax2.set_xlabel(" $Re(\Gamma)$ ")
	ax2.set_ylabel(" $Im(\Gamma)$ ")
	ax2.set_aspect(1)

	# boundary, gamma=0.5,0.1
	theta = np.arange(0, 2.5*np.pi, 0.1)
	ax2.plot( np.cos(theta), np.sin(theta) , 'black', linewidth=3 )
	ax2.plot( 0.5np.cos(theta), 0.5np.sin(theta) , 'blue', linewidth=0.5 )
	ax2.plot( 0.1np.cos(theta), 0.1np.sin(theta) , 'green', linewidth=0.5 )
	ax2.plot( 0.01np.cos(theta), 0.01np.sin(theta) , 'red', linewidth=0.5 )

	# R, G =1
	x = np.logspace(-5,5,num=300)
	x = np.concatenate((x,-x[::-1]))
	z = 1+1j*x
	g = (z-1)/(z+1)
	ax2.plot( g.real, g.imag, 'black')
	y = 1-1j*x
	g = (1-y)/(1+y)
	ax2.plot( g.real, g.imag, 'black', linestyle=':')

	# \|X, B\|=1
	r = np.logspace(-5,5,num=300)
	z = r+1j*1.
	g = (z-1)/(z+1)
	ax2.plot( g.real, g.imag, 'black' )

	z = r-1j*1.
	g = (z-1)/(z+1)
	ax2.plot( g.real, g.imag, 'black' )

	y = (r+1j*1.)
	g = (1-y)/(1+y)
	ax2.plot( g.real, g.imag, 'black', linestyle=':' )

	y = (r-1j*1.)
	g = (1-y)/(1+y)
	ax2.plot( g.real, g.imag, 'black', linestyle=':' )


	interval = int(len(freqs)/numPoints)

	f1 = freqs[::interval]
	f1.sort()
	Gamma11 = np.array([Gamma_i[0,0] for Gamma_i in Gamma])
	G1 = Gamma11[::interval]
	lSmith = ax2.scatter( G1.real, G1.imag, c=(f1[::-1]), cmap = plt.cm.spectral)

	###### Widgets #########

	axcolor = 'lightgoldenrodyellow'
	ax = []
	slider = []
	for i in range(len(tunableParameters)):
	p = tunableParameters[i]
	ax.append( plt.axes([0.25, 0.1+0.05*i, 0.65, 0.03], axisbg=axcolor) )
	slider.append( Slider(ax[i],p.label + ', '+p.symbol,p.limit[0],p.limit[1], valinit=x_opt[i]) )

	def update(val):
	x_new = [s.val for s in slider]
	Gamma = [S(y0,x_new,fi) for fi in freqs]
	S11_dB = np.array([10 * np.log10(np.abs(Gamma_i[0,0])**2) for Gamma_i in Gamma])
	S12_dB = np.array([10 * np.log10(np.abs(Gamma_i[0,1])**2) for Gamma_i in Gamma])
	S22_dB = np.array([10 * np.log10(np.abs(Gamma_i[1,1])**2) for Gamma_i in Gamma])
	phase11 = np.array( [ np.angle(Gamma_i[0,0]) for Gamma_i in Gamma] )
	phase12 = np.array( [ np.angle(Gamma_i[0,1]) for Gamma_i in Gamma] )

	l1.set_ydata(S11_dB)
	l2.set_ydata(S12_dB)
	#l3.set_ydata(S22_dB)
	l4.set_ydata(phase11 * 180/np.pi)
	l5.set_ydata(phase12 * 180/np.pi)

	Gamma11 = np.array([Gamma_i[0,0] for Gamma_i in Gamma])
	G1 = Gamma11[::interval]
	lSmith.set_offsets( [(Gi.real,Gi.imag) for Gi in G1] )
	plt.draw()

	for s in slider: s.on_changed(update)

	resetax = plt.axes([0.8, 0.025, 0.1, 0.04])
	button = Button(resetax, 'Optimize', color=axcolor, hovercolor='0.975')
	def reset(event):
	for s in slider: s.reset()
	button.on_clicked(reset)
	return figure, button

	def sensitivityPlot(fixedParameters, tunableParameters, f0, i=0,j=1):
	h = 1e-6
	allParameters = fixedParameters + tunableParameters
	y0,x0 = np.array([p.value for p in fixedParameters]), np.array([p.value for p in tunableParameters])
	Nf, Nt = len(fixedParameters), len(tunableParameters)
	If, It = np.eye(Nf), np.eye(Nt)

	DeltaPhi_f,DeltaPhi_t = [],[]
	Phi0 = np.angle( S(y0,x0,f0)[i,j] )

	for k in range( Nf ): DeltaPhi_f.append( ( np.angle(S(y0+h*If[k],x0,f0)[i,j]) - Phi0 )/h )
	for k in range( Nt ): DeltaPhi_t.append( ( np.angle(S(y0,x0+h*It[k],f0)[i,j]) - Phi0 )/h )
	DeltaPhi = np.array(DeltaPhi_f + DeltaPhi_t)

	DeltaG_f,DeltaG_t = [],[]
	G0 = np.abs( S(y0,x0,f0)[i,j] )
	for k in range( Nf ): DeltaG_f.append( ( np.abs(S(y0+h*If[k],x0,f0)[i,j]) - G0 )/h )
	for k in range( Nt ): DeltaG_t.append( ( np.abs(S(y0,x0+h*It[k],f0)[i,j]) - G0 )/h )
	DeltaG = np.array(DeltaG_f + DeltaG_t)

	sensitivityFigure = plt.figure(figsize=(14,6))
	ax1 = plt.subplot(121)
	ax2 = plt.subplot(122)

	plt.subplots_adjust(left=0.07)
	plt.subplots_adjust(right=0.98)

	index = np.arange(Nf+Nt)
	width = 0.35

	ax1.set_ylabel(" $S_{12}$ amplitude sensitivity @ "+`f0`+ " MHz [ $10^{-6}$ /mm]")
	ax1.set_xticks(index+width/2.)
	ax1.set_xticklabels([p.symbol for p in allParameters],rotation=17,fontsize=16)
	ax1.bar(index,DeltaG*1e6/10,width,color='b')

	ax2.set_ylabel(" $S_{12}$ phase sensitivity @ "+`f0`+ " MHz [mrad/mm]")
	ax2.set_xticks(index+width/2.)
	ax2.set_xticklabels([p.symbol for p in allParameters],rotation=17,fontsize=16)
	ax2.bar(index,DeltaPhi*1000/10.,width,color='r')

	return sensitivityFigure

	####################### end of setup ################################


	############# Problem definition #########################

	a = Parameter(24.,'width (cm)',' $a$ ',(23,25))
	#l_probe = Parameter(3.,'probe length (cm)',' $l_{probe}$ ',(3,3))
	b = Parameter(3.,'height (cm)',' $b$ ',(3,3))
	probeDiameter = Parameter(0.62,'probe dia. (cm)',' $d_{probe}$ ',(0.4,0.7))

	c = Parameter(14,'backshort distance (cm)', ' $c$ ',(10,22))
	l_seriesStub = Parameter(0.03,'series stub length (cm)',' $l_{seriesStub}$ ',(0.01,0.8))
	z_shuntStub = Parameter(3,'stub position (cm)',' $z_{shunt}$ ',(0.01,8.25))
	l_shuntStub = Parameter(3,'stub length (cm)',' $l_{shunt}$ ',(0.01,8.25))
	l_WG = Parameter(27,'distance to WG symmetry plane (cm)',' $l_{WG}$ ',(10,45))

	fixedParams = [a,b,probeDiameter]
	tunableParams = [c,l_seriesStub,z_shuntStub,l_shuntStub,l_WG]

	y0 = [p.value for p in fixedParams]
	x0 = [p.value for p in tunableParams]
	limits = [p.limit for p in tunableParams]


	################ Optimization #####################

	res = optimize.minimize(partial(targetFunction,y0),x0,method='SLSQP',jac=False,bounds=limits)
	# available methods: 'TNC' (truncated Newton), 'COBYLA' (constrained optimization by lin. approx.), 'SLSQP' (sequential least squares programming)
	# http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html#scipy.optimize.minimize
	for p in tunableParams: p.value = res.x[tunableParams.index(p)]
	freqs = np.arange(800, 1000, 1)

	print res
	launcher = makeLauncher(y0,res.x,910,verbose=True)
	figure,button = tuner(S,fixedParams,tunableParams,freqs)
	sense = sensitivityPlot(fixedParams,tunableParams,910,0,1)
	plt.show()

view raw twiddle-v7.py hosted with ❤ by GitHub

It's a lot of fun watching things move around on a Smith chart as you twiddle parameters.

I've been wondering why the curves always curve clockwise (as frequencies range from low to high). For purely reactive load, it makes sense in terms of Foster's theorem. But what happens with loads that have a resistive component?

A sketch of a proof/reason:

$z(f)$ , with passive resistances, lives on the right half plane of

$r + jx$ . The y-axis,

$r=0$ , ends up as the

$\Gamma=1$ circle. The vertical line

$r=1$ gets mapped to the resistance circle. Everything to the right of the line

$r=1$ ends up on the inside of the resistance circle.

The map

$\Gamma(z(f)) = (z-1)/(z+1)$ is a Mobius transformation of the form

$\Big(\frac{az+b}{cz+d}\Big)$ with determinant

$ad-bc > 0$ , so it should preserve handedness. Foster's theorem implies that the curve

$z(f)$ always moves upwards monotonically, therefore the curve

$\Gamma(f)$ always circulates clockwise.

p.s.

$\textrm{LaTeX}$ -t in here is thanks to this tip about MathJax, and code snippets thanks to Gist.

Pythagon

Sunday, April 14, 2013

Smith chart tuner