2910.sagews

var ('F,D,Vd,ka_high,ka_low,ke,t,_a,_b')

class Bateman:
    def __init__(self, ka=2, ke=0.2, F=1, D=50, Vd=7):
        self.ka = ka
        self.ke = ke
        self.F = F
        self.D = D
        self.Vd = Vd
    def __call__(self, t):
        return self.F*self.D/self.Vd * (e**(-self.ke*t) - e**(-self.ka*t)) * self.ka / (self.ka - self.ke)
    def getplot(self, min=0, max=24, color='blue', mark_tmax=False):
        p = plot(self(t), min, max, color=color)
        if mark_tmax:
            self.calculate_cmax()
            p += point((self.tmax, self.cmax), color=color)
        return p
    def getlnplot(self, min=0, max=24, color='blue'):
        return plot(ln(self(t)), min, max, color=color)
    def calculate_tmax(self):
        self.tmax = ln(self.ka/self.ke)/(self.ka-self.ke)
        return self.tmax
    def calculate_cmax(self):
        self.cmax = self(self.calculate_tmax())
        return self.cmax
    def f(self, x, _a, _b):
        return -_a*x+_b
    def getshaveplot(self, min=0, max=24, color='blue'):
        self.calculate_tmax()
        data_1 = [(t, ln(self(t))) for t in xsrange(10*self.tmax, 11*self.tmax, self.tmax/5)]
        #print data_1
        fit_1 = find_fit(data_1, self.f, parameters = [_a, _b], variables = [x], solution_dict = True)
        #print fit_1
        _ke = fit_1[_a]
        _keb = fit_1[_b]
        data_2 = [(t, ln(e**(-_ke*t+_keb) - self(t))) for t in xsrange(0, self.tmax, self.tmax/5)]
        #print data_2
        fit_2 = find_fit(data_2, self.f, parameters = [_a, _b], variables = [x], solution_dict = True)
        _ka = fit_2[_a]
        _kab = fit_2[_b]
        #print fit_2
        try:
            _kamax = find_root(-_ka*x+_kab, 0, max)
        except Exception as a:
            _kamax = max/5
        if _kab < 0:
            _kaby = -_ka * _kamax + _kab
            _kax = _kamax
        else:
            _kaby = _kab / 2
            _kax = _kamax / 2 + max / 50
        p = plot(-_ke*x+_keb, x, min, max, color=color, linestyle='-.')
        p += plot(-_ka*x+_kab, t, min, _kamax , color=color, linestyle='-.')
        p += self.getlnplot(min, max, color)
        p += text("ka: " + str(n(_ka, digits=2)), (_kax, _kaby), color=color, horizontal_alignment='left')
        p += text("ke: " + str(n(_ke, digits=2)), (max/2, _keb/2), color=color, horizontal_alignment='left')
        return p
        

normal = Bateman()
flipflop = Bateman(ka=0.05)

p = normal.getplot(max=24)
p += flipflop.getplot(max=24, color='red')
p.show()

p = normal.getshaveplot(max=12)
p.show()
p = flipflop.getshaveplot(max=12, color='red')
p.show()