%% Code for Figure 5.9
%% Carol Lucas

clear CT Cray
jimax=1.5;init=[.5;1];
jimax=input('Enter jin (1.5 in F5.9) ');
init(1,1)=input('Enter starting Cai (range=0 to 3 in F5.9) ');
init(2,1)=input('Enter starting CT (range=1.2 to 3.2 in F5.9) ');
fi=0.01;kserca=0.2;vserca=100;vRyR=5;vleak=.2;sig=.02;
vpmca=5.; kpmca=.6;
ka=0.4;kb=0.6;kc=.1;
jistart=0;jiend=1000;kcm=.1;

Cray=[.01:.01:3];

for k=1:length(Cray)
    C=Cray(k);
    winf=(1+(ka/C)^4+(C/kb)^3)/(1+(1/kc)+(ka/C)^4+(C/kb)^3);
    Po=winf*(1+(C/kb)^3)/(1+(ka/C)^4 + (C/kb)^3);
    CT(k)=(C+(C^2*vserca/(C^2+kserca^2))/(vRyR*Po+vleak))*sig+C;
end

figure(1);clf;plot(CT,Cray);xlabel('CT');ylabel('[Ca2+]');
title(['Figure 5.9, jin = ' num2str(jimax)]);
axis([0 4 0 3]);hold on;
mess=' Hit return to move to do open cell simulation',pause;


[t,x]=ode45('klopenfun_red',[0:1000],[init],[],jimax,jistart,jiend,kcm);
figure(2);clf;plot(t,x);title('kl reduced open model');
mess=' Hit return to see phase plane',pause;

figure(1);plot(x(:,2),x(:,1),'r',x(1,2),x(1,1),'ro');
CTnull=sqrt((jimax*kpmca^2)/(vpmca-jimax));
plot([0 4],[CTnull,CTnull],'m:');