void medfit(double x[], double y[], int ndata, double& a, double& b, double& abdev)
{
	int j;
    double bb,del,chisq,f,sx = 0.0;
    double sigb,b1,b2,f1,f2,sy = 0.0;
    double sxy = 0.0;
    double sxx = 0.0;
    for (j = 1; j<=ndata; j++)
	{
        xt[j] = x[j];
        yt[j] = y[j];
        sx = sx + x[j];
        sy = sy + y[j];
        sxy = sxy + x[j] * y[j];
        sxx = sxx + x[j] * x[j];
    }
    ndatat = ndata;
    del = ndata * sxx - sx * sx;
    aa = (sxx * sy - sx * sxy) / del;
    bb = (ndata * sxy - sx * sy) / del;
    chisq = 0.0;
    for (j = 1; j<=ndata; j++)
	{
        chisq = chisq + pow((y[j] - (aa + bb * x[j])) , 2);
    }
    sigb = sqrt(chisq / del);
    b1 = bb;
    f1 = rofunc(b1);
    b2 = bb + fabs(3 * sigb) * sgn(f1);
    f2 = rofunc(b2);
    while (f1 * f2 > 0.0)
	{
        bb = 2.0 * b2 - b1;
        b1 = b2;
        f1 = f2;
        b2 = bb;
        f2 = rofunc(b2);
	}
    sigb = 0.01 * sigb;
    while (fabs(b2 - b1) > sigb)
	{
        bb = 0.5 * (b1 + b2);
        if (bb == b1 || bb == b2)
		{
			break;
		}
        f = rofunc(bb);
        if (f * f1 >= 0.0)
		{
            f1 = f;
            b1 = bb;
		}
        else
		{
            f2 = f;
            b2 = bb;
        }
    }
    a = aa;
    b = bb;
    abdev = abdevt / ndata;
}