%% Testing of the B-micro model
% Verify against MuPad "analysis" document:  B_family_model_analysis.html
%%

close all 
clear all

%%
% Here you will find all figures
figures_folder='B_micro_testing_figures';


%%
% *TEST COMPUTATION OF EQUILIBRIUM POPULATIONS*
%
% Set some meaningful parameters matching settings of MuPad analysis
% document
Rtotal=1e-3; % Receptor concentration, M
LRratio_array=[0 : 0.02 : 3]; % Array of L/R 
K_A_1_I   = 5000;
K_A_1_I_I = 2500;
K_A_2_I_I = 20000;

% 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*
concentrations_array=[];
for counter=1:length(LRratio_array)
    % compute
    [concentrations species_names] = equilibrium_thermodynamic_equations.B_micro_model(...
        Rtotal, LRratio_array(counter), K_A_1_I, K_A_1_I_I, K_A_2_I_I,...
        model_numeric_solver, model_numeric_options);
    % collect
    concentrations_array = [concentrations_array ; concentrations];
end

%% 
% Plot
Figure_title= 'B-micro monodentate ligand model';
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(...
    LRratio_array, concentrations_array, species_names, Figure_title, X_range, Y_range);
% save it
results_output.output_figure(figure_handle, figures_folder, 'Concentrations_plot');

%% 
% *Observations*
%  The simulation matches results from the "analysis" MuPad notebook. The
%  code of B-micro model is correct.















