clear all
global G1 G2 G_Cl1 G_Cl2 V1 V2 r1 r2 i1 i2 E_Na E_K E_Cl R el C Iin eps

%
% Parameter values for the Spike Burster IV-characteristics
%  1/r1<-G2
%
G1=1; G2=-0.5; V1=.5; V2=3;

G_Cl1=.01; G_Cl2=0; 

r1=3; r2=-5; i1=.5; i2=.75;

E_Na = 0.5; E_K = - 0.75; E_Cl = - 0.6; 


%
% Circuit simulation
%


R=0.01; el=.05; L=6; C=.01; Iin=0; eps=.001; % L=.5, 1.5, 2, 3 => 1, 2, 3, 4 spikes
% R=0.1; L=.5; C=.01; Iin=-7; eps=.001; % L=.3, .5 => 2, 3 spikes
%%%% R=0.1; L=.1; C=.01; Iin=-7; eps=.001; % close to Shilnikov orbit
% R--resistance, L--inductance, C--capacitance, 
% E--battery(Nernst potential)
% Increasing L to increase the number of spikes
% Increase Iin to suppress the spike bursts. 

%E_K=1;r2=-1;L=1;C=.001;eps=.0001;



T=1; % transcient time
Tend=8;% Numerical solutions
xinit=[-0.2965;    2.8680;    3.0566;   -0.2655];
[t,x]=ode15s('pKsNaModel',[0,Tend],xinit);


figure(20)
hold off
N=find(t>T); % Skip the first T interval solution 
plot(t(N:end),x(N:end,1),'r-',t(N:end),x(N:end,4),'g-',t(N:end),x(N:end,2),'m-',...
    t(N:end),x(N:end,3),'k-',t(N:end),x(N:end,2)-x(N:end,3),'b-')
hold on
legend('V_C','I_{Na}','A_{Na}','A_K','I_L = A_{Na} - A_K')

figure(30)
hold off
plot3(x(N:end,1),x(N:end,2)-x(N:end,3),x(N:end,4),'r-')
hold on

figure(50)
hold off
j_K=K(x(N:end,1)-E_K);
plot(t(N:end),j_K, 'b-', t(N:end),x(N:end,4),'r-')
hold on


figure(40)
hold off
V=-1:.05:3;
i=-.5:.05:1.5;
plot(V,K(V-E_K),'r-',Na(i)+E_Na,i,'b-',[-1 3], [0 0], [0 0], [-1 1.5])

figure(40)
hold on
plot(x(N:end,1),x(N:end,4),'g-')
