Skip to content

Instantly share code, notes, and snippets.

@amarvutha

amarvutha/twiddle-v7.py

Last active December 16, 2015 05:19
Show Gist options
  • Save amarvutha/5383239 to your computer and use it in GitHub Desktop.
Save amarvutha/5383239 to your computer and use it in GitHub Desktop.
Download ZIP
Microwave match tuner, with interactive Smith chart plot
Raw
twiddle-v7.py
# 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(-1j*beta*l)],[np.exp(-1j*beta*l),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-Gamma**2)],[np.sqrt(1-Gamma**2),-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(pi*d/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 = 2*pi*f/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/(2*pi)) * (pi*probeDiameter/a)**2 /( 1+ 11./24 * (pi*d/a)**2 ) # 5.11 (4b)
s = 0
for n in range(3,100,2): s+= (1/np.sqrt(n**2 - (beta0*a/pi)**2) - 1./n)
S0 = np.log(4*a/(pi*probeDiameter)) - 2. + 2*s
X_a = 0.5*X_b + Z_eq * beta_WG*a/(4*pi) * ( S0 -(beta0*probeDiameter/4)**2 ) # 5.11 (3)
#X_a = 0.55 * Z_eq * (a*beta_WG/pi) /np.sin(pi*d/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 = alpha*Z_eq * np.tan(beta_WG * c) # prescription from Slater
# launcher
Z_12 = -1j*X_b + 1j*X_backshort
Z_s1 = 1j*X_a + 1j*X_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 = (alpha*Z_eq - Z0)/(alpha*Z_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 = Z0*np.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/(2*pi*f*1e6) *1e9,"nH"
print "L_seriesStub = ",X_series/(2*pi*f*1e6) *1e9,"nH"
#print "L_shuntStub = ",X_shuntStub/(2*pi*f*1e6) *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.5*np.cos(theta), 0.5*np.sin(theta) , 'blue', linewidth=0.5 )
ax2.plot( 0.1*np.cos(theta), 0.1*np.sin(theta) , 'green', linewidth=0.5 )
ax2.plot( 0.01*np.cos(theta), 0.01*np.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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment