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

void main()
{
    //program d9r8a
    //driver for routine mrqmin
    int i,j,npt = 100;
    int mfit,ma = 6;
    double ochisq,chisq,spread = 0.001;
    double x[101], y[101], sig[101], a[7];
	int lista[7];
    double covar[37], alpha[37], gues[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.5; gues[2] = 2.2; gues[3] = 2.8;
    gues[4] = 2.5; gues[5] = 4.9; gues[6] = 2.8;
    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;
    }
    double alamda = -1;
    for (i = 1; i<=ma; i++)
	{
        a[i] = gues[i];
    }
    mrqmin(x, y, sig, npt, a, ma, lista, mfit, covar, alpha, ma, chisq, alamda);
    int k = 1;
    int itst = 0;
	cout.setf(ios::left);
	do
	{
	    cout<<"Iteration #"<<k<<"   Chi-squared: ";
		cout<<chisq;
		cout<<"   alamda: "<<alamda;
		cout<<endl;
		cout<<endl;
		cout<<" a(1)     a(2)    a(3)    a(4)    a(5)    a(6)"<<endl;
		cout<<endl;
		for (i = 1; i<=ma; i++)
		{
			cout.precision(5);
			cout<<"  ";
			cout.width(9);
			cout<<a[i];
		}
		cout<<endl;
		k = k + 1;
		ochisq = chisq;
		mrqmin(x, y, sig, npt, a, ma, lista, mfit, covar, alpha, ma, chisq, alamda);
		if (chisq > ochisq)
		{
			itst = 0;
		}
		else
		{
			if (fabs(ochisq - chisq) < 0.1)
			{
				itst = itst + 1;
			}
		}
		cout<<"press any key to continue...."<<endl;
		getch();
	}while(itst < 2);
    alamda = 0.0;
    mrqmin(x, y, sig, npt, a, ma, lista, mfit, covar, alpha, ma, chisq, alamda);
    cout<<endl;
	cout<<endl;
	cout<<"Uncertainties:"<<endl;
	cout<<endl;
    for (i = 1; i<=ma; i++)
	{
		cout<<"  ";
		cout.width(9);
		cout.precision(3);
        cout<<sqrt(covar[(i-1)*ma+i]);
    }
	cout<<endl;
	cout<<endl;
    cout<<"Expected results:"<<endl;
	cout<<endl;
    cout<<" 5.0      2.0      3.0      2.0      5.0      3.0"<<endl;
	cout<<endl;
}