%-------------------------------------------------------------------
% Example for the PH-LIBRARY to test the solver: advancment from:
%   Roger Luff and Matthias Haeckel 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 aqueous systems", Computers & Geosciences.
%
% This program will create Figure 1 in Luff et al. using the 
% advancement approach
%-------------------------------------------------------------------
clear;
clf;
%-------------------------------------------------------------------
% Parameters of the environment
%-------------------------------------------------------------------
t=3.0;
s=34.84;
p=190;
%-------------------------------------------------------------------
% Calculate the density of sewater
%-------------------------------------------------------------------
rsw= ph_density(t,s,p);
%-------------------------------------------------------------------
% Calculate the total borate concentration from seawater
%-------------------------------------------------------------------
[ssum,tboh,caw,so4w] = ph_swconst(s,rsw);
%-------------------------------------------------------------------
% Calculate the thermodynamic constants
%-------------------------------------------------------------------
k=ph_tdconst(t,s,p);
%-------------------------------------------------------------------
% Initialize parameters
%-------------------------------------------------------------------
alk=2.42;
tco2=2.8;
th2s=0.0;
co2=tco2*0.1;
co3=tco2*0.1;
hco3=tco2*0.8;
%-------------------------------------------------------------------
% Calculate the borate species distribution with the alkalinity approach
%-------------------------------------------------------------------
[hp,co2i,co3i,hco3i,hsi,h2si,boh3i,boh4i]=ph_eqalk(k,alk,tco2,th2s,tboh);
subplot(1,2,1)
hold on
disp('Start random test for carbonate species using advancement approach')
%-------------------------------------------------------------------
% Calculate 5000 pH-values with randomly generated carbonate species
% distributions
%-------------------------------------------------------------------
for i=1:5000
  co2i=rand(1)*100+0.0001;
  co3i=rand(1)*100+0.0001;
  hco3i=rand(1)*100+0.0001;
[hp,co2,co3,hco3,hs,h2s,boh3,boh4] = ph_eqadvan(k,co2i,co3i,hco3i,hsi,h2si,boh3i,boh4i);
  ph=-log10(hp/1.0E+03);
  tco2 = co2+co3+hco3;
  alk  = 2.0*co3+hco3+hs+boh4;
  plot (alk/tco2,ph,'p');
end;

%-------------------------------------------------------------------
% Calculate 5000 pH-values with randomly generated carbonate and
% sulfide species distributions
%-------------------------------------------------------------------
disp('Start random test for carbonate and sulfate species using advancement approach')
subplot(1,2,2)
hold on
for i=1:5000
  co2i=rand(1)*100+0.0001;
  co3i=rand(1)*100+0.0001;
  hco3i=rand(1)*100+0.0001;
  hsi=rand(1)*10;
  h2si=rand(1)*10;
[hp,co2,co3,hco3,hs,h2s,boh3,boh4] = ph_eqadvan(k,co2i,co3i,hco3i,hsi,h2si,boh3i,boh4i);
  ph=-log10(hp/1.0E+03);
  tco2 = co2+co3+hco3;
  alk  = 2.0*co3+hco3+hs+boh4;
  plot (alk/tco2,ph,'p');
end;
%-------------------------------------------------------------------
% Finish plots
%-------------------------------------------------------------------
for i=1:2
  subplot(1,2,i)
  ylabel ('pH')
  xlabel ('Alkalinity/TCO_2')
  ph_font(12,12,12,'times');
  xlim([0 3]);
  ylim([2 12]);
end;
orient portrait
print -deps ph_eqadvan.eps
