/*  Fourth order Runge-Kutta numerical integration */

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

/* general numerical integration parameters */

#define NUMVARS		2
#define MAXTIME		40.0
#define	OUTTIMEINC	0.1
#define STEPSIZE	0.01


/* specific model parameters */

#define	INITHERM	0.2
#define	INITFEM		0.01
#define	RH		40.0
#define	RF		50.0
#define	H0		3.5

double	state1[NUMVARS], state2[NUMVARS];
double	stepsize = STEPSIZE;
void OneStep(double, double *, double *, double);
void functions(double, double *, double *, double);

int main(void)
{
	int	flag=1;
	double	time = 0.0, printtime = 0.0;
	double	*oldstate, *newstate, *temp;

	oldstate = state1; newstate = state2;

	newstate[0] = INITHERM; /* model dependent initial conditions */
	newstate[1] = INITFEM;

	while(time<=MAXTIME)
	{						/* switch pointers */
		temp = oldstate; oldstate = newstate; newstate = temp;

		OneStep(time, oldstate, newstate, stepsize);	/* take a step */
		time += stepsize;

		if(time>printtime)
		{				/* printout every OUTTIMEINC time units */
			printtime += OUTTIMEINC;
			printf( " %3.6f  %3.6f  %3.6f\n", time, newstate[0], newstate[1]);
		}
	}
	return (0);
}

void OneStep(double time, double *in, double *out, double step) /* Runge-Kutte numerical integration method. */
{
	int 	i;
	double	v1[NUMVARS], v2[NUMVARS], v3[NUMVARS], v4[NUMVARS];
	double	vals[NUMVARS];

	for(i=0;i<NUMVARS;i++) vals[i] = in[i];
	functions(time, vals, v1, step);

	for(i=0;i<NUMVARS;i++) vals[i] = in[i] + v1[i]/2.0;
	functions(time+step/2.0, vals, v2, step);

	for(i=0;i<NUMVARS;i++) vals[i] = in[i] + v2[i]/2.0;
	functions(time+step/2.0, vals, v3, step);

	for(i=0;i<NUMVARS;i++) vals[i] = in[i] + v3[i];
	functions(time+step, vals, v4, step);

	for(i=0;i<NUMVARS;i++)
		out[i] = in[i] + (v1[i] + 2.0*v2[i] + 2.0*v3[i] + v4[i])/6.0;
}


void functions(double time, double *state, double *derivs, double step)  /* model-specific equations */
{
	int i;
	double hh, ff;

	hh = state[0]; ff = state[1];		/* modify derivs[] for specific problem */

	derivs[0] = RH*(1.0-hh-ff)*hh*hh/(H0+hh) - hh;
	derivs[1] = RF*(1.0-hh-ff)*ff*hh/(H0+hh) - ff;

	for(i=0;i<NUMVARS;i++) derivs[i] *= step;
}