#include <iostream.h>
#include <fstream.h>
#include <cfortran/cfortran.h>
#include <cfortran/minuit.h>
#include "fit.hh"

void get_data (int* num_points, double* x, double* y, double* dy)
{

  ifstream dataIn ("galileo.dat");
  int n;
  dataIn >> n;
  *num_points = n;
  for(int i=0; i<n; i++)
  {
    dataIn >> x[i] >> y[i];
  }

  cout << "enter error in punti" << endl;
  double sigma;
  cin >> sigma;
  for(int i=0; i<n; i++)
  {
    dy[i] = sigma;
  }

}

extern "C" {void *fcn_();}

void do_fit ()
{

  extern int num_points;
  int n_dof;
  double chi2;

  const int NPAR = 2;
  char chnam[10];
 
  int  ierr;
  int  istat;
  int  ivarbl;
  int  npari, nparx;

  int error_flag = 0;
  double f_null = 0.;
  const int MAX_ARGS = 10;
  double arglis[MAX_ARGS];
  double bnd1, bnd2;
  double deriv[NPAR];
  double dpar[NPAR];
  double fmin, fedm, errdef;
  double par[NPAR];

  int ird = 5;                 // unit number for input to Minuit (keyboard)
  int iwr = 6;                 // unit number for output from Minuit (screen)
  int isav = 7;                // unit number for use of the SAVE command
  MNINIT (ird, iwr, isav);

//  Define parameters alpha and beta, give initial values and step sizes.

  MNPARM (1, "alpha", 50., 1., f_null, f_null, error_flag);
  MNPARM (2, "beta", 0.5, 0.1, f_null, f_null, error_flag);

//  Minimize chi2 using SIMPLEX and MIGRAD, get cov. mat. with HESSE

  MNEXCM (fcn_, "SIMPLEX", 0, 0, error_flag, 0);
  MNEXCM (fcn_, "MIGRAD", 0, 0, error_flag, 0);
  MNEXCM (fcn_, "HESSE", 0, 0, error_flag, 0);

//  Get results of fit

  MNSTAT (fmin, fedm, errdef, npari, nparx, istat);
  chi2 = fmin;
  n_dof = num_points - npari;
  MNPOUT (1, chnam, par[0], dpar[0], bnd1, bnd2, ivarbl);
  MNPOUT (2, chnam, par[1], dpar[1], bnd1, bnd2, ivarbl);

  cout << endl;
  cout << "Fit results:" << endl;
  cout << endl;
  cout << "alpha           = " << par[0] << " +- " << dpar[0] << endl;
  cout << "beta            = " << par[1] << " +- " << dpar[1] << endl;
  cout << "chi2            = " << chi2 << endl;
  cout << "n_dof           = " << n_dof << endl;
  cout << "chi2/n_dof      = " << chi2/double(n_dof) << endl;

}