// plot pH as a function of concentration
// for various pKa values (including weak and strong acids)
//
// writes:
//    pa-pc.csv     plotted versus concentration for 8 different Ka values
//    ph-pc.csv     ditto
//    pha-pc.csv    ditto
//    ph-pc-allroots.csv
// see also ./ph-scan.gnumeric
// see also ./acid-pa-scan.c

// References:
// http://www.av8n.com/physics/ph-versus-concentration.htm
// http://laude.cm.utexas.edu/courses/ch302/lecture/ln15s09.pdf
// http://www.chem1.com/acad/pdf/c1xacid2.pdf

// http://www2.csudh.edu/oliver/chemdata/data-ka.htm
// Name	                Formula	Ka1             pKa1	Ka2             pKa2
// Acetic acid          CH3CO2H	1.75 × 10−5	4.756
// Sulfuric acid	H2SO4	Strong          Strong	1.0 × 10−2	1.99

//https://chemistry.stackexchange.com/questions/50880/...
// combining-acid-dissociation-constants-to-determine-ph-of-diprotic-acid

using namespace std;
#include <iostream>
#include <fstream>
#include "brent.hpp"
#include <cmath>
#include <stdlib.h>       /* for exit() */
#include <vector>

void usage() {
  string msg = R"EoF(
Calculates pH as a function of Ta (total acid), for a range of pKa values,
or optionally for one or more given pKa values.

Works for strong or weak acids.
Works for high or low concentrations.
Works even for diprotic acids.

Options include:
  -h | -help    # print this messaged (and immediately exit
  acetic        # calculate the curve for acetic acid
  sulfuric      # calculate the curve for sulfuric acid
  -details      # add 8 columns of additional info in output file


The resulting data can be plotted and analyzed using:
  ~/hobbies/tap-water-titration-acid-base.gnumeric

In a standard car battery, the electrolyte is a mixture of around 35% sulfuric acid and 65% water by weight. This leads to an approximate molarity of about 4.2 M and a density of 1.28 g/cm³

HDX #2 is phenol red < 0.1%
https://www.manualshelf.com/manual/hdx/62364/sds-english/page-15.html

HDX #5 contains Bromophenol Blue 0.1% and Sodium Thiosulfate 0.24%
https://www.manualshelf.com/manual/hdx/62364/sds-english/page-47.html
In pH testing of bleach substances, sodium thiosulfate neutralizes the color-removing effects of bleach and allows one to test the pH of bleach solutions with liquid indicators. The relevant reaction is akin to the iodine reaction: thiosulfate reduces the hypochlorite (the active ingredient in bleach) and in so doing becomes oxidized to sulfate.

BPB solutions vs. PH without (a) and with (b) CTApTs (0.0123 M), with increasing concentrations of indicator: (1) 1.42 × 10−5 M; (2) 2.91 × 10−5 M; (3) 5.72 × 10−5 M.

)EoF";

  cout << msg.substr(1);
}

template <typename T>
inline int sgn(const T& value) {
  return value < 0 ? -1 : value > 0;
}

//template <typename T>
// inline T max(T a, T b) { return a >= b ? a : b; }

double const NaN(nan(""));
double bracket(double const val, double const top, double const bot) {
  if (val == top) return NaN;
  if (val == bot) return NaN;
  return val;
}

struct Job {
  double pKa1;
  double pKa2;
  double pKa3;
};

int main(int const argc, char * const * argv){
  double NS = 4.;
  int details(0);
  vector<Job> todo;
  for (int ii = 1; ii < argc; ii++) {
    string arg(argv[ii]);
    if (arg[0] == '-') arg = arg.substr(1);
    if (0) {}
    else if (arg == "h" || arg == "help") {
      usage();
      exit(0);
    }
    else if (arg == "NS" || arg == "ns") {
      ii++;
      NS = atof(argv[ii]);
    }
    else if (arg == "pKa") {
      ii++;
      string val = argv[ii];
      todo.push_back({atof(val.c_str()), INFINITY, INFINITY});
    }
    else if (arg == "acetic") {
      todo.push_back({4.745, INFINITY, INFINITY});
    }
    else if (arg == "sulfuric") {
      todo.push_back({-3, 1.99, INFINITY});
    }
    else if (arg == "details") {
      details++;
    }
    else {
      cerr << "Unrecognized option: " << arg << endl;
      exit(1);
    }
  }

  if (todo.size() == 0) {
    for (double pKa = -4; pKa <= 10; pKa++){
      todo.push_back({pKa, INFINITY, INFINITY});
    }
  }

  for (Job &job : todo) {
    cout << "pKa1: " << job.pKa1
       << "  pKa2: " << job.pKa2
        << endl;
  }

  string roots_fn = "ph-pc-allroots.csv";
  ofstream roots(roots_fn.c_str());

  string ph_pc_fn = "ph-pc.csv";
  ofstream ph_pc(ph_pc_fn.c_str());

  double Kw = pow(10., -14.);
    for (Job job : todo) {
    ph_pc  << ",pKa";
    if (details) ph_pc << ",,,,,,,,";
  }
  ph_pc << endl;        // end of first line

  for (Job job : todo) {
    ph_pc << ", " << job.pKa1;
    if (details) ph_pc << ",,,,,,,,";
  }
  ph_pc << endl;        // end of pKa1 values line

  for (Job job : todo) {
    ph_pc << ", " << job.pKa2;
    if (details) ph_pc << ",,,,,,,,";
  }
  ph_pc << endl;        // end of pKa2 values line

  for (Job job : todo) {
    ph_pc << ", " << job.pKa3;
    if (details) ph_pc << ",,,,,,,,";
  }
  ph_pc << endl;        // end of pKa3 values line

  ph_pc << "pTa";
  for (Job &job : todo) {
    ph_pc << ", pH";
    if (details) {
      ph_pc << ",weak,strong,a0,a1,a2,OH,acid,charge";
    }
  }
  ph_pc  << endl;       // end of main body headers line

  brent::Poly Hpoly(5);
  brent::Poly Hderiv(4);

  double const sigma(0);
// loop over lines in output file,
// i.e. scan over concentrations.
  for (double pTax = -1*NS; pTax <= 10*NS; pTax += 1.){
    double pTa = pTax / NS;
    double Ta = pow(10., -pTa);
    int didTa(0);
// loop over columns in output file,
// i.e. loop ovber pKa values:
    for (Job &job : todo) {
      if (!didTa) {
        ph_pc << -log10(Ta);
        didTa++;
      }
      double Ka1 = pow(10., -job.pKa1);
      double Ka2 = pow(10., -job.pKa2);
      Hpoly.coeff[4] = 1;                     // [H]^4 term.
      Hpoly.coeff[3] = Ka1;                   // [H]^3 term.
      Hpoly.coeff[2] = -(Ka1*Ta + Kw - Ka1*Ka2);
      Hpoly.coeff[1] = -(Ka1*Kw + 2*Ka1*Ka2*Ta);
      Hpoly.coeff[0] = -Ka1*Ka2*Kw;               // constant term.

      for (int ii = 0; ii < Hderiv.coeff.size(); ii++) {
        Hderiv.coeff[ii] = (1 + ii) * Hpoly.coeff[1 + ii];
      }

// reference:
// http://people.sc.fsu.edu/~jburkardt/cpp_src/brent/brent.html
      double xstrong = Ta;
      double fancy_strong = sqrt(Kw + Ta*Ta);
//>>> Even fancier:
//>>> WS = (weak^-1.5 + strong^-1.5)^(-1/1.5)
//>>>> Ku = sqrt(Kw)
//>>> WSD = (WS^1.5 + Ku^1.5)^(1/1.5)
//>>> see   ~/hobbies/tap-water-titration-acid-base.gnumeric
//
// FIXME factor of 2 for diprotic acid KLUDGE
      double Htop = fancy_strong * 2;
      double weak = sqrt(Ka1 * Ta);
      double bot(min(weak, fancy_strong) / 2.);
      bot = 1e-20;
      if (0) cerr << "pKa1: " << job.pKa1
           << "  pTa: " << pTa
           << "  Htop: " << Htop
           << "  Hpoly(Htop): " << Hpoly(Htop)
           << "  Hpoly(0): " << Hpoly(0)
           << endl;

      if (Hpoly(Htop) < -1e-20) {
        cerr << "****** oops: Hpoly(Htop): " << Hpoly(Htop) << endl;
        cerr << "      Hpoly(0): " << Hpoly(0) << endl;
        exit(1);
      }
      double Htolerance = min(sqrt(Kw), Ka1)/1e10;
      double Hroot = zero(bot, Htop, Htolerance, Hpoly);
      double check_a0 = Ta * Hroot*Hroot / (Ka1*Ka2 + Ka1*Hroot + Hroot*Hroot);
      double check_a1 = check_a0 * Ka1 / Hroot;
      double check_a2 = check_a1 * Ka2 / Hroot;
      double check_OH = Kw / Hroot;
      double conserve_acid = (check_a2 + check_a1 + check_a0 - Ta) / Ta;
      double conserve_charge = (Hroot - 2*check_a2 - check_a1 - check_OH)
                                 / fancy_strong;

      cout      << "  Ta: " << Ta
                << "  Ka1: " << Ka1
                << "  root: " << Hroot
                << "  miss: " << Hpoly(Hroot)
                << "  deriv: " << Hderiv(Hroot)
                << endl;

      if (abs(conserve_charge) > 1e-10) {
        cerr << "bad conserve_charge " << conserve_charge << endl;
        cerr << (Hroot - 2*check_a2 - check_a1 - check_OH)
                << " versus " << fancy_strong << endl;
        exit(1);
      }
      if (abs(conserve_acid) > 1e-10) {
        cerr << "bad conserve_acid " << conserve_acid << endl;
        exit(1);
      }

      ph_pc << ", " << -log10(Hroot);
// details:
      if (details) {
        ph_pc << ", " << -log10(weak);
        ph_pc << ", " << -log10(xstrong);
        ph_pc << ", " << -log10(check_a0);
        ph_pc << ", " << -log10(check_a1);
        ph_pc << ", " << -log10(check_a2);
        ph_pc << ", " << -log10(check_OH);
        ph_pc << ", " << conserve_acid;
        ph_pc << ", " << conserve_charge;
      }

//xxxx      double A0 = Ta*Hroot*Hroot / (Ka1*Ka2 + Ka1*Hroot + Hroot*Hroot);
//xxxx      double A1 = A0 * Ka1 / Hroot;

#if 0
/////////////////////////////////
// find the other two roots
// synthetic division is cute but not sufficiently accurate:
#ifdef synth
          double k2 = 1.;           // not really needed
          double k1 = a + k2 * Hroot;
          double k0 = b + k1 * Hroot;
          double remainder = c + k0 * Hroot;   // should be very small

    // now use quadratic formula to find roots of what's left:
          double bigroot = -0.5 * (k1 + sgn(k1)*sqrt(k1*k1 - 4*k0));
          if (1) roots << bigroot << ", ";
          if (0) roots << Hpoly(bigroot) << ", ";
          double smallroot = k0/bigroot;
          if (1) roots << smallroot << ", ";
          if (0) roots << Hpoly(smallroot) << ", ";
#endif
// find estimates that bracket the root positions:
      double midlow = -Ka*Kw / (Kw + Ka*Ta);
      if (Hpoly(midlow) < 0) {
        cerr << "Hpoly(midlow): " << Hpoly(midlow) << endl;
        cerr << "Hpoly(0): " << Hpoly(0) << endl;
        exit(1);
      }
      double verylow = -Ka - Ta - 1.;
      if (Hpoly(verylow) > 0) {
        cerr << "Hpoly(verylow): " << Hpoly(verylow) << endl;
        cerr << "verylow: " << verylow
                << endl;
        exit(1);
      }
// find actual roots by iteration:
      double smallroot2 = zero(midlow, 0, abs(midlow)/1e9, Hpoly);
      double negroot2 = zero(verylow, midlow, abs(verylow)/1e9, Hpoly);
      roots << Hroot << ", ";
      roots << negroot2 << ", ";
      roots << smallroot2 << ", ";
      roots << ", ";
#endif
    }
    roots << endl;
    ph_pc << endl;
  }
}

/*
@ #include <gsl/gsl_poly.h>
@ int bad_main(){
@   gsl_complex z0, z1, z2;
@   double Kw = pow(10., -14.);
@   cout << "pKa:, ";
@   for (double pKa = -4; pKa <= 10; pKa++){
@     cout << pKa << ", ";
@   }
@   cout << endl;
@   cout << "v pTa v" << endl;    // mostly blank line
@
@   for (double pTa = -1; pTa <= 10; pTa++){
@     cout << pTa << ", ";
@     double Ta = pow(10., -pTa);
@     for (double pKa = -4; pKa <= 10; pKa++) {}
@     for (double pKa = -4; pKa <= -4; pKa++){
@       double Ka = pow(10., -pKa);
@       double a = Ka;
@       double b = -Kw - Ka*Ta;
@       double c = -Ka*Kw;
@       gsl_poly_complex_solve_cubic(a, b, c, &z0, &z1, &z2);
@       cout << (GSL_REAL(z2)) << ", ";
@       cout << (GSL_IMAG(z2)) << ", ";
@
@       cout << (GSL_REAL(z1)) << ", ";
@       cout << (GSL_IMAG(z1)) << ", ";
@
@       cout << (GSL_REAL(z0)) << ", ";
@       cout << (GSL_IMAG(z0)) << ", ";
@
@       double Hroot = GSL_REAL(z2);
@       cout << Hroot*Hroot*Hroot + a*Hroot*Hroot + b*Hroot + c << ", ";
@       cout << a << ", ";
@       cout << b << ", ";
@       cout << c << ", ";
@       cout << Hroot*Hroot*Hroot + a*Hroot*Hroot + b*Hroot + c << ", ";
@
@  // int gsl_poly_complex_solve_cubic
@  //     (double a, double b, double c,
@  //     gsl_complex * z0, gsl_complex * z1, gsl_complex * z2)
@     }
@     cout << endl;
@   }
@   return 0;
@ }
*/