#include <stdio.h>
#include "nrutil.h"

void lfit(double x[], double y[], double sig[], int ndata, double a[],
	int ma, int lista[], int mfit, double **covar, double *chisq,
	void (*funcs)(double,double *,int))
{
	int k,kk,j,ihit,i;
	double ym,wt,sum,sig2i,**beta,*afunc;

	beta=matrix(ma,1);
	afunc=vector(ma);
	kk=mfit;
	for (j=0;j<ma;j++) {
		ihit=0;
		for (k=0;k<mfit;k++)
			if (lista[k] == j) ihit++;	// lista[]に順列が正しく入っているかどうかチェックし，入っていればihitをインクリメント
		if (ihit == 0)	// lista[]に順列が正しく入っているかどうかチェックし，不足があれば補う．
			lista[kk++]=j;	// mfit個のlista[]の設定がひとつもなければ，残りのma-mfit個のパラメータに無理矢理lista[]を設定する．
		else if (ihit > 1) nrerror("Bad lista permutation in LFIT-1");
	}
	if (kk != (ma)) nrerror("Bad lista permutation in LFIT-2");	// lista[]のパラメータ数が合うか最終確認．
	for (j=0;j<mfit;j++) {
		for (k=0;k<mfit;k++) 
			covar[j][k]=0.0;	// covar[][]行列の初期設定
		beta[j][0]=0.0;
	}

	for (i=0;i<ndata;i++) {
		(*funcs)(x[i],afunc,ma);	// x[]:変数，afunc: x,x^2,x^3, ma:ParaNum
		ym=y[i];
		if (mfit < ma)
			for (j=mfit;j<ma;j++)
				ym -= a[lista[j]]*afunc[lista[j]];
		sig2i=1.0/POW(sig[i]);	// POWは2乗．よってsig2iは分散
		for (j=0;j<mfit;j++) {
			wt=afunc[lista[j]]*sig2i;
			for (k=0;k<=j;k++)
				covar[j][k] += wt*afunc[lista[k]];
			beta[j][0] += ym*wt;
		}
	}
	if (mfit > 1)
		for (j=1;j<mfit;j++)
			for (k=0;k<=j-1;k++)
				covar[k][j]=covar[j][k];	// covar[][]は対称だから，下三角要素を上三角要素に代入している．
	gaussj(covar,mfit,beta,1);				// 行列方程式を解く, beta[j][0]が行列方程式の解！
	for (j=0;j<mfit;j++) a[lista[j]]=beta[j][0];	// 解を適当な係数a[]に代入．
	*chisq=0.0;
	for (i=0;i<ndata;i++) {
		(*funcs)(x[i],afunc,ma);
		for (sum=0.0,j=0;j<ma;j++) sum += a[j]*afunc[j];
		*chisq += POW((y[i]-sum)/sig[i]);
	}
	covsrt(covar,ma,lista,mfit);
	free_vector(afunc);
	free_matrix(beta);
}