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

void main()
{
    //program d9r8b
    //driver for routine mrqcof
    int i,j,mfit,npt = 100;
    int ma = 6;
    double chisq,spread = 0.1;
    double x[101], y[101], sig[101], a[7];
	int lista[7];
    double alpha[37], gues[7], beta[7];
    a[1] = 5.0; a[2] = 2.0; a[3] = 3.0; a[4] = 2.0; a[5] = 5.0; a[6] = 3.0;
    gues[1] = 4.9; gues[2] = 2.1; gues[3] = 2.9;
    gues[4] = 2.1; gues[5] = 4.9; gues[6] = 3.1;
    long idum = -911;

    //first try a sum of two gaussians
    for (i = 1; i<=100; i++)
	{
        x[i] = 0.1 * i;
        y[i] = 0.0;
        for (j = 1; j<=4; j=j+3)
		{
            y[i] = y[i] + a[j] * exp(-pow(((x[i] - a[j + 1]) / a[j + 2]), 2));
        }
        y[i] = y[i] * (1.0 + spread * gasdev(idum));
        sig[i] = spread * y[i];
    }
    mfit = ma;
    for (i = 1; i<=mfit; i++)
	{
        lista[i] = i;
    }
    for (i = 1; i<=ma; i++)
	{
        a[i] = gues[i];
    }
    mrqcof(x, y, sig, npt, a, ma, lista, mfit, alpha, beta, ma, chisq);
    cout<<endl;
    cout<<"matrix alpha"<<endl;
    cout<<endl;
	cout.setf(ios::fixed|ios::left);
	cout.precision(4);
    for (i = 1; i<=ma; i++)
	{
        for (j = 1; j<=ma; j++)
		{
			cout.width(12);
            cout<<alpha[(i-1)*ma+j];
        }
        cout<<endl;
    }
    cout<<endl;
    cout<<"Vector beta"<<endl;
    cout<<endl;
    for (i = 1; i<=ma; i++)
	{
		cout.width(12);
        cout<<beta[i];
    }
    cout<<endl;
	cout<<endl;
    cout<<"Chi-squared: "<<chisq;
	cout<<endl;
    cout<<endl;
    //next fix one line and improve the other
    for (i = 1; i<=3; i++)
	{
        lista[i] = i + 3;
    }
    mfit = 3;
    for (i = 1; i<=ma; i++)
	{
        a[i] = gues[i];
    }
    mrqcof(x, y, sig, npt, a, ma, lista, mfit, alpha, beta, ma, chisq);
    cout<<"Matrix alpha"<<endl;
    cout<<endl;
    for (i = 1; i<=mfit; i++)
	{
        for (j = 1; j<=mfit; j++)
		{
			cout.width(12);
            cout<<alpha[(i-1)*ma+j];
        }
        cout<<endl;
    }
    cout<<endl;
    cout<<"Vector beta"<<endl;
    cout<<endl;
    for (i = 1; i<=mfit; i++)
	{
		cout.width(12);
        cout<<beta[i];
    }
    cout<<endl;
	cout<<endl;
    cout<<"Chi-squared: "<<chisq<<endl;
	cout<<endl;
}