%
% heasper.m
%
% hybrid euler on the active sphere with periodic input
%
% usage:  heasper(dt,Tfin,I0)
%

function heas(dt,Tfin,I0)

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

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

B = [-0.6773   76.3634    2.2602 -27.6994
    0.0264   -4.2236         0         0
   -0.0041         0   -0.1174         0
    0.0028         0         0   -0.1832];

B = eye(4)-dt*B;

f = zeros(4,1);

figure(1)

for k=1:25,

j = 2;
t(1) = 0;
v(1) = 0;
qv(1) = 0;
pv(1) = 0;
y = zeros(4,1);
n = an(0)/(an(0)+bn(0));  % 0.3177;
m = am(0)/(am(0)+bm(0));  % 0.0529;
h = ah(0)/(ah(0)+bh(0));  % 0.5961;
vCl = -(GNa*m^3*h*vNa + GK*n^4*vK)/GCl

omega(k) = 4*k/1000;

while j*dt < Tfin,

      t(j) = j*dt;

      Istim = I0*sin(2*pi*t(j)*omega(k))*(1e-6)/A;

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

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

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

      v(j) = top/bot; 

      f(1) = Istim;

      y = B\(y+dt*f);

      qv(j) = y(1);

      pv(j) = (pv(j-1)+dt*Istim)/(1+dt*GCl);

      j = j + 1;

end

plot(t,v)
hold on
plot(t,qv,'r')
plot(t,pv,'k')
drawnow
hold off
disp('pause'); pause

vamp(k) = max(abs(v(ceil(end/2):end)));
qvamp(k) = max(abs(qv(ceil(end/2):end)));
pvamp(k) = max(abs(pv(ceil(end/2):end)));

%hold on

end

%hold off

%figure(2)
plot(omega*1000,vamp/(I0*1e-6)/1000,omega*1000,qvamp/(I0*1e-6)/1000)
hold on
plot(omega*1000,pvamp/(I0*1e-6)/1000,'--')
xlabel('\omega  (Hz)','fontsize',16)
ylabel('R_{in}  (M\Omega)','fontsize',16)

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