class Excitable():
    def __init__(self):
        self.timestep=0.01
        self.maxtime=160
        self.meshpoints=int(self.maxtime/self.timestep)
        self.v=self.meshpoints*[0]
        self.w=self.meshpoints*[0]
        self.wdisp=self.meshpoints*[0]
        self.Iext=self.meshpoints*[0]
        self.time=self.meshpoints*[0]
        self.damp=0.0
        self.a=0.08;self.b=0.7;self.c=0.8;self.tau=1.0
    def dvdtstep(self,v,w,Iext):
        return self.timestep*(Iext+(1-self.damp)*v-v**3-w)
    def dvdt(self,point):
        return self.dvdtstep(self.v[point],self.w[point],self.Iext[point])
    def dwdtstep(self,v,w,Iext):
        return self.timestep*self.a*(v+self.b-self.c*w)/self.tau
    def dwdt(self,point):
        return self.dwdtstep(self.v[point],self.w[point],self.Iext[point])
    def StepForward(self,point):
        if point<7000: self.Iext[point]=point*0.0001
        self.time[point+1]=self.time[point]+self.timestep
        self.v[point+1]=self.v[point]+self.dvdt(point)
        self.w[point+1]=self.w[point]+self.dwdt(point)
        self.wdisp[point+1]=(self.w[point+1]+0.215)
        
if __name__ == "__main__":
    import matplotlib
    matplotlib.use('MacOSX')
    from pylab import *
    excited=Excitable()
    print 'Timestep:',excited.timestep,'  Max time:',excited.maxtime,'   Meshpoints:',excited.meshpoints
    for i in range(excited.meshpoints-1):
        excited.StepForward(i)
    plot(excited.time,excited.v,linewidth=1.0)
    plot(excited.time,excited.w,linewidth=1.0)
    plot(excited.time,excited.wdisp,linewidth=1.0)
#    plot(excited.v,excited.w,linewidth=1.0)
    xlabel('time (s)')
    ylabel('ODE')
    title('Excitable system')
    grid(True)
    show()