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

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

#define	NumVars		2

#ifndef M_PI
#define	M_PI		3.14159265358979323846
#endif

static double	omega;

static void differential(int n, double x[], double t, double value[])
{
	value[0] = x[1];
	value[1] = -(omega*omega) * x[0];
}

main()
{
	int	kcnt;
	double	t, dt;
	double	x[NumVars];
	char	*filename = "rungekutta-t1.dat";
	FILE	*fd;
	
	if ((fd = fopen(filename,"w")) == NULL) {
		fprintf(stderr, "cannot open %s\n", filename);
		exit(1);
	}
	
	omega = 2 * M_PI / 10.0;
	
	kcnt = 5000;
	dt = 0.1;
	
	/* initial values */
	t = 0.0;
	x[0] = 10.0;
	x[1] = 0.0;
	RungeKuttafprint(fd, NumVars, x, t);	// print initial values
	
	RungeKutta(fd, NumVars, x, differential, &t, dt, kcnt);
	RungeKutta(fd, NumVars, x, differential, &t, dt, kcnt);
	
	fclose(fd);
}