%-------------------------------------------------------------------
% Example for the PH-LIBRARY to test the solver: alkalinity 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 Figure 1 in Luff et al. using the 
% alkalinity approach
%-------------------------------------------------------------------
clear;
clf;
%-------------------------------------------------------------------
% Parameters of the environment
%-------------------------------------------------------------------
t=3.0;
s=34.84;
p=190;
%-------------------------------------------------------------------
% Initialize parameters
%-------------------------------------------------------------------
alk=2.42;
tco2=2.2;
th2s=0.0;
%-------------------------------------------------------------------
% 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);
[hp,co2,co3,hco3,hs,h2s,boh3,boh4] = ph_eqalk(k,alk,tco2,th2s,tboh);
subplot(1,2,1)
hold on
disp('Start random test for carbonate species using alkalinity conservation')
%-------------------------------------------------------------------
% Calculate 5000 pH-values with randomly generated carbonate species
% distributions
%-------------------------------------------------------------------
for i=1:5000
  alk=rand(1)*100+0.00001;
  tco2=rand(1)*100+0.00001;
  [hp,co2,co3,hco3,hs,h2s,boh3,boh4] = ph_eqalk(k,alk,tco2,th2s,tboh);
  ph=-log10(hp/1.0E+03);
  plot (alk/tco2,ph,'p');
end;
ylabel ('pH')
xlabel ('Alkalinity/TCO_2')
xlim([0 3]);
ylim([2 12]);
ph_font(12,12,12,'times')
subplot(1,2,2)
th2s = 0.0;
hold on
disp('Start random test for carbonate and sulfide species using alkalinity conservation')
%-------------------------------------------------------------------
% Calculate 5000 pH-values with randomly generated carbonate and
% sulfide species distributions
%-------------------------------------------------------------------
for i=1:5000
  alk=rand(1)*100+0.00001;
  tco2=rand(1)*100+0.00001;
  th2s=rand(1)*15;
  [hp,co2,co3,hco3,hs,h2s,boh3,boh4] = ph_eqalk(k,alk,tco2,th2s,tboh);
  ph=-log10(hp/1.0E+03);
  plot (alk/tco2,ph,'p');
end;
ylabel ('pH')
xlabel ('Alkalinity/TCO_2')
xlim([0 3]);
ylim([2 12]);
ph_font(12,12,12,'times')
orient portrait
print -deps ph_eqalk.eps
