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

/*
	solve an ODE

		d^2 x/dt^2 + damper dx/dt + spring x = u(t)
		x(0)     = 0
		dx/dt(0) = 0
	
	let x0 = x and x1 = dx/dt.
	the above can be converted as
	
	(d/dt) x0 = x1
	(d/dt) x1 = - damper x1 - spring x0 + u(t)
*/

#define	NumVars		2

#ifndef M_PI
#define	M_PI		3.14159265358979323846
#endif

static double pulse(double t)
{
	static double	period = 2.00;
	static double	duty = 0.50;
	static double	amplitude = 4.00;
	static double	offset = 0.00;
	
	if ( t - floor(t/period)*period <= period * duty ) {
		return amplitude + offset;
	} else {
		return offset;
	}
}


static void differential(int n, double x[], double t, double value[])
{
	static double	damper = 1.00;
	static double	spring = 4.00;
	
	value[0] = x[1];
	value[1] = - damper * x[1] - spring * x[0] + pulse(t);
}

main()
{
	int	kcnt;
	double	t, dt;
	double	x[NumVars];
	char	*filename = "twoorder.dat";
	FILE	*fd;
	
	if ((fd = fopen(filename,"w")) == NULL) {
		fprintf(stderr, "cannot open %s\n", filename);
		exit(1);
	}
	
	kcnt = 2000;
	dt = 0.01;
	/* kcnt * dt = 20 [sec] */

	/* initial values */
	t = 0.0;
	x[0] = 0.00;	/* position */
	x[1] = 0.00;	/* velocity */
	RungeKuttafprint(fd, NumVars, x, t);	// print initial values
	
	RungeKutta(fd, NumVars, x, differential, &t, dt, kcnt);
	
	fclose(fd);
}