mrqcof.cpp

来自「工程算法 这是一个很有用的工程数值算法集锦」· C++ 代码 · 共 41 行

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