%
% hybrid Euler on the active fork with soma
%

function hefork(dt,Tfin,t1,t2,I0,pinc)

a = 1e-4*[1 1 1]; % 1e-3*[1 1/2 2];
ell = 0.02*[1 1 1];  %[1 1 2];
dx = .01;
%dx = 0.2;
N = ell/dx;
A3 = 2*pi*a(3)*dx;
As = 4*pi*1e-6;
rho = A3/As;
R2 = 0.15;  %  0.034;
gL = 0.3;
Cm = 1;
tau = Cm/gL;
lam = a/(2*R2*gL)/dx^2;   % lammda^2
r = a/a(3);
Hd = [2*lam(1)*ones(1,N(1)) 2*lam(2)*ones(1,N(2)) 2*lam(3)*ones(1,N(3)+1)];
Hd(1) = lam(1);
Hd(N(1)+1) = lam(2);
Hd(N(1)+N(2)+1) =  lam*r';
Hd(end) = rho*lam(3);
Hlen = length(Hd);

Hu = [-lam(1)*ones(1,N(1)-1) 0 -lam(2)*ones(1,N(2)) -lam(3)*ones(1,N(3))];
Hl = [-lam(1)*ones(1,N(1)-1) 0 -lam(2)*ones(1,N(2)-1) -r(2)*lam(2) -lam(3)*ones(1,N(3))];
Hl(end) = rho*Hl(end);

H = spdiags( [[Hl 0]' Hd' [0 Hu]'], -1:1, Hlen, Hlen);

H(N(1)+N(2)+1,N(1)) = -r(1)*lam(1);
H(N(1),N(1)+N(2)+1) = -lam(1);

%I = speye(Hlen);

%H = dt*H/tau;
dH = diag(H);

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

GK = 36;        % mS/(cm)^2
GNa = 120;      % mS/(cm)^2
GL = 0.3;      % mS/(cm)^2
GK = GK/GL;
GNa = GNa/GL;

%Bb = I+(I+H)*(dt/tau); 
%[L,U] = lu(Bb); 

%keyboard

x3 = 0:dx:ell(3);
x1 = ell(3):dx:ell(3)+ell(1);
x2 = ell(3):dx:ell(3)+ell(2);

v = zeros(Hlen,1);         % initial conditions
rhs = zeros(Hlen,1);
e = ones(Hlen,1);
n = 0.3177*e;
m = 0.0529*e;
h = 0.5961*e;

t = 0;
tcnt = 0;
x = linspace(0,ell(3)+max(ell(1:2)),Hlen);
%plot3(x,t*ones(Hlen,1),v)
%hold on

%stimloc = N(1)+N(2)+5;
stimloc = round(N(1)/2); 

recs = round([N(1)/10 stimloc N(1)+N(2)/2 N(1)+1]);
%vrecs(:,1) = v(recs);

while (t <= Tfin)

      t = t + dt;

      Istim = I0*(t>t1)*(t<t2);

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

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

      ed = tau/dt + 1 + cNa + cK;

      H(1:Hlen+1:end) = dH + ed;         % update the diagonal

      rhs = v*tau/dt + cNa*vNa + cK*vK + vCl;
%      rhs(Hlen) = rhs(Hlen) + dt*Istim/(Cm*As);
      rhs(stimloc) = rhs(stimloc) + Istim/(2*pi*a(1)*dx*gL);

      v = H \ rhs;

      tcnt = tcnt + 1;
      
      vrecs(:,tcnt) = v(recs);

      if mod(tcnt,pinc) == 0 & 1==0
         v1 = fliplr(v(1:N(1))');
         v2 = fliplr(v(N(1)+1:N(1)+N(2))');
         v3 = fliplr(v(N(1)+N(2)+1:end)');
         plot3(x3,t*ones(size(x3)),v3,'k')
         plot3(x2,t*ones(size(x2)),[v3(end) v2],'b')
         plot3(x1,t*ones(size(x1)),[v3(end) v1],'r')
      end

end

%xlabel('x (cm)','fontsize',16)
%ylabel('t (ms)','fontsize',16)
%zlabel('v (mV)','fontsize',16)
%hold off

t = linspace(0,Tfin,tcnt);
plot(t,vrecs)

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