%
% iaf cell of D&A with gsra (spike rate adaptation)
%
%  taum V' = EL - V + Rm*I - rm*gsra*(V-EK)
%
%  tausra gsra' = -gsra
%
% where I is clamped to Ie after 1 ms
%
%  usage   iaf2(Ie,Vth,dt,T)
%
%  where dt is the time step and T is the final time
%
%  example   iaf2(10,-50,.01,40)
%

function iaf2(Ie,Vth,dt,T)

EL = -65;
Vreset = -65;
taum = 10;
rm = 1;		% M Ohm/ mm^2
A = .1;         % mm^2;
Rm = rm/A;	%M Ohm
Ie = 10;		%n A
EK = -70;
tausra = 10;
ginc = 50*0.06/rm;

s = taum/dt;
st = tausra/dt;
v = EL;
vshow(1) = v;
gsra(1) = 0;
t = 0;
j = 1;

while t < T

    j = j + 1;
    t = (j-1)*dt;
    I = 0;
    if t > 1,
       I = Ie;
    end

    if v < Vth

       gsra(j) = st*gsra(j-1)/(1+st);
 
       v = (s*v+EL+Rm*I+rm*gsra(j)*EK)/(s+1+rm*gsra(j));

       vshow(j) = v;

    else

       v = Vreset;

       gsra(j) = gsra(j-1) + ginc;

       vshow(j) = 10;

    end

end

tim = 0:dt:(j-1)*dt;
subplot(1,2,1)
plot(tim,vshow)
xlabel('t (ms)','fontsize',16)
ylabel('V (mV)','fontsize',16)
subplot(1,2,2)
plot(tim,gsra)
xlabel('t (ms)','fontsize',16)
ylabel('g_{sra} (mS/cm^2)','fontsize',16)

head = ['iaf2(' num2str(Ie) ',' num2str(Vth) ','];
head = [head num2str(dt) ',' num2str(T) ')'];
title(head,'fontsize',16)