/* Mackey-Glass time series */
/* Jyh-Shing Roger Jang, 1992 */

#include <stdio.h> 
#include <strings.h> 
#include <math.h> 

main(argc, argv)
int argc;
char **argv;
{
	double runge();
	double x, x_next;
	int tau;         /* delay constant */
	double x_before; /* x(t-tau) */
	double step_size;
	double time, init_time = 0;
	int sample_n;
	int i, l, index = 0;
	double *x_history; /* array to store delay information */
	int history_length;
	int out_data_n; /* number of output data entries is 
			   (out_data_n+1), one of which is initial
			   contidtion */

	if (argc != 6) {
		printf("Usage: %s step_size sample_number init_x tau output_data_number\n", argv[0]);
		exit(1);
	}
	step_size = atof(argv[1]);
	sample_n = atoi(argv[2]);
	x = atof(argv[3]);
	tau = atoi(argv[4]);
	out_data_n = atoi(argv[5]);

	history_length = (int)(tau/step_size);
	if (history_length != 0) {
		x_history = (double *)calloc(history_length, sizeof(double));
		for (i = 0; i < history_length; i++)
			x_history[i] = 0;
	}

	l = (int)(sample_n/out_data_n);

	time = init_time;
	for (i = 0; i <= sample_n; i++) {
		if (i%l == 0)
			printf("%g %f\n", i*step_size, x);
		x_before = tau == 0 ? x:x_history[index];
		x_next = runge(x, time, step_size, x_before);
		if (tau != 0) {
			x_history[index] = x_next;
			index = (index + 1)%history_length; 
		}
		time += step_size;
		x = x_next;
	}
}
 
/* 4-th order runge-kutta integration formula for one variable */
double 
runge(x_now, time, step, x_before)
double x_now, time, step, x_before;
{
	double equation();
	double a, b, c, d;

	a = step*equation(x_now,         time,          x_before);
	b = step*equation(x_now + 0.5*a, time+step/2.0, x_before);
	c = step*equation(x_now + 0.5*b, time+step/2.0, x_before);
	d = step*equation(x_now + c,     time+step,     x_before);
	return(x_now + a/6 + b/3 + c/3 + d/6);
}

/* return dx/dt */
double
equation(x, time, x_before)
double x;
double time;
double x_before;
{
	double a = 0.2;
	double b = 0.1;
	double x_dot;
	x_dot = -b*x + a*x_before/(1 + pow(x_before, 10.0));
	return(x_dot);
}
