clear all;

global vNa vK vL vCa vH ...
    gNa gK gL gCa gSK gNap gKNa gH gA ...
    tauNbar tauHPbar tauE tauH Ca_ex RTF ...
    thetaM sigmaM thetaN sigmaN thetaS sigmaS ...
    thetaMP sigmaMP thetaHP sigmaHP thetaRF sigmaRF thetaRS sigmaRS ...
    kr prf taurs f ...
    thetaA sigmaA thetaE sigmaE ...
    C iapp eps kca ks pmax ec nh alphaNa rpump kp neq ...
    gT thetaaT sigmaaT thetabT sigmabT thetarT sigmarT taur0 taur1 thrT sgmrt phirT ...
    from to
  
  
% Nernest potentials in mV
vNa=50; 
vK=-90; 
vL=-70; 
vH=-30;

% Conductances in nS
gNa=450;
gK=50;
gL=2;
gCa=19;
gNap=1;
gSK=6;
gKNa=40;
gH=4;
gA=5;
gT=2.65;

% Time constants in mS
tauNbar=10;
tauHPbar=1000;
tauE=20;
tauH=1;
taurs=1500;
RTF=26.7;
Ca_ex=2.5;

% Function constants in mV
thetaM=-35;
thetaN=-30;
sigmaM=-5;
sigmaN=-5;
thetaS=-20;
sigmaS=-0.05;
thetaMP=-40;
sigmaMP=-6;
thetaHP=-48;
sigmaHP=6;
thetaA=-20;
sigmaA=-10;
thetaE=-60;
sigmaE=5;


thetaaT = -65;
sigmaaT = -7.8;
thetabT = 0.4;
sigmabT = -0.1;
thetarT = -67;
sigmarT = 2; 
taur0 = 40;
taur1 = 17.5;
thrT = 68;
sgmrt = 2.2;
phirT=0.2;

thetaRF=-105;
sigmaRF=5;
thetaRS=-105;
sigmaRS=25;

% Capacitance in pF
C=100;

current=[-280 -240 -200 -160 -120 -80 -40 0];

from=160;
to=600;

time=0:.01:800;

% Other constants
kr=0.3;
prf=100;

f=0.1;
eps=0.0015;
kca=0.3;
ks=0.5;
pmax=0.37;
ec=38.7;
nh=3.5;
alphaNa=0.0001;
rpump=0.0008;
kp=15;
neq=8;

% Initial Conditons
v_i=-72.14; 
n_i=.0002;
ca_i=0.103;
h_i=0.01;
na_i=7.09;
rf_i=0.03;
rs_i=0.03;
e_i=0.05;
rT_i=0.01;
hp_i=0.1;

init=[n_i;h_i;rf_i;rs_i;ca_i;v_i;na_i;e_i;rT_i;hp_i];

vEnd = [];
vMin = [];

for i=1:size(current,2)
    iapp=current(i);
    
    [time,output]=ode113('DiffEquations',time,init);
    
    v = output(:,6);
        
    fromIndex=1;
    s=size(time,1);    
    while (time(fromIndex) < from) && (fromIndex <= s)
        fromIndex = fromIndex + 1;
    end
    
    toIndex=1;
    while (time(toIndex) < to) && (toIndex <= s)
        toIndex = toIndex + 1;
    end

    minVolt = min(v(fromIndex:toIndex));
    
    if(iapp < 0)
        vMin = [vMin minVolt];
        vEnd = [vEnd v(toIndex-1)];
    else
        vEnd = [vEnd v(toIndex)];
        vMin = [vMin v(toIndex)];
    end    
end

figure();
plot(current, vMin, 'b^-', 'MarkerSize', 8, 'LineWidth', 2);
hold on
plot(current, vEnd, 'bo-', 'MarkerSize', 8, 'LineWidth', 2);
plot(current, vMin, 'b-', 'LineWidth', 2);
plot(current, vEnd, 'b-', 'LineWidth', 2);

