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

void get_data (int* num_points, double* z, double* n)
{

  ifstream dataIn ("perrin.dat");
  int num;
  dataIn >> num;
  *num_points = num;
  for(int i=0; i<num; i++)
  {
    dataIn >> z[i] >> n[i];
  }

}

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

void do_fit ()
{

  extern int num_points;

  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 nu0 and k, give initial values and step sizes.

  MNPARM (1, "nu0", 1000., 100., f_null, f_null, error_flag);
  MNPARM (2, "k", 1.4, 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);
  MNPOUT (1, chnam, par[0], dpar[0], bnd1, bnd2, ivarbl);
  MNPOUT (2, chnam, par[1], dpar[1], bnd1, bnd2, ivarbl);

  double kB = par[1]*1.0e-23;           // J/K
  double dkB = dpar[1]*1.0e-23;
  const double R_GAS = 8.32;            // J/K/mole
  double NAv = R_GAS/kB;
  double dNAv = NAv * dpar[1]/par[1];

  cout << endl;
  cout << "Fit results:" << endl;
  cout << endl;
  cout << "nu0             = " << par[0] << " +- " << dpar[0] << endl;
  cout << "kB              = " << kB << " +- " << dkB << endl;
  cout << "logL            = " << - fmin / 2. << endl;;
  cout << "NAvagadro     = " << NAv << " +- " << dNAv << endl;

}