/*
 *   Juan Velasquez - 26
 *   Numerical Analysis - cs407
 *   2006-04-09 - Spring 2006
 * 
 *   Analysis #A - To compute a set of functions that compute properties for the earth.
 */

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


#define MAXN 10

// Least Squares vars
double PressureCurveFitted1[5], PressureCurveFitted2[5];

//  uses polynomial coefficients generated on the fly
double Pressure(double alt);
//  uses a hard coded x array from prior calculation
double Pressure26(double alt);
double Temperature26(double alt);
double Celerity26(double temp);
double Gravity26(double alt);

double convertMilliBarsToKSM(double millibars);



// Lagrange Interpolation functions
struct pair {
	double x;
	double y;
};
int lagrange_interpolate( int n, struct pair table[], int k, double x, double *y, int debug );






int main()
{	
	double alt=0, temp=0;
	
	header();
	
	
	calculatePressurePolyValues();
	return 0;
	
	
	while(1)
	{
		printf("Altitude (double): ");
		scanf("%lf",&alt);
		if (alt<=0)
		{
			return 0;
		}
		
		printf("Pressure:   %lf\n", Pressure26(alt));
		printf("    Temp:   %lf\n", temp=Temperature26(alt));
		printf("Celerity:   %lf\n", Celerity26(temp));
		printf(" Gravity:   %lf\n", Gravity26(alt));
		
		printf("\n");
	}
	
	
	return 0;
}






double Pressure26(double alt)
{
	//double PressureCurveFit1[] = { 0.012217, -0.835692, 19.939858, -207.500562, 997.557278};
	//double PressureCurveFit2[] = { 0.000052, -0.010532, 0.801042, -28.366402, 422.738095};
	double PressureCurveFit1[] = { -9.6666667e-04, 7.8333333e-02, -2.3791667e+00, 3.4291667e+01, -2.5116667e+02};
	double PressureCurveFit2[] = { -2.0833333e-06, 5.2083333e-04, -5.1041667e-02, 2.4729167e+00, -6.1125000e+01};
	
	alt = alt/1000;
	
	if (alt < 25) {
		return convertMilliBarsToKSM( func(PressureCurveFit1, 5, alt) );
		//return func(PressureCurveFit1, 5, alt);
	} else if (alt < 65) {
		return convertMilliBarsToKSM( func(PressureCurveFit2, 5, alt) );
		//return func(PressureCurveFit2, 5, alt);
	} else if (alt >= 65) {
		return 0;
	} else {
		printf("Error: Pressure(%lf): Can't process invalid altitude.\n", alt);
	}
}

double Temperature26(double alt)
{
	int n=15, i=0;
	double y=0;
	struct pair table[15];
	double TemperatureAry[2][15] = { {0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0},
						{290.0, 240.0, 220.0, 225.0, 230.0, 245.0, 255.0, 265.0, 270.0, 275.0, 270.0, 255.0, 235.0, 215.0, 200.0} };
	
	
	// since my lagrange function used a different format for data, I have to convert
	for (i=0; i<n; i++)
	{
		table[i].x = TemperatureAry[0][i] * 1000.0;
		table[i].y = TemperatureAry[1][i];
	}
	
	lagrange_interpolate(n, table, 4, alt, &y, 0);
	
	return y;
}

double Celerity26(double temp, double alt)
{
	if (alt >= 65000.0)
	{
		return 1000.0;
	}
	else
	{
		return ( pow(1.4 * 8314.32 / 28.9644 * temp, 0.5) );
	}
}

double Gravity26(double alt)
{
	return (398622000000000.0/pow(6378178.0+alt,2.0));
}

double convertMilliBarsToKSM(double millibars)
{
	return (millibars * 1000.0 * 0.010197162);
}




double func(double x[10], int degree, double xVal)
{
	int i, power;
	double sum = 0.0;
	power = degree-1;
	for (i=0;i<degree;i++)
	{
		sum = sum + pow(xVal, power) * x[i];
		power--;
	}
	return sum;

}



int lagrange_interpolate( int n, struct pair table[], int k, double x, double *y, int debug )
{
	int i=0, j=0;
	int halfk=0, halfk_top=0;
	int start_point=0, points_less_x=0, points_greater_x=0;
	double l=0.0, ans=0.0, denom=0.0;
	
	
	// if k > n,  then fail.
	if (k>n) {
		return 0;
	}
	
	// if less than floor(k/2) data points to left or right
	halfk = k/2;  // using an int here is on purpose
	halfk_top = k - halfk;
	
	
	// count the number of data points above and below x.
	if (debug) { printf("lagrange_interpolate(): calculating location of x in table.\n"); }
	for (i=0; i<n; i++) {
		if (table[i].x < x) {
			points_less_x++;
		} else {
			points_greater_x++;
		}
	}
	if (debug>2) { printf("lagrange_interpolate(): points_less_x=%i\tpoints_greater_x=%i.\n", points_less_x, points_greater_x); }
	
	
	// figure out the starting data point for K points
	//  start_point is 0 indexed
	if (debug) { printf("lagrange_interpolate(): calculating starting data point.\n"); }
	if ( halfk > points_less_x ) {
		// then use all of them
		start_point = 0;
	} else if ( halfk_top > points_greater_x ) {
		// then use all of those
		start_point = n - k;
	} else {
		start_point = points_less_x - halfk;
	}
	
	
	if (debug) { printf("lagrange_interpolate(): starting interpolation.\n"); }
	for (i=0; i<k; i++) {
		
		// calc l
		l = 1;
		for (j=0; j<k; j++) {
			if ( i!=j ) {
				
				if (debug>2) { printf("lagrange_interpolate(): i=%i\tj=%i.\n", (start_point+i), (start_point+j)); }
				
				// double check demoinator is not zero
				//   even though this should never happen
				denom = table[start_point+i].x - table[start_point+j].x;
				if ( denom ) {
					l = l * (  (x - table[start_point+j].x)/denom );
					
					if (debug>3) { printf("lagrange_interpolate(): calculating l%i\tl=%lf.\n", (i+1), l); }
				} else {
					if (debug>2) { printf("lagrange_interpolate(): strange, demoninator in l%i is zero.\n", (i+1)); }
				}
			}
		}
		
		ans += table[start_point+i].y * l;
		
		if (debug>2) { printf("lagrange_interpolate(): using data point: y%i=%lf.\n", (start_point+i), table[start_point+i].y ); }
		if (debug>2) { printf("lagrange_interpolate(): calculating y\ty=%lf.\n", ans); }
	}
	
	*y = ans;	
	
	return 1;	
}