%plot the calcium release function for CICR
formatSpecF = '%5.2f\n';

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.0, ...
   'defaultlinelinewidth', 1.2, ...
   'defaultpatchlinewidth', 0.7); 

cmykflg=1;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;

%parameters
Vp=0.9;
kf=1.1;
ct=2;
K1=0.1;
K2=0.08;
gm=5.5;
Kp=0.1;
c=[0:0.01:ct];
plist=[0.12,0.15,0.2];
xplace=[0.6,0.6,0.6];
yplace=[-0.21,0.34,0.78];
for j=1:3
p=plist(j);
f = kf*(p*c./((p+K1)*(c+K2))).^3.*(gm*(ct-c)-c)-Vp*c.^2./(Kp^2+c.^2);

figure(1)
plot(c,f,'linewidth',2)
xlabel('Ca^{++}  (\mu M)','fontsize',20)
ylabel('J_{tot}','fontsize',20)
set(gca, 'XTick', 0 : 0.25 : 1.5);
text(xplace(j),yplace(j),strcat('p = ',sprintf(formatSpecF,p)),'fontsize',20)

hold on
end
axis([0 1.5 -0.4 1])
plot(c,zeros(length(c),1),'--','linewidth',2)
hold off

if(cmykflg==1) 
print('../../figs_c/chapt_7/calc_flux','-deps','-cmyk')
end

k2=K2/ct; 
k1 =(Kp/ct)^2;
u = [0.005:0.01:gm/(gm+2)];
p =(u+k2).^3./(u.*(u.^2+k1).*(gm*(1-u)-u));

figure(2)
plot(p,u,[0,max(p)],[0,0],'linewidth',2)
xlabel('\rho','fontsize',20)
ylabel('u','fontsize',20)
set(gca, 'XTick', 0 : 0.25 : 1.5);
axis([0 1.5 0 0.8])

ptmp = [0:.001:1];
ctmp=[0:0.001:3]';
ct= ctmp*ones(1,length(ptmp));
 p = ones(length(ctmp),1)*ptmp;
res = 0.531e-4*p-0.174e-3*ct-0.922e-3*ct.^2 ...
-0.937e-1*ct.^4.*p.^2+0.128e-2*ct.*p+0.934e-2*ct.^2.*p ...
    +0.16e-1*ct.^3.*p+0.3857e-1*ct.^4.*p.^3 ...
+0.191e-4*ct.^2.*p.^3+0.18e-4*ct.*p.^3+0.672e-2*ct.^3.*p.^3 ...
-0.499e-2*ct.^2.*p.^2-0.722e-4*ct.*p.^2 ...
-0.368e-1*ct.^3.*p.^2+0.211*ct.^5.*p.^3-0.3162e-4*ct.^2.*p.^4 ...
-0.4528e-2*ct.^4.*p.^4-0.1621*ct.^6.*p.^4-0.8245e-5 ...
    +0.1138e-5*p.^3+0.5308e-5*p.^2;
 
figure(3)
 contour(p,ct, res,[0 0],  'linewidth',2)
 xlabel('p','fontsize',20)
 ylabel('c_t','fontsize',20)
