function [p, K, R2_vector, v_vector, p_other, Model]=B_R2_R2L2_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_R2_R2L2 model was developed in 
% [ek_mac]: /Users/kovrigin/Documents/Workspace/Data/Data.XV/EKM16.Analysis_of_multistep_kinetic_mechanisms/Equilibria/B_R2_R2L2_model/B_R2_R2L2_derivation.mn 
%
% Bugs and fixes: 
% 2/25/2012 - 'LevenbergMarquardt option is no longer valid' issue
%                     addressed.


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

 

%--- assign model parameters corresponding to the order 
%    in 'model_setup.m' 
 
K_A_1     =  Model.Ka( 1 );   
K_A_2_m =  Model.Ka( 2 );   
K_A_3_m =  Model.Ka( 3 );   
K_B_1      =  Model.Ka( 4 );   
 
k_2_A_1      = Model.k2( 1 );
k_2_A_2_m = Model.k2( 2 );
k_2_A_3_m = Model.k2( 3 );
k_2_B_1      = Model.k2( 4 );
k_2_B_2_m = Model.k2( 5 );
k_2_B_3      = Model.k2( 6 );




%-- 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 
            
            error('B_R2_R2L2 model does not have analytical solution. Use numeric mode!');
             
            
            
    
    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_R2_R2L2_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_B_2 = (2*K_A_2_m*K_B_1)/K_A_1;                  %   eq4_9
 K_B_2_m = K_B_2/2;                                           % eq3_4a;
 
 K_B_3 = (K_A_2_m*K_A_3_m*K_B_1)/K_A_1^2 ; %  eq4_10

% save dependent constants for reporting
Model.Ka_dependent_names={ 'B2'   'B3' };
Model.Ka_dependent = [ K_B_2  K_B_3 ]; 

%-- compute concentrations of  other species

%  eq4_7
Req = -(K_A_1*Leq^(1/2) - (K_A_1^2*Leq - 8*K_A_2_m*K_B_1*Leq + 8*K_A_2_m*K_B_1*Ltot - 8*K_A_2_m*K_A_3_m*K_B_1*Leq^2 + 8*K_A_2_m*K_A_3_m*K_B_1*Leq*Ltot)^(1/2))/(4*K_A_2_m*K_B_1*Leq^(1/2) + 4*K_A_2_m*K_A_3_m*K_B_1*Leq^(3/2))  ;
 

RLeq = K_A_1*Leq*Req    ; %  eq4_1


R2eq = K_B_1*Req^2  ;  %  eq4_2


R2Leq = 2*K_A_2_m*K_B_1*Leq*Req^2  ; %  eq4_3
% switch to microscopic species in the symmetrical  system
RaLRbeq=R2Leq/2;
RaRbLeq=RaLRbeq;


R2L2eq = K_A_2_m*K_A_3_m*K_B_1*Leq^2*Req^2  ; %  eq4_4



%% Line shape equations

% Equilibrium mole fractions of species containting P 
% NOTE: make sure to count dimers twice because we are counting nuclei!
%           'R'                'R2'              'RL'                'RaLRb'                'RaRbL'                'R2L2' 
p=[        Req/Rtot,   2*R2eq/Rtot,    RLeq/Rtot,  2*RaLRbeq/Rtot,  2*RaRbLeq/Rtot,  2*R2L2eq/Rtot ];
%p=[        Req/Rtot,   R2eq/Rtot,    RLeq/Rtot,  RaLRbeq/Rtot,  RaRbLeq/Rtot,  R2L2eq/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      =k_2_A_1     *K_A_1;
k_1_A_2_m =k_2_A_2_m*K_A_2_m;
k_1_A_3_m =k_2_A_3_m*K_A_3_m;
k_1_B_1     =k_2_B_1      *K_B_1;
k_1_B_2_m=k_2_B_2_m *K_B_2_m;
k_1_B_3     =k_2_B_3      *K_B_3;


% save them in the order of Model.k_names !!!
Model.k1 = [ k_1_A_1   k_1_A_2_m   k_1_A_3_m   k_1_B_1    k_1_B_2_m  k_1_B_3];
 

% Write matrix elements for the rate equations for the system
%                                  'R'                                                                                 'R2'                                             'RL'                                                                                                          'RaLRb'                                                             'RaRbL'                                                'R2L2' 
%  d[R]/dt     =   -(k_1_A_1*Leq+k_1_B_1*Req)                                                   k_2_B_1            k_2_A_1                                                                                                                           0                                                                     0                                                         0
                       
%  d[R2]/dt   =                              k_1_B_1*Req                 -(2*k_1_A_2_m*Leq+k_2_B_1)                                  0                                                                                                                         k_2_A_2_m                                                  k_2_A_2_m                                 0
                                 
%  d[RL]/dt   =      k_1_A_1*Leq                                                                           0                      -(k_2_A_1+k_1_B_3*RLeq+2*k_1_B_2_m*Req)                                                                           k_2_B_2_m                                                  k_2_B_2_m                                       k_2_B_3     

%  d[RaLRb]/dt =               0                                                     k_1_A_2_m*Leq                                                                            k_1_B_2_m*Req                        -(k_1_A_3_m*Leq+k_2_A_2_m+k_2_B_2_m)                                       0                                      k_2_A_3_m  

%  d[RaRbL]/dt =               0                                                     k_1_A_2_m*Leq                                                                            k_1_B_2_m*Req                                                                 0                           -(k_1_A_3_m*Leq+k_2_A_2_m+k_2_B_2_m)         k_2_A_3_m                                                                                                       

%  d[R2L2]/dt =                 0                                                                                  0                                         k_1_B_3*RLeq                                                              k_1_A_3_m*Leq                                         k_1_A_3_m*Leq                                               -(2*k_2_A_3_m+k_2_B_3)                   



% Copy constants to a K-matrix
K=[    
    -(k_1_A_1*Leq+k_1_B_1*Req)                                                   k_2_B_1            k_2_A_1                                                                                                                           0                                                                     0                                                         0
                       
                               k_1_B_1*Req                 -(2*k_1_A_2_m*Leq+k_2_B_1)                                  0                                                                                                                         k_2_A_2_m                                                  k_2_A_2_m                                 0
                                 
      k_1_A_1*Leq                                                                           0                      -(k_2_A_1+k_1_B_3*RLeq+2*k_1_B_2_m*Req)                                                                           k_2_B_2_m                                                  k_2_B_2_m                                       k_2_B_3     

                0                                                     k_1_A_2_m*Leq                                                                            k_1_B_2_m*Req                        -(k_1_A_3_m*Leq+k_2_A_2_m+k_2_B_2_m)                                       0                                      k_2_A_3_m  

                0                                                     k_1_A_2_m*Leq                                                                            k_1_B_2_m*Req                                                                 0                           -(k_1_A_3_m*Leq+k_2_A_2_m+k_2_B_2_m)         k_2_A_3_m                                                                                                       

                  0                                                                                  0                                         k_1_B_3*RLeq                                                              k_1_A_3_m*Leq                                         k_1_A_3_m*Leq                                               -(2*k_2_A_3_m+k_2_B_3)                   
 
    ];


 


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