%here we estimate the critical sodium conductance need
% to guarantee propagation

function HH_speed
global gnabar

dv = 0.0001;
V = [-12.000016:0.0001:115];
gnabar = 120;  %nominal value
gstep = 0.01;
for j = 1:20
    gnabar = 1.4+j*gstep;
    gbar(j) =gnabar;
Ion_current = HH_rhs(V);
%find 3 zeros

ndx = find(Ion_current(1:end-1).*Ion_current(2:end)<=0);
 
sion=  Ion_current(ndx(1):ndx(3));
 
figure(6)
plot(sion)
I_int = Ion_current(ndx(1))/2  +sum(sion)+...
     Ion_current(ndx(3))/2;
 
INT(j) = I_int*dv;
end
 idx = find(INT(1:end-1).*INT(2:end)<=0);
 gbstar =gbar(idx)- gstep*INT(idx)/(INT(idx+1)-INT(idx))
figure(2)
plot(gbar,INT)
Veq = -65;
figure(1)
plot(V,Ion_current,V(ndx(1)),0,'*',V(ndx(2)),0,'*',...
    V(ndx(3)),0,'*','linewidth',2);
xlabel('V','fontsize',20)
ylabel('-I_{ion}(V)','fontsize',20)

function Ion_current = HH_rhs(V)
global gnabar
gkbar = 36;

h = 0.5961;
n = 0.3177;

AM=.1*(25.-V)./(exp(.1*(25.-V))-1.);
BM=4.*exp(-V/18);
m = AM./(AM+BM);

gna=gnabar*m.^3.*h;
 
gk=gkbar*n.^4;

Vna = 115 ;
Vk =  -12 ;
Vl =  10.5988;
Icl=.3*(V-Vl);

Ion_current = -gna.*(V-Vna)-gk.*(V-Vk)-Icl;
