%-----------------------------------------------------------------------
% Plot the fractions of three forms of a chemical species with pH.
% Amino Acid Double Dissociation.  An amino acid becomes a
% positively charged cation (A+) in an acid solution; it becomes
% a negatively charged anion (A-) in an alkaline solution.  The
% dissociation constants for these steps are K1 and K2, which are
% defined as:
%   A+ --> A  + H+    K1 = 10-pK1 = A*H/A+   ... eqn (1)
%   A  --> A- + H+    K2 = 10-pK2 = A-*H/A   ... eqn (2)
%   A+ +  A + A- = 1                         ... eqn (3)
%  For glutamic acid pK1=2.10 and pK2=9.47.
% Progamming Note:  Linear algebraic equations
%-----------------------------------------------------------------------
% Instructor: Nam Sun Wang
%-----------------------------------------------------------------------

% Start fresh ----------------------------------------------------------
      clear, clc

% Program header -------------------------------------------------------
      disp('This program finds the fractions of the three different forms')
      disp('when the chemical species, e.g., amino acid, can double dissociate.')

% Dissociation constants for glutamic acid -----------------------------
% Definition of pH & pK
      pK1 = 2.10;  K1 = 10^(-pK1);
      pK2 = 9.47;  K2 = 10^(-pK2);
%                  H  = 10^(-pH);

% Solve a set of linear algebraic equations of the form A*x=b
%     K1*Acation - H*A            = 0  % Definition for K1
%                 K2*A - H*Aanion = 0  % Definition for K2
%        Acation   + A   + Aanion = 1  % Conservation of mass
      A=@(pH) [ K1 -10^(-pH)  0
                0    K2      -10^(-pH)
                1    1        1        ];
      b = [ 0; 0; 1];
      fraction=@(pH) A(pH)\b;
      Acation =@(pH) [1 0 0]*fraction(pH);
      Aneutral=@(pH) [0 1 0]*fraction(pH);
      Aanion  =@(pH) [0 0 1]*fraction(pH);
% The following xi function returns ith element of vector x
%     xi=@(x,i) x(i);  Aneutral=@(pH) xi(fraction(pH),2);

% Range pH and find the fractions of the three different forms ---------
      pH = (0:0.2:12)';    % Finer divisions for plotting
      Acationv =arrayfun(Acation,  pH);
      Aneutralv=arrayfun(Aneutral, pH);
      Aanionv  =arrayfun(Aanion,   pH);

% Print a table --------------------------------------------------------
      disp(' ')
      disp('            ----------Fractions---------')
      disp('     pH      Acation      A       Aanion')
      disp('  --------------------------------------')
      disp([pH Acationv Aneutralv Aanionv])
      
% Plot and add legends -------------------------------------------------
      plot(pH,Acationv, pH,Aneutralv, pH,Aanionv)
      xlabel('pH')
      ylabel('Fraction')
      legend('Cation +', 'Neutral', 'Anion -')
      title('Double Dissociation of Glutaminic Acid')
      text(pK1-1., 0.05, ['pK_1=' num2str(pK1)])
      text(pK2-1., 0.05, ['pK_2=' num2str(pK2)])