function [hp,co2,co3,hco3,hs,h2s,boh3,boh4] = ph_eqadvan(k,co2i,co3i,hco3i,hsi,h2si,boh3i,boh4i)
%-------------------------------------------------------------------
% PH_EQADVAN $Revision: 1.0 $ $Date: 1998/12/07 $
%            Copyright (C) GEOMAR, Roger Luff, Matthias Haeckel 1998.
%
% USAGE:[hp,co2,co3,hco3,hs,h2s,boh3,boh4]=ph_eqadvan(k,co2i,co3i,hco3i,hsi,h2si,boh3i,boh4i);
%
% DESCRIPTION:
%   This function calculates the H+ concentration and the concentration
%   of the dissociated species to the closest equilibrium of co2,co3,hco3,hs,
%   h2s,boh3,boh4 from the dissociated species not in equilibrium co2,co3,
%   hco3,hs,h2s,boh3,boh4, and the equilibrium constants. To calculate 
%   the pH from H+ use: ph=-log10(hp/1.0E+03)
%
% INPUT:
%   k    = equilibrium constants from ph_thermo
%   hpi  = Initial concentration H+ [mmol/l]
%   co2i = Initial concentration of CO2 [mmol/l]
%   co3i = Initial concentration of CO3[mmol/l]
%   hco3i= Initial concentration of HCO3 [mmol/l]
%   hsi  = Initial concentration of HS [mmol/l]
%   h2si = Initial concentration of H2S [mmol/l]
%   boh3i= Initial concentration of BOH3 [mmol/l]
%   boh4i= Initial concentration of BOH4 [mmol/l]
%
% OUTPUT:
%   hp   = Equilibrium concentration H+ [mmol/l]
%   co2  = Equilibrium concentration of CO2 [mmol/l]
%   co3  = Equilibrium concentration of CO3[mmol/l]
%   hco3 = Equilibrium concentration of HCO3 [mmol/l]
%   hs   = Equilibrium concentration of HS [mmol/l]
%   h2s  = Equilibrium concentration of H2S [mmol/l]
%   boh3 = Equilibrium concentration of BOH3 [mmol/l]
%   boh4 = Equilibrium concentration of BOH4 [mmol/l]
%
% AUTHOR: 
%   Roger Luff 98-12-07  (rluff@gmx.de)
%   GEOMAR Research Centre for Marine Geosciences
%   Wischhofstr. 1-3, D-24148 Kiel, Germany
%   visumod@technologist.com
%
% DISCLAIMER:
%   This software is provided "as it is" without warranty of any kind,
%   see ph_info for details.
%
% REFERENCES:
%   Luff, Haeckel and Wallmann (submitted) "Robust and fast FORTRAN
%   and MATLAB libraries to calculate pH distributions in
%   marine systems", Computers & Geosciences.
%-------------------------------------------------------------------

%-------------------------------------------------------------------
% CHECK INPUT ARGUMENTS
%-------------------------------------------------------------------
if nargin ~=8
   error('ph_eqtherm.m: Must pass 8 parameters')
end
%-------------------------------------------------------------------
% COPY THE INPUT VALUES AS A FIRST GUESS TO THE FINAL EQUILIBRIUM 
% VALUES
%-------------------------------------------------------------------
if co2i<0.1 
  hpi   = hco3i*k(3)/co3i;
else
  hpi   = co2i*k(2)/hco3i;
end
hp   = hpi;
co2=co2i;
co3=co3i;
hco3=hco3i;
hs=hsi;
h2s=h2si;
boh3=boh3i;
boh4=boh4i;
%-------------------------------------------------------------------
% SET UP LOCAL VARIABLES
%-------------------------------------------------------------------
iter = 0;
delta1 = 0;
delta2 = 0;
delta3 = 0;
delta4 = 0;
atol=100;
tol = 1E-11;
maxiter = 500;
%-------------------------------------------------------------------
% ITERATION LOOP STARTS HERE 
%-------------------------------------------------------------------
while atol>tol & iter<maxiter
  iter   = iter + 1;
%-------------------------------------------------------------------
% SET UP THE JACOBI
%-------------------------------------------------------------------
  MATRIX=zeros(12);
% line 1
  MATRIX(1,1) = -1.0;
  MATRIX(1,7) = -k(2) ;
% line 2
  MATRIX(2,2) = -1.0;
  MATRIX(2,7) = hp;
  MATRIX(2,8) = -k(3);
% line 3
  MATRIX(3,3) = -1.0;
  MATRIX(3,8) = hp;
% line 4
  MATRIX(4,4) = -1.0;
  MATRIX(4,9) = -k(10);
% line 5
  MATRIX(5,5) = -1.0 ;
  MATRIX(5,9) = hp;
% line 6
  MATRIX(6,6) = -1.0;
  MATRIX(6,7) = hco3;
  MATRIX(6,8) = co3;
  MATRIX(6,9) = hs;
  MATRIX(6,12) = boh4;
% line 7
  MATRIX(7,1) = -1.0;
  MATRIX(7,2) = 1.0;
  MATRIX(7,6) = 1.0;
% line 8
  MATRIX(8,2) = -1.0;
  MATRIX(8,3) = 1.0;
  MATRIX(8,6) = 1.0;
% line 9
  MATRIX(9,4) = -1.0;
  MATRIX(9,5) = 1.0;
  MATRIX(9,6) = 1.0;
% line 10
  MATRIX(10,6) = 1.0;
  MATRIX(10,10) = -1.0 ;
  MATRIX(10,11) = 1.0 ;
% line 11
  MATRIX(11,10) = -1.0 ;
  MATRIX(11,12) = -k(4);
% line 12
  MATRIX(12,11) = -1.0;
  MATRIX(12,12) = hp;
%-------------------------------------------------------------------
% INVERT THE JACOBIAN
%-------------------------------------------------------------------
  MATRIX=inv(MATRIX);
%-------------------------------------------------------------------
% SET UP THE FUNCTIONS WHICH HAVE TO BE MINIMIZED
%-------------------------------------------------------------------
  f(1) = co2i-delta1-co2;
  f(2) = hco3i+delta1-delta2-hco3;
  f(3) = co3i+delta2-co3;
  f(4) = h2si-delta3-h2s;
  f(5) = hsi+delta3-hs;
  f(6) = hpi+delta1+delta2+delta3+delta4-hp;
  f(7) = hco3*hp-k(2)*co2;
  f(8) = co3*hp-k(3)*hco3;
  f(9) = hs*hp-k(10)*h2s;
  f(10)= boh3i-delta4-boh3;
  f(11)= boh4i+delta4-boh4;
  f(12)= boh4*hp-k(4)*boh3;
  dd(1:12)=0;
%-------------------------------------------------------------------
% CALCULATE THE DISTANCE TO THE NEXT INTERATION STEP
%-------------------------------------------------------------------
  for kk=1:12
    for j=1:12
      dd(kk) =dd(kk) +MATRIX(j,kk)*f(j);
    end
  end;
%-------------------------------------------------------------------
% CALCULATE THE NEW DISTRIBUTION
%-------------------------------------------------------------------
  co2=co2-dd(1);
  hco3=hco3-dd(2);
  co3=co3-dd(3);
  h2s=h2s-dd(4);
  hs=hs-dd(5);
  hp =hp-dd(6);
  delta1=delta1-dd(7);
  delta2=delta2-dd(8);
  delta3=delta3-dd(9);
  delta4=delta4-dd(10);
  boh3=boh3-dd(11);
  boh4=boh4-dd(12);
  if (hp <= 0)
    hp =rand(1);
    dd(1)=2.0;
  end;
%-------------------------------------------------------------------
% CHECK THE CONVERGENCE
%-------------------------------------------------------------------
  atol = sum(abs(dd));
end;
%-------------------------------------------------------------------
% IN THE CASE THAT THE FUNCTION DOES NOT CONVERGE, END WITH AND
% ERROR MESSAGE
%-------------------------------------------------------------------
if (iter==maxiter)
  error('pH_eqadvan: too much iterations');
end
