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


double product(double x[], double bigx, int i, int upperbound,int lowerbound) ;
double summation(double x[], double y[], int numberofpairs, double bigx, int upperbound, int lowerbound);
double interpolate(int numberofpairs, double x[], double y[], double bigx, int k);
double pressure23(double altitude);
double temperature23(double altitude);
double gravity23(double alt);
double celerity23(double temperature);

double product(double x[], double bigx, int i, int upperbound,int lowerbound) 
{ 
     
	double result = 1;  
	int j = 0;
     
	for(j=lowerbound; j<=upperbound; j++) 
     	{ 
		if (j!=i)
		{ 
			result = (result) * ((bigx - x[j]) / (x[i]- x[j]));  
		} 
	} 
	return result; 
}




double summation(double x[], double y[], int numberofpairs, double bigx, int upperbound, int lowerbound)
{

	double sum = 0;
	int i=lowerbound;
	numberofpairs=upperbound;
	
	/*
	printf("inside summation\n");
	printf("bigx is %lf\n", bigx);
	printf("upperbound is is %d lowerbound is %d\n", upperbound, lowerbound);
	printf("npair is %d\n", numberofpairs);
	*/
	

	for(i=lowerbound; i<=upperbound; i++) 
	{ 
		sum+=(y[i] * product(x,bigx,i,upperbound,lowerbound)); 
	} 
	return sum; 
}





double interpolate(int numberofpairs, double x[], double y[], double bigx, int k)
{
	int i=0,j=0,width=0;
	int upperbound=0,lowerbound = 0;
	double fx=0;
	
/*
	printf("bigx is %lf\n", bigx);
	printf("k is %d\n", k);
	printf("npair is %d\n", numberofpairs);
*/	



	for (i=0;i<numberofpairs;i++) 
	{
		if (x[i]>bigx)
		{
			lowerbound=i-1;
			upperbound=i;
			break;
		}
	}
	width = (k/2) - 1;
	upperbound+=width;
	lowerbound-=width;

	if (lowerbound < 0) 
	{
		upperbound = (abs(lowerbound)+upperbound);
		lowerbound = 0; 	
	}
	else 
	{

		if(upperbound > numberofpairs)  
		{
			lowerbound=lowerbound-(abs(upperbound-numberofpairs));
			upperbound=numberofpairs;
		}
	}

	fx = summation(x,y,numberofpairs,bigx,upperbound,lowerbound);
	//printf("fx is %lf\n", fx);
	return fx;
}


double pressure23(double altitude)
{
	double pressure =0.0;
	
	double X[6];
	//X = (0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
	double results[6];	
	
	int i =0;
	double tempsum=0;

	if (altitude > 65)
	{
		pressure =0;
	}
	else	
	{
		if (altitude > 25)
		{

			X[0] = -0.000002;
			X[1] = 0.000521;
			X[2] = -0.051042;
			X[3] = 2.472917;
			X[4] = -61.125000;
			X[5] = 665.000000;

			
			
			
		}
		else
		{
	
			X[0] = -0.000967;
			X[1] = 0.078333;
			X[2] = -2.379167;
			X[3] = 34.291670;
			X[4] = -251.166700;
			X[5] = 1000.000000;



		}
		
		
		results[0] = (X[0] * (pow ( altitude, 5)) );
		results[1] = (X[1] * (pow ( altitude, 4)) );
		results[2] = (X[2] * (pow ( altitude, 3)) );
		results[3] = (X[3] * (pow ( altitude, 2)) );
		results[4] = (X[4] * (pow ( altitude, 1)) );
		results[5] = (X[5]);
		
		/*
		printf("%lf * %lf ^ 5  = %lf\n", X[0], altitude, results[0]);
		printf("%lf * %lf ^ 4  = %lf\n", X[1], altitude, results[1]);
		printf("%lf * %lf ^ 3  = %lf\n", X[2], altitude, results[2]);
		printf("%lf * %lf ^ 2  = %lf\n", X[3], altitude, results[3]);
		printf("%lf * %lf ^ 1  = %lf\n", X[4], altitude, results[4]);
		printf("%lf * %lf ^ 0  = %lf\n", X[5], altitude, results[5]);
		*/

		pressure = results[0] + results[1] + results[2] + results[3] + results[4] + results[5];   



		/*
		pressure =(X[0] * (pow ( altitude, 5)) ) 
			 + (X[1] * (pow ( altitude, 4)) )
			 + (X[2] * (pow ( altitude, 3)) )
			 + (X[3] * (pow ( altitude, 2)) )
			 + (X[4] * (pow ( altitude, 1)) )
			 + (X[5]);
		*/

	}
	
	return pressure * 1000 * .010197162;
}


double temperature23(double altitude)
{

	double X[20];
	double Y[20];
	double temp = 0;


	X[0]=0.0, Y[0]=290.0;      
	X[1]=5000.0, Y[1]=240.0;      
	X[2]=10000.0, Y[2]=220.0;      
	X[3]=15000.0, Y[3]=225.0;      
	X[4]=20000.0, Y[4]=230.0;      
	X[5]=25000.0, Y[5]=245.0;      
	X[6]=30000.0, Y[6]=255.0;      
	X[7]=35000.0, Y[7]=265.0;      
	X[8]=40000.0, Y[8]=270.0;      
	X[9]=45000.0, Y[9]=275.0;      
	X[10]=50000.0, Y[10]=270.0;    
	X[11]=55000.0, Y[11]=255.0;  
	X[12]=60000.0, Y[12]=235.0; 
	X[13]=65000.0, Y[13]=215.0;      
	X[14]=70000.0, Y[14]=200.0;

	
	temp= interpolate(15,X,Y,altitude,4);

	return temp;


}


double gravity23(double alt)
{
	return (3.98622e+14) / (pow(((6.378178e+06) + (alt)),2));
}

double celerity23(double temperature)
{
	return sqrt(((1.4)*((8.31432e+03) / (28.9644)) * temperature));
}

int main()
{
	double pressureval=0.0;
	double altval =0.0;
	double altvalkm=0.0;
	double tempval =0.0;
	double celery =0.0;
	double gravity =0.0;

	printf ("\n Devin Scott NN: 23 Numerical Analysis. Analysis A: reprise\n");  //header
	do
	{
		printf("enter value of altitude in meters");
		scanf("%lf",&altval);
		altvalkm = altval / 1000; //convert to KM
	
		pressureval=pressure23(altvalkm);
		printf("the pressure is %lf kg/m^2\n",pressureval);
	
		tempval = temperature23(altval);
		printf("the temperature is %lf degrees Kelvin\n",tempval);

		celery = celerity23(tempval);
		printf("the celerity is %lf m/s\n",celery); 

		gravity = gravity23(altval);
		printf("the gravity is %lf m/s^2\n",gravity); 
	
	}
	while(altval != 0);


	return 0;

}