/* James Cherry
 * EGR 345
 * September 11, 2007
 * Lab 2 - C Proggy to solve given state equation
 */
#include <stdio.h>
#include <math.h>
void multiply(double, double[],double[]);
void add(double[], double[], double[]);
void step(double, double, double[]);
void derivative(double, double[], double[]);

#define SIZE 4	//Global variable 

FILE *fp;

double  ks1=10,	//The coefficients for the variables were arbitraily
		ks2=10, //chosen and can be modified here in the list.  These are
		ks3=10, //placed in this fashion for when i made them #define for global
		kd1=15, //variables, the results went crazy
		kd2=15, 
		kd3=15, 
		F=35, 
		m1=15, 
		m2=15;
int main(void)
{
	double h=0.001;		//step height
	double t;		//time
	int j=0;		//loop variable
	double X[SIZE]; //create state variable list
	X[0]=0; 		// initial condition for x1 
	X[1]=0; 		// initial condition for v1
	X[2]=0; 		// initial condition for x2 
	X[3]=0; 		// initial condition for v2
	//Below is the main loop.  It creates the file and prints in the results.
	if((fp = fopen("EGR345_Lab2_b.txt", "w")) !=NULL)
	{
		fprintf(fp, "   t(s)	 	  x1		   v1		   x2		   v2\n\n");
		for( t = 0.0; t < 10.0; t += h) // Here we specify the total time 
		{
			step(t, h, X);
			if (j == 0) fprintf(fp, "%9.5f	%9.5f	%9.5f	%9.5f	%9.5f\n", 
					t, X[0], X[1], X[2], X[3]);
			j++; if(j >= 1/h) j=0; // 1/h prints only when 't' is a whole integer like 1, 2, etc.
		}
	}
	fclose(fp);
	printf("\n\nIt appears everything has worked so far!!\n\n");
	return 0;
}

void step(double t, double h, double X[])
{
	double dX[SIZE], F1[SIZE];

	derivative(t, X, dX);
	multiply(h, dX, F1);
	add(F1, X, X);
}
//Below are the derived state functions.  They are the same as used in SciLab.
void derivative(double t, double X[], double dX[])
{
	dX[0] = X[1];
	dX[1] = ((-ks1-ks2)/m1)*X[0] + ((-kd1-kd2)/m1)*X[1] + (ks1/m1)*X[2] 
				+ (kd2/m1)*X[3];
	dX[2] = X[3];
	dX[3] = ((-ks2-ks3)/m2)*X[2] + ((-kd2-kd3)/m2)*X[3] + (ks2/m2)*X[0] 
				+ (kd2/m2)*X[1] + F/m2;
}
void add(double X1[], double X2[], double R[])
{
	int i = 0;
	for(i = 0; i < SIZE; i++) 
		R[i] = X1[i] + X2[i];
}
void multiply(double X, double V[], double R[])
{
	int i = 0;
	for(i = 0; i < SIZE; i++)
		R[i] = X*V[i];
}