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

/*
	Runge-Kutta method
	
	solve an ordinary differential equation dx/dt = f(x,t)
	where x is an n-dimensional vector consisting of n variables.

	you can find how to use this function in
		rungekutta-t0.c
		rungekutta-t1.c
		rungekutta-t2.c
*/

/*
	solve an ordinary differential equation dx/dt = f(x,t)
	from the k-th time step to the (k+1)-th time step
	
	input:
		n	dimension of vector x
		x[]	value of vector x at the k-th time step
		diff()	procedure to compute f(x,t)
		*pt	time at the k-th time step
		dt	time step (interval)
		
		procedure diff() should be in the following form:

			diff(int n, double x[], double t, double v[])
				n	dimension of vector x
				x[]	vector (n-dimensional array)
				t	time
				v[]	array where value of f(x,t) is stored
					(n-dimensional array)

	output:
		x[]	value of vector x at the (k+1)-th time step
		*pt	time at the (k+1)-th time step
*/

extern void RungeKuttaOneStep(int n, double x[], void diff(int, double *, double, double *), double *pt, double dt)
{
	int	i;	
	double	x1[MaxVars], x2[MaxVars], x3[MaxVars];
	double	d1[MaxVars], d2[MaxVars], d3[MaxVars], d4[MaxVars];
	
	diff(n, x, *pt, d1);				/* comp. d1 */
	for(i=0;i<n;i++) {
		x1[i] = x[i] + (dt/2.00)*d1[i];
	}
	diff(n, x1, (*pt)+(dt/2.00), d2);		/* comp. d2 */
	for(i=0;i<n;i++) {
		x2[i] = x[i] + (dt/2.00)*d2[i];
	}
	diff(n, x2, (*pt)+(dt/2.00), d3);		/* comp. d3 */
	for(i=0;i<n;i++) {
		x3[i] = x[i] + dt*d3[i];
	}
	diff(n, x3, (*pt)+dt, d4);			/* comp. d4 */
#if DEBUG
	printf("--- d1 ---\n");
	for (i=0; i<n; i++) printf("%5.2lf ", d1[i]);
	printf("\n");
	printf("--- d2 ---\n");
	for (i=0; i<n; i++) printf("%5.2lf ", d2[i]);
	printf("\n");
	printf("--- d3 ---\n");
	for (i=0; i<n; i++) printf("%5.2lf ", d3[i]);
	printf("\n");
	printf("--- d4 ---\n");
	for (i=0; i<n; i++) printf("%5.2lf ", d4[i]);
	printf("\n");
#endif
	/* update x[] */
	for(i=0;i<n;i++) {
		x[i] += (dt/6.00)*(d1[i] + 2*d2[i] + 2*d3[i] + d4[i]);
	}
	/* update *t */
	*pt += dt;
}

/*
	solve an ordinary differential equation dx/dt = f(x,t)
	computation results are stored in a file

	input:
		fd	file discripter where results are stored
		n	dimension of vector x
		x[]	initial value of vector x
		diff()	procedure to compute f(x,t)
		&t	pointer to time
		dt	time step
		kcnt	number of iteration
		
		procedure diff() should be in the following form:

			diff(int n, double x[], double t, double v[])
				n	dimension of vector x
				x[]	vector (n-dimensional array)
				t	time
				v[]	array where value of f(x,t) is stored
					(n-dimensional array)

	output:
		the following values are printed to FILE *fd:
		time x[0] x[1] ... x[n-1]
*/

static char format[] = "%lf\t";
extern void RungeKuttafprint(FILE *fd, int n, double x[], double t)
{
	int	i;
	
	fprintf(fd, format, t);
	for(i=0;i<n;i++) {
		fprintf(fd, format, x[i]);
	}
	fprintf(fd, "\n");
}

extern void RungeKutta(FILE *fd, int n, double x[], void diff(int, double *, double, double *), double *pt, double dt, int kcnt)
{
	int	k, i;
	//double	t;
	
	if (n > MaxVars) {
		fprintf(stderr, "%d exceeds MaxVars (%d)\n", n, MaxVars);
		fprintf(stderr, "increase MaxVars in opt.h and recompile all\n");
		exit(1);
	}
	
	for(k=1; k<=kcnt; k++) {
#if DEBUG
		printf("<<<<< count = %d >>>>>\n", k);
#endif
		RungeKuttaOneStep(n, x, diff, pt, dt);
		RungeKuttafprint(fd, n, x, *pt);
		//fprintf(fd, format, *pt);
		//for(i=0;i<n;i++) {
		//	fprintf(fd, format, x[i]);
		//}
		//fprintf(fd, "\n");
	}
}