function [hp,co2,co3,hco3,hs,h2s,boh3,boh4] = ph_eqalk(k,alk,tco2,th2s,tboh)
%-------------------------------------------------------------------
% PH_EQALK $Revision: 1.0 $ $Date: 1998/12/07 $
%          Copyright (C) GEOMAR, Roger Luff 1998.
%
% USAGE:[hp,co2,co3,hco3,hs,h2s,boh3,boh4]=ph_eqalk(k,alk,tco2,th2s,tboh);
%
% DESCRIPTION:
%   This function calculates the H+ concentration and the concentration
%   of the dissociated species co2,co3,hco3,hs,h2s,boh3,boh4 from
%   alkalinity, total carbonate, total sulfide and total borate
%   concentration and the equilibrium constants. To calculate the pH from H+ 
%   use: ph=-log10(hp/1.0E+03)
%
% INPUT:
%   k    = equilibrium constants from ph_thermo
%   alk  = Alkalinity [meq/l]
%   tco2 = Total carbonate concentration [mmol/l]
%   th2s = Total suffate concentration [mmol/l]
%   tboh = Total borate concentration [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.
%
% 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 ~=5
   error('ph_eqalk.m: Must pass 5 parameters')
end
%-------------------------------------------------------------------
% INITIALIZE THE CARBONAT SPECIES
%-------------------------------------------------------------------
co2=tco2*0.1;
co3=tco2*0.1;
hco3=tco2*0.8;
%-------------------------------------------------------------------
% SET UP LOCAL VARIABLES
%-------------------------------------------------------------------
maxiter=1500;
iter = 0;
atol = 100;
tol = 1E-12;
k3k2 = k(3)/k(2);
alktco2=alk/tco2;
th2stco2=th2s/tco2;
relaxfac=0.5;
%-------------------------------------------------------------------
% alktco2 > 1 : CALCULATE EQUILIBRIUM USING CO3 AND HCO3 AS DEP. VAR.
% alktco2 < 1 : CALCULATE EQUILIBRIUM USING CO2 AND HCO3 AS DEP. VAR.
%-------------------------------------------------------------------
if alktco2 > 1
%-------------------------------------------------------------------
% CHECK IF FUNCTION CAN CALCULATE A SULUTION, SEE REFERENCE
%-------------------------------------------------------------------
  if alktco2 < 2+(th2stco2)
%-------------------------------------------------------------------
% ITERATION LOOP STARTS HERE FOR CO3 AND HCO3 AS DEPENDENT VARIABLES
%-------------------------------------------------------------------
    while atol>tol & iter<maxiter
      iter = iter + 1;
      co2 = k3k2*hco3^2/co3;
      helptb = k(3)/k(4)*hco3/co3;
      helpts = k(3)/k(10)*hco3/co3;
%-------------------------------------------------------------------
% SET UP THE JACOBI
%-------------------------------------------------------------------
% line 1
      MATRIX(1,1) = 1.0+2.0*k3k2*hco3/co3;
      MATRIX(1,2) = 1.0-tboh/(1.0+helptb)^2*helptb/hco3-th2s/(1.0+helpts)^2*helpts/hco3;
% line 2
      MATRIX(2,1) = 1.0-k3k2*hco3^2/co3^2;
      MATRIX(2,2) = 2.0+tboh/(1.0+helptb)^2*helptb/co3+th2s/(1.0+helpts)^2*helpts/co3;
%-------------------------------------------------------------------
% INVERT THE JACOBI
%-------------------------------------------------------------------
      MATRIX=inv(MATRIX);
%-------------------------------------------------------------------
% SET UP THE FUNCTIONS WHICH HAVE TO BE MINIMIZED
%-------------------------------------------------------------------
      f(1)=co2+hco3+co3-tco2;
      f(2)=hco3+2.0*co3+tboh/(1.0+helptb)+th2s/(1.0+helpts)-alk;
      dhco3 = 0;
      dco3 = 0;
      for j=1:2
        dhco3 = dhco3 + MATRIX(j,1)*f(j);
        dco3 = dco3 + MATRIX(j,2)*f(j);
      end
      hco3=hco3-dhco3*relaxfac;
      co3=co3-dco3*relaxfac;
%-------------------------------------------------------------------
%     if Newton-Raphson is going in the wrong direction try it
%     again with a new choice of first guess concentration.
%-------------------------------------------------------------------
      if (hco3 <= 0) | (co3 <= 0)
        co3=rand(1)*10;
        hco3=rand(1)*10;
        dhco3=2.0;
      end;
      atol = abs(dhco3)+abs(dco3);
    end;
%-------------------------------------------------------------------
% SET RESULT NOT NaN, WHEN NO SOLUTION CAN BE CALCULATED
%-------------------------------------------------------------------
  else
    hco3 = NaN;
    co3  = NaN; 
  end;
%-------------------------------------------------------------------
% CALCULATE THE DIRECT DEPENDING VARIABLES
%-------------------------------------------------------------------
  co2 = k3k2*hco3^2/co3;
  hp   = hco3*k(3)/co3;
else;
%-------------------------------------------------------------------
% ITERATION LOOP STARTS HERE FOR CO2 AND HCO3 AS DEPENDENT VARIABLES
%-------------------------------------------------------------------
  while atol>tol & iter<maxiter
    iter = iter + 1;
    co3 = k3k2*hco3^2/co2;
    helptb = k(2)/k(4)*co2/hco3;
    helpts = k(2)/k(10)*co2/hco3;
%-------------------------------------------------------------------
% SET UP THE JACOBI
%-------------------------------------------------------------------
% line 1
    MATRIX(1,1) = 1.0+2.0*k3k2*hco3/co2;
    MATRIX(1,2) = 1.0+4.0*k3k2*hco3/co2+tboh/(1.0+helptb)^2*helptb/hco3+th2s/(1.0+helpts)^2*helpts/hco3;
% line 2
    MATRIX(2,1) = 1.0-k3k2*hco3^2/co2^2;
    MATRIX(2,2) = -2.0*k3k2*hco3^2/co2^2-tboh/(1.0*helptb)^2*helptb/co2-th2s/(1.0*helpts)^2*helpts/co2;
%-------------------------------------------------------------------
% INVERT THE JACOBI
%-------------------------------------------------------------------
    MATRIX=inv(MATRIX);
%-------------------------------------------------------------------
% SET UP THE FUNCTIONS WHICH HAVE TO BE MINIMIZED
%-------------------------------------------------------------------
    f(1)=co2+hco3+co3-tco2;
    f(2)=hco3+2.0*co3+tboh/(1.0+helptb)+th2s/(1.0+helpts)-alk;
    dhco3 = 0;
    dco2 = 0;
    for j=1:2
      dhco3 = dhco3 + MATRIX(j,1)*f(j);
      dco2 = dco2 + MATRIX(j,2)*f(j);
    end
    hco3=hco3-dhco3*relaxfac;
    co2=co2-dco2*relaxfac;
%-------------------------------------------------------------------
% CHECK THE CONVERGENCE
%-------------------------------------------------------------------
    atol = abs(dhco3)+abs(dco2);
  end;
%-------------------------------------------------------------------
% CALCULATE THE DIRECT DEPENDING VARIABLES
%-------------------------------------------------------------------
  co3 = k3k2*hco3^2/co2;
  hp   = co2*k(2)/hco3;
end;
%-------------------------------------------------------------------
% IN THE CASE THAT THE FUNCTION DOES NOT CONVERGENTE, END WITH AND
% ERROR MESSAGE
%-------------------------------------------------------------------
if (iter==maxiter)
  error('pH_eqalk: too much iterations');
end
%-------------------------------------------------------------------
% IF THE FUNCTION CONVERGENTES, CALCULATE THE OTHER DISSOSIATED SPECIES
%-------------------------------------------------------------------
helptb = k(2)/k(4)*co2/hco3;
helpts = k(2)/k(10)*co2/hco3;
boh4 = tboh/(1.0+helptb);
hs   = th2s/(1.0+helpts);
boh3 = tboh-boh4;
h2s  = th2s-hs;
%-------------------------------------------------------------------
% END OF PH_EQALK
%-------------------------------------------------------------------
