function [numSSPR,diffProjFullEq,fval,exitflag,output,Jacob,eigJacob,nzeig] = GetFzeroSSPR_db(aPR,VdsOut,quasiSS)
%% What does it do?
% This serves as another check for the steady state by using fsolve on
% the full 8-dimensional problem. It uses as a starting point the
% equilibrium state found by projecting the problem to the Vs Vd plane in
% ProjGateEquilVsVd_db. The differenxe between the two should be minimal
% One advantage is we get the full 8-d linearization with eigenvalues.
%options=optimset('TolFun',1e-14,'Display','iter','MaxIter',1000,'MaxFunEvals',100000);
options=optimset('TolFun',1e-14,'MaxIter',1000,'MaxFunEvals',100000,'Display','none');
[numSSPR,fval,exitflag,output,Jacob] = fsolve(@PR94rhs,quasiSS,options);
eigJacob=eig(Jacob);
nzeig=0;
for i=1:size(eigJacob,1)
if eigJacob(i) > 0
% display('Warning one of the eigenvalues of the Jacobian is positive')
nzeig=nzeig+1;
end
end
diffProjFullEq=numSSPR-quasiSS;
function dY = PR94rhs(Y)
dY = zeros(8,1); % a column vector
dY(1) = (1/aPR.Cm)*(-aPR.gL*(Y(1)-aPR.EL)-aPR.gNa*MInfPR94(Y(1),aPR.WRT)*Y(4)*(Y(1)-aPR.ENa)-aPR.gKDR*Y(5)*(Y(1)-aPR.Ek) ...
+(aPR.gc/aPR.p)*(Y(2)+VdsOut-Y(1))+aPR.Isinj/aPR.p);
dY(2) = (1/aPR.Cm)*(-aPR.gL*(Y(2)-aPR.EL)-aPR.gCa*Y(6)^2*(Y(2)-aPR.ECa)-aPR.gKAHP*Y(8)*(Y(2)-aPR.Ek)-aPR.gKC*Y(7)*Chi(Y(3))*(Y(2)-aPR.Ek)...
+ aPR.gc/(1-aPR.p)*(Y(1)-VdsOut-Y(2))+aPR.Idinj/(1-aPR.p));
dY(3) = -0.13*aPR.gCa*Y(6)^2*(Y(2)-aPR.ECa)-0.075*Y(3);
dY(4) = (GateEquil_db(alphah_db(Y(1),aPR.WRT),betah_db(Y(1),aPR.WRT))-Y(4))/GateTimeCnst_db(alphah_db(Y(1),aPR.WRT),betah_db(Y(1),aPR.WRT));
dY(5) = (GateEquil_db(alphan_db(Y(1),aPR.WRT),betan_db(Y(1),aPR.WRT))-Y(5))/GateTimeCnst_db(alphan_db(Y(1),aPR.WRT),betan_db(Y(1),aPR.WRT));
dY(6) = (GateEquil_db(alphas_db(Y(2),aPR.WRT),betas_db(Y(2),aPR.WRT))-Y(6))/GateTimeCnst_db(alphas_db(Y(2),aPR.WRT),betas_db(Y(2),aPR.WRT));
dY(7) = (GateEquil_db(alphac_db(Y(2),aPR.WRT),betac_db(Y(2),aPR.WRT))-Y(7))/GateTimeCnst_db(alphac_db(Y(2),aPR.WRT),betac_db(Y(2),aPR.WRT));
dY(8) = (GateEquil_db(alphaq_db(Y(3)),betaq_db)-Y(8))/GateTimeCnst_db(alphaq_db(Y(3)),betaq_db);
end
function CaSatChi= Chi(Ca)
CaSatChi = min(Ca/250,1);
end
function MInfsqr = MInfPR94(Vs,WRT)
MInfsqr = GateEquil_db(alpham_db(Vs,WRT),betam_db(Vs,WRT))^2;
end
end