%
%  he_afib.m
%
%  Steve Cox, 1/20/03
%
% solve the active fiber problem via sparse hybrid euler
%
% usage:	he_afib(a,ell,N,dt,t1,t2,T,Iapp)
%
%		a = fiber radius (cm)
%		ell = fiber length (cm)
%             	N = number of patches
%               dt = timestep (ms)
%		t1 = start of current pulse (ms)
%		t2 = end of current pulse (ms)
%               T = stopping time (ms)
%		Iapp = size of current pulse (micro amps)
%               pinc = number of time steps between plots
%                      if pinc > 0 then 3D full
%                      if pinc < 0 then 2D movie
%
% examples:	he_afib(.0001,1,200,.01,1,2,20,.1)
%

function he_afib(a,ell,N,dt,t1,t2,T,Iapp)

C_M = 1;		% micro F / cm^2
R_2 = 2*0.034;		% k Ohm cm
G_i = a/2/R_2;		% scaled axial conductance
dx = ell/N;		% patch length
A = 2*pi*a*dx;		% patch surface area
x = dx/2:dx:ell-dx/2;	% vector of patch midpoints

vK = -6;        % mV
vNa = 127;      % mV
vCl = 3;        % mV

GK = 36;        % mS/(cm)^2
GNa = 40;      % mS/(cm)^2
GCl = 0.3;      % mS/(cm)^2

e = ones(N,1);

v = 0*e;		% initial conditions
n = 0.3177*e;
m = 0.0529*e;
h = 0.5961*e;

rhs = zeros(N,1);

Bd = 2*e;
Bd(1) = 1; Bd(N) = 1;
Bd = G_i*Bd/dx/dx + C_M/dt + GCl;
B = spdiags([-e 2*e -e], -1:1, N, N)/dx/dx;
B(1,1) = 1/dx/dx;
B(N,N) = 1/dx/dx;
B = G_i*B;

t = 0;
tcnt = 0;
vcnt = 1;

while (t <= T)

      t = t + dt;

      Istim = 0;
      if (t>t1 & t<t2)
         Istim = Iapp;
      end

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

      cNa = GNa * m.^3 .* h;
      cK = GK * n.^4;

      B(1:N+1:end) = (Bd + cNa + cK)';   % update diagonal of B

      rhs = v*C_M/dt + cNa*vNa + cK*vK + GCl*vCl;
      rhs(1) = rhs(1) + Istim/A;

      v = B \ rhs;

      tcnt = tcnt + 1;
   
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);
