Maxima : intro computation ,base and acid, 'appendix 2,version 1.0'
peter.vlasschaert@gmail.com, 16/12/2018

kill(all)$

Load 'packages'
==============
rem :factorout,factorsum,lratsubst

load(lrats)$
load(facexp)$

convert base : 10
--------------------------

p2:log10(x):= log(x)/log(10);

definition : pH,pOH
================

p1: pH(x)   :=-log10(x)$
p1a:pOH(x):=-log10(x)$

Distribution Diagrams of Polyprotic Acids
===================================

acetic acid              :  pK_a1 = 4.74    , mono
===================================
hydrogen sulfide  :  pK_a1 = 7.00    ,poly
                                      pK_a2 = 14.00
phosphoric acid  :   pK_a1 = 2.12   , poly
                                      pK_a2 = 7.21
                                      pK_a3 = 12.00
H4-EDTA                    pK_a1 = 1.99    ,poly
                                      pK_a2 = 2.67
                                      pK_a3 = 6.16
                                      pK_a4 =10.26
===================================

 ==================================

Weak acid: H(n)A    (example :  n=2)
===============

reaction : H(2)A+H(2)O        ⇔ H(3)O(+) + HA(-)           ....  '1'
reaction : HA(-)+H(2)O         ⇔ H(3)O(+) + A(2-)            ....  '1a'
reaction : H(2)O+H(2)O        ⇔ H(3)O(+) + OH(-)           .... '1b'
--------------------------------------------------------------------------------------
rem:H(3)O(+) =H, A(2-) = A,OH(-)=OH,HA(-)=HA H(2)A=H2A
--------------------------------------------------------------------------------------

definition : pKai    i=1,2,...,n
======================

p2:pKa(x)   :=-log10(x)$

definition : K_a i (i=1,2,..n)     .... '2'
============================

p2a:K_a1 = ((H)*(HA))/(H2A);
p2b:K_a2 = ((H)*(A))/(HA);

p2c:p2a*p2b;

=>     : H2A ⇔ 2*H(+) + A(2-)  for 'p2c'

consecutive dissociation constants 
==============================

general:   K_a1 > K_a2> ......> K_an 

C_H2A = polyprotic  Acid,(n=2)  ... '3'
======================================

p3:C_H2A=(H2A)+(HA)+(A);

general charge balance
====================

a : charge balance :  '1' and '1a'
==========================

[H(+)] = [HA(-)]+[A(2-)]

p4:H = HA+A;

b : charge balance :  '1' and '1a' and '1b'
==================================

[H(+)] = [HA(-)]+[A(2-)]+[OH(-)]

p4a:H = HA+A+OH;

distribution : polyprotic acids :
=========================

n:2;

p5:sum(α(k), k, 0, n)=1;

p6a:α(0) = H2A/rhs(p3);
p6b:α(1) = HA/rhs(p3);
p6c:α(2) = A/rhs(p3);

check : for 'p5'
============

p6:rhs(p6a)+rhs(p6b)+rhs(p6c);

p7:ratsimp(p6);

computation: from 'p2a,p2b,p2c'
===========================

find α(0):
=======

p8a:solve(p2a,HA)[1];
p8b:solve(p2c,A)[1];

p8c:ratsimp(subst([p8a,p8b],p6a));

from 'p8c'
========

p8c1:num(rhs(p8c));
p8c2:denom(rhs(p8c));

how to find 'p8c2' : from 'p3,p2a,p2b'
===============================

p8c3:eliminate([rhs(p3),p2a,p2b],[HA,A]);

p8c4:part(p8c3[1],1,2);

find α(1):
=======

p9a:solve(p2a,H2A)[1];
p9b:solve(p2b,A)[1];

p9c:ratsimp(subst([p9a,p9b],p6b));

find α(2):
=======

p10a:solve(p2b,HA)[1];
p10b:solve(p2c,H2A)[1];

p10c:ratsimp(subst([p10a,p10b],p6c));

check : sum ' p5'  ? = 1
==============

p11:rhs(p8c+p9c+p10c);

p11a:ratsimp(p11);

rem :   log10(100) = 2  ⇔  100 = 10^2
analog :                               H = 10^(-ph)
-------------------------------------------------------
ph = hydrogen pressure
-------------------------------------------------------

p12:[H = 10^(-ph)];
p12a:subst(p12,p8c);
p12a1:rhs(p12a);
p12b:subst(p12,p9c);
p12b1:rhs(p12b);
p12c:subst(p12,p10c);
p12c1:rhs(p12c);

α(0)  => f1(ph,K_a1,K_a2)
α(1)  => f2(ph,K_a1,K_a2)
α(2)  => f3(ph,K_a1,K_a2)
===========================
How to create functions : f1,f2,f3
===========================

f1(ph,K_a1,K_a2):=''p12a1;

f2(ph,K_a1,K_a2):=''p12b1;

f3(ph,K_a1,K_a2):=''p12c1;

example : f1(ph =5,K_a1=7,K_a2=14), hydrogen sulfide.

f1(5,4.76,6.12);

example : graphics
================
============================
succinic acid  :  K_a1 = 6.2*10^(-5)
                               K_a2 = 2.3*10^(-6) 
=============================
see : above : K_a1  >  K_a2

How to make a plot : 
=================

Menu : Plot: Plot2d
-----------------------------
-------------------------------------------------------------------------------------------------------------------------------------------------------
Expression(s) :f1(ph,6.2*10^(-5),2.3*10^(-6)),f2(ph,6.2*10^(-5),2.3*10^(-6)),f3(ph,6.2*10^(-5),2.3*10^(-6))
variable : x=> ph
range     :  -5 => 0
range     :   5 =>  14
variable : y => y
to put inside document : after excute  , plot2d => wxplot2d
-------------------------------------------------------------------------------------------------------------------------------------------------------

q:[f1(ph,6.2*10^(-5),2.3*10^(-6)),f2(ph,6.2*10^(-5),2.3*10^(-6)),f3(ph,6.2*10^(-5),2.3*10^(-6))];

wxplot2d(q, [ph,0,14], [y,0,1],
 [plot_format, xmaxima])$