%
% heas.m    Steve Cox 1/11/03
%
% Hybrid Euler on the Active Sphere
%
% usage:  heas(dt,Tfin,I0)   
%
% where    dt   = timestep (ms)
%          Tfin = duration of simulation (ms)
%          I0   = amplitude of current step (pA)
%
% e.g.,   heas(.01,40,20)
%

function heas(dt,Tfin,I0)

A = 4*pi*(1e-6);% cm^2
vK = -6;	% mV
vNa = 127;	% mV
vCl = 2.8417; 	% mV

Cm = 1;		% micro F/cm^2
GK =  36;	% mS/(cm)^2
GNa = 120;	% mS/(cm)^2
GCl = 0.3;	% mS/(cm)^2

N = ceil(Tfin/dt);

t = zeros(1,N);   % allocate space for long vectors
v = t;
n = t;
m = t;
h = t;
ICl = t;
IK = t;
INa = t;

n(1) = an(0)/(an(0)+bn(0));  
m(1) = am(0)/(am(0)+bm(0)); 
h(1) = ah(0)/(ah(0)+bh(0));

ICl(1) = GCl*(v(1)-vCl);
IK(1) = GK*n(1)^4*(v(1)-vK);
INa(1) = GNa*m(1)^3*h(1)*(v(1)-vNa);

for j=2:N,

    t(j) = (j-1)*dt;

    Istim = I0*(t(j)>2)*(t(j)<22)*(1e-6)/A;

    n(j) = ( n(j-1) + dt*an(v(j-1)) )/(1 + dt*(an(v(j-1))+bn(v(j-1))) );
    m(j) = ( m(j-1) + dt*am(v(j-1)) )/(1 + dt*(am(v(j-1))+bm(v(j-1))) );
    h(j) = ( h(j-1) + dt*ah(v(j-1)) )/(1 + dt*(ah(v(j-1))+bh(v(j-1))) );

    top = Cm*v(j-1)+dt*(GNa*m(j)^3*h(j)*vNa + GK*n(j)^4*vK + GCl*vCl + Istim);

    bot = Cm + dt*(GNa*m(j)^3*h(j) + GK*n(j)^4 + GCl);

    v(j) = top/bot; 

    ICl(j) = GCl*(v(j)-vCl);

    IK(j) = GK*n(j)^4*(v(j)-vK);

    INa(j) = GNa*m(j)^3*h(j)*(v(j)-vNa);
 
end

close all

figure(1)
subplot(2,2,1)
plot(t,v)
ylabel('v (mV)','fontsize',16)
subplot(2,2,2)
plot(t,n)
ylabel('n','fontsize',16)
subplot(2,2,3)
plot(t,m)
xlabel('t (ms)','fontsize',16)
ylabel('m','fontsize',16)
subplot(2,2,4)
plot(t,h)
ylabel('h','fontsize',16)
xlabel('t (ms)','fontsize',16)

figure(2)
plot(t,ICl,t,IK,t,INa)
xlabel('t (ms)','fontsize',16)
ylabel('I (\muA/cm^2)','fontsize',16)
legend('I_{Cl}','I_K','I_{Na}')

% rate functions from page 519 of HH
% (with polarity switched to agree with modern usage)

function val = an(v)
val = .01*(10-v)./(exp(1-v/10)-1);

function val = bn(v)
val = .125*exp(-v/80);

function val = am(v)
val = .1*(25-v)./(exp(2.5-v/10)-1);

function val = bm(v)
val = 4*exp(-v/18);

function val = ah(v)
val = 0.07*exp(-v/20);

function val = bh(v)
val = 1./(exp(3-v/10)+1);