#
# Gabbiani & Cox,  Mathematics for Neuroscientists
# 
# Solve the Passive Membrane Equation subject to a
# current pulse, via forward Euler
#
#  usage   feps(dt,Tfin)
#
#  e.g.    feps(0.01,40)
#

from scipy import *
from matplotlib.pyplot import *

close('all')

def feps(dt,Tfin):

    VCl = -68	    # mV
    A = 4*pi*1e-6	# cm^2 patch area
    Cm = 1         # micro F/cm^2
    gCl = 0.3      # mS/cm^2
    tau = Cm/gCl   # ms
    
    Nt = round(1+Tfin/dt)                  # number of time steps
    v = zeros(Nt)  
    t = zeros(Nt) 
    
    v[0] = VCl                 # initialize v
    
    for j in range(1,Nt):
    
       t[j] = j*dt
    
       Istim = (t[j]>2)*(t[j]<22)*1e-5    # 10 pA 20 ms pulse
    
       v[j] = v[j-1]*(1-dt/tau) + dt*(VCl/tau + Istim/A/Cm)  # forward euler
    
    plot(t,v,linewidth=2)
    xlabel('t  (ms)',fontsize=16)
    ylabel('V  (mV)',fontsize=16)