#include <cmath>
#include <iostream>
#include <fstream>

using namespace std;

double alpha_M,
       ca, // analit derişimi 
       ct; // titrant derişimi

void cunh3(const double c_cu, const double v_cu, const double c_nh3, 
     char *da);
void f_ve_df(double x, double &f, double &df);

int main()
{
   cunh3(0.01, 25.0, 0.01, "cu0.01amin0.01");
   cunh3(0.05, 25.0, 0.05, "cu0.05amin0.05");
   cunh3(0.1, 25.0, 0.1, "cu0.1amin0.1");
   return(0);
}

void cunh3(const double c_cu, const double v_cu, const double c_nh3, 
     char *da)
{
   int N = 250;
   double v_nh3, dv_nh3=v_cu/N, c, x, x_eski, err, eps=1.0e-4, f, df;

   ofstream fp(da);
   for(v_nh3 = 0.0; v_nh3 <= 6.0*v_cu; v_nh3 += dv_nh3){
      ca = c_cu*v_cu/(v_cu+v_nh3);
      ct = c_nh3*v_nh3/(v_cu+v_nh3);
      x=1.0e-14;
      do{
         x_eski = x;
         f_ve_df(x_eski, f, df);
         x = x_eski - f/df;
         err = fabs(x - x_eski)/fabs(x_eski);
      } while(err > eps);
      c = alpha_M*ca;
      fp << v_nh3 << " " << -log10(c) << " " << -log10(x)  << endl;
   }
   fp.close();
}

void f_ve_df(double x, double &f, double &df)
{
   double Kw = 1.0e-14, Kb = 1.75e-5, beta1 = pow(10.0, 4.04),
          beta2 = pow(10.0,7.47), beta3 = pow(10.0,10.27),
          beta4 = pow(10.0,11.75);
   double y, dy, n, d, dn, dd;

   y = (Kw/Kb)*( (Kw/(x*x)) - 1.0);
   dy = -2.0*Kw*Kw/(Kb*x*x*x);
   n = beta1*y + 2.0*beta2*y*y + 3.0*beta3*y*y*y + 4.0*beta4*y*y*y*y;
   dn = beta1 + 4.0*beta2*y + 9.0*beta3*y*y + 16.0*beta4*y*y*y;
   d = 1.0 + beta1*y + beta2*y*y + beta3*y*y*y + beta4*y*y*y*y;
   dd = n/y;
   f = y + (Kb/Kw)*x*y + (n/d)*ca - ct;
   df = dy + (Kb/Kw)*(y+x*dy) + ca*dy*(dn*d -n*dd)/(d*d);
   alpha_M = 1.0/d;
}