function fraction=aminofr(pHx,pK1x,pK2x)
%-----------------------------------------------------------------------
% Programming Note:
%   Call fsolve to find the roots to a coupled set of nonlinear algebraic equations.
%-----------------------------------------------------------------------
% Instructor: Nam Sun Wang
%-----------------------------------------------------------------------

% Pass pH, pK1, and pK2 values
      global pH pK1 pK2

% Place the function argument onto "global" variables to pass them to where these numbers are needed.
      pH  = pHx;
      pK1 = pK1x;
      pK2 = pK2x;

% Provide an initial guess of the three fractions of A (It seems reasonable to divide evenly among the three.)
      fraction_guess=[1/3 1/3 1/3]';

% Call a routine to solve for the root of f(x)=0 -----------------------
      fraction = fsolve(@aminofcn, fraction_guess);
end


%-----------------------------------------------------------------------
function f = aminofcn(x)
%-----------------------------------------------------------------------
% Specify the three coupled algebraic equations arising from the zwitterion
% behavior of an amino acid or other chemical species that can present in
% an aqueous solution in several forms.
%-----------------------------------------------------------------------
% Zwitterion behavior of an amino acid.
% An amino acid acts as a zwitterion, i.e., it can be either:
%   1) a positively charged cation (A+) in an acid solution,
%   2) a negatively charged anion (A-) in an alkaline solution, or
%   3) a neutral molecule (A) at the isoelectric point.
% Zwitterion Reaction:
%   NH3+-CHR-COOH --> NH3+-CHR-COO- + H+
%   NH3+-CHR-COO- --> NH2-CHR-COO-  + H+
% which can be abbreviated as:
%   A+ --> A  + H+     K1
%   A  --> A- + H+     K2
% The dissociation constants for these steps are K1 and K2:
%   K1=10^(-pK1)=[A]*[H+]/[A+]     (1)
%   K2=10^(-pK2)=[A-]*[H+]/[A]     (2)
% Conservation of mass (All three possible fractions must add up to unity):
%   [A+]+[A]+[A-]=1                (3)
% This function specifies eqns (1)-(3)
%-----------------------------------------------------------------------
% Instructor: Nam Sun Wang
%-----------------------------------------------------------------------

% Pass pH, pK1, and pK2 values
      global pH pK1 pK2

% Create a vector f having the same shape as x
      f    = zeros(size(x));

% Use my own variables
      Acation = x(1);
      A       = x(2);
      Aanion  = x(3);

% Definition of pH & pK
      H  = 10^(-pH);
      K1 = 10^(-pK1);
      K2 = 10^(-pK2);

% Specify three equations f (good for pH1--12)
%     f(1) = A - K1/H*Acation;     % Definition for K1 -- o.k.
      f(1) = A*H - K1*Acation;     % Definition for K1 -- o.k.
      f(2) = Aanion/A  - K2/H;     % Definition for K2 -- o.k.
      f(3) = Acation + A + Aanion - 1.;   % Conservation of mass

%-----------------------------------------------------------------------
% Substituting one of the functions above with one of following is no good.
%     f(1) = A/Acation - K1/H;     % Definition for K1 -- no good, divide by almost 0 at high pH
%     f(2) = Aanion*H/A  - K2;     % Definition for K2 -- no good at high pH
%     f(2) = Aanion*H  - K2*A;     % Definition for K2 -- no good at high pH
%-----------------------------------------------------------------------
% Note: The following combinations are good for only a subrange
% (Lession: avoid dividing by a very small number.)
% o.k. for pH=1--6
%     f(1) = A*H/Acation - K1;     % Definition for K1
%     f(2) = Aanion*H/A  - K2;     % Definition for K2
% o.k. for pH=1--7
%     f(1) = A*H - K1*Acation;     % Definition for K1
%     f(2) = Aanion*H  - K2*A;     % Definition for K2
% o.k. for pH=9--12
%     f(1) = Acation/(A*H) - 1/K1; % Definition for K1
%     f(2) = A/(Aanion*H)  - 1/K2; % Definition for K2
end