/* Two-species competition simulation */

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

#define	MAXTIME		100
#define	OUTTIME		5
#define	PSTIME		1
#define	CELLS		1000
#define	PSSKIP		2

#define	ALPHA1	0.25
#define	ALPHA2	0.25
#define	BETA12	0.075 /* b12+b21<1 */
#define	BETA21	0.1
#define	MU1	0.05
#define	MU2	0.05
#define	MIXED	1

int	state[CELLS], dist1, dist2;
double	n1cnt, n2cnt, alpha1, alpha2, beta12, beta21, mu1, mu2;

void InitStuff(void);
void Competition(void);
void Reproduction(void);
void NewBabe(int, int, int);
void Mortality(void);
#if(MIXED)
void Mix(void);
#endif
void Measurements(void);
void Spit_EPS_Header(FILE *, int, int);
void Spit_EPS_Row(FILE *, int);
void Spit_EPS_Trailer(FILE *);

int main(void)
{
	int i;
	FILE	*fileid1, *fileid2;

	fileid1 = fopen("compete1.eps","w");
	fileid2 = fopen("compete2.eps","w");
	InitStuff();

	printf("Enter dispersal distances dist1 and dist2:\n");
	scanf("%d %d",&dist1,&dist2);

	Spit_EPS_Header(fileid1,CELLS/PSSKIP,MAXTIME/PSTIME);
	Spit_EPS_Header(fileid2,CELLS/PSSKIP,MAXTIME/PSTIME);

	for(i=0;i<MAXTIME;i++)
	{
		Measurements();
		if(i%PSTIME==0)
		{
			Spit_EPS_Row(fileid1,1);
			Spit_EPS_Row(fileid2,2);
		}
		if(i%OUTTIME==0) printf("%d %f %f\n", i, n1cnt, n2cnt);

		Mortality();
		Competition();
		Reproduction();
#if(MIXED)
		Mix();
#endif
		Measurements();
	}
	Spit_EPS_Trailer(fileid1);
	Spit_EPS_Trailer(fileid2);

	return (0);
}

void InitStuff(void)
{
	int i;
	int seed=145739853;
	double xxx;

	srand48(seed);
	alpha1 = ALPHA1; alpha2 = ALPHA2;
	beta12 = BETA12; beta21 = BETA21;
	mu1 = MU1; mu2 = MU2;

	for(i=0;i<CELLS;i++)
	{
		xxx = drand48();
		if(xxx<0.33)      state[i] = 1;
		else if(xxx<0.66) state[i] = 2;
		else              state[i] = 0;
	}
}

void Competition(void)
{
	int i, j, k, pairs[CELLS], npair, totpairs;
	double xxx;

	totpairs = 0;
	for(i=0;i<CELLS-1;i++)
		if(state[i]*state[i+1]==2)
			pairs[totpairs++] = i;

	k = totpairs;
	for(i=0;i<totpairs;i++)
	{
		j = k*drand48();
		npair = pairs[j];
		if(state[npair]*state[npair+1]==2) /* states might've changed */
		{
			xxx = drand48();
			if(xxx<beta12)
			{
				if(state[npair]==1) state[npair] = 0;
				else state[npair+1] = 0;
			}
			else if(xxx<beta12+beta21)
			{
				if(state[npair]==2) state[npair] = 0;
				else state[npair+1] = 0;
			}
		}
		pairs[j] = pairs[--k];
	}
}



void Reproduction(void)
{
	int i, j, k, sites[CELLS], nsite, totpop=0;

	for(i=0;i<CELLS;i++) /* determine reproducing individuals */
	{
		if(state[i]==1)
		{
			if(drand48()<alpha1)
				sites[totpop++] = i;
		}
		else if(state[i]==2)
		{
			if(drand48()<alpha2)
				sites[totpop++] = i;
		}
	}

	k = totpop;
	for(i=0;i<totpop;i++)
	{
		j = k*drand48();
		nsite = sites[j];
		if(state[nsite]==1) NewBabe(nsite,1,dist1);
		else		    NewBabe(nsite,2,dist2);
		sites[j] = sites[--k];
	}
}

void NewBabe(int isite, int species, int dist)
{
	int newsite;

	newsite = dist*drand48() + 1;
	if(drand48()<0.5) newsite = -newsite;
	newsite += isite;
	if(newsite>=0&&newsite<CELLS)
		if(state[newsite]==0) state[newsite] = species;
}

void Mortality(void)
{
	int i;
	double xxx;

	for(i=0;i<CELLS;i++)
	{
		xxx = drand48();
		if(state[i]==1 && xxx<mu1) state[i] = 0;
		if(state[i]==2 && xxx<mu2) state[i] = 0;
	}
}

#if(MIXED)
void Mix(void)
{
	int i, j, k, temp;

	for(i=0;i<CELLS;i++)
	{
		j = CELLS*drand48();
		k = CELLS*drand48();
		temp = state[j];
		state[j] = state[k];
		state[k] = temp;
	}
}
#endif


void Measurements(void)
{
	int i;

	n1cnt = n2cnt = 0.0;

	for(i=0;i<CELLS;i++)
	{
		if(state[i]==1) 	n1cnt += 1.0;
		else if(state[i]==2) 	n2cnt += 1.0;
	}
	n1cnt /= (double)(CELLS);
	n2cnt /= (double)(CELLS);
}



void Spit_EPS_Header(FILE *fileid, int width, int length)
{
	double neww, newl, rat, xorigin, yorigin;

	neww = (double)width; newl = (double)length;

/* scale size */
	rat = 500.0/neww; neww = 500.0; newl = rat*newl;
	if(newl>700.0)
	{
		rat = 700.0/newl; newl = 700.0; neww = rat*neww;
	}

/* 8.5x11 page with no margins is 612x792 points */

	xorigin = (612 - neww)/2.0;
	yorigin = (792 - newl)/2.0;

/* defines compliance */
	fprintf(fileid,"%%!PS-Adobe-2.0 EPSF-2.0\n");

/* defines size of image */
	fprintf(fileid,"%%%%BoundingBox: %f %f %f %f\n",
			xorigin-5, yorigin-5,xorigin+neww+5,yorigin+newl+5);

/* save present graphics state */
	fprintf(fileid,"gsave\n");

/* put box around image */
	fprintf(fileid,"%f %f moveto\n",xorigin-1,yorigin-1);
	fprintf(fileid,"%f %f lineto\n",xorigin-1,yorigin+newl+1);
	fprintf(fileid,"%f %f lineto\n",xorigin+neww+1,yorigin+newl+1);
	fprintf(fileid,"%f %f lineto\n",xorigin+neww+1,yorigin-1);
	fprintf(fileid,"%f %f lineto\n",xorigin-1,yorigin-1);
	fprintf(fileid,"stroke\n\n",xorigin-1,yorigin-1);

/* define procedure "bufstr" to read strings of "width" characters */
	fprintf(fileid,"/bufstr %d string def\n\n",width);

/* move to origin of image */
	fprintf(fileid,"%f %f translate\n",xorigin,yorigin);

/* define the scale (default is 1 pt I guess) */
	fprintf(fileid,"%f %f scale\n\n",neww,newl);

/* width, height, and number of bits per pixel */
	fprintf(fileid,"%d %d 4\n",width,length);

/* a matrix definition */
	fprintf(fileid,"[%d 0 0 %d 0 %d]\n",width,-length,length);

/* guidelines on how to read the image data */
	fprintf(fileid,"{currentfile bufstr readhexstring pop} bind image\n");
}

void Spit_EPS_Row(FILE *fileid, int species)
{
	int	i;

	for(i=0;i<CELLS;i+=PSSKIP)
	{
		if(state[i]==species) 	fprintf(fileid,"0");
		else 		 	fprintf(fileid,"F");
		if((i+1)%60==0) fprintf(fileid,"\n");
	}
	fprintf(fileid,"\n");
}

void Spit_EPS_Trailer(FILE *fileid)
{
	fprintf(fileid,"grestore\n showpage\n");
}