option nolet
! serie : analytical chemistry
print "************************************************"
print "* calculation : species present in solutions   *"
print "* of polyprotic acid of the conc CZ ? pH       *"
print "* version 1.0                                  *"
print "************************************************"
print "------------------------------------------------"
print " [email protected],09/09/2016          "
print "------------------------------------------------"
print "************************************************"
print " polyprotic acid: H(2)S = H(+) + HS(-) => ka(1) "
print "                  HS(-) = H(+) + S(2-) => ka(2) "
print "************************************************"
print " general : polyprotic acid ,m=n-1 => H(m)A ,m=2 "
print "************************************************"
print " (di)protic acid : n=3,m=2 (see code)           "
print " sum = (H[+])^2+ka(1)*H[+] + Ka(1)*ka(2)        "
print "  CZ= [H(2)S]+[HS(-)]+[S(2-)]                   "
print " alpha(0) = (H[+])^2/ sum   = [H(2)S]/CZ        "
print " alpha(1) = ka(1)*H[+]/sum  = [HS(-)]/CZ        "
print " alpha(2) = Ka(1)*ka(2)/sum = [S(2-)]/CZ        "
print " check = alpha(0)+alpha(1)+alpha(2)=1           "
print "************************************************"
!file :FRACTIONCALANALYTICALSPECIESusePHweakacidfinal.tru
print
!*******************************************************
! use : charge balance ( general use S=A)
! Positive ion : [H(+)]
! Negative ion : [HS(-)]+[S(2-)]+[OH(-)] 
! [H(+)] = [HS(-)]+[S(2-)]+[OH(-)] 
! build : second_approx_h (see below sub)
! [H(+)] = alpha(1)*CZ+2*alpha(1)*CZ+[OH(-)]
!*******************************************************
! to use above : first approximation [H+]
!   [H(+)] = sqr(ka(1)*CZ) to CVG
!*******************************************************
n=3     ! diprotic acid ( n = 4, triprotic acid H(3)PO4)
dim ka(1)
dim p(1)
m=n-1
mat redim ka(m) ! acid equilibruim constant (see above)
! p(i) numerator
mat redim p(n) 
print "*********"
print " given = "
print "*********"
ka(1) = 1E-7
ka(2) = 1E-14
!*******************************************
! 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
CZ=0.01     ! Molair concentration diprotic acid
print "******************************"
print " first approximation [H(+)] = "
print "******************************"
print
print " concentration weak acid(= CZ) = ";CZ;" M"
x=sqr(ka(1)*CZ)
print " [H(+)] = ";x
print " pH     = ";-log(x)/log(10)
call frac_denominator(x,ka(),m,som,p())
!**********************
! build the fractions 
!**********************
print
print "**********"
print " result = "
print "**********"
print
dim alpha(0 to 1)
mat redim alpha(0 to m)
call distr_alpha(m,p(),som,alpha())
call sec_approx_h(x,m,CZ,alpha(),x_new)
call check_alpha(m,alpha())
print "******************************"
print "* Second approximation        "
print "******************************"
print " [H(+)]_new = ";x_new
print " pH([H(+)]_new) = ";-log(x_new)/log(10)
print
call frac_denominator(x_new,ka(),m,som,p())
call distr_alpha(m,p(),som,alpha())
call sec_approx_h(x_new,m,CZ,alpha(),x_new)
call check_alpha(m,alpha())
print "******************************"
print "* Third approximation         "
print "******************************"
print " [H(+)]_new = ";x_new
print " pH([H(+)]_new) = ";-log(x_new)/log(10)
end

sub sec_approx_h(x,m,CZ,alpha(),x_new)
kw=1E-14
!------------------------------------------
! [OH(-)] = kw/[H(+)]
!------------------------------------------
som2 = kw/x
for i=1 to m
  som2 = som2 + i*alpha(i)*CZ
next i
!------------------------------------------
![H(+)] = alpha(1)*CZ+2*alpha(1)*CZ+[OH(-)]
! x_new = [H(+)]
!------------------------------------------
x_new = som2
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
!--------------------------
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