/* density dependent birth-death simulation (made avedatF.txt)*/

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

#define	DELTAT		1.0
#define	NPERT		((int)(1.0/DELTAT))
#define	NUMRUNS		50
#define	SKIPTIME 	(100)
#define	AVETIME 	2000
#define	MEASSTEPS	20
#define	KK		40.0
#define	BETA		0.1

/*
Model 1 ---   dN/dt = alpha N (1 - N / K) - beta N
Model 2 ---   dN/dt = alpha N (1 - N / 2K) - beta N (1 + alpha N / 2K beta)
Model 3 ---   dN/dt = alpha N - beta N (1 + alpha N / K beta)
*/

double	alpha, kk, beta;
int	PopChange(int, int);
void Update_Ave(double, double *, double *, double);

int main(void)
{
	int	irun, ttt, n1, n2, n3;
	long	seed;

	double	mean1, var1, pnts1;
	double	mean2, var2, pnts2;
	double	mean3, var3, pnts3;
	double	a1, a2, b1, b2;
	double	mean, rr, ss;

	beta = BETA*DELTAT;
	kk = KK;

	seed = 1456739853;
	srand48(seed);

	for(alpha=0.15;alpha<=0.35;alpha+=0.01)
	{
		mean1 = var1 = pnts1 = 0.0;
		mean2 = var2 = pnts2 = 0.0;
		mean3 = var3 = pnts3 = 0.0;

		for(irun=0;irun<NUMRUNS;irun++)
		{
			n1 = n2 = n3 = (int)((alpha-beta)*kk/alpha) + 1;
			for(ttt=0;ttt<SKIPTIME+AVETIME;ttt++)
			{
				n1 += PopChange(n1,1);
				n2 += PopChange(n2,2);
				n3 += PopChange(n3,3);

				if(ttt%MEASSTEPS==0 && ttt>SKIPTIME)
				{
					if(n1>0) Update_Ave((double)n1, &mean1, &var1, ++pnts1);
					if(n2>0) Update_Ave((double)n2, &mean2, &var2, ++pnts2);
					if(n3>0) Update_Ave((double)n3, &mean3, &var3, ++pnts3);
				}
			}

		}
		if(pnts1>0.0) var1 /= pnts1;
		if(pnts2>0.0) var2 /= pnts2;
		if(pnts3>0.0) var3 /= pnts3;

		mean = (alpha-beta)*kk/alpha;	/* deterministic mean */

/* print simulation measurements */

		printf("%f %f %f %f %f %f %f %f  ",alpha,
			mean, mean-mean1, mean-mean2, mean-mean3,
			var1, var2, var3);

/* calculate and print stochastic model predictions */

		a1 = alpha; a2 = beta;
		rr = a1-a2; ss = alpha/kk;

		b1 = alpha/kk;		/* Model 1 */
		b2 = 0.0;
		var1 = (mean/2.0/ss)*((a1+a2)/mean + (b2-b1));
		mean1 = mean - var1/mean;

		b1 = alpha/2.0/kk;
		b2 = alpha/2.0/kk;
		var2 = (mean/2.0/ss)*((a1+a2)/mean + (b2-b1));
		mean2 = mean - var2/mean;

		b1 = 0.0;
		b2 = alpha/kk;
		var3 = (mean/2.0/ss)*((a1+a2)/mean + (b2-b1));
		mean3 = mean - var3/mean;

		printf(" %f %f %f %f %f %f\n",
			mean-mean1, mean-mean2, mean-mean3,
			var1, var2, var3);

	}
	return (0);
}


int PopChange(int pop, int modnum)
{
	int temp=0, i;
	double xpop;

	xpop = (double)pop;
	if(modnum==1)
	{
		for(i=0;i<pop;i++)
			if(drand48()<alpha*(1.0-xpop/kk)) temp++;
		for(i=0;i<pop;i++)
			if(drand48()<beta) temp--;
	}

	if(modnum==2)
	{
		for(i=0;i<pop;i++)
			if(drand48()<alpha*(1.0-xpop/kk/2.0)) temp++;
		for(i=0;i<pop;i++)
			if(drand48()<beta*(1.0+alpha*xpop/kk/beta/2.0)) temp--;
	}

	if(modnum==3)
	{
		for(i=0;i<pop;i++)
			if(drand48()<alpha) temp++;
		for(i=0;i<pop;i++)
			if(drand48()<beta*(1.0+alpha*xpop/kk/beta)) temp--;
	}
	return(temp);
}


void Update_Ave(double datum, double *mean, double *sum_sq, double n)
{
	double	deviation;

	deviation = datum - *mean;
	*mean += (deviation/n);
	*sum_sq += (deviation*(datum-*mean));
}