%Example for the PH-LIBRARY from:
%   Roger Luff 98-12-07
%   GEOMAR Research Centre for Marine Geosciences
%   Wischhofstr. 1-3, D-24148 Kiel, Germany
%   visumod@technologist.com
%
% REFERENCES:
%   Luff R., Haeckel M. and Wallmann K. (submitted) "Robust and fast 
%   FORTRAN and MATLAB libraries to calculate pH distributions in  
%   marine systems", Computers & Geosciences.
%
% This program will create a Table using the 
% different solvers to test the internal consistency
%
clear;
clf;
% Parameters of the environment
t=3.0;
s=34.84;
p=190;
% Initialization parameters
alk=2.42;
tco2=2.2;
th2s=0.0;
co3  = 0.11*tco2;
co2  = 0.09*tco2;
hco3 = 0.80*tco2;
% Calculate the density of sewater
rsw=ph_density(t,s,p);
% Calculate the total borate concentration from seawater
[ssum,tboh,ca,so4] = ph_swconst(s,rsw);
% Calculate the thermodynamic constants
k=ph_tdconst(t,s,p);
% Initialization parameters
hs = th2s/2;
h2s=th2s/2;
boh3=tboh/2;
boh4=tboh/2;
disp('Recalculate the results from Table 1 in Luff et al.:')

% First calculation, using the advancement approach
[hp,co2,co3,hco3,hs,h2s,boh3,boh4] = ph_eqadvan(k,co2,co3,hco3,hs,h2s,boh3,boh4);
ph=-log10(hp/1.0E+03);
[me,e]=sprintf('ph_eqadvan : %10.10f %10.10f %10.10f %10.10f %10.10f %10.10f\n',ph,co2,co3,hco3);
disp(me)
alk  = 2.0*co3+hco3+hs+boh4;

% Second calculation, using the alkalinity conservation approach
[hp,co2,co3,hco3,hs,h2s,boh3,boh4] = ph_eqalk(k,alk,tco2,th2s,tboh);
ph=-log10(hp/1.0E+03);
[me,e]=sprintf('ph_eqalk   : %10.10f %10.10f %10.10f %10.10f %10.10f %10.10f\n',ph,co2,co3,hco3);
disp(me)

% Thrid calculation, using the advancement approach again
[hp,co2,co3,hco3,hs,h2s,boh3,boh4] = ph_eqadvan(k,co2,co3,hco3,hs,h2s,boh3,boh4);
hp   = hco3*k(3)/co3;
[me,e]=sprintf('ph_eqadvan : %10.10f %10.10f %10.10f %10.10f %10.10f %10.10f\n',ph,co2,co3,hco3);
disp(me)
tnh4=0.0;
tpo4=0.0;
no3=0.0;
mnii = 0;
feii = 0;
alk  = 2.0*co3+hco3+hs+boh4;

% Forth calculation, using the charge conservation approach
[hp,co2,co3,hco3,hs,h2s,h3po4,h2po4,hpo4,po4,nh3,nh4,boh3,boh4]=ph_eqcharge(k,ssum,tco2,th2s,tboh,tnh4,tpo4,mnii,feii,ca,so4,no3);
ph=-log10(hp/1.0E+03);
[me,e]=sprintf('ph_eqcharge: %10.10f %10.10f %10.10f %10.10f %10.10f %10.10f\n',ph,co2,co3,hco3);
disp(me)

% Fifth calculation, using the advancement approach again
[hp,co2,co3,hco3,hs,h2s,boh3,boh4] = ph_eqadvan(k,co2,co3,hco3,hs,h2s,boh3,boh4);
ph=-log10(hp/1.0E+03);
[me,e]=sprintf('ph_eqadvan : %10.10f %10.10f %10.10f %10.10f %10.10f %10.10f\n',ph,co2,co3,hco3);
disp(me)

alk  = 2.0*co3+hco3+hs+boh4;
[hp,co2,co3,hco3,hs,h2s,boh3,boh4] = ph_eqalk(k,alk,tco2,th2s,tboh);
ph=-log10(hp/1.0E+03);
% Last  calculation, using the alkalinity conservation approach
[me,e]=sprintf('ph_eqalk   : %10.10f %10.10f %10.10f %10.10f %10.10f %10.10f\n',ph,co2,co3,hco3);
disp(me)
