%plot the hh right hand side of the ``fast" upstroke  HH equations

function HH_plot
global gnabar

 set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7); 

cmykflg=0;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;

dv=0.01;
V = [-20.000016:dv:120];  %This is the relevant range of V
gnabar = 120;  %nominal value

%find 3 zeros of I_ion
Ion_current = HH_rhs(V);
%ndx has three numbers
ndx = find(Ion_current(1:end-1).*Ion_current(2:end)<=0)
   
%use bisection to determine the three zeros (fairly accurately)
a = V(ndx);
b = V(ndx+1);
r = bisect(a,b,@HH_rhs)

% Estimate the slope of I_ion at the rest state 
 % slope at the rest point
 RM_inv = (Ion_current(ndx(1)+1)-Ion_current(ndx(1)))/dv;
 %Rm is the membrane resistance
 Rm = 1/RM_inv

figure(1)
plot(V,Ion_current ,'linewidth',2);
hold on
plot( r(1),0,'*',r(2),0,'*',...
    r(3),0,'*','linewidth',4);
hold off

xlabel('V (mV)','fontsize',20)
ylabel('-I_{ion}(V) (\mu A/cm^2)','fontsize',20)
set(gca, 'XTick', 0 : 20 : 120);
axis([-20 120 -100 4000])

if(cmykflg==1) 
print('../../figs_c/chapt_7/hh_I_ion','-deps','-cmyk')
end


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;

 function root = bisect(a,b,feval)

ul = a;
fl = feval(ul);
uu = b;
fu = feval(uu);

% we make the assumption, without checking, that 
% fu*fl<0

% if not, the algoritm fails to find a root.

N = 20;  % number of iterates
% the main bisection algorithm
for j = 1:N
uc = (ul+uu)/2;
fc = feval(uc);
ftest = (fc.*fl>=0);
ul = ftest.*uc+(1-ftest).*ul;
fl = ftest.*fc + (1-ftest).*fl;

uu = (1-ftest).*uc+ ftest.*uu;
fu = (1-ftest).*fc + ftest.*fu;
end
root = uc;