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


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

  
%--- assign model parameters
 
K_a_A =  Model.Ka(1);   
K_a_B =  Model.Ka(2); 

k_2_A = Model.k2(1);
k_2_B = Model.k2(2);




%-- 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('Analytical solution for U-L available but not implemented!');
            
    
    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);
            
            % this is Y-array 
            y_local=Rtot;

             %-- physically meaningful bounds
             lower_bounds=[0];  
              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);
             
             
             %  --- 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_L2_model_numeric(B, [X  Model.Ka(2) Ltot]) );
             
             % --- 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 concentrations of other other species

RLeq = Ltot - Leq - 2*K_a_B*Leq^2; % eq2.7

Req = RLeq/(K_a_A*Leq);  % eq2.6

LLeq = K_a_B*Leq^2;  % eq2.5





%% Line shape equations

% Equilibrium mole fractions of species containting P 
%                     R         RL
p=[      Req/Rtot,     RLeq/Rtot];
p_other=[Leq/Rtot, 2*LLeq/Rtot];   % normalized concentrations for plotting

% Notation for transitions
% A: R + L <=> RL 
% B: R <=> R*
%-- express forward rate constants through defined reverse rate constants
k_1_A=k_2_A*K_a_A;
k_1_B=k_2_B*K_a_B;

% write rate equations to generate a K-matrix
%                               R                             RL
%  dReq/dt =     - k_1_A*Leq                    k_2_A
%  dRLeq/dt =      k_1_A*Leq                   -k_2_A

% Copy constants to a K-matrix
K=[  - k_1_A*Leq                    k_2_A;...
           k_1_A*Leq                   -k_2_A
    ];



%-- 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={   };
Model.Ka_dependent = [   ];
% save forward rate constant values in the order of Model.k_names !!!
Model.k1 = [ k_1_A  k_1_B ]; 
