function [p, K, R2_vector, v_vector, p_other, Model]=U_R_RL_RM_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
%
% Current model  U-R-RL-RM   developed in 
%  /Users/kovrigin/Documents/Workspace/Data/Data.XV/
% EKM16.Analysis_of_multistep_kinetic_mechanisms/Equilibria/U_R_RL_RM_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_3 =  Model.Ka( 2);   
K_B_1 =  Model.Ka( 3); 
K_B_2 =  Model.Ka( 4); 
K_B_3 =  Model.Ka( 5); 

k_2_A_1 = Model.k2( 1);
k_2_A_2 = Model.k2( 2);
k_2_A_3 = Model.k2( 3);
k_2_A_4 = Model.k2( 4);
k_2_B_1 = Model.k2( 5);
k_2_B_2 = Model.k2( 6);
k_2_B_3 = Model.k2( 7);




%-- assign experimental parameters:
Rtot = Rtotal;    
Ltot = Rtot*LRratio;  
Mtot = Rtot*Model.MRratio;

 

%% Determine equilibrium concentrations of species

switch mode
    
    %---- traditional explicit calculation of [L]
    case 'analytical'
            %% Analytical solution for [L] out of 
            
            error('Model U_R_RL_RM does not have analytical solution! Use numeric mode instead.'); 
            
            
    
    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) ( U_R_RL_RM_model_numeric(B, [X  Model.Ka(2:end) Ltot Mtot]) );
            
            % 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 = (K_A_1*K_B_2)/K_B_1 ;   % eq3.12  

K_A_4 = (K_A_3*K_B_3)/K_B_1 ;   % eq  3.11



%-- compute concentrations of  other species

% eq 3.9   
Meq = (K_A_1*Leq*Mtot + K_A_1*K_B_2*Leq*Mtot)/(K_A_1*Leq - K_A_3*Leq + K_A_3*Ltot + K_A_1*K_B_2*Leq - K_A_3*K_B_3*Leq + K_A_3*K_B_3*Ltot);
  

% eq   3.7
Req = -(Leq - Ltot)/(K_A_1*Leq + K_A_1*K_B_2*Leq) ;  

% eq  3.5
RMeq = K_A_3*Meq*Req  ; 

% eq  3.4
RLeq = K_A_1*Leq*Req ; 

% eq  3.3
Rstareq = K_B_1*Req ;  

% eq 3.2
RstarMeq = K_B_3*RMeq; 

% eq 3.1
RstarLeq = K_B_2*RLeq ; 



%% Line shape equations

% Equilibrium mole fractions of species containting P 
%                     R        R*             RL               R*L             RM                  R*M
p=[      Req/Rtot,   Rstareq/Rtot,  RLeq/Rtot,  RstarLeq/Rtot,  RMeq/Rtot,  RstarMeq/Rtot];

% Normalized concentration of the ligand for plotting
p_other=[Leq/Rtot,  Meq/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=k_2_A_2*K_A_2;
k_1_A_3=k_2_A_3*K_A_3;
k_1_A_4=k_2_A_4*K_A_4;

k_1_B_1=k_2_B_1*K_B_1;
k_1_B_2=k_2_B_2*K_B_2;
k_1_B_3=k_2_B_3*K_B_3;

% Write matrix elements for the rate equations for the system
%                                        R                                                            R*                                                                   RL                                   R*L                                     RM                                    R*M 
%  dReq/dt     =    -(k_1_B_1+k_1_A_1*Leq+k_1_A_3*Meq)       k_2_B_1                                                        k_2_A_1                               0                              k_2_A_3                                        0
%  dR*eq/dt   =       k_1_B_1                                                     -(k_2_B_1+k_1_A_2*Leq+k_1_A_4*Meq)             0                                k_2_A_2                                    0                        k_2_A_4      
%  dRLeq/dt   =                       k_1_A_1*Leq                                     0                                                            -(k_2_A_1+k_1_B_2)                          k_2_B_2                    0                                       0   
%  dR*Leq/dt =                    0                                                                       k_1_A_2*Leq                                                k_1_B_2        -(k_2_A_2+k_2_B_2)                   0                                       0                      
%  dRMeq/dt   =                                             k_1_A_3*Meq             0                                                                      0                                      0                          -(k_2_A_3+k_1_B_3)                         k_2_B_3
%  dR*Meq/dt =                    0                                                                                              k_1_A_4*Meq              0                                      0                                             k_1_B_3       -(k_2_A_4+k_2_B_3)                          

% Consistency checking: 
% 1. each row must have complete pairs of constants so you for every forward constant 
%     you should have reverse one as well (complete cycle rule). 
% 2. each column must sum to 0 (mass conservation law). 

% Copy constants to a K-matrix
K=[ ...   
  -(k_1_B_1+k_1_A_1*Leq+k_1_A_3*Meq)       k_2_B_1                                                        k_2_A_1                               0                              k_2_A_3                                        0                  ;...
     k_1_B_1                                                     -(k_2_B_1+k_1_A_2*Leq+k_1_A_4*Meq)             0                                k_2_A_2                                    0                        k_2_A_4                      ;...
                        k_1_A_1*Leq                                     0                                                            -(k_2_A_1+k_1_B_2)                          k_2_B_2                    0                                       0               ;...                  
                    0                                                                       k_1_A_2*Leq                                                k_1_B_2        -(k_2_A_2+k_2_B_2)                   0                                       0                ;...
                                              k_1_A_3*Meq             0                                                                      0                                      0                          -(k_2_A_3+k_1_B_3)                         k_2_B_3  ;...                       
                     0                                                                                              k_1_A_4*Meq              0                                      0                                             k_1_B_3       -(k_2_A_4+k_2_B_3) ;...                       
      ];



%-- prepare relaxation rates and frequencies for the individual species
 R2_vector= Model.R2;  
 v_vector=  Model.w0;  
 
 
% Output of dependent constants
Model.Ka_dependent_names={ 'A2'  'A4'  };
Model.Ka_dependent = [  K_A_2   K_A_4 ];
% save forward rate constant values in the order of Model.k_names !!!
Model.k1 = [ k_1_A_1 k_1_A_2 k_1_A_3 k_1_A_4  k_1_B_1 k_1_B_2 k_1_B_3 ]; 
 
