kill(all)$

Maxima : intro computation ,base and acid, 'part 15c,version 1.0'
peter.vlasschaert@gmail.com,31/12/2018

load packages :
------------------------

load(lrats)$
load(facexp)$

rem : load(lrats), use for command ' lratsubst'

convert:
------------

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

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

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

definition : pKa
=============

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

log plot for polyprotic weak acid. 'HnA'
================================
================================
x: ph , y: log(concentration)

example : n=2  => H2A
===================
see: ( n=1 , HA)
=============
Maxima : intro computation ,base and acid, 'part :15,15a,15b,version 1.0',
(weak acid,distribution,log10(ph))

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 : K_a i (i=1,2,..n)     .... '2'
============================

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

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

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

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

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

distribution : polyprotic acids :  α(i) , i=0,1,..,n
==========================

see:Maxima : intro computation ,base and acid, 'appendix 2,version 1.0'

p4a:α(0) = H2A/lhs(p3);
p4b:α(1) = HA/lhs(p3);
p4c:α(2) = A/lhs(p3);

α(i) =>   f(ph,C_H2A,K_ai)
=====================

p5:H^2+K_a1*H+K_a1*K_a2;

p5a:args(p5);

p6a:lhs(p4a) = p5a[3]/p5;
p6b:lhs(p4b) = p5a[2]/p5;
p6c:lhs(p4c) = p5a[1]/p5;

combination : 'p4a,b,c' and 'p6a,b,c'
==============================

p7a: rhs(p4a) = rhs(p6a);
p7b: rhs(p4b) = rhs(p6b);
p7c: rhs(p4c) = rhs(p6c);

solve : {H2A,HA,A} from above
==========================

p8a:solve(p7a,H2A)[1];
p8b:solve(p7b,HA)[1];
p8c:solve(p7c,A)[1];

 n = 1,  2 :regions  (Maxima : intro computation ,base and acid, 'part :15,version 1.0',(weak acid,mono,log10(ph)))
 n = 2,  3 :regions
 n = 3,  4 :regions
..............................

example (n=2) :  C_H2A,K_a1,K_a2
=============================

1e) first        region :   ph < pK_a1                  , log10 plot , log10(H2A)  = constant
2e) second region :   pK_a1 < ph < pK_a3 , log10 plot , log10(HA)    = constant
3e) thrid       region :   ph > pK_a2                  ,log10 plot , log10(A)       = constant
-----------------------------------------------------------------------------------------------------------------
constant = log10(C_H2A)

conversion : H => ph

p8d:H= 10^(-ph);

p8e:ph = -log10(H);

p8f:rhs(p8e) =lhs(p8e);

p9:p5a;

first region :
==========

p9a:[p9[1]=0,p9[2]=0];

p9a1:subst(p9a,p5);

second region :
=============

p9b:[p9[1]=0,p9[3]=0];;

p9b1:subst(p9b,p5);

thrid region :
===========

p9c:[p9[2]=0,p9[3]=0];;

p9c1:subst(p9c,p5);

[H2A] => H2A ,  'p10a3,p10b6,p10c6'
===========

1e) region
---------------

p10a:lhs(p8a) =num(rhs(p8a))/p9a1;

p10a1:'log10(lhs(p10a)) = log10(rhs(p10a));

p10a2:'log10(lhs(p10a)) =rhs(p10a1),logexpand=super;

p10a3:lhs(p10a2) = subst([p8f],rhs(p10a2));

p10a4:p10a3,numer;

2e) region
---------------

p10b:lhs(p8a) =num(rhs(p8a))/p9b1;

p10b1:'log10(lhs(p10b)) = log10(rhs(p10b));

p10b2:'log10(lhs(p10b)) =rhs(p10b1),logexpand=super;

p10b3:lhs(p10b2) = expand(rhs(p10b2));

p10b4:lhs(p10b3) = lratsubst([p8f],rhs(p10b3));

p10b5:lhs(p10b4) = logcontract(rhs(p10b4));

p10b6:lhs(p10b5) = expand(rhs(p10b5));

p10b7:p10b6,numer;

3e) region
---------------

p10c:lhs(p8a) =num(rhs(p8a))/p9c1;;

p10c1:'log10(lhs(p10c)) = log10(rhs(p10c));

p10c2:'log10(lhs(p10c)) =rhs(p10c1),logexpand=super;

p10c3:lhs(p10c2) = expand(rhs(p10c2));

p10c4:lhs(p10c3) = lratsubst([p8f],rhs(p10c3));

p10c5:lhs(p10c4) = logcontract(rhs(p10c4));

p10c6:lhs(p10c5) = expand(rhs(p10c5));

p10c7:p10c6,numer;

[HA(-)] => HA ,,  'p11a6,p11b6,p11c6'
===========

1e) region
---------------

p11a:lhs(p8b) =num(rhs(p8b))/p9a1;

p11a1:'log10(lhs(p11a)) = log10(rhs(p11a));

p11a2:'log10(lhs(p11a)) =rhs(p11a1),logexpand=super;

p11a3:lhs(p11a2) = expand(rhs(p11a2));

p11a4:lhs(p11a3) = lratsubst([p8f],rhs(p11a3));

p11a5:lhs(p11a4) = logcontract(rhs(p11a4));

p11a6:lhs(p11a5) = expand(rhs(p11a5));

2e) region
---------------

p11b:lhs(p8b) =num(rhs(p8b))/p9b1;

p11b1:'log10(lhs(p11b)) = log10(rhs(p11b));

p11b2:'log10(lhs(p11b)) =rhs(p11b1),logexpand=super;

p11b3:lhs(p11b2) = expand(rhs(p11b2));

p11b4:lhs(p11b3) = lratsubst([p8f],rhs(p11b3));

p11b5:lhs(p11b4) = logcontract(rhs(p11b4));

p11b6:lhs(p11b5) = expand(rhs(p11b5));

3e) region
---------------

p11c:lhs(p8b) =num(rhs(p8b))/p9c1;;

p11c1:'log10(lhs(p11c)) = log10(rhs(p11c));

p11c2:'log10(lhs(p11c)) =rhs(p11c1),logexpand=super;

p11c3:lhs(p11c2) = expand(rhs(p11c2));

p11c4:lhs(p11c3) = lratsubst([p8f],rhs(p11c3));

p11c5:lhs(p11c4) = logcontract(rhs(p11c4));

p11c6:lhs(p11c5) = expand(rhs(p11c5));

[A(2-)] => A     ,  'p12a6,p12b6,p12c5'
===========

1e) region
---------------

p12a:lhs(p8c) =num(rhs(p8c))/p9a1;

p12a1:'log10(lhs(p12a)) = log10(rhs(p12a));

p12a2:'log10(lhs(p12a)) =rhs(p12a1),logexpand=super;

p12a3:lhs(p12a2) = expand(rhs(p12a2));

p12a4:lhs(p12a3) = lratsubst([p8f],rhs(p12a3));

p12a5:lhs(p12a4) = logcontract(rhs(p12a4));

p12a6:lhs(p12a5) = expand(rhs(p12a5));

2e) region
---------------

p12b:lhs(p8c) =num(rhs(p8c))/p9b1;

p12b1:'log10(lhs(p12b)) = log10(rhs(p12b));

p12b2:'log10(lhs(p12b)) =rhs(p12b1),logexpand=super;

p12b3:lhs(p12b2) = expand(rhs(p12b2));

p12b4:lhs(p12b3) = lratsubst([p8f],rhs(p12b3));

p12b5:lhs(p12b4) = logcontract(rhs(p12b4));

p12b6:lhs(p12b5) = expand(rhs(p12b5));

3e) region
---------------

p12c:lhs(p8c) =num(rhs(p8c))/p9c1;;

p12c1:'log10(lhs(p12c)) = log10(rhs(p12c));

p12c2:'log10(lhs(p12c)) =rhs(p12c1),logexpand=super;

p12c3:lhs(p12c2) = expand(rhs(p12c2));

p12c4:lhs(p12c3) = lratsubst([p8f],rhs(p12c3));

p12c5:lhs(p12c4) = logcontract(rhs(p12c4));

How draw  log10(ph)  ,polyprotic acid: H2A
=====================================

given : C_H2A , 'oxalic acid' = 0.01M solution
              pK_a1 = 1.25
              pK_a2 = 4.28

solution :
========

start:
====

K_ai i=1,2,..n
------------------

p13a:pKa1 =1.25;
p13b:pKa2 =4.28;

p14a:K_a1 = subst([p13a],10^(-pKa1));
p14b:K_a2 = subst([p13b],10^(-pKa2));

concentration acid :
------------------------------

p15:C_H2A = 0.01;

how to construct piecewise function.
===============================

f(x):= block(
if x < 0
then 0
else if x≥ 0 and  x < 3
then x^2
else if x≥3
then   
x^3);
       

wxplot2d([f(x)], [x,-5,5],
 [plot_format, xmaxima])$

general : piecewise function with three regions(r1,r2,r3)

g(x,r1,r2,r3,conc,pKa1,pKa2,K_a1,K_a2) := 
block (
if x< pKa1
then r1
else if  x ≥ pKa1 and x < pKa2
then r2
else if  x ≥ pKa2 and x ≤ 14
then r3);

part 1:
---------

x -axes : ph
y-axes  : log10 (concentration)

xmin : 0;
xmax :14;
ymin : -14;
ymax:0;

part 2: log10(H),log10(OH)
---------------------------------------

definition : line trough two points
============================
-------------------------------------------------------
p1(x1,y1) , p2(x2,y2)
-------------------------------------------------------
m= (y2-y1)/(x2-x1)
y - y1 = m*(x-x1)  => y = y1+m*(x-x1)
-------------------------------------------------------
rem : y = m *x +n,first order equation
-------------------------------------------------------

twoptline(x1,y1,x2,y2):= y1+((y2-y1)/(x2-x1))*(x-x1);

1e)  log10(H) ,blue,points
-------------------

line : point (14,0) to (0,-14)

s1:twoptline(xmax,0,0,ymin);

2e)  log10(OH),red,points
---------------------

line : point (0,0) to (14,-14)

s2:twoptline(0,0,xmax,ymin);

wxplot2d( [s1,s2],
[x,xmin,xmax],[y,ymin,ymax],[same_xy],[style,[points,0.01,red],[points,0.01,blue]]
    ,grid2d,[xlabel,"ph"],[ylabel,"log(conc)"]);

p16:[ph=x,C_H2A =rhs(p15),pKa1 = rhs(p13a),pKa2=rhs(p13b),K_a1=rhs(p14a),K_a2=rhs(p14b)];

function : [H2A]   (3 regions,see above)  , 'p10a3,p10b6,p10c6'
==================================

q1:subst(p16,rhs(p10a4));
q2:subst(p16,rhs(p10b7));
q3:subst(p16,rhs(p10c7));

wxplot2d([s1,s2,g(x,q1,q2,q3,rhs(p15),rhs(p13a),rhs(p13b),rhs(p14a),rhs(p14b))],
    [x,xmin,xmax],[y,ymin,ymax],grid2d,[title,"log10([H2A]),green line"],
[legend,"H(+)","OH(-)","q1,q2,q3"])$

function : [HA(-)]   (3 regions,see above) , 'p11a6,p11b6,p11c6'
==================================

q1a:subst(p16,rhs(p11a6));
q2a:subst(p16,rhs(p11b6));
q3a:subst(p16,rhs(p11c6));

wxplot2d([s1,s2,g(x,q1a,q2a,q3a,rhs(p15),rhs(p13a),rhs(p13b),rhs(p14a),rhs(p14b))],
    [x,xmin,xmax],[y,ymin,ymax],grid2d,[title,"log10([HA(-)]),green line"],
[legend,"H(+)","OH(-)","q1a,q2a,q3a"])$

function : [A(2-)]   (3 regions,see above)  , 'p12a6,p12b6,p12c5'
==================================

q1b:subst(p16,rhs(p12a6));
q2b:subst(p16,rhs(p12b6));
q3b:subst(p16,rhs(p12c5));

wxplot2d([s1,s2,g(x,q1b,q2b,q3b,rhs(p15),rhs(p13a),rhs(p13b),rhs(p14a),rhs(p14b))],
    [x,xmin,xmax],[y,ymin,ymax],grid2d,[title,"log10([A2-]),green line"],
[legend,"H(+)","OH(-)","q1b,q2b,q3b"])$

combination : [H2A] and [HA(-)] and [A(2-)]
===================================

z1:[style,[points,0.01,red],[points,0.01,blue],[lines,0.5,red],[lines,0.5,green],[lines,0.5,blue]];
z2:[same_xy];
z3:[x,xmin,xmax];
z4:[y,ymin,ymax];
z5:[xlabel,"ph"];
z6:[ylabel,"log10(conc)"];
z7:[title,"log10 (conc ),[H2A],[HA(-)],[A(2-)]"] ;
z8:[legend,"H(+),red points","OH(-),blue points","H2A","HA(-)","A(2-)"];

wxplot2d([s1,s2,g(x,q1,q2,q3,rhs(p15),rhs(p13a),rhs(p13b),rhs(p14a),rhs(p14b))
                ,g(x,q1a,q2a,q3a,rhs(p15),rhs(p13a),rhs(p13b),rhs(p14a),rhs(p14b)),
                g(x,q1b,q2b,q3b,rhs(p15),rhs(p13a),rhs(p13b),rhs(p14a),rhs(p14b)) ]
                 ,z3,z4,z2,z1,z5,z6,z7,z8)$
               

rem : v1:line @ -2 = log10(C_H2A). ' [C_H2A] = 0.01 M'
============================================
----------------------------------------------------------------------------------------
level (v1) and intersection(red line and green line)    = pKa1
level (v1) and intersection(blue line and green line)  = pKa2
----------------------------------------------------------------------------------------
regio 1 : 'below' ,pKa1
regio 2 : 'between' , pKa1 < ph < pKa2
regio 3 : 'above',pKa2
---------------------------------------------------------------------------------------