/*
cc detspat2b.c -o bdetspat -lm
*/
#include <stdlib.h>
#include <stdio.h>

#define	MAXTIME		5000
#define	OUTTIME		100
#define	PSTIME		50
#define	CELLS		100
#define	PSSKIP		1

#define	DELTAT		0.1
#define	STEPSPERT	10 /* 1/DELTAT */
#define	ALPHA1	0.25
#define	ALPHA2	0.25
#define	BETA12	0.075 /* b12+b21<0 */
#define	BETA21	0.1
#define	MU1	0.05
#define	MU2	0.05

double	n1[CELLS], n2[CELLS];
double	babes1[CELLS], babes2[CELLS];
double	disp1[CELLS],  disp2[CELLS];
int	dist1, dist2;
double	n1cnt, n2cnt;
double	alpha1, alpha2, beta12, beta21, mu1, mu2;

void InitStuff(void);
void MakeBabes(void);
void CompeteAndDie(void);
void DisperseBabes(double *, double *, double *, double *, int);
void AddBabes(void);
void Measurements(void);
void Spit_EPS_Header(FILE *, int, int);
void Spit_EPS_Row(FILE *, double *);
void Spit_EPS_Trailer(FILE *);

int main(void)
{
	int i, j;
	FILE	*fopen(const char *, const char *), *fileid1, *fileid2;

	fileid1 = fopen("detspat1.eps","w");
	fileid2 = fopen("detspat2.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++)
	{
		if(i%PSTIME==0)
		{
			Spit_EPS_Row(fileid1,n1);
			Spit_EPS_Row(fileid2,n2);
		}

		if(i%OUTTIME==0)
		{
			Measurements();
			printf("%d %f %f\n", i, n1cnt, n2cnt);
		}

		for(j=0;j<STEPSPERT;j++)
		{
			MakeBabes();
			DisperseBabes(disp1,babes1,n1,n2,dist1);
			DisperseBabes(disp2,babes2,n1,n2,dist2);
			CompeteAndDie();
			AddBabes();
		}
	}
	Spit_EPS_Trailer(fileid1);
	Spit_EPS_Trailer(fileid2);

	for(j=0;j<CELLS;j++) printf("%3d %f %f\n", j, n1[j], n2[j]);

	return (0);
}

void InitStuff(void)
{
	int i;

	alpha1 = ALPHA1*DELTAT;
	alpha2 = ALPHA2*DELTAT;
	beta12 = BETA12*DELTAT;
	beta21 = BETA21*DELTAT;
	mu1 = MU1*DELTAT;
	mu2 = MU2*DELTAT;

	for(i=0;i<CELLS/2;i++)
	{
		n1[i] = 1.0;
		n2[i] = 0.0;
	}
	for(i=CELLS/2;i<CELLS;i++)
	{
		n1[i] = 0.0;
		n2[i] = 1.0;
	}
}

void MakeBabes(void)
{
	int j;

	for(j=0;j<CELLS;j++)
	{
		babes1[j] = alpha1*n1[j];
		babes2[j] = alpha2*n2[j];
	}
}

void CompeteAndDie(void)
{
	int j;
	double oldpop1, oldpop2;

	for(j=0;j<CELLS;j++)
	{
		oldpop1 = n1[j]; oldpop2 = n2[j];

		n1[j] -= (beta12*oldpop2 + mu1)*oldpop1;
		n2[j] -= (beta21*oldpop1 + mu2)*oldpop2;
		if(n1[j]<0.0) n1[j] = 0.0;
		if(n2[j]<0.0) n2[j] = 0.0;
	}
}


void DisperseBabes(double *disp, double *babes, double *pop1, double *pop2, int dist)
{
	int j, kmin, k, kmax;
	double babes_per_cell;

	for(j=0;j<CELLS;j++) disp[j] = 0.0;
	for(j=0;j<CELLS;j++)
	{
		babes_per_cell = babes[j]/(2.0*dist+1.0);
		kmin = j - dist;
		if(kmin<0)
		{
			disp[0] += (-kmin)*babes_per_cell*(1.0-pop1[0]-pop2[0]);
			kmin = 0;
		}
		kmax = j + dist + 1;
		if(kmax>CELLS)
		{
			disp[CELLS-1] += (kmax-CELLS)*babes_per_cell
				*(1.0-pop1[CELLS-1]-pop2[CELLS-1]);
			kmax = CELLS;
		}
		for(k=kmin;k<kmax;k++)
			disp[k] += babes_per_cell*(1.0-pop1[k]-pop2[k]);
	}
}


void AddBabes(void)
{
	int j;

	for(j=0;j<CELLS;j++)
	{
		n1[j] += disp1[j];
		n2[j] += disp2[j];
	}
}

void Measurements(void)
{
	int i;

	n1cnt = n2cnt = 0.0;
	for(i=0;i<CELLS;i++)
	{
		n1cnt += n1[i];
		n2cnt += n2[i];
	}
	n1cnt /= (double)(CELLS);
	n2cnt /= (double)(CELLS);
}

void Spit_EPS_Header(FILE *fileid, int width, int length)
{
	double neww, newl, rat;
	double 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, yorigin,xorigin+neww,yorigin+newl);

/* 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 that I don't understand */
	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, double *pop)
{
	int	i;

	for(i=0;i<CELLS;i+=PSSKIP)
	{
		if(pop[i]>0.5) 	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");
	fprintf(fileid,"showpage\n");
}