%-------------------------------------------------------------------
% Example for the PH-LIBRARY to test both solvers: advancment and 
%   alkalinity conservation 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 2 in Luff et al. using the 
% alkalinity conservation approach with results from the advancement
% approach
%-------------------------------------------------------------------
close;
close;
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
%-------------------------------------------------------------------
anz=5000;
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:anz
  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);
  pht1(i)=-log10(hp/1.0E+03);
  alkt1(i)= 2.0*co3+hco3+hs+boh4;
  tco2t1(i)=co2+co3+hco3;
  th2s1(i)=hs+h2s;
end;
plot (alkt1./tco2t1,pht1,'p');

%-------------------------------------------------------------------
% 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:anz
  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);
  pht2(i)=-log10(hp/1.0E+03);
  alkt2(i)= 2.0*co3+hco3+hs+boh4;
  tco2t2(i)=co2+co3+hco3;
  th2st2(i)=hs+h2s;
end;
plot (alkt2./tco2t2,pht2,'p');

%-------------------------------------------------------------------
% Finish plots
%-------------------------------------------------------------------
for i=1:2
  subplot(1,2,i)
  ylabel ('pH')
  xlabel ('TA/TCO_2')
  ph_font(12,12,12,'times');
  xlim([0 3]);
  ylim([2 12]);
end;
orient portrait
print -deps ph_eqadvan.eps

figure;
subplot(1,2,1)
hold on
th2s=0.0;
disp('Test again for carbonate species using alkalinity conservation')
for i=1:anz
  [hp,co2,co3,hco3,hs,h2s,boh3,boh4] = ph_eqalk(k,alkt1(i),tco2t1(i),th2s,tboh);
  ph=-log10(hp/1.0E+03);
  plot (pht1(i),ph,'p');
end;
subplot(1,2,2)
hold on
disp('Test again for carbonate and sulfate species using alkalinity conservation')
for i=1:anz
  [hp,co2,co3,hco3,hs,h2s,boh3,boh4] = ph_eqalk(k,alkt2(i),tco2t2(i),th2st2(i),tboh);
  ph=-log10(hp/1.0E+03);
  plot (pht2(i),ph,'p');
end;
x=2:0.01:12;
for i=1:2
  subplot(1,2,i)
  ylabel ('pH (Advancement)')
  xlabel ('pH (TA-conservation)')
  ph_font(12,12,12,'times');
  xlim([2 12]);
  ylim([2 12]);
  plot(x,x);
end;
orient portrait
print -deps ph_diff.eps
