%
% heas_sp.m    Steve Cox 
%
% Count HH spikes in time interval [0,Tfin]
%
% this a stripped down copy of heas.m, for it does not
% accumulate any info (other than spikes) and so does
% not plot anything
%
% usage:  heas_sp(dt,Tfin,I0,vth)   
%
% where    dt   = timestep (ms)
%          Tfin = duration of simulation (ms)
%          I0   = amplitude of current step (pA)
%          vth  = voltage threshold
%
% e.g.,   heas_sp(.01,1000,100,20)
%

function spcnt = heas_sp(dt,Tfin,I0,vth)

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

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

up = 0;
spcnt = 0;
Iapp = I0*(1e-6)/A;

while t < Tfin+2

    t = t + dt;

    Istim = Iapp*(t>2);

    ant = an(v);
    n = ( n + dt*ant )/(1 + dt*(ant+bn(v)) );
    amt = am(v);
    m = ( m + dt*amt )/(1 + dt*(amt+bm(v)) );
    aht = ah(v);
    h = ( h + dt*aht )/(1 + dt*(aht+bh(v)) );

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

    bot = Cm + dt*(iNa + iK + GCl);

    v = top/bot; 

    if v > vth
       up = 1;  	     % mark upward crossing
    else
       spcnt = spcnt + up;   % increment spcnt on downward crossing
       up = 0;
    end

end

return

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);