option nolet
! split screen two parts
!***************************************
!***************************************
! 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
!***************************************
! serie : analytical chemistry 'graphics'
!***************************************
window #1 ! 'open left part of the screen
print "**********************************************"
print "* draw distibution curves,pH versus fraction *"
print "* in the solution for H(m)A. 'version 1.0' *"
print "**********************************************"
print "* [email protected],02/10/2016 *"
print "**********************************************"
print "*species solution H3PO4 => species *"
print "* H(3)PO4,H(2)PO4(2-),HPO4(-),PO4(3-) *"
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(0) ,alpha(2),........,alpha(m) "
print "************************************************"
n=3 ! diprotic acid ( n = 4, triprotic acid H(3)PO4)
dim ka(1),alpha(0)
dim p(1),pp(0 to 1)
m=n-1
mat redim ka(m) ! acid equilibruim constant (see above)
! p(i) numerator
mat redim p(n),alpha( 0 to m),pp(0 to m)
print "*********"
print " given = "
print "*********"
print
ka(1) = 1.51E-3
ka(2) = 2.2E-6
print " polyprotic acid = H(";m;")A"
print
for i=1 to m
print " ka(";i;") = "; ka(i)
print " pka(";i;") = ";-log(ka(i))/log(10)
next i
print
print "**********"
print " result = "
print "**********"
dim tp(1)
mat redim tp(3) ! vertical lines,'distribution diagram'
print
for i= 1 to 2
print " intersection = alpha(";i-1;") and alpha(";i;") = "
print
tp(i)= -log(ka(i))/log(10)
print " pH = ";tp(i)
print
next i
print " intersection = alpha(";0;") and alpha(";m;") = "
print
!**************************************************
prodx=1
for i=1 to m
prodx=prodx*ka(i)
next i
hh =(prodx)^(1/m)
print " [H(+)] =(prod(ka(i),i=1..";m;")^(1/";m;")=";hh
tp(n)= -log(hh)/log(10)
print
print " pH = ";-log(hh)/log(10)
print
! sum of pka's divide by m
call som_mean(n,tp(),som)
print
print " sum(pka(i),i=1..";m;")/";m;" = ";som
!**************************************************
!**********************************************
window #2 ! 'open right part to draw graphics
!**********************************************
dim xp(2),yp(2)
ntx = 14 ! tick marks x-axis,ph values
nty = 10 ! tick marks y-axis
xp(1)=0 ! xmin
xp(2)=14 ! xmax
yp(1)=0 ! ymin
yp(2)=1 ! ymax
ask pixels px,py
asp_r = 1! px/py !aspect ratio correction(rectangular window)
set window xp(1)*asp_r,xp(2)*asp_r,yp(1),yp(2)
call draw_axes(xp(1)*asp_r,xp(2)*asp_r,yp(1),yp(2),ntx,nty)
! output distribution curve
dim ll$(0 to 1)
mat redim ll$(0 to m)
for i1=0 to 14 step 0.001
x=10^(-i1)
call frac_denominator(x,ka(),m,som1,p())
for i=0 to m
! color selector
if i=0 then
set color "red"
lt$ = "alpha("& str$(i)&")"
plot text,at 1,(i+1)*0.1:lt$
end if
if i=1 then
set color "green"
lt$ = "alpha("& str$(i)&")"
plot text,at 1,(i+1)*0.1:lt$
end if
if i=2 then
set color "blue"
lt$ = "alpha("& str$(i)&")"
plot text,at 1,(i+1)*0.1:lt$
end if
if i=3 then
set color "black"
lt$ = "alpha("& str$(i)&")"
plot text,at 1,(i+1)*0.1:lt$
end if
if i=4 then
set color "Magenta"
lt$ = "alpha("& str$(i)&")"
plot text,at 1,(i+1)*0.1:lt$
end if
if i=5 then
set color "yellow"
lt$ = "alpha("& str$(i)&")"
plot text,at 1,(i+1)*0.1:lt$
end if
if i=6 then
set color "cyan"
lt$ = "alpha("& str$(i)&")"
plot text,at 1,(i+1)*0.1:lt$
end if
! end selector
if i=0 then
pp(i) = x^m
plot points: i1,pp(i)/som1
else
prod1 = 1
for j=1 to i
prod1=prod1*ka(j)*x^(m-i)
next j
pp(i) = prod1
plot points: i1,pp(i)/som1
end if
next i
!**** with out color **************************
! plot points: i1,x^m/som1
! plot points: i1,(ka(1)*(x)^(m-1))/som1
! plot points: i1,(ka(1)*ka(2)*(x)^(m-2))/som1
!**********************************************
next i1
set color "black"
title$ = " distribution curve H("& str$(m) &")A"
plot text,at 1,.9: title$
!*****************************************
!**** range x axis 0 to 14 ****
!*****************************************
for i=0 to 14
plot text,at i,0.2/10:str$(i)
next i
plot text,at 13.5,0.5/10:"pH"
!*****************************************
!**** range y axis 0 to 1 (fraction)****
!*****************************************
for i=1 to 10
plot text,at 0.2,i/10:str$(i/10)
next i
plot text,at 0.3,9.5/10:"fraction"
!*************************************
! make vertical lines tp matrix values
!*************************************
for i=1 to n
plot lines :tp(i),0;tp(i),1
next i
end
sub som_mean(m,tp(),som)
som2=0
for i=1 to m
som2 = som2 + tp(i)
next i
! arithmic mean
som = som2/m
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
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 draw_axes(xmin,xmax,ymin,ymax,ntx,nty)
! axis :for distribution function
! see extra,website
! distance between tick marks on x-axis
LET dx = (xmax - xmin)/ntx
! distance between tick marks on y-axis
LET dy = (ymax - ymin)/nty
SET WINDOW xmin ,xmax ,ymin ,ymax
LET x0 = max(0,xmin)
LET y0 = max(0,ymin)
IF ymin*ymax < 0 then
LET y0 = 0
ELSE
LET y0 = ymin
END IF
PLOT LINES: xmin,y0;xmax,y0 ! horizontal axis
PLOT LINES: x0,ymin;x0,ymax ! vertical axis
LET tx = 0.1*dy ! size of tick mark
LET ty = 0.1*dx
FOR itick = 0 to ntx
LET x = xmin + itick*dx
PLOT LINES: x,y0 - tx; x,y0 + tx ! draw ticks on x axis
next itick
FOR itick = 0 to nty
LET y = ymin + itick*dy
PLOT LINES: x0 - ty,y; x0 + ty,y ! draw ticks on y axis
NEXT itick
END SUB