/*
 * Example program for a runge kutta integration for two accelerating masses
 * with the diff eqn
 */

#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 x2
	X2[1] = 0;	// set initial condition for v2

        if( ( fp = fopen("outjason.txt", "w")) != NULL){
		fprintf(fp, "   t(s)      X1        V1      X2        V2 \n\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] = X1[1] * ((-Kd1 - Kd2)/M1) + X1[0] * ((-Ks1 - Ks2)/M1) + X2[1] * (Kd2/M1) + X2[0] * (Ks2/M1);
	dX2[0] = X2[1];
        dX2[1] = X2[1] * ((-Kd2 - Kd3)/M2) + X2[0] * ((-Ks2 - Ks3)/M2) + X1[1] * (Kd2/M2) + X1[0] * (Ks2/M2) + (F/M2);
}