Pentode Modeling

  • Model Parameter Extraction
  • Model Parameter Verification

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_{g1k} + D_{g2}V_{g2k} + D_aV_{ak})^\frac{3}{2}$$

Now we're adding the van der Veen K modifier,

$$\alpha = \alpha_0\left(\frac{2}{\pi}\arctan \left(\frac{V_{ak}}{V_{g2k}}\right)\right)^\frac{1}{n}$$$$I_a = \alpha K (V_{g1k} + D_{g2}V_{g2k} + D_aV_{ak})^\frac{3}{2}$$$$I_a = \alpha_0\left(\frac{2}{\pi}\arctan\left(\frac{V_{ak}}{V_{g2k}}\right)\right)^\frac{1}{n} K \left(V_{g1k} + D_{g2}V_{g2k} + D_aV_{ak}\right)^\frac{3}{2}$$

We are going to use curve fitting to determine $$K, D_a, D_{g2},\alpha_0, \text{ and } n$$

Then we can use Leach's pentode SPICE model with the van der Veen modifier.


In [1]:
import scipy
from scipy.optimize import curve_fit
import numpy as np
import matplotlib.pyplot as plt
from math import pi,atan,log,pow,exp


/home/holla/anaconda2/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')

Starting with the Philips EL34 data sheet, create a PNG of the import this image into engauge

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 export a csv file


In [2]:
%cat data/el34-philips-1958-360V.csv


x,Curve1
9.66,0.17962
22.99,0.26382
41.49,0.3227
70.55,0.37863
116.61,0.42335
167.87,0.44981
208.47,0.45979

x,Curve2
10.69,0.14899
31.64,0.22789
55.38,0.27793
80.22,0.30912
107.97,0.3356
142.14,0.35266
183.36,0.36911
233.38,0.38262
270.44,0.38907

x,Curve3
10.54,0.11601
26.63,0.16252
49.13,0.19902
79.27,0.23256
112.89,0.25786
159.43,0.27962
208.86,0.29312
259.46,0.30309
315.33,0.31012
358.25,0.31244

x,Curve4
10.97,0.08244
35.85,0.12306
66.58,0.15896
102.56,0.18602
159.69,0.2113
217.95,0.22657
280.3,0.23653
338.51,0.24061
419.65,0.24643
489.03,0.24991

x,Curve5
10.89,0.06595
47.46,0.09301
83.43,0.11831
125.86,0.14183
195.91,0.16298
295.88,0.17409
379.38,0.18109
469.93,0.1875
644.55,0.20032

x,Curve6
12.53,0.04004
29.04,0.05121
49.61,0.05002
73.19,0.0659
115.62,0.08884
176.84,0.10822
243.3,0.11642
443.79,0.12922
740.72,0.15077

x,Curve7
10.73,0.03179
36.62,0.03707
66,0.03646
117.81,0.05468
194.3,0.0717
275.44,0.07693
381.26,0.08156
498.26,0.08913
745.19,0.10189

x,Curve8
11.26,0.01884
37.72,0.02235
74.77,0.02644
144.18,0.03817
231.2,0.04458
364.06,0.0486
549.83,0.05493
741.47,0.06126

x,Curve9
12.39,0.01059
35.33,0.01293
64.68,0.00525
117.06,0.01934
194.07,0.02282
319.29,0.02684
439.79,0.0291
584.99,0.0337
739.6,0.03829

Need to create scipy array like this

x = scipy.array( [[360, -0.0, 9.66], [360, -0.0, 22.99], ...

y = scipy.array( [0.17962, 0.26382, 0.3227, 0.37863, ...

Vaks = scipy.array( [9.66, 22.99, 41.49, 70.55, 116.61, ...

from the extracted curves


In [3]:
fname = "data/el34-philips-1958-360V.csv"
f = open(fname,'r').readlines()

deltaVgk = -4.0
n = 1.50


VgkVak = []
Iak    = []
Vaks   = []

f = open(fname,'r').readlines()

vg2k = 360

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:plt.xkcd()
            (Vak,i) = l.split(',')
        VgkVak.append([vg2k,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 [4]:
%matplotlib inline
def func(x,K,Da,Dg2,a0,n):
    rv = []
    for VV in x:
        Vg2k = VV[0]
        Vg1k = VV[1]
        Vak  = VV[2]
        t = Vg1k + Dg2 * Vg2k + Da * Vak

        if t > 0:
            a = a0 * ((2/pi) * atan(Vak/Vg2k))**(1/n)
            Ia = a * K * t**n
        else:
            Ia = 0
        # print "func",Vg2k,Vg1k,Vak,t,K,Da,Dg2,a0,n
        rv.append(Ia)

    return rv

popt, pcov = curve_fit(func, x, y,p0=[0.5,0.05,0.05,0.02,5])
#print popt,pcov

(K,Da,Dg2,a0,n) = popt
print "K   =",K
print "Da  =",Da
print "Dg2 =",Dg2
print "a0  =",a0
print "n   =",n


K   = 0.0055663327717
Da  = 9.93851673822e-06
Dg2 = 0.153150838142
a0  = 0.000192943868262
n   = 3.33116708071

In [5]:
Vg2k = x[0][0]

def IaCalc(Vg1k,Vak):
    t = Vg1k + Dg2 * Vg2k + Da * Vak
    if t > 0:
        a = a0 * ((2/pi) * atan(Vak/Vg2k))**(1/n)
        Ia = a * K * t**n
    else:
        Ia = 0
    # print "IaCalc",Vgk,Vak,t,Ia
    return Ia

Vgk = np.linspace(0,-32,9)
Vak = np.linspace(0,400,201)

vIaCalc = np.vectorize(IaCalc,otypes=[np.float])

Iavdv = vIaCalc(Vgk[:,None],Vak[None,:])

plt.figure(figsize=(14,6))

for i in range(len(Vgk)):
    plt.plot(Vak,Iavdv[i],label=Vgk[i])

plt.scatter(Vaks,y,marker="+")

plt.legend(loc='upper left')
plt.suptitle('EL34@%dV Child-Langmuir-Compton-VanDerVeen Curve-Fit K/Da/Dg2 Model (Philips 1949)'%Vg2k, fontsize=14, fontweight='bold')
plt.grid()
plt.ylim((0,0.5))
plt.xlim((0,400))
plt.show()


Trying the Koren's triode phenomenological model.

$$E_1 = \frac{E_{G2}}{k_P} log\left(1 + exp^{k_P (\frac{1}{u} + \frac{E_{G1}}{E_{G2}})}\right)$$$$I_P = \left(\frac{{E_1}^X}{k_{G1}}\right) \left(1+sgn(E_1)\right)atan\left(\frac{E_P}{k_{VB}}\right)$$

Need to fit $X, k_{G1}, k_P, k_{VB}$


In [6]:
mu = 11.0

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:
        EG2 = VV[0]
        EG1 = VV[1]
        EP  = VV[2]

        if kP < 0:
            kP = 0
        #print EG2,EG1,EP,kG1,kP,kVB,exp(kP*(1/mu + EG1/EG2))
        
        E1 = (EG2/kP) * log(1 + exp(kP*(1/mu + EG1/EG2))) 
        if E1 > 0:
            IP = (pow(E1,X)/kG1)*(1 + sgn(E1))*atan(EP/kVB)
        else:
            IP = 0
        rv.append(IP)

    return rv

popt, pcov = curve_fit(funcKoren,x,y,p0=[1.3,1000,40,20])
#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


X=1.68255534   kG1=2266.00960590   kP=44.79909890   kVB=23.92582971
SPICE model
see http://www.normankoren.com/Audio/Tubemodspice_article_2.html#Appendix_A
.SUBCKT 6550 1 2 3 4 ; P G1 C G2 (PENTODE) 
+ PARAMS: MU=7.9 EX=1.35 KG1=890 KG2=4200 KP=60 KVB=24 
E1 7 0 VALUE={V(4,3)/KP*LOG(1+EXP((1/MU+V(2,3)/V(4,3))*KP))} 
G1 1 3 VALUE={(PWR(V(7),EX)+PWRS(V(7),EX))/KG1*ATAN(V(1,3)/KVB)} 
G2 4 3 VALUE={(EXP(EX*(LOG((V(4,3)/MU)+V(2,3)))))/KG2} 

In [7]:
EG2 = x[0][0]

def IaCalcKoren(EG1,EP):
    global X,kG1,kP,kVB,mu
    E1 = (EG2/kP) * log(1 + exp(kP*(1/mu + EG1/EG2))) 
    if E1 > 0:
        IP = (pow(E1,X)/kG1)*(1 + sgn(E1))*atan(EP/kVB)
    else:
        IP = 0
    return IP

Vgk = np.linspace(0,-32,9)
Vak = np.linspace(0,400,201)

vIaCalcKoren = np.vectorize(IaCalcKoren,otypes=[np.float])

Iakoren = vIaCalcKoren(Vgk[:,None],Vak[None,:])

plt.figure(figsize=(14,6))

for i in range(len(Vgk)):
    plt.plot(Vak,Iakoren[i],label=Vgk[i])

plt.scatter(Vaks,y,marker="+")

plt.legend(loc='upper left')
plt.suptitle('EL34@%dV Child-Langmuir-Compton-Koren Curve-Fit Model (Philips 1949)'%Vg2k, fontsize=14, fontweight='bold')
plt.grid()
plt.ylim((0,0.5))
plt.xlim((0,400))
plt.show()



In [8]:
plt.figure(figsize=(14,6))

for i in range(len(Vgk)):
    plt.plot(Vak,Iavdv[i],label=Vgk[i],color='red')
    plt.plot(Vak,Iakoren[i],label=Vgk[i],color='blue')

plt.scatter(Vaks,y,marker="+")

plt.legend(loc='upper left')
plt.suptitle('EL34@%dV CLCVDV & CLCK Curve-Fit Model (Philips 1949)'%Vg2k, fontsize=14, fontweight='bold')
plt.grid()
plt.ylim((0,0.5))
plt.xlim((0,400))
plt.show()



In [ ]: