function [p, K, R2_vector, v_vector, p_other, Model]=U_R_RL_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
%
% U-R-RL model developed in 
% [ek_mac]: /Users/kovrigin/Documents/Workspace/Data/Data.XV/
% EKM16.Analysis_of_multistep_kinetic_mechanisms/Equilibria/U_R_RL_derivation.mn 



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

% Current model:  U-R-RL 

%--- assign model parameters corresponding to the order 
%    in 'model_setup.m' 
 
K_A_1 =  Model.Ka(1);   
K_B_1 =  Model.Ka(2); 
K_B_2 =  Model.Ka(3); 

k_2_A_1 = Model.k2(1);
k_2_A_2 = Model.k2(2);
k_2_B_1 = Model.k2(3);
k_2_B_2 = 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 
            
            % eq 3.4
            Leq = -(K_B_2 - (K_A_1^2*K_B_2^4*Ltot^2 - 2*K_A_1^2*K_B_2^4*Ltot*Rtot + K_A_1^2*K_B_2^4*Rtot^2 + 2*K_A_1^2*K_B_2^3*Ltot^2 - 4*K_A_1^2*K_B_2^3*Ltot*Rtot + 2*K_A_1^2*K_B_2^3*Rtot^2 + K_A_1^2*K_B_2^2*Ltot^2 - 2*K_A_1^2*K_B_2^2*Ltot*Rtot + K_A_1^2*K_B_2^2*Rtot^2 + 2*K_A_1*K_B_1*K_B_2^3*Ltot + 2*K_A_1*K_B_1*K_B_2^3*Rtot + 2*K_A_1*K_B_1*K_B_2^2*Ltot + 2*K_A_1*K_B_1*K_B_2^2*Rtot + 2*K_A_1*K_B_2^3*Ltot + 2*K_A_1*K_B_2^3*Rtot + 2*K_A_1*K_B_2^2*Ltot + 2*K_A_1*K_B_2^2*Rtot + K_B_1^2*K_B_2^2 + 2*K_B_1*K_B_2^2 + K_B_2^2)^(1/2) + K_B_1*K_B_2 - K_A_1*K_B_2*Ltot + K_A_1*K_B_2*Rtot - K_A_1*K_B_2^2*Ltot + K_A_1*K_B_2^2*Rtot)/(2*K_A_1*K_B_2^2 + 2*K_A_1*K_B_2);
            
             
            
            
    
    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_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',...
                  'LevenbergMarquardt', '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;   % eq 3.2



%-- compute concentrations of  other species


Req  = Rtot/(K_B_1 + K_A_2*K_B_1*Leq + (K_A_2*K_B_1*Leq)/K_B_2 + 1);  % eq  2.11

RLeq = (K_A_2*K_B_1*Leq*Req)/K_B_2  ; % eq 2.9


Rstareq = K_B_1*Req ;  % eq 2.8 | eq 2.9


RstarLeq = K_A_2*K_B_1*Leq*Req ; % eq2_7 | eq2_9



%% Line shape equations

% Equilibrium mole fractions of species containting P 
%                R             R*               RL              R*L
p=[      Req/Rtot,   Rstareq/Rtot,  RLeq/Rtot,  RstarLeq/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=k_2_A_2*K_A_2;
k_1_B_1=k_2_B_1*K_B_1;
k_1_B_2=k_2_B_2*K_B_2;

% Write matrix elements for the rate equations for the system
%                                        R                                                  R*                                  RL                                   R*L 
%  dReq/dt     =    -(k_1_B_1+k_1_A_1*Leq)                              +k_2_B_1         +k_2_A_1                                      0
%  dR*eq/dt   =       k_1_B_1                               -(k_1_A_2*Leq+k_2_B_1)                       0                      k_2_A_2
%  dRLeq/dt   =                   +k_1_A_1*Leq                                 0                     -(k_2_A_1+k_1_B_2)                        k_2_B_2      
%  dR*Leq/dt =                   0                                k_1_A_2*Leq                                          k_1_B_2        -(k_2_A_2+k_2_B_2)      

% Copy constants to a K-matrix
K=[     -(k_1_B_1+k_1_A_1*Leq)                              +k_2_B_1         +k_2_A_1                                      0;...
             k_1_B_1                               -(k_1_A_2*Leq+k_2_B_1)                       0                      k_2_A_2     ;...
                        +k_1_A_1*Leq                                 0                     -(k_2_A_1+k_1_B_2)                        k_2_B_2 ;...                  
                        0                                k_1_A_2*Leq                                          k_1_B_2        -(k_2_A_2+k_2_B_2)
     ];


 


%-- 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'  };
Model.Ka_dependent = [ K_A_2  ];
% save forward rate constant values in the order of Model.k_names !!!
Model.k1 = [    k_1_A_1  k_1_A_2   k_1_B_1   k_1_B_2  ]; 

