This experiment uses data extracted from a vacuum tube datasheet and scipy.optimize to calculate the Child-Langmuir parameters used for circuit simulation.
$$I_a = K (V_{gk} + D_aV_{ak})^{3/2}$$We are going to use curve fitting to determine $$K \text{ and } D_A$$
Then we can use Leach's triode SPICE model.
In [14]:
import scipy
from scipy.optimize import curve_fit
import numpy as np
import matplotlib.pyplot as plt
import math
Starting with the Philips ECC83 data sheet, create a PNG of the
Create 9 curves then use 'curve point tool' to add points to each curve
Change export options, "Raw Xs and Ys" and "One curve on each line", otherwise engauge will do some interrupting of your points
In [15]:
%cat data/ecc83-philips-1954.csv
Need to create scipy array like this
x = scipy.array([[0,150],[-3.0,300],[-1.0,200],[-2.0,250]])
y = scipy.array([3.2e-3,0.6e-3,2.2e-3,1.2e-3])
from the extracted curves
In [16]:
fname = "data/ecc83-philips-1954.csv"
f = open(fname,'r').readlines()
deltaVgk = -0.5
VgkVak = []
Iak = []
Vaks = []
for l in f:
l = l.strip()
if len(l): # skip blank lines
if l[0] == 'x':
vn = float(l.split("Curve")[1]) - 1.0
Vgk = vn * deltaVgk
continue
else:
(Vak,i) = l.split(',')
VgkVak.append([float(Vgk),float(Vak)])
Iak.append(float(i))
Vaks.append(float(Vak))
x = scipy.array(VgkVak)
y = scipy.array(Iak)
Vaks = scipy.array(Vaks)
In [17]:
%matplotlib inline
exp = 1.5
def func(x,K,Da):
rv = []
for VV in x:
Vgk = VV[0]
Vak = VV[1]
t = Vgk + Da * Vak
if t > 0:
Ia = K * t**exp
else:
Ia = 0
# print "IaCalc",Vgk,Vak,t,Ia
rv.append(Ia)
return rv
popt, pcov = curve_fit(func, x, y,p0=[0.001,0.01])
#print popt,pcov
(K,Da) = popt
print "K=%.8f Da=%.8f"%(K,Da)
In [18]:
def IaCalc(Vgk,Vak):
t = Vgk + Da * Vak
if t > 0:
Ia = K * t**exp
else:
Ia = 0
# print "IaCalc",Vgk,Vak,t,Ia
return Ia
Vgk = np.linspace(0,-5,11)
Vak = np.linspace(0,400,21)
vIaCalc = np.vectorize(IaCalc,otypes=[np.float])
Ia = vIaCalc(Vgk[:,None],Vak[None,:])
plt.figure(figsize=(12,6))
for i in range(len(Vgk)):
plt.plot(Vak,Ia[i],label=Vgk[i])
plt.scatter(Vaks,y,marker="+")
plt.legend(loc='upper left')
plt.suptitle('ECC83 Child-Langmuir Curve-Fit K/Da Model (Philips 1954)', fontsize=14, fontweight='bold')
plt.grid()
plt.ylim((0,0.005))
plt.xlim((0,400))
#plt.savefig("ecc83-philips-1954-curvefit.png")
plt.show()
Trying the Koren's triode phenomenological model.
$$E_1 = \frac{E_P}{k_P} log\left(1 + exp^{k_P (\frac{1}{u} + \frac{E_G}{\sqrt{k_{VB} + {E_P}^2}})}\right)$$$$I_P = \frac{{E_1}^X}{k_{G1}} \left(1+sgn(E_1)\right)$$Need to fit $X, k_{G1}, k_P, k_{VB}$
SPICE Model
see http://www.normankoren.com/Audio/Tubemodspice_article_2.html#Appendix_A
.SUBCKT 12AX7 1 2 3 ; P G C; NEW MODEL (TRIODE)
+ PARAMS: MU=100 EX=1.4 KG1=1060 KP=600 KVB=300 RGI=2000
E1 7 0 VALUE={V(1,3)/KP*LOG(1+EXP(KP*(1/MU+V(2,3)/SQRT(KVB+V(1,3)*V(1,3)))))}
G1 1 3 VALUE={(PWR(V(7),EX)+PWRS(V(7),EX))/KG1}
In [19]:
mu = 100
def sgn(val):
if val >= 0:
return 1
if val < 0:
return -1
def funcKoren(x,X,kG1,kP,kVB):
rv = []
for VV in x:
EG = VV[0]
EP = VV[1]
E1 = (EP/kP) * math.log(1 + math.exp(kP*(1/mu + EG / math.sqrt(kVB + EP*EP))))
if E1 > 0:
IP = (math.pow(E1,X)/kG1)*(1 + sgn(E1))
else:
IP = 0
rv.append(IP)
return rv
popt, pcov = curve_fit(funcKoren,x,y,p0=[1.3,1000,610,305])
#print popt,pcov
(X,kG1,kP,kVB) = popt
print "X=%.8f kG1=%.8f kP=%.8f kVB=%.8f"%(X,kG1,kP,kVB)
# koren's values 12AX7 mu=100 X=1.4 kG1=1060 kP=600 kVB=300
In [20]:
'''
X=1.4
kG1=1060
kP=600
kVB=300
mu=100
'''
def IaCalcKoren(Vgk,Vak):
global X,kG1,kP,kVB,mu
E1 = (Vak/kP) * math.log(1 + math.exp(kP*(1/mu + Vgk / math.sqrt(kVB + Vak*Vak))))
Ia = (math.pow(E1,X)/kG1)*(1 + sgn(E1))
return Ia
Vgk = np.linspace(0,-5,11)
Vak = np.linspace(0,400,21)
vIaCalcKoren = np.vectorize(IaCalcKoren,otypes=[np.float])
Ia = vIaCalcKoren(Vgk[:,None],Vak[None,:])
plt.figure(figsize=(12,6))
for i in range(len(Vgk)):
plt.plot(Vak,Ia[i],label=Vgk[i])
plt.scatter(Vaks,y,marker="+")
plt.legend(loc='upper left')
plt.suptitle('ECC83 Child-Langmuir Curve-Fit Koren Model (Philips 1954)', fontsize=14, fontweight='bold')
plt.grid()
plt.ylim((0,0.005))
plt.xlim((0,400))
#plt.savefig("ecc83-philips-1954-curvefit.png")
plt.show()
In [ ]: