/* density dependent birth-death simulation */

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

#define	DELTAT		1.0
#define	NPERT		((int)(1.0/DELTAT))
#define	MAXTIME 	(1000*NPERT)
#define	PRNTTIME	(10*NPERT)
#define	PNTS		(MAXTIME/PRNTTIME)
#define	ALPHA		0.2
#define	KK			20.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);

int main(void)
{
	int	irun, ttt, n1, n2, n3, cnter, seed, numruns;
	double	pops1[PNTS], pops2[PNTS], pops3[PNTS];

	alpha = ALPHA*DELTAT;
	beta = BETA*DELTAT;
	kk = KK;
	seed = 1456739853;
	srand48(seed);

	printf("Enter number of runs:\n");
	scanf("%d",&numruns);

	for(ttt=0;ttt<PNTS;ttt++)
		pops1[ttt] = pops2[ttt] = pops3[ttt] = 0.0;

	for(irun=0;irun<numruns;irun++)
	{
		cnter = 0;
		n1 = n2 = n3 = 5;
		for(ttt=0;ttt<MAXTIME;ttt++)
		{
			n1 += PopChange(n1,1);
			n2 += PopChange(n2,2);
			n3 += PopChange(n3,3);
			if(ttt%PRNTTIME==0)
			{
				pops1[cnter]   += n1;
				pops2[cnter]   += n2;
				pops3[cnter++] += n3;
			}
			if(n1==0 && n2==0 && n3==0) break;
		}
	}

	for(ttt=0;ttt<PNTS;ttt++)
		printf("%d %f %f %f\n",ttt*PRNTTIME,
				pops1[ttt]/(double)numruns,
				pops2[ttt]/(double)numruns,
				pops3[ttt]/(double)numruns);
	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);
}