function Rtot = B_R2_R2L2_model_numeric(params, constants )
%------------------- LineshapeKin Simulation -----------------------------
% Function for numeric solving of equilibrium thermodynamic equations
%
% 5/12/09
% Evgenii Kovrigin
% Medical College of Wisconsin
%
% 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 


%--- this is what I need for K-matrix
Leq=params(1);


%-- these are constants named as used in the equation below
K_A_1      =constants( 1 );
K_A_2_m =constants( 2 );
K_A_3_m =constants( 3 );
K_B_1     =constants( 4 );
Ltot         =constants( 5 );



%-- equation  4.8

Rtot = (2*K_B_1*(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))^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))^2 - (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)) - (K_A_1*Leq*(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)) + (4*K_A_2_m*K_B_1*Leq*(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))^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))^2 + (2*K_A_2_m*K_A_3_m*K_B_1*Leq^2*(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))^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))^2;










