function [BufferCa0, BufferCa2, Calcium, r] = Get_RBA_AnyRestingCa(nu2, eps, gamma, CaInfty, lineType)
% 
% Compute RBA approximation for concentrations of Ca2+ 
%   and 2:1 buffer state concentrations near a point channel / source
%        Case of arbitrary background / resting [Ca2+] at infinity, CaInfty
%
% Input variables:
%  (1-4) Dimensionless parameters: lambda2, nu2, eps, gamma, CaInfty
%  (7)   lineType: MATLAB plot linetype string, e.g. 'b--'
%
% Return arrays: BufferCa0 = [Free buffer]        vs. distance from channel 
%                BufferCa2 = [Fully bound buffer] vs. distance from channel 
%                Calcium   = [Ca2+]               vs. distanc(see Eq. 28)
%                r         = Array of distances from channel see Eq. 28)
%
%       Victor Matveev and Yinbo Chen (NJIT) September 2020
%  (GPL) Free to use, share, modify, with appropriate attribution
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

r    = logspace(-3, 2, 100);
wInf = 2 * eps * CaInfty;
vInf =     eps * CaInfty^2;

bT   = 1 + wInf + vInf;
cT   = CaInfty + nu2 * (wInf/2 + vInf);

%%%%%%%%%%%%  RBA BLOCK :  See Equation S31  %%%%%%%%%%%%%%%%

F =  ( 1./r - 2 - nu2*bT +   cT ) / 3;           
G = ( -2./r +     nu2*bT - 2*cT + 1/eps ) / 3;      
S =  F.^2 - G;                                    
R =  F .* (F.^2 - 3/2*G) + (1./r + cT) / (2*eps);
A =  ( sqrt(R.^2 - S.^3) + R ).^(1/3);           

Calcium   = real( A + S ./ A  + F );
BufferCa0 =                    bT ./ (1 + eps * Calcium .* (Calcium + 2));
BufferCa2 = eps * Calcium.^2 * bT ./ (1 + eps * Calcium .* (Calcium + 2));
BufferCa1 = bT - BufferCa0 - BufferCa2;

if ~isempty(lineType)
    
    subplot(411); hold on;
    set(gca, 'fontsize', 14);
    title(['[Free buffer] b (\nu_2=', num2str(nu2), ' \epsilon=', num2str(eps), ' \gamma=', num2str(gamma),' c_\infty=', num2str(CaInfty), ')']);
    plot(r, BufferCa0, lineType, 'linewidth', 2);
    set(gca, 'xscale', 'log');
    xlabel('r'); ylabel('b');
    axis tight;

    subplot(412); hold on;
    set(gca, 'fontsize', 14);
    title(['[Single-bound buffer] b*']);
    plot(r, BufferCa1, lineType, 'linewidth', 2);
    set(gca, 'xscale', 'log');
    xlabel('r'); ylabel('b*');
    axis tight;
    
    subplot(413); hold on;
    set(gca, 'fontsize', 14);
    title(['[Fully-bound buffer] b**']);
    plot(r, BufferCa2, lineType, 'linewidth', 2);
    set(gca, 'xscale', 'log');
    xlabel('r'); ylabel('b**');
    axis tight;
    
    subplot(414); hold on;
    set(gca, 'fontsize', 14);
    title(['[Ca2+] Calcium']);
    plot(r, Calcium, lineType, 'linewidth', 2);
    set(gca, 'xscale', 'log', 'yscale', 'log');
    xlabel('r'); ylabel('Calcium');
    axis tight;
end