%------------------------------------
% Zwitterion behavior of an amino acid.
% Matlab version of amino.mcd
% See amino.mcd for instructions and comments.
%-----------------------------------------------------------------------
% Instructor: Nam Sun Wang
%-----------------------------------------------------------------------

% Start fresh ----------------------------------------------------------
      clear all

% Symbolic solution of 3 simultaneous algebraic equations
% Matlab Bug: The orders returned by "solve" are switched around (lexicographic-alphabetically ordered)!!
[Aanion Acation Aneutral] = solve(...
  'Aneutral*10^(-pH)/Acation=10^(-pK1)', ...  % 1st equation; definition for 1st dissociation constant K1: A+ --> A + H+
  'Aanion*10^(-pH)/Aneutral=10^(-pK2)' , ...  % 2nd equation; definition for 2nd dissociation constant K2: A --> A- + H+
  'Acation+Aneutral+Aanion=1'          , ...  % 3rd equation; Conservation of mass (i.e., all fractions add up to unity)
  'Acation,Aneutral,Aanion')                  % solve for 3 variables; however, Matlab seems to return variables in lexicographic (alphabetical) order
A=inline([Acation Aneutral Aanion])           % turn into a vector function A(pH,pK1,pK2)
%Acation=inline([Acation])                    % turn into a scalar functionn Acation(pH,pK1,pK2)

% Example: Alanine has the following dissociation constants:
      pK1=2.34; pK2=9.69;
      disp('In a neutral solution, the composition is: ')
        pH=7, A(pH,pK1,pK2)
      disp('In an acid solution, the composition is: ')
        pH=1, A(pH,pK1,pK2)
      disp('In an alkaline solution, the composition is: ')
        pH=11, A(pH,pK1,pK2)
%plot
      pH=[0.1:0.1:12]';
      plot(pH, A(pH,pK1,pK2))
      xlabel('pH')
      ylabel('Amino Acid Fraction')
      legend('A+', 'A', 'A-')