// Routine for fitting to an exponential (for loading curves)
// Using Numerical Recipes in C, 15.5 Nonlinear Models Example

#include <fstream.h>
#include <time.h>
#include <stdlib.h>
#include <math.h>

// NUMERICAL RECIPES ROUTINES - CHECK PATHNAMES BEFORE USING

#include "/home/gensemer/NumRec/Other/NRUTIL.C"
#include "/home/gensemer/NumRec/Other/nr.h"
#include "/home/gensemer/NumRec/ANSI/MRQCOF.C"
#include "/home/gensemer/NumRec/ANSI/MRQMIN.C"
#include "/home/gensemer/NumRec/ANSI/GAUSSJ.C"
#include "/home/gensemer/NumRec/ANSI/COVSRT.C"

// USER-SUPPLIED FUNCTION HERE

#include "/home/gensemer/datanalysis/analysis/betafit/beta.cpp" 

void
betafit (char *filename, char *outfile, int verbose)
{
  
    
  // OPEN INPUT AND OUTPUT FILES--------------------------------- 
    
    ifstream file (filename);
    if (!file)
    cerr << "Cannot open infile for input";

    ofstream file2 (outfile);
    if (!file2)
    cerr << "Cannot open outfile for output";
  
    int i, xmax, j;
    float *y,*x;
    y = vector (1, 25000);
    x = vector (1, 25000);
    
  // LOAD DATA FROM FILE----------------------------------------
    
    i = 0;
    while (file.good () == 1 && i <= 25000)
    {
	i++;
	file >> y[i];
    }
  
    const int ndata = i-1;
    if (i == 25000)
      cout << "Data Overflow (file contains > 25,000 points)\n";
  
    
  // DECLARE VARIABLES AND ARRAYS---------------------------------
    
    const int one = 1;
    const int ma = 3;		// number of fitting params in func


    float *yfit, *sig, *a, *dyda;

    yfit = vector (one, ndata);
    sig = vector (one, ndata);
    a = vector (one, ma);
    dyda = vector (one, ma);
    
    int *ia;
    ia = ivector (one, ma);
    
    float **covar, **alpha;
    covar = matrix (one, ma, one, ma);
    alpha = matrix (one, ma, one, ma);
    
    void (*funcs) ();
    
    float chi, oldchi, lamda, oldlamda, *chisq, *alamda;
    alamda = &lamda;
    chisq = &chi;
    
  // SET x[], sig[], AND ia[] VALUES---------------------------------
    
    for (i = 1; i <= ndata; i++)
    x[i] = i/10;				// x axis values (optional)
    
    for (i = 1; i <= ndata; i++)
      sig[i] = 0.02;		// estimated experimental uncertainty 
    
    for (i = 1; i <= ma; i++)
      ia[i] = 1;			// weighting
    
  // GUESS FIT PARAMETERS---------------------------------
 
  double v1,v2;
  v1 = 0;
  v2 = 0;

  a[1] = y[1];					// Initial value

  for (i=1;i<=5;i++)
	v1 += y[4*ndata/5-i]/5;

  for (i=1;i<=5;i++)
	v2 += y[ndata-i]/5;

  a[2] = log(v1/v2)/(ndata/50);
				// Gamma (exponential time constant)
  a[3] = a[2]*(1/v2 - exp(a[2]*ndata/10)/a[1])/(exp(a[2]*ndata/10) - 1);
				//Beta (non-exponential correction)	

  // INITIALIZE ROUTINE---------------------------------
    *alamda = -1;
  
    if (verbose == 1)
    {
	cout << "i	alamda	chisq	amatrix\n";
	cout << -1 << "   " << *alamda << "   " << *chisq << "   ";
	for (j = 1; j <= ma; j++)
	  cout << a[j] << "   ";
	cout << "\n";
    }
    
    i = 0;
    mrqmin (x, y, sig, ndata, a, ia, ma, covar, alpha, chisq, beta, 
	    alamda);
    if (verbose == 1)
    {
	cout << i << "   " << *alamda << "   " << *chisq << "   ";
	for (j = 1; j <= ma; j++)
	  cout << a[j] << "   ";
	cout << "\n";
    }
    i++;
    
  // FIT WHILE CONDITION---------------------------------
    
    oldchi = chi * 2;
  
    while (((oldchi - chi) / oldchi > .001 || chi > oldchi) && i < 40)
    {
	oldchi = chi;
	mrqmin (x, y, sig, ndata, a, ia, ma, covar, alpha, chisq, 
		beta, alamda);
	if (verbose == 1)
	{
	    cout << i << "   " << *alamda << " " << *chisq << "   ";
	    for (j = 1; j <= ma; j++)
	      cout << a[j] << "   ";
	    cout << "\n";
	}
	i++;
    }
  
    if (i == 40)
      cout << "No Convergence\n";
  
    
  // FINAL CALL TO DETERMINE CONVERGED VALUES----------------------------
    
    *alamda = 0;
    mrqmin (x, y, sig, ndata, a, ia, ma, covar, alpha, chisq, beta, 
	    alamda);
  
    
  // OUTPUT RESULTS TO STDOUT (pipe to file if needed)
    
    if (verbose == 1)
      file2 << "Results:" << filename << "\n i	a[i]	error[i]\n";

    for (i = 1; i <= ma; i++)      
       file2 << i << "   " << a[i] << "   " << sqrt (3.53 * covar[i][i]) << "\n";
 
}