%%  Program for Figure 5.12, Li and Rinzel Model
%%  Carol Lucas
%%  calls:  function xdot=lrclosefun(t,x,flag,I,L,Ve)

clear all

fi=0.01;Vi=4;P=26640;
Ki=1.0;Ka=0.4;
Ke=0.2;A=0.5;Kd=0.4;
sigma=0.185;L=.37;I=.6;Ct=2;Ve=400;

C=[.015:.01:1];Ce=(Ct-C)/sigma;
hnull=Kd./(C+Kd);

Iray=[.6:.01:.75 .8 .9 1 1.1:.01:1.3 1.4 1.5 1.6];
Lray=length(Iray);
cnull=zeros(3,length(C));
for k=1:Lray
    I=Iray(k);
    vars=['[IP3]= ' num2str(I) ' L= ' num2str(L) '  Vserca= ' num2str(Ve)];
    [t,x]=ode45('lrclosefun',[0 200],[.2;.8],[],I,L,Ve);
    [beat,tb,xbeat]=fanal(t,x(:,1));
    period(k)=beat(80);
    maxim(k)=max(x([length(t)-80:length(t)],1));
    minim(k)=min(x([length(t)-80:length(t)],1));
    Css(k)=lrclose_intersect(I,L,Ve);
    figure(1);clf;plot(t,x,t,beat);
    xlabel('t');ylabel('C and h'); 
    title([vars ' Period (red line) = ' num2str(period(k))]);pause(1);
    cnull(k,:)=(((Ve*C.*C./((Ke*Ke+C.*C).*(Ce - C))-L)/P).^(1/3)).*(I+Ki).*(C+Ka)./(I*C);
end
figure(4);plot(C,hnull,C,real(cnull));axis([0 1 .2 1.4]);
    xlabel('[C]');ylabel('h'); 
    title(['F5.13A: [IP3] Variations = '  num2str(Iray(1)) ' to ' num2str(Iray(Lray))]);
    mess=' Hit return to see bifurcation plot',pause;
figure(6);clf;plot(Iray,Css,Iray,maxim,Iray,minim);
xlabel('[IP3]');ylabel('[Ca2+i]');
title('F5.12A');
mess=' Hit return to see change in period with [IP3]',pause;
figure(7);clf;plot(Iray,period);xlabel('[IP3]');ylabel('Period (s)');
title('F5.12B');axis([.7 1.2 0 100]);