%------------------- LineshapeKin Simulation -----------------------------
%-- version 4.0 
%
% 5/8/09
% Written by Evgenii Kovrigin
% Medical College of Wisconsin
%
% 
function R2 = R2_Alltimescale_Twosite(CPMG_frequency_array, Model, species_indices, p, kex, field_strength)
%--- Two-site All timescale equation for R2 from {Palmer, 2001 #1779} 
%     page 214
%  CPMG_frequency_array - array of frequencies for a smooth curve
%  Model - model object
%  species_indices - array of two indices corresponding to species in the Model
%                           with major populated species listed first in this array
% p - vector of populations for all species in the system
% kex - exchange rate computed for  the current conditions
% field_strength - scaling factor for dw used for simulation of data at different
%                          field strengths. Use 1 for 500, 1.2 for 600 etc.


% calculate tau_cp
tcp=0.5./CPMG_frequency_array;

% by convention major conformer is which index is listed first.
pa=p(species_indices(1));


% assign relaxation rates of major and minor species
R20_a=Model.R2(species_indices(1)); 
R20_b=Model.R2(species_indices(2)); 


% set array of dWs
dw=(Model.w0(species_indices(2)) - Model.w0(species_indices(1)))*field_strength;






pb = 1 - pa;




zeta = 2 .* dw .* ( -pa .* kex + pb .* kex );

psi = ( -pa .* kex + pb .* kex ).^2 - dw.^2  +  4 .* pa .* pb .* kex.^2;

eta_plus =  tcp .* sqrt( + psi + sqrt( psi.^2 + zeta.^2 ) ) ./ sqrt(2.0);
eta_minus = tcp .* sqrt( - psi + sqrt( psi.^2 + zeta.^2 ) ) ./ sqrt(2.0);

D_plus  = 0.5 .* ( + 1 + ( psi + 2.*dw.^2) ./ sqrt( psi.^2 + zeta.^2 ) );
D_minus = 0.5 .* ( - 1 + ( psi + 2.*dw.^2) ./ sqrt( psi.^2 + zeta.^2 ) );


R2 =  0.5 .* (R20_a + R20_b + kex - acosh( D_plus .* cosh(eta_plus) - D_minus .* cos(eta_minus)  )./tcp);



