
%
% the reduced mec model from Koppel and co (sec 3.2.1)
%
%
% C v' = -Gp*pinf(v)(v-ENa) - GL*(v-EL) - Gh*(.65rf+.35rs)(v-Eh)
%
%  rf' = (rfinf(v)-rf)/taurf(v)
%
%  rs' = (rsinf(v)-rs)/taurs(v)
%

%syms v

v = -72;

pinf = 1/(1+exp(-(v+38)/6.5));

rfinf = 1/(1+exp((v+79.2)/9.78));
rfinf = rfinf*(1+0.2*(v+70)*exp(-(1/80)*(v+70)^2));

drfinf = -50/489/(1+exp(50/489*v+1320/163))^2*(1+(1/10*v+7)*exp(-1/80*(v+70)^2))*exp(50/489*v+1320/163)+1/(1+exp(50/489*v+1320/163))*(1/10*exp(-1/80*(v+70)^2)+(1/10*v+7)*(-1/40*v-7/4)*exp(-1/80*(v+70)^2));

taurf = 1 + 0.51/(exp((v-1.7)/10)+exp(-(v+340)/52));

rsinf = 1/(1+exp((v+2.83)/15.9))^58;

taurs = 1 + 5.6/(exp((v-1.7)/14)+exp(-(v+260)/43));

ENa = 50;
EL = -70;
Eh = -80;
GL = .05;
Gh = .75;
%Gp = 2;

% restfun = Gp*pinf*(v-ENa) + GL*(v-EL) + Gh*rfinf*(v-Eh);

Gp =  - (GL*(v-EL) + Gh*rfinf*(v-Eh) )/(pinf*(v-ENa))

%ezplot(restfun,[-100 0])

%tv = ginput(1);

%bv = tv(1)

%dpinf = diff(pinf);

dpinf = 2/13/(1+exp(-2/13*v-76/13))^2*exp(-2/13*v-76/13);

%syms tau

for tau = 1:.1:15,

   A = [-Gp*dpinf*(v-ENa)-GL-Gh*rfinf  -Gh*(v-Eh)
       drfinf/tau -1/tau];

   e = eig(A);

   plot(real(e),imag(e),'*')
   hold on

end

hold off

%hopf = Gp*( subs(dpinf,'v',bv)*(bv-ENa) + subs(pinf,'v',bv) ) + GL + ...
%       Gh*subs(rfinf,'v',bv)
