/*

	Prelab 3 
	This is a C program that solves the state equations for a 
	spring-mass-damper system

	By John Fountain

	EGR 345 

*/


#include <stdio.h>

void  step(double, double, double[]);
void  derivative(double, double[], double[]);

int main()
{

	FILE *fp;

	double h = 0.01;
	double t;

	double X[4];  /*create state variable list*/
	X[0] = 0;     /*set initial condition for x1*/
	X[1] = 0;     /*set initial condition for v1*/
	X[2] = 0;     /*set initial condition for x2*/
	X[3] = 0;     /*set initial condition for v2*/

	if( ( fp = fopen("out.txt", "w")) != NULL){
		fprintf(fp, "   t(s)      x1         v1            x2            v2        \n\n"); 
		for( t = 0.0; t < 10.0; t += h ){
			step(t, h, X);
			fprintf(fp, "%8.3f %8.3f %8.3f %8.3f %8.3f\n", t, X[0], X[1], X[2], X[3]);
		}
	}
	fclose(fp);
return 0;

}

/*First order integration done here*/

void step(double t, double h, double X[])
{
	double 	dX[4];
	derivative(t, X, dX);
	X[0] = X[0] + h * dX[0];
	X[1] = X[1] + h * dX[1];
	X[2] = X[2] + h * dX[2];
	X[3] = X[3] + h * dX[3];
}

/*State equations developed here*/

void derivative(double t, double X[], double dX[])
{

	double F = 100;
	double m1 = 300;
	double m2 = 200;
	double kd1 = 50;
	double kd2 = 60;
	double ks1 = 15;
	double ks2 = 20;
	double ks3 = 10;
	double kd3 = 70;

	if(t > 0 ){

	dX[0] = X[1];
	dX[1] = ((-ks1-ks2)/m1)*X[0] + ((-kd1-kd2)/m1)*X[1] + ((ks2)/m1)*X[2] + ((kd2)/m1)*X[3];
	dX[2] = X[3];
	dX[3] = (ks2/m2)*X[0] + (kd2/m2)*X[1] + ((-ks2-ks3)/m2)*X[2] + ((-kd2-kd3)/m2)*X[3] + F/m2;
	}
}