function R2 = R2_2st_alltmscl_2fields( params, X)
%--- Two-site All timescale equation for R2 using two-field data.  ----
%
% 3/21/08 - swapped order of parameters
%                - uses X with column-oriented structure
global X_column



%-- extract taucp data and field dependence of dw



tcp=X(:,X_column);
field_dependence=X(:,2);


R20=15; %-- base R2 rate

kex=params.param(1); 
pa=params.param(2);  


%conversion coefficient
dw=params.param(3)*field_dependence;  



pb = 1 - pa;


% ---- IMPORTANT!!! For correct multi-field operation one MUST have all *
% and / operators replaced by .*, ./ and .^

%--- In this equation we consider R20a=R20b so R20a-R20b=0 ---
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 = R20 + 0.5 .* ( kex - acosh( D_plus .* cosh(eta_plus) - D_minus .* cos(eta_minus)  )./tcp);



