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

/*
	d^2 x0/dt^2 = -omega^2 (      2 x0 - x1)
	d^2 x1/dt^2 = -omega^2 (-x0 + 2 x1 - x2)
	d^2 x2/dt^2 = -omega^2 (-x1 + 2 x2 - x3)
	d^2 x3/dt^2 = -omega^2 (-x2 + 2 x3     )
	
	let y0=x0, y1=x1, y2=x2, y3=x3,
	    y4=(d/dt)x0, y5=(d/dt)x1, y6=(d/dt)x2, y7=(d/dt)x3 
	
	the above can be converted as
	
	(d/dt) y0 = y4
	(d/dt) y1 = y5
	(d/dt) y2 = y6
	(d/dt) y3 = y7
	(d/dt) y4 = -omega^2 (      2 y0 - y1)
	(d/dt) y5 = -omega^2 (-y0 + 2 y1 - y2)
	(d/dt) y6 = -omega^2 (-y1 + 2 y2 - y3)
	(d/dt) y7 = -omega^2 (-y2 + 2 y3)
*/

#define	NumVars		8

static double	omega, mass, cmp;

void diff(int n, double y[], double t, double value[])
{
	int	i;
	int	h;

#if DEBUG	
	printf("comp. diff.\n");
	for (i=0; i<n; i++) {
		printf("%5.2lf ", y[i]);
	}
	printf("\n");
#endif
	h = n/2;
	
	for (i=0; i<n/2; i++) {
		value[i] = y[h+i];
	}
	
	value[h+0] = -omega*omega*( 2*y[0] - y[1]);
	for (i=1; i<h-1; i++) {
		value[h+i] = -omega*omega*( -y[i-1] + 2*y[i] -y[i+1]);
	}
	value[h+(h-1)] = -omega*omega*( -y[h-2] + 2*y[h-1] );

#if DEBUG	
	for (i=0; i<n; i++) {
		printf("%5.2lf ", value[i]);
	}
	printf("\n--------------------\n");
#endif
}

main()
{
	int	i, kcnt;
	double	t, dt;
	double	y[NumVars];
	char	*filename = "rungekutta-t2.dat";
	FILE	*fd;
	
	if ((fd = fopen(filename,"w")) == NULL) {
		fprintf(stderr, "cannot open %s\n", filename);
		exit(1);
	}
	
	mass = 1.0;
	cmp = 4.0;
	omega = cmp/mass;
	
	kcnt = 500;
	dt = 0.01;
	
	/* initial values */
	t = 0.0;
	for (i=0; i<NumVars; i++) {
		y[i] = 0.00;
	}
	y[0] = -2.0;
/*	y[4] = 2.0; */
/*	y[0] = y[1] = y[2] = y[3] 2.0; */
	
	RungeKutta(fd, NumVars, y, diff, &t, dt, kcnt/5);
	RungeKutta(fd, NumVars, y, diff, &t, dt, kcnt/5);
	RungeKutta(fd, NumVars, y, diff, &t, dt, kcnt/5);
	RungeKutta(fd, NumVars, y, diff, &t, dt, kcnt/5);
	RungeKutta(fd, NumVars, y, diff, &t, dt, kcnt/5);
	
	fclose(fd);
}