#!/usr/bin/python
"""
Optimize parameters for the triode model developed by Norman Koren
http://www.normankoren.com/Audio/Tube_params.html

Triode equations
Parameters:

ep: plate voltage 
ip: plate current 
eg: grid voltage 

m         - amplification factor (mu). 
kp        - emperical 
kvb       - emperical 
ex        - 3/2 from the Child-Langmuir law. Start with 1.5
kgi       - emperical 

"""
from numpy import sqrt,exp, log, sign, arange
from pylab import plot, show, ylim
from scipy.optimize import minimize

# file with measured data 
filename = ""

# function for the plate current
def ip(mu, ex, kg1, kp, kvb, eg, ep):
    e1 = ep / kp * log(1.0 + exp( kp*( (1.0/mu) + (eg / sqrt(kvb + pow(ep,2.0)) ))))
    return  pow(e1, ex) / kg1 * (1.0 + sign(e1))

# parse measured data 
Vp, Vg, Idata = [], [], []
for line in open(filename):
    if line.count("#"):
        continue
    point =  line.split()
    if (len(point) == 3):
        # only include measured data for plate currents under 200 mA in this case.
        if (float(point[2]) > 200.0):
            continue
        else:
            Vg.append(float(point[0]))
            Vp.append(float(point[1]))
            Idata.append(float(point[2]) * 1e-3)

# initial values 
mu_,ex_,kg1_,kp_,kvb_ = 3.0, 1.0, 1000.0, 10.0, 1.0
x0 = [mu_,ex_,kg1_,kp_,kvb_]

# function to minimize
def cost_function(x):
    mu, ex, kg1, kp, kvb = x[0],x[1],x[2],x[3],x[4]
    err = 0
    for idata, ep, eg in zip(Idata, Vp, Vg):
        err += pow( (idata - ip(mu, ex, kg1, kp, kvb, eg , ep)), 2)
    return err

# minimize cost function function
res = minimize(cost_function, x0, method='Nelder-Mead')
print res.x
print cost_function(res.x)

# plot model  
mu, ex, kg1, kp, kvb = res.x[0],res.x[1],res.x[2],res.x[3],res.x[4]
for eg in range(0, -70, -10):
    ipp = []
    for ep in range(0, 500):
        ipp.append(ip(mu, ex, kg1, kp, kvb, eg, ep))
    plot(ipp)

# plot measured data
plot(Vp,Idata,'o')
ylim(0,0.1)

show()