#include "iostream.h"
#include "math.h"
#include "stdlib.h"

void funcs(double x, double afunc[], int ma)
{
	int i;
    afunc[1] = 1.0;
    for (i = 2; i<=ma; i++)
	{
        afunc[i] = x * afunc[i - 1];
    }
}

void main()
{
    //program d9r2
    //driver for routine lfit
    int i,j,npt = 100;
    double chisq,spread = 0.1;
    int mfit,nterm = 3;
    double x[101], y[101], sig[101], a[4], covar[10];
	int lista[4];
    long idum = -911;
    for (i = 1; i<=npt; i++)
	{
        x[i] = 0.1 * i;
        y[i] = nterm;
        for (j = nterm - 1; j>=1; j--)
		{
            y[i] = j + y[i] * x[i];
        }
        y[i] = y[i] + spread * gasdev(idum);
        sig[i] = spread;
    }
    mfit = nterm;
    for (i = 1; i<=mfit; i++)
	{
        lista[i] = i;
    }
    lfit(x, y, sig, npt, a, nterm, lista, mfit, covar, nterm, chisq);
    cout<<"Parameter               Uncertainty"<<endl;
	cout.setf(ios::scientific|ios::left);
	cout.precision(4);
    for (i = 1; i<=nterm; i++)
	{
        cout<<"a("<<i<<")= ";
		cout.width(18);
		cout<<a[i];
		cout.width(18);
        cout<<sqrt(covar[(i-1)*nterm+i]);
		cout<<endl;
    }
    cout<<endl;
    cout<<"Chi-squared =  "<<chisq<<endl;
    cout<<endl;
    cout<<"Full covariance matrix"<<endl;
    for (i = 1; i<=nterm; i++)
	{
        for (j = 1; j<=nterm; j++)
		{
			cout.width(18);
            cout<<covar[(i-1)*nterm+j];
        }
		cout<<endl;
    }
    //now test the lista feature
    for (i = 1; i<=nterm; i++)
	{
        lista[i] = nterm + 1 - i;
    }
    lfit(x, y, sig, npt, a, nterm, lista, mfit, covar, nterm, chisq);
    cout<<endl;
    cout<<"Parameter               Uncertainty"<<endl;
    for (i = 1; i<=nterm; i++)
	{
        cout<<"a("<<i<<")= ";
		cout.width(18);
		cout<<a[i];
		cout.width(18);
		cout<<sqrt(covar[(i-1)*nterm+i])<<endl;
    }
    cout<<endl;
    cout<<"Chi-squared =  "<<chisq<<endl;
    cout<<endl;
    cout<<"Full covariance matrix"<<endl;
    for (i = 1; i<=nterm; i++)
	{
        for (j = 1; j<=nterm; j++)
		{
			cout.width(18);
            cout<<covar[(i-1)*nterm+j];
        }
		cout<<endl;
    }
    //now check results of restricting fit parameters
    int ii = 1;
	int aaa;
    for (i = 1; i<=nterm; i++)
	{
        aaa = i - int(i / 2) * 2;
        if (aaa == 1)
		{
            lista[ii] = i;
            ii = ii + 1;
        }
    }
    mfit = ii - 1;
    lfit(x, y, sig, npt, a, nterm, lista, mfit, covar, nterm, chisq);
    cout<<endl;
    cout<<"Parameter               Uncertainty"<<endl;
    for (i = 1; i<=nterm; i++)
	{
        cout<<"a("<<i<<")= ";
		cout.width(18);
		cout<<a[i];
		cout.width(18);
		cout<<sqrt(covar[(i-1)*nterm+i])<<endl;
    }
    cout<<endl;
    cout<<"Chi-squared =  "<<chisq<<endl;
    cout<<endl;
    cout<<"Full covariance matrix"<<endl;
    for (i = 1; i<=nterm; i++)
	{
        for (j = 1; j<=nterm; j++)
		{
			cout.width(18);
            cout<<covar[(i-1)*nterm+j];
        }
		cout<<endl;
    }
}