function [p, K, R2_vector, v_vector, p_other, Model]=B_model_matrices(Model, Rtotal, LRratio)
%------------------- LineshapeKin Simulation -----------------------------
% 5/12/09
% Evgenii Kovrigin
% Medical College of Wisconsin
%
%
% This function builds  model-specific matrices to supply to Spectra:
% rate (K) matrix for the spin exchange,
% equilibrium populations
% vectors of pure frequencies and R2
%
%
% B model was developed in 
% /Users/kovrigin/Documents/Workspace/Data/Data.XV/EKM16.Analysis_of_multistep_kinetic_mechanisms/Equilibria/B_model/B_model_derivation.mn
%
% Bugs and fixes: 
% 2/25/2012 - 'LevenbergMarquardt option is no longer valid' issue
%                     addressed.

%mode='analytical';
mode='numeric';

% Current model:  B 

%--- assign model parameters corresponding to the order 
%    in 'model_setup.m' 
 
K_A_1_a =  Model.Ka( 1);   
K_A_1_b =  Model.Ka( 2);   
K_A_2_a =  Model.Ka( 3);   
 
k_2_A_1_a = Model.k2( 1);
k_2_A_1_b = Model.k2( 2);
k_2_A_2_a = Model.k2( 3);
k_2_A_2_b = Model.k2( 4);
 



%-- assign experimental parameters:
Rtot = Rtotal;    
Ltot = Rtot*LRratio;  


 

%% Determine equilibrium concentrations of species

switch mode
    
    %---- traditional explicit calculation of [L]
    case 'analytical'
            %% Analytical solution for [L] out of 
            
            % B model does not have analytical solution
             
            
            
    
    case 'numeric'
            %%  minimization to obtain [L]
            

            Leq=0;

            if Ltot>0   %-- to skip 0 concentration point
            %   obtain equilibrium concentration of the ligand by numeric
            %   minimization.

            % --STARTING CONDITIONS. It is essential for numeric stability
            % of the solution to start from just below full concentration.
            % Sometimes it may be unstable - try using different starting
            % conditions.
           
            Leq0=Ltot*0.9;
            params=[Leq0];
            
            % -- I formally define X. The first equilibrium constant
            % goes here. The other constants are supplied directly to
            % the function below.
            x_local=Model.Ka(1);
            
            %  --- we define the target function using 'anonymous' synthax
            %  and pass additional constants to it here. Order of the
            %  constants in defined in model_numeric.m
            myfun = @(B,X) ( B_model_numeric(B, [X  Model.Ka(2:end) Ltot]) );
            
            % this is Y-array
            y_local=Rtot;

             %-- physically meaningful bounds
             lower_bounds=[1e-32];
             upper_bounds=[Ltot];
             
             %-- default fitting options. Sometimes TolFun and TolX need
             %    to be adjusted for better convergence.
             options=optimset('Diagnostics','off', ...
                  'Display','off',...
                  'LargeScale','on',...
                  'TolFun',1e-18,...
                  'TolX',1e-18,...
                 'MaxFunEvals', 1e9);
                                     
             % --- run least-squares minimization protocol to extract
             % equilibrium concentrations we need for K-matrix ('params'
             % vector)
             [Leq, resnorm, residual, exitflag, output,lambda, Jacobian]=...
                  lsqcurvefit(myfun, params, x_local, y_local,...
                  lower_bounds,upper_bounds,options);

  
            end
            
end


%-- compute dependent equilibrium constant

 K_A_2_b = (K_A_1_a*K_A_2_a)/K_A_1_b;   % eq 4_10  

% save dependent constants for reporting
Model.Ka_dependent_names={ 'A2b' };
Model.Ka_dependent = [ K_A_2_b ]; 

%-- compute concentrations of  other species


Req = -(Leq - Ltot)/(K_A_1_a*Leq + K_A_1_b*Leq + 2*K_A_1_a*K_A_2_a*Leq^2)  ;  % eq   4.7

LaReq = K_A_1_a*Leq*Req   ; % eq  4.4


RLbeq = K_A_1_b*Leq*Req  ;  % eq  4.5


RL2eq = K_A_2_a*LaReq*Leq  ; % eq  4.3



%% Line shape equations

% Equilibrium mole fractions of species containting P 
%                R               LaR               RLb             RL2
p=[        Req/Rtot,   LaReq/Rtot,    RLbeq/Rtot,  RL2eq/Rtot];

% Normalized concentration of the ligand for plotting
p_other=[Leq/Rtot]; 

%-- express forward rate constants through defined reverse rate constants
k_1_A_1_a=k_2_A_1_a*K_A_1_a;
k_1_A_1_b=k_2_A_1_b*K_A_1_b;
k_1_A_2_a=k_2_A_2_a*K_A_2_a;
k_1_A_2_b=k_2_A_2_b*K_A_2_b;

% save them in the order of Model.k_names !!!
Model.k1 = [    k_1_A_1_a,    k_1_A_1_b,    k_1_A_2_a,     k_1_A_2_b   ];
 

% Write matrix elements for the rate equations for the system
%                                        [R]                                            [LaR]                                        [RLb]                                       [RL2] 
%  d[R]/dt     =    -(k_1_A_1_a+k_1_A_1_b)*Leq      k_2_A_1_a                                k_2_A_1_b                                               0
%  d[LaR]/dt   =     k_1_A_1_a*Leq                        -(k_2_A_1_a+k_1_A_2_a*Leq)                       0                            k_2_A_2_a 
%  d[RLb]/dt   =                       k_1_A_1_b*Leq                        0                            -(k_2_A_1_b+k_1_A_2_b*Leq)                       k_2_A_2_b
%  d[RL2]/dt =                       0                                                 k_1_A_2_a*Leq                          k_1_A_2_b*Leq   -(k_2_A_2_a+k_2_A_2_b)                    

% Copy constants to a K-matrix
K=[    ...
-(k_1_A_1_a+k_1_A_1_b)*Leq      k_2_A_1_a                                k_2_A_1_b                                               0                        ;...
  k_1_A_1_a*Leq                        -(k_2_A_1_a+k_1_A_2_a*Leq)                       0                            k_2_A_2_a                         ;...
                      k_1_A_1_b*Leq                        0                            -(k_2_A_1_b+k_1_A_2_b*Leq)                       k_2_A_2_b    ;...
                    0                                                 k_1_A_2_a*Leq                          k_1_A_2_b*Leq   -(k_2_A_2_a+k_2_A_2_b)   ;...
    ];


 


%-- prepare relaxation rates and frequencies for the individual species
 R2_vector= Model.R2;  
 v_vector=  Model.w0; 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
