/* a least-squares fitter to a straight line.  given two arrays of
	   doubles, x and y, and the number of points in each array 'num',
	   it calculates a and b, such that
	                       y = a + bx
	   and returns them in the final two parameters.  the function returns
	   0 if all goes well, or -1 if there's a problem. 	
	*/

#include <stdio.h>
#include <math.h>

#define TINY  0.000001

static double get_t( /* int n */ );

int 
fitline(x, y, num, a, b)
double x[], y[], *a, *b;
int num;
{
	int i;
	double sumx, sumy, sumxy, sumxx, sumyy, meanx, meany;
	double bconf, syx, r;

#ifdef DEBUG
printf("FITLINE: num = %d    data follows\n", num);
for (i = 0; i < num; i++)
printf("%lf %lf \n", x[i], y[i]);
printf("----------------------\n");
#endif

	sumx = 0;
	sumy = 0;
	sumxy = 0;
	sumxx = 0;
	sumyy = 0;
	for (i = 0; i < num; i++) {
		sumx += x[i];
		sumy += y[i];
		sumxy += x[i]*y[i];	
		sumxx += x[i]*x[i];
		sumyy += y[i]*y[i];
	}
	meanx = sumx/num;
	meany = sumy/num;

	if ((num*sumxx - sumx*sumx) < TINY) {
#ifdef DEBUG
		fprintf(stderr, "fitline: indeterminate - zero denominator\n");
#endif
		return(-1);
	}

	*b = ((num*sumxy) - (sumx*sumy))/((num*sumxx - sumx*sumx));
	*a = meany - (*b)*meanx;

	return(0);

	/* go on to calculate some other things */
	/* standard error of the estimate */
	if (num > 2) {
		syx = sqrt((sumyy - ((*a)*sumy + (*b)*sumxy))/(num - 2.0));
	}
	else {
		syx = 0.0;
	}
	/* 95% confidence interval in the slope 'b' */
	if ((num > 2) && ((sumxx - ((sumx*sumx)/(double) num)) > TINY)) {
		bconf = get_t(num - 2)*syx/(sumxx - ((sumx*sumx)/(double) num));
	}
	else {
		bconf = 0.0;
	}
	/* correlation coefficient r */
	if (((sumxx - ((sumx*sumx)/(double) num)) < TINY) ||
	    ((sumyy - ((sumy*sumy)/(double) num)) < TINY)) {
		r = -99.0;
	}
	else {
		r = (sumxy - ((sumx*sumy)/(double) num));
		r /= sqrt((sumxx - ((sumx*sumx)/(double) num))*((sumyy - ((sumy*sumy)/(double) num))));
	}

	return (0);
}


	/* return the student t statistic with upper-tail area 0.025 (that will
	   give us a test for 95% confidence) for the given number of degrees
	   of freedom */

double tt[30] = { 12.706, 4.303, 3.182, 2.776, 2.571, 2.447, 2.365, 2.306,
                   2.262, 2.228, 2.201, 2.179, 2.160, 2.145, 2.131, 2.120,
                   2.110, 2.101, 2.093, 2.086, 2.080, 2.074, 2.069, 2.064,
                   2.060, 2.056, 2.052, 2.048, 2.045, 2.042 };
double tt40 = 2.021;
double tt60 = 2.000;
double tt120 = 1.980;
double ttinf = 1.960;

static double
get_t(n)
int n;
{
	if (n <= 30)
		return(tt[n - 1]);
	else if (n < 40) {
		return(tt[29] + 0.1*(tt40 - tt[29])*(n - 30.0));
	}
	else if (n == 40) {
		return(tt40);
	}
	else if (n < 60) {
		return(tt40 + (((tt60 - tt40)*(n - 40.0))/20.0));
	}
	else if (n == 60) {
		return(tt60);
	}
	else if (n < 120) {
		return(tt60 + (((tt120 - tt60)*(n - 60))/60.0));
	}
	else if (n == 120) {
		return(tt120);
	}
	else {
		return(ttinf);
	}
}