// Irma Montigny		
// Martin Przekadzinski
// Chris Hough
// Final Progam
// Program 4 - Write a program in C to solve the state equations of the diagram shown.
// September 16, 2003

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


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


int main(){
	FILE	*fp;

	double	h = 0.1;
	double	t;

	double	X1[2];	// create state variable list
	double	X2[2];	// create state variable list
	X1[0] = 0;	// set initial condition for x1
        X1[1] = 0;      // set initial condition for v1
	X2[0] = 0;	// set initial condition for x1
        X2[1] = 0;      // set initial condition for v1

	if( ( fp = fopen("out.txt", "w")) != NULL){
		fprintf(fp, "   t(s)		x1		v1		x2		v2 \n\n");
                fprintf(fp, "__________________________________________________________________________\n");
                for( t = 0.0; t <= 10.0; t += h ){
                        fprintf(fp, "%8.3lf        %8.3lf      %8.3lf      %8.3lf       %8.3lf\n", t, X1[0], X1[1], X2[0], X2[1]);
			step(h, X1, X2);

        	}
	}
	fclose(fp);
	return 0;
}

//
// First order integration done here (could be replaced with runge kutta)
//
void step(double h, double X1[], double X2[]){
	double 	dX1[2];
	double	dX2[2];

	derivative(X1, dX1, X2, dX2);
	
	X1[0] = X1[0] + h * dX1[0];
	X1[1] = X1[1] + h * dX1[1];
	
	X2[0] = X2[0] + h * dX2[0];
	X2[1] = X2[1] + h * dX2[1];

}

//
// State Equations Calculated Here
//
void derivative(double X1[], double dX1[], double X2[], double dX2[]){
        double  F = 1;
	double	M1 = 1;
	double	M2 = 1;
        double  kd1 = 1;
        double  kd2 = 1;
        double  kd3 = 1;
        double  ks1 = 1;
        double  ks2 = 1;
        double  ks3 = 1;

	dX1[0] = X1[1];
	dX1[1] = ((kd1+kd2)/-M1)*X1[1]+((ks1+ks2)/-M1)*X1[0]+(kd2/M1)*X2[1]+(ks2/M1)*X2[0];
	dX2[0] = X2[1];
	dX2[1] = ((kd2+kd3)/-M2)*X2[1]+((ks2+ks3)/-M2)*X2[0]+(kd2/M2)*X1[1]+(ks2/M2)*X1[0]+(F/M2);
}