/*  Fourth order Runge-Kutta numerical integration */
/*
gcc LVadapt.c -o xlv -lm

xmgrace -hardcopy -printfile many.ps -param many.par -nxy pops1.dat

*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>

/* model parameters and stuff */

#define NUMSPECS	(20)

#define NUMVARS		(NUMSPECS)
#define MAXTIME		100.0
#define	OUTTIMEINC	1.0
#define STEPSIZE	0.001

#define MAX_ERROR	1.0e-10
#define MIN_ERROR	(MAX_ERROR/20.0)   	
#define SCALER		1.2			
#define MIN_STEP	0.000001
#define MAX_STEP	0.1	
#define	ADAPT_TESTS	5

#define RRR			((double)1.0)
#define	KAVE		((double)1.0)
#define	KWID		((double)0.1)
#define ALPHAWIDTH	((double)0.25)
#define	INITDENS	(KAVE/(double)NUMSPECS)
#define ZERO_DENS	(0.0)

#define RRR			((double)1.0)
#define MATSENS		((double)0.001)


double	state1[NUMVARS], state2[NUMVARS];
double	rrr[NUMVARS], alphas[NUMVARS][NUMVARS];
int		kkkk[NUMVARS], atesting=0;
double	stepsize = STEPSIZE;
double	exp(double), pow(double,double);
double	zero_dens;

int main(int argc, char **argv)
{
	int		i, seed;
	double	time = 0.0, printtime = 0.2, nbar;
	double	*oldstate, *newstate, *temp;
	FILE	*fopen(), *fileid, *fileid2, *fileid3, *fileid4, *fileid5;
	char fname[30];


	InitSeed(&seed);
	srand48(seed);
	if(argc>1)
	{
 		sscanf(argv[1],"%lf",&zero_dens);
	}
	
	InitRAlphs();
	InitKs();
	InitDens(state1);
	fileid3 = fopen("matrix.dat","w");
	OutputMatrix(state1, fileid3);
	fflush(fileid3);
		
	fileid = fopen("pops1.dat","w");
	time = printtime = 0.0;

	oldstate = state1; newstate = state2;
	while(time<=MAXTIME)
	{
		if(time>=printtime)
		{
			printtime += OUTTIMEINC;
			printf( " %3.3f (%7.5f)", time, stepsize);
			fprintf(fileid," %f ", time);
			nbar=0.0;
			for(i=0;i<NUMSPECS;i++)
			{	
				nbar += oldstate[i];
				if(i<10) printf( " %3.3f ", oldstate[i]);
				fprintf(fileid," %f ", oldstate[i]);
			}
			fprintf(fileid,"\n");
			printf("------ %3.3f  \n", nbar);
	   }
		
/* take a step */
		OneAdaptiveStep(&time, oldstate, newstate);
/* switch pointers */
		temp = oldstate; oldstate = newstate; newstate = temp;

	}
	fclose(fileid);
	OutputMatrix(oldstate, fileid3);
}

InitSeed(newseed)
long unsigned *newseed;
{
	int	j;
	struct	timeval	tp;
	struct	timezone tzp;

	j = gettimeofday(&tp,&tzp);
	*newseed = (long)fabs((double)(tp.tv_usec * tp.tv_sec));
}

InitRAlphs()
{
	int i, j;
	double xxx;
	
	for(i=0;i<NUMSPECS;i++)
	{
		rrr[i] = RRR;	
		for(j=0;j<NUMSPECS;j++)
		{
			alphas[i][j] = (0.5-ALPHAWIDTH/2.0) + 
							drand48()*(ALPHAWIDTH);
				
			if(drand48()<zero_dens)
				alphas[i][j] = (double)0.0;
		   if(i==j) alphas[i][j] = (double)1.0;
	   	}
	}
}


InitDens(state)
double *state;
{
	int i, j;
	
	for(i=0;i<NUMSPECS;i++)
	{
		state[i] = drand48()*INITDENS;
	}
}

InitKs()
{
	int i, j;
	
	for(i=0;i<NUMSPECS;i++)
	{
		kkkk[i] = KAVE-KWID/2.0+(drand48()*KWID);
	}
}

OutputMatrix(in, fileid)
double *in;
FILE *fileid;
{
	int i, j, oks[NUMSPECS];
	
	fprintf(fileid,"        ");
	for(i=0;i<NUMSPECS;i++)
	{
		oks[i] = 0;
		if(in[i]>MATSENS) oks[i] = 1;
		if(oks[i]==1) fprintf(fileid,"%-5d  ",i);
	}
	fprintf(fileid,"\n");
	
	
	for(i=0;i<NUMSPECS;i++)
	{
		if(oks[i]==1)
		{
			fprintf(fileid,"%-5d",i);
			
			for(j=0;j<NUMSPECS;j++)
			{
				if(oks[j]==1)
				{
					fprintf(fileid," %6.3f",alphas[i][j]);
				}
			}
			fprintf(fileid,"\n");
	   	}
	}
	fprintf(fileid,"\n\n\n");

}


OneAdaptiveStep(double *time, double *in, double *out) 
{
	int 	i;
	double	v1[NUMVARS], v2[NUMVARS];
	double	vals[NUMVARS], diff=0.0, xxx;
	
	if(atesting==0)
	{
		OneStep(*time, in, v1, stepsize/2.0); 	/* half steps */
		OneStep(*time+stepsize/2.0, v1, v2, stepsize/2.0); 	
	}
	OneStep(*time, in, v1, stepsize); 	/* full step */
	*time += stepsize;

	
	if(atesting==0)
	{
		for(i=0;i<NUMVARS;i++)
		{
			xxx = fabs((v1[i]-v2[i])/v2[i]);
			if(xxx>diff && in[i]>MAX_ERROR) diff = xxx;
		}
		if(diff>MAX_ERROR && stepsize>MIN_STEP) 
			stepsize /= SCALER;	/* reduce */
		else if(diff<MIN_ERROR && stepsize<MAX_STEP) 
			stepsize *= SCALER;	/* increase */
		atesting = ADAPT_TESTS+1;
	}
	atesting--;
	
	for(i=0;i<NUMVARS;i++)
		out[i] = v1[i];
}



OneStep(time, in, out, step) /* Runge-Kutte numerical integration method. */
double time, *in, *out, step;
{
	int 	i;
	double	v1[NUMVARS], v2[NUMVARS];
	double	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;
}


functions(time, state, derivs, step)
double time, *state, *derivs, step;
{
	int	i, j;
	double chis[NUMSPECS];
	
	for(i=0;i<NUMSPECS;i++)
	{
		chis[i] = 0.0;
		for(j=0;j<NUMSPECS;j++) 
			chis[i] += alphas[i][j]*state[j];

		derivs[i] = rrr[i]*(kkkk[i]-chis[i])*state[i];
	}
	
	for(i=0;i<NUMVARS;i++) derivs[i] *= step;
}