%Attempt to find the integral of f for calcium rhs
formatSpecF = '%5.2f\n';

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

%parameters
Vp=0.9;
kf=1.1;
ct=2;
K1=0.1;
K2=0.08;
gm=5.5;
Kp=0.1;
Jer = 0.0;

cbig = gm*ct/(gm+1);

c=[0.001:cbig/200:cbig];
 
% make a plot of the zeros of F as a function of P

%find P as a function of c
Jserca = Vp*c.^2./(Kp^2+c.^2);
ce =  gm*(ct-c);
P=((c+K2)./c).^3.*(Jserca./(ce-c)-Jer)/kf;

%f = kf*P(c./( c+K2)).^3.*(gm*(ct-c)-c)-Vp*c.^2./(Kp^2+c.^2);

figure(1)
plot(P,c,'linewidth',2)
axis([0 1, 0 ct])
xlabel('P')
ylabel('c, the zero of f(c)')

% Now for each value of c, hence P, find the integral 0 to c of f

for j = 1:length(c)
    C  = c(j);  %This is the zero of the function f
    p = Vp*C.^2./(Kp^2+C.^2).*((C+K2)./C).^3./ (gm*(ct-C)-C)/kf;  %this is P
    
    u=[0:C/100:C];  % discretize the interval from 0 to C
    fu = kf*p.*(u./(u+K2)).^3.*(gm*(ct-u)-u)-Vp*u.^2./(Kp^2+u.^2);
    intf(j)=sum(fu)*C/100; % estimate the integral with trapezoidal rule
end
 % now estimate where this changes sign
 fndx = find(intf>=0);
 ndx = min(fndx);
 Pcr = P(ndx)  % this is the critical value of P
 P(ndx-1)
 Cst = c(ndx) %the largest zero
 intf(ndx)
 intf(ndx-1)

 %For this value of P:
 p =-K2/(1-Pcr^(-1/3))
 figure(2)
plot(P,intf,[0,max(P)],[0,0],'--')
axis([0.14 0.2 -.1 .1])
xlabel('P')
ylabel('The integral of f(c)')

figure(4)
plot(P,c,P,intf,[0,max(P)],[0,0],'--')
axis([0.14  0.2 -.1 .8])
xlabel('P')

 
ce = gm*(ct-c);
Jserca = Vp*c.^2./(Kp^2+c.^2);
JIPR = Pcr*(c./(c+K2)).^3;
Jtot = (kf*JIPR+Jer).*(ce-c)-Jserca;

c_flux = Jtot;
 
figure(3)
plot(c ,c_flux ,c ,zeros(1,length(c)),'--','linewidth',2)
xlabel('c','fontsize',20)
ylabel('Total Flux J_{tot}','fontsize',20)
axis([ 0 1 -0.2 0.2])
 