%% Testing of the U-2R-RL model
% Verify against MuPad "analysis" document:  U-2R-RL_analysis.html
% Run in "Publish" mode to have a report.
%
%%

close all 
clear all



%% Set input parameters
%  
model_name='U-2R-RL';

% Set some meaningful parameters matching settings of MuPad analysis
% document
LRratio_array=[0.0001 : 0.05 : 1.5]'; % Array of L/R 

Rtotal_array=1e-3*ones(length(LRratio_array),1); % Receptor concentration, M


K_A   = 1e5;
K_B_1_s_1 = 2;
K_B_1_s_2 = 4;
K_B_2 = 2;

% Set appropriate options for the model (see model file for details)
% model_numeric_solver='fminbnd'  ; 
% 
% model_numeric_options=optimset('Diagnostics','off', ...
%                     'Display','off',...
%                     'TolX',1e-9,...
%                     'MaxFunEvals', 1e9);
%%
% The important option here is "TolX" that sets termination tolerance on free ligand
% concentation in molar units. With our solution concentrations in 1e-3
% range TolX should be set to some 1e-9.



%% Compute arrays for populations and plot
% Call as non-vectorized:
concentrations_array=[];
tic
for counter=1:length(LRratio_array)
    % compute
    [concentrations species_names] = equilibrium_thermodynamic_equations.U_2R_RL_model(...
        Rtotal_array(counter), LRratio_array(counter), K_A, K_B_1_s_1, K_B_1_s_2, K_B_2, ...
        'analytical', 'none');
    % collect
    concentrations_array = [concentrations_array ; concentrations];
end
toc

%% 
% Plot
Figure_title= model_name;
X_range=[0 max(LRratio_array)+0.1 ]; % extend X just a bit past last point
Y_range=[ ]; % keep automatic scaling for Y
% display figure
figure_handle=equilibrium_thermodynamic_equations.plot_populations(...
    Rtotal_array, LRratio_array, concentrations_array, species_names, Figure_title, X_range, Y_range);

%% 
% *Observations*
%  The result matches simulations in the 'analysis' MuPad notebook. The
%  code is correct.

% 
% %% Test vectorization (if analytical model was vectorized!)
% % Call as vectorized:
% tic
% [concentrations_array species_names] = equilibrium_thermodynamic_equations.U_5R_model(...
%         Rtotal_array, LRratio_array, K_A,...
%         K_B_s_1, K_B_s_2, K_B_s_3, K_B_s_4, K_B_s_5,...
%         'analytical', 'none');
% toc
% 
% Figure_title= sprintf(' %s - vectorized version', model_name);
% X_range=[0 max(LRratio_array)+0.1 ]; % extend X just a bit past last point
% Y_range=[ ]; % keep automatic scaling for Y
% % display figure
% figure_handle=equilibrium_thermodynamic_equations.plot_populations(...
%     Rtotal_array, LRratio_array, concentrations_array, species_names, Figure_title, X_range, Y_range);
% 
% 
% %% Conclusions
% % Both vectorized and non-vectorized calls work well.








