#include "lib.h"


// Minimization wrapper.
//
double chi2(double p_max, void* params)
  {
  // Pars parameters.
  struct parameters p = *((struct parameters*) params);
  struct bin* histo = p.histo;
  size_t n = p.n;
  double step = p.step;

  // Compute χ².
  double chi = 0;
  double expecto = 0; 
  for (size_t i = 0; i < n; i++)
    {
    expecto = expected((i + 0.5) * step, p_max);
    chi += pow(histo[i].sum - expecto, 2)/expecto;
    };
  return chi;
  }


// Expected function.
//
double expected (double x, double p_max)
  {
  // When p_max < x, the argument under sqrt is negative. Return great number
  // to lead away the minimization.
  if (p_max < x)
    {
    double num = 500;
    return num;
    }
  return x * log(p_max/x)/atan(sqrt(fabs(pow(p_max, 2)/pow(x,2) - 1)));
  }


// First derivative of the expected function.
//
double exp1d (double x, double p_max)
  {
  // When p_max < x, the argument under sqrt is negative. Return great number
  // to lead away the minimization.
  if (p_max < x)
    {
    double num = 500;
    return num;
    }
  double a = p_max;
  double A = sqrt(fabs(pow(a, 2)/pow(x,2) - 1));
  double B = atan(A);
  double C = log(a/x);
  return x/(a*B) - x*C/(a*A*pow(B, 2));
  }


//
// Second derivative of the expected function.
//
double exp2d (double x, double p_max)
  {
  // When p_max < x, the argument under sqrt is negative. Return great number
  // to lead away the minimization.
  if (p_max < x)
    {
    double num = 500;
    return num;
    }
  double a = p_max;
  double a2 = pow(a, 2);

  double A = sqrt(fabs(a2/pow(x,2) - 1));
  double B = atan(A);
  double C = log(a/x);
  double D = 1 - a2/pow(x,2);

  double E = -2/(a2*D*B);
  double F = 1/(a2*A);
  double G = 1/(pow(x, 2)*pow(A, 3));
  double H = 1/a2 + 2/(a2 * A * B);

  return x*((E + F + G) * C/B - H)/B;
  }