%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.
%   Example for the usage of the solver, based on the charge balance
%   equation
%
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;
tnh4=0.0;
tpo4=3.0000000000000D-03;
no3=3.7000000000000D-02;
mnii = 0;
feii = 1;
%-------------------------------------------------------------------
% Calculate the density of sewater
%-------------------------------------------------------------------
rsw=ph_density(t,s,p);
%-------------------------------------------------------------------
% Calculate the total borate and sulfate concentration and the 
% concentration of the conservative species (ssum) from seawater
%-------------------------------------------------------------------
[ssum,tboh,ca,so4] = ph_swconst(s,rsw);
%-------------------------------------------------------------------
% Calculate the thermodynamic constants
%-------------------------------------------------------------------
k=ph_tdconst(t,s,p);
[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);
plot (alk/tco2,ph,'p');
hold on;
for i=1:500
  tco2=rand(1)*10+0.5;
  [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);
  alk  = 2.0*co3+hco3+hs+boh4;
  plot (alk/tco2,ph,'p');
end;
ylabel ('pH')
xlabel ('Alkalinity/TCO_2')
xlim([0 3]);
ylim([2 12]);

