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

double ran1(long& idum)
{
    int j,iff=-1;
	static long ix1,ix2,ix3;
	static double r[98];
    long m1 = 259200; long m2 = 134456; long m3 = 243000;
	long ia1 = 7141; long ia2 = 8121; long ia3 = 4561;
	long ic1 = 54773; long ic2 = 28411; long ic3 = 51349;
	double rm1 = 0.0000038580247; double rm2 = 0.0000074373773;    
    if (idum < 0 || iff == 0)
	{
        iff = 1;
        ix1 = (ic1 - idum) % m1;
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix2 = ix1 % m2;
        ix1 = (ia1 * ix1 + ic1) % m1;
        ix3 = ix1 % m3;
        for (j = 1; j<=97; j++)
		{
            ix1 = (ia1 * ix1 + ic1) % m1;
            ix2 = (ia2 * ix2 + ic2) % m2;
            r[j] = (double(ix1) + double(ix2) * rm2) * rm1;
        }
        idum = 1;
    }
    ix1 = (ia1 * ix1 + ic1) % m1;
    ix2 = (ia2 * ix2 + ic2) % m2;
    ix3 = (ia3 * ix3 + ic3) % m3;
    j = 1 + int((97 * ix3) / m3);
    if (j > 97 || j < 1)
	{
		cout<<"abnormal exit in ran1"<<endl;
		exit(1);
	}
	double temp=r[j];
    r[j] = (double(ix1) + double(ix2) * rm2) * rm1;
	return temp;
}

double gasdev(long& idum)
{
    static int iset;
	static double gset;
	double v1,v2,r,fac;
    if (iset == 0)
	{
		do
		{
			v1 = 2.0 * ran1(idum) - 1.0;
			v2 = 2.0 * ran1(idum) - 1.0;
			r = v1 * v1 + v2 * v2;
		}while (r >= 1.0 || r == 0);
		fac = sqrt(-2.0 * log(r) / r);
		gset = v1 * fac;
		iset = 1;
		return v2 * fac;
	}
    else
	{
		iset = 0;
		return gset;
    }
}
void fgauss(double x, double a[], double& y, double dyda[], int na)
{
	y = 0.0;
	for (int i = 1; i<=na - 1; i=i+3)
	{
		double arg = (x - a[i + 1]) / a[i + 2];
		double ex = exp(-(arg * arg));
		double fac = a[i] * ex * 2.0 * arg;
		y = y + a[i] * ex;
		dyda[i] = ex;
		dyda[i + 1] = fac / a[i + 2];
		dyda[i + 2] = fac * arg / a[i + 2];
	}
}

void mrqcof(double x[], double y[], double sig[], int ndata, double a[],
			int ma, int lista[], int mfit, double alpha[], double beta[],
			int nalp, double& chisq)
{
	int i,j,k;
	double wt,ymod,sig2i,dy;
    double  dyda[21];
    for (j = 1; j<=mfit; j++)
	{
        for (k = 1; k<=j; k++)
		{
            alpha[(j-1)*nalp+k] = 0.0;
        }
        beta[j] = 0.0;
    }
    chisq = 0.0;
    for (i = 1; i<=ndata; i++)
	{
        fgauss(x[i], a, ymod, dyda, ma);
        sig2i = 1.0 / (sig[i] * sig[i]);
        dy = y[i] - ymod;
		for (j = 1; j<=mfit; j++)
		{
            wt = dyda[lista[j]] * sig2i;
			for (k = 1; k<=j; k++)
			{
                alpha[(j-1)*nalp+k] = alpha[(j-1)*nalp+k] + wt * dyda[lista[k]];
            }
            beta[j] = beta[j] + dy * wt;
        }
        chisq = chisq + dy * dy * sig2i;
    }
    for (j = 2; j<=mfit; j++)
	{
        for (k = 1; k<=j - 1; k++)
		{
            alpha[(k-1)*nalp+j] = alpha[(j-1)*nalp+k];
        }
    }
}

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;
}