## Copyright (C) 2023 Jens
##
## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <https://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {} {@var{retval} =} brewwater_optimizer (@var{input1}, @var{input2})
##
## @seealso{}
## @end deftypefn

## Author: Jens <jens@orion>
## Created: 2023-06-08

function brewwater_optimizer (input1, input2)

# Load constants
water_constants()

k_ve = 0.25;

m_Ca_0   = 94.0;	% mg/l
m_Mg_0   = 18.1;	% mg/l
m_Na_0   =  3.6;	% mg/l
m_S_0    = 44.1;	% mg/l
m_Cl_0   =  4.6;	% mg/l
k_HCO3_0 =  5.5;	% mmol/l

M_CaSO4 = create_mol([i_Ca i_S i_O], [1 1 4], mol_masses)
M_CaCl2 = create_mol([i_Ca i_Cl], [1 2], mol_masses)
M_NaCl = create_mol([i_Na i_Cl], [1 1], mol_masses)
M_MgSO4  = create_mol([i_Mg i_S i_O], [1 1 4], mol_masses)
M_NaHCO3 = create_mol([i_Na i_H i_C i_O], [1 1 1 3], mol_masses)
M_H20 = create_mol([i_H i_O], [2 1], mol_masses)

# ignore oxygen
mask = create_mol([i_Ca i_Mg i_S i_Cl i_Na], [1 1 1 1 1], mol_masses) > 0

T = create_mol([i_Ca i_Mg i_Na i_S i_Cl], [0.05 0.01 0.02 0.05 0.10], mol_masses)
A = create_mol([i_Ca i_Mg i_Na i_S i_Cl], [0.094 0.018 0.004 0.044 0.005], mol_masses)
A = 0*T

Pa = [0 0 0 0]';
Pp = [M_CaSO4; M_CaCl2; M_MgSO4; M_NaCl];

CONV_COUNTER_RELOAD = 10;
Dn0 = 1;
conv_counter = CONV_COUNTER_RELOAD;
while true,
  [A1, Pa, Dn1] = vecfit(Pp,mask,T,A,Pa, 0.0001.*(1-exp(-Dn0)))
  Dc = round_n(Dn1,3);
  if Dn0 == Dc,
    if conv_counter == 0,
      break;
    else
      conv_counter -= 1;
    end
  else
    conv_counter = CONV_COUNTER_RELOAD;
  end
  Dn0 = Dc;
end
T=T
A1=A1
Pa=Pa
endfunction

function zr = round_n(z, n)
   zr = round(z*10^n)/10^n;
endfunction