ercheng.h

实现二次曲线的拟合.不是很完善

#include "math.h"
//appendix functions:
float gammln(float xx)
{
	int j;
	float temp;
    float cof[6],stp,half,one,fpf,x,tmp,ser;
    cof[1] = 76.18009173;
    cof[2] = -86.50532033;
    cof[3] = 24.01409822;
    cof[4] = -1.231739516;
    cof[5] = 0.00120858003;
    cof[6] = -0.00000536382;
    stp = 2.50662827465;
    half = 0.5;
    one = 1.0;
    fpf = 5.5;
    x = xx - one;
    tmp = x + fpf;
    tmp = (x + half) * log(tmp) - tmp;
    ser = one;
    for (j = 1;j<=6;j++)
	{
        x = x + one;
        ser = ser + cof[j] / x;
    }
    temp = (float)tmp + log(stp * ser);
	return temp;
}


void gser(float& gamser, float& a, float& x, float& gln)
{
	int itmax,n;
	float ap,sum,del,eps;
    itmax = 100;
    eps = 0.0000003;
    gln = gammln(a);
    if (x <= 0.0)
	{
        if (x < 0.0)
		{
            //cout<<" pause";
           // exit(1);
		}
        gamser = 0.0;
        exit(1);
    }
    ap = a;
    sum = 1.0 / a;
    del = sum;
    for( n = 1;n<=itmax;n++)
	{
        ap = ap + 1.0;
        del = del * x / ap;
        sum = sum + del;
        if (fabs(del) < (fabs(sum) * eps))
		{//	goto loop;
		gamser = sum * exp(-x + a * log(x) - gln);
		break;
		}
    }
}


void gcf(float& gammcf, float& a, float& x, float& gln)
{
	int itmax,n;
	float eps,a0,a1,b0,b1,fac,an,ana,anf,gold,g;
    itmax = 100;
    eps = 0.0000003;
    gln = gammln(a);
    gold = 0.0;
    a0 = 1.0;
    a1 = x;
    b0 = 0.0;
    b1 = 1.0;
    fac = 1.0;
    for (n = 1;n<=itmax;n++)
	{
        an = n;
        ana = an - a;
        a0 = (a1 + a0 * ana) * fac;
        b0 = (b1 + b0 * ana) * fac;
        anf = an * fac;
        a1 = x * a0 + anf * a1;
        b1 = x * b0 + anf * b1;
        if (a1 != 0.0)
		{
            fac = 1.0 / a1;
            g = b1 * fac;
            if (fabs((g - gold) / g) < eps )
			{
				gammcf = exp(-x + a * log(x) - gln) * g;
				break;
			}
            gold = g;
        }
    }
}



float gammq(float &a, float &x)
{
	float temp,gamser,gammcf,gln;
    if (x < 0.0 || a <= 0.0)
	{
        //cout<< "pause";
       // exit(1);
    }
    if (x < a + 1.0)
	{
        gser(gamser, a, x, gln);
        temp= 1.0 - gamser;
	}
    else
	{
        gcf(gammcf, a, x, gln);
        temp = gammcf;
	}
	return temp;   
}









//

//****************
//线性最小二乘法函数
//****************************

//usage:
//ndata: int,input,数据点的个数
//MA:int, input,参数的总个数
//mfit: int, input,欲求参数的总个数
//ncvm, int ,input,存储协放查方差矩阵的物理为数
//x[]: input ,ndata个元素的一维数组存放自变量x
//y[]:input,ndata个元素的一维数组存放自变量y
//SIG:inppt,ndata个元素的一维实型数组，输入全为1，表示同等权重
//lista:ma个一维整型数组，input,对参数重编号，使前面的MFIT个元素相应于要确定的参数，其余ma-mfit
	//个相应于已知参数
//a[]: ma个元素数组，某些为input,某些为output
//convar[]:二维实型叔祖，output
//chisq:output, to save the minimam tao	
void fit(float x[], float y[], int ndata, float sig[], int mwt, float& a,
		float& b, float& siga, float& sigb, float& chi2, float& q)
{
	int i;
    float sigdat,t,sxoss,wt,ss,sx = 0.0;
    float sy = 0.0;
    float st2 = 0.0;
    b = 0.0;
    if (mwt != 0)
	{
        ss = 0.0;
        for (i = 1; i<=ndata; i++)
		{
            wt = 1.0 / (sig[i] * sig[i]);
            ss = ss + wt;
            sx = sx + x[i] * wt;
            sy = sy + y[i] * wt;
        }
	}
    else
	{
        for (i = 1; i<=ndata; i++)
		{
            sx = sx + x[i];
            sy = sy + y[i];
        }
        ss = ndata;
    }
    sxoss = sx / ss;
    if (mwt != 0)
	{
        for (i = 1; i<=ndata; i++)
		{
            t = (x[i] - sxoss) / sig[i];
            st2 = st2 + t * t;
            b = b + t * y[i] / sig[i];
        }
	}
    else
	{
        for (i = 1; i<=ndata; i++)
		{
            t = x[i] - sxoss;
            st2 = st2 + t * t;
            b = b + t * y[i];
		}
    }
    b = b / st2;
    a = (sy - sx * b) / ss;
    siga = sqrt((1.0 + sx * sx / (ss * st2)) / ss);
    sigb = sqrt(1.0 / st2);
    chi2 = 0.0;
    if (mwt == 0)
	{
        for (i = 1; i<=ndata; i++)
		{
            chi2 = chi2 + pow((y[i] - a - b * x[i]) , 2);
        }
        q = 1.0;
        sigdat = sqrt(chi2 / (ndata - 2));
        siga = siga * sigdat;
        sigb = sigb * sigdat;
	}
    else
	{
        for (i = 1; i<=ndata; i++)
		{
            chi2 = chi2 + pow(((y[i] - a - b * x[i]) / sig[i]) , 2);
        }
       float tmpa,tmpx;
	   tmpa=0.5 * (ndata - 2);
	   tmpx=0.5 * chi2;
		q = gammq(tmpa, tmpx);
    }
}