option nolet
! serie : analytical chemistry !***************************************
!***************************************
! graphics (split the screen in two )
! open #n :screen xl,xh,yl,yh
!***************************************
open #1:screen 0,0.5,0,1
open #2:screen .5,1,.5,1
!***************************************
!***************************************
!*******************************************
window #1 ! ' open left part of the screen '
!*******************************************
declare function CA,CB
print "*************************************************"
print "*calculation : general equation H(m)A titration *"
print "*of polyprotic acid 'weak acid' with strong base*"
print "* version 1.0a ( example : diprotic acid,m=2) *"
print "*************************************************"
print "-------------------------------------------------"
print " [email protected],14/09/2016 "
print "-------------------------------------------------"
print "************************************************"
print " polyprotic acid: H(2)A = H(+) + HA(-) => ka(1) "
print " HA(-) = H(+) + A(2-) => ka(2) "
print "************************************************"
print " general : polyprotic acid ,m=n-1 => H(m)A ,m=2 "
print "************************************************"
print " use : alph(1) ,alpha(2),........,alpha(m) "
print "************************************************"
!file :TITRATIONGENERALpolyPROTICACIDWITHSTRONGBASEversion1a.tru
print
!********************************************************
!********************************************************
! tot_vol = v + vx
!********************************************************
! v = initial volume H(m)A ! vx = volume of base added from buret
!--------------------------------------------------------
! CA = CA0(v/tot_vol)
!--------------------------------------------------------
!CA0 = molarity initial concentration of H(m)A
!--------------------------------------------------------
! CB = CB0(vx/tot_vol)
!--------------------------------------------------------
!CB0 = molarity of strong base in the buret
!********************************************************
!********************************************************
!********************************************************
! charge balance
! [H(+)] = [HA(-)]+2*[A(2-)]+[OH(-)] (1)
!********************************************************
! [H(+)] = [HA(-)]+2*[A(2-)]+[OH(-)]-CB
! = (H(+))! CB negativ sign use up [H(+)]
!********************************************************
! [H(+)] = alpha(1)*CA+2*alpha(1)*CA+kw/[H(+)]-CB (2)
!********************************************************
! use above : definition CA,CB , solve (2) for !********************************************************
!********************************************************
! vx = ratio = num/denum
!********************************************************
!num = (alpha(1)*CA0+2*alpha(1)*CA0+kw/[H(+)]-[H(+)])*v
!denum = (CB0-kw/[H(+)]+[H(+)])
!********************************************************
!********************************************************
n=3 ! diprotic acid ( n = 4, triprotic acid H(3)PO4)
dim ka(1),alpha(0)
dim p(1)
m=n-1
mat redim ka(m) ! acid equilibruim constant (see above)
! p(i) numerator
mat redim p(n),alpha( 0 to m)
print " given = "
print
ka(1) = 6.46E-5
ka(2) = 3.31E-6
v = 50.2
CA0 = 0.0617
CB0 = 0.0954
!*******************************************
! ka(3) = ? value
!*******************************************
! change this n=4,m=4-1 =3,extra ka(3) .....
!*******************************************
for i=1 to m
print " ka(";i;") = ";ka(i)
next i
print " CA0 initial conc in vessel H(";m;")A = ";CA0;" M"
print " v initial volume in vessel = ";v;" ml"
print " CB0 initial conc in buret 'strongbase' = ";CB0;" M"
print
print " result = "
print
print " vx = ' added strong base ( ml) ' = num/denum "
print
print "num = (alpha(1)*CA0+2*alpha(1)*CA0+kw/[H(+)]-[H(+)])*v"
print "denum = (CB0-kw/[H(+)]+[H(+)]) "
!**********************************************************
!**********************************************************
! definition pH = - log10[H(+)],[H(+)] = 10^(-pH)
!**********************************************************
! x = 10^(-pH)
!**********************************************************
!**********************************************************
print
print "nu","pH","alpha(1)","2*alpha(2)","vx ml NaOH","vx_approx ml"
print
dim xp(2),yp(2)
nn =0
tb = 2.5
te = 12.5
st = .45
nt = int((te-tb)/st)
! calculation limit values (=> draw graphics )
! use : set window
for pH = tb to te step st
nn = nn+1
x = 10^(-pH)
call frac_denominator(x,ka(),m,som,p())
call distr_alpha(m,p(),som,alpha())
call num(x,m,CA0,alpha(),v,num_ratio)
call denum(CB0,x,denum_ratio)
vx =(num_ratio)/denum_ratio
call approx_vx(alpha(),CA0,CB0,ratio_approx,v)
call limit_draw(nn,nt,pH,ratio_approx,xp(),yp())
print nn,pH,alpha(1),2*alpha(2),vx,ratio_approx
next pH
!**********************************************
!**********************************************
! graphics part =>
!**********************************************
window #2 ! open right part to draw graphics
!**********************************************
!**********************************************
! plot the graph ::x_value = vx ,y_value = pH
!**********************************************
!**********************************************
! problem :aspect_ratio
! reason:number pixels => px > py => (px/py>1)
!**********************************************
ask pixels px,py
asp_r = px/py
!**********************************************
! set window limits
!**********************************************
set window xp(1)*asp_r,xp(2)*asp_r,yp(1),yp(2)
fact_x = 1/10
fact_y = 5
ii = 0
for pH = tb to te step 1
ii = ii+1
x = 10^(-pH)
call frac_denominator(x,ka(),m,som,p())
call distr_alpha(m,p(),som,alpha())
call num(x,m,CA0,alpha(),v,num_ratio)
call denum(CB0,x,denum_ratio)
vx =(num_ratio)/denum_ratio
call approx_vx(alpha(),CA0,CB0,ratio_approx,v)
plot points:ratio_approx*fact_x,ph*fact_y
plot text,at ratio_approx*fact_x,ph*fact_y :str$(ii)
next pH
plot text,at 10,10*fact_y: " vx:'ml added NaOH' versus pH "
end
sub limit_draw(nn,nt,pH,ratio_approx,xp(),yp())
if nn = 1 then
xp(1) = pH
yp(1) = ratio_approx
elseif ( nt = nn) then
xp(2) = pH
yp(2) = ratio_approx
end if
end sub
FUNCTION CA(v,vx,CA0)
tot_vol = v+vx
CA = CA0*(v/tot_vol)
END FUNCTION
FUNCTION CB(v,vx,CB0)
tot_vol = v+vx
CB = CB0*(vx/tot_vol)
END FUNCTION
sub approx_vx(alpha(),CA0,CB0,ratio_approx,v)
num=(alpha(1)*CA0+2*alpha(2)*CA0)*v
denum = CB0
! vx = ratio_approx
ratio_approx = num/denum
end sub
sub num(x,m,CZ,alpha(),v,num_ratio)
kw=1E-14
!------------------------------------------
! [OH(-)] = kw/[H(+)],v
!------------------------------------------
som2 = kw/x- x
for i=1 to m
som2 = som2 + i*alpha(i)*CZ
next i
!------------------------------------------------
!num=(alpha(1)*CZ+2*alpha(1)*CZ+[OH(-)]-[H(+)])*v
!------------------------------------------------
! num_ratio = som2*v
!-----------------------------------------------
num_ratio = som2*v
end sub
sub denum(C_b,x,denum_ratio)
! C_b = CB0
kw = 1E-14
denum_ratio = (C_b - kw/x+x)
end sub
sub distr_alpha(m,p(),som,alpha())
!---------------------------------
! alpha(0), alpha(1)....,alpha(m)
!---------------------------------
for i=0 to m
alpha(i) = p(i+1)/som
!print " alpha (";i;") = "; alpha(i)
next i
end sub
sub check_alpha(m,alpha())
!------------------------
! always one for som1
!------------------------
som1 = 0
for i=0 to m
som1 = som1 + alpha(i)
next i
print " sum ( alpha(i),i=0..";m;" ) = ";som1
end sub
sub frac_denominator(x,k(),m,som,p())
!********************************************************
!--------------------------
! see above : diprotic acid
!--------------------------
!print " (di)protic acid : n=3,m=2 (see code) "
!print " sum = (H[+])^2+ka(1)*H[+] + Ka(1)*ka(2) "
!********************************************************
som = x^m
p(1) = som
for i=1 to m
prod = 1
for j=1 to i
prod = prod * x^(m-i)*k(j)
next j
p(i+1) = prod
!*** denominator ***
som = som + p(i+1)
!*******************
next i
end sub