%
% Plots Pade Approximations produced by PadeApproximation(N, lambda, mu) 
%
%                    Compare to Figure 3 in paper:
%  "Padé approximation of a stationary single-channel Ca2+ nanodomain"
%     Victor Matveev, Sep 10 2016 * Biophysical Journal Supp Info
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

syms Buff Cal r;  % Reserve symbolic names returned by PadeApproximation()

lambda  = 0.05;   % Set to any value > 0
mu      = 0.05;   % Set to any value > 0

Rleft  = 1e-3;    % Minimal distance to channel
Rright = 100;     % Maximal distance to channel

Rarray = exp( linspace(log(Rleft), log(Rright), 80) ); % Log-distance grid

subplot(121);                        % Set properties of [Buffer] subplot
hold off; set(gca, 'fontsize', 10);
xlabel('\rho'); ylabel('b(\rho)');
set(gca, 'Xscale', 'log');
title('A       Buffer concentration', 'fontsize', 14);

subplot(122);                        % Set properties of [Ca2+] subplot
hold off; set(gca, 'fontsize', 10); 
xlabel('\rho'); ylabel('c(\rho)');
set(gca, 'Xscale', 'log', 'Yscale', 'log');
title('B       Ca^{2+} concentration', 'fontsize', 14);

for N = 1 : 5
    [Buff, Cal] = PadeApproximation(N, lambda, mu);  % Get the Pade approdimation
    Buffer  = subs(Buff, r, Rarray);                 % Compute buffer  on the r-grid
    Calcium = subs(Cal,  r, Rarray);                 % Compute calcium on the r-grid
    shade = 0.2*(5-N);
    
    subplot(121); hold on;
    plot(Rarray, Buffer,  '-', 'color', [1, shade, shade], 'linewidth', 3-N/2);
    subplot(122); hold on;
    plot(Rarray, Calcium, '-', 'color', [1, shade, shade], 'linewidth', 3-N/2);
    
    drawnow;
end;

subplot(121);
legend('Pade 1^{st} order', 'Pade 2^{nd} order', 'Pade 3^{rd} order', 'Pade 4^{th} order', 'Pade 5^{th} order'); 
axis([Rarray(1) Rarray(end) min(eval(Buffer))  max(eval(Buffer)) ]);

subplot(122);
legend('Pade 1^{st} order', 'Pade 2^{nd} order', 'Pade 3^{rd} order', 'Pade 4^{th} order', 'Pade 5^{th} order'); 
axis([Rarray(1) Rarray(end) min(eval(Calcium)) max(eval(Calcium)) ]);



