function [p, K, R2_vector, v_vector, p_other, Model]=U_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: Two-site exchange -------------------
%   
%      P+L=PL, Kd=k2/k1, where k1=kon*[L], k2=Koff
%
%  Here 'numeric' mode is 4x slower than 'analytical'



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

  
%--- extract model parameters
 
Kd = 1/Model.Ka(1);   
koff=Model.k2(1);




%-- extract experimental parameters:


Ptot = Rtotal;    
Ltot = Ptot*LRratio;  


 
%% Analytical solution for [L] out of 
%
%   Kd = [P][L]/[PL]


switch mode
    
    %---- traditional explicit calculation of [L]
    case 'analytical'
            a=1;
            b=-(Ptot+Ltot+Kd);
            c=Ptot*Ltot;
            sola=-b/(2*a);
            solb=sqrt(b^2-(4*a*c))/(2*a);
            sol1=[(sola+solb) (sola-solb)];
            Peq=max(Ptot-sol1); %-- get positive solution
            PLeq=Ptot-Peq;
            Leq=Ltot-PLeq;
            
            
            
            
    %--- minimization to obtain [L]
    case 'numeric'

            

            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. This is just one of the constants in
            % the equation. The other constants are supplied directly to
            % the function below.
            x_local=Kd;
            
            % this is Y-array 
            y_local=Ptot;

             %-- 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',...
                  'LevenbergMarquardt', 'on',...
                  'TolFun',1e-18,...
                  'TolX',1e-18,...
                 'MaxFunEvals', 1e9);
             
             
             %  --- we define the target function using 'anonymous' synthax
             %  and pass additional parameters to it here
             %myfun = @(B,X) ( lnKd1fun(B, [Ptot,X]) );
             myfun = @(B,X) ( U_model_numeric(B, [X,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
            
            %-- compute concentrations of other other species
            Peq= Ptot-Ltot+Leq;
            PLeq=Ptot-Peq;

end



% Equilibrium mole fractions of species containting P 
p=[Peq/Ptot,PLeq/Ptot];
p_other=[Leq/Ptot];   % normalized concentration for plotting

%-- relationships of true constants
kon=koff/Kd;

% K matrix elements
kab=kon*Leq;
kba=koff;

K=[-kab,kba;kab,-kba];



%-- 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 them in the order of Model.k_names !!!
Model.k1 = [ kon ];
 
