/* gynodioecy simulation
cc gyno.c -o bgyno -lm
*/

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

#define	MAXGENS	 	1600
#define	PRNTTIME	5
#define	DIMEN		100
#define	NPLANTS		(DIMEN*DIMEN)
#define	PSOUT		1		/* toggle postscript on/off */
#define	PSOUTTIME	100
#define	MIX		0
#define	DELTAT		(1.0)		/* allows an arbitrary step size */
#define	STEPSPERT	((int)(1.0/DELTAT))
#define	POLLPROD	(1.0)
#define	POLLZERO	(3.5)
#define	HSEEDS		(4)
#define	FSEEDS		(5)
#define	POLLDIST	(4)
#define	HSEEDDIST	(2)
#define	FSEEDDIST	(2)
#define	DEATHRATE	(0.1)
#define	INITHERMS	(0.2)
#define	INITFEMS	(0.025)

typedef struct plant_def
{
	int 	alive;
	int	herm;
	double	inpollen;
	int	seeds;
	int	inseeds[2];  /*1-herms, 0-fems */
} Plant;
Plant plants[DIMEN][DIMEN];

double deathrate;
int	numh, numf, num_pnabes, num_hnabes, num_fnabes;
int	polldist, hsdist, fsdist;
int	pnabes[4*POLLDIST*POLLDIST][2];   /* neighbor index offsets for pollen, */
int	hnabes[4*HSEEDDIST*HSEEDDIST][2]; /*      hermaphrodite's seeds,  */
int	fnabes[4*FSEEDDIST*FSEEDDIST][2]; /*        and female's seeds.   */


void ArgumentControl(int argc, char **argv);
void InitStuff(void);
void InitSeed(void);
void InitPlants(Plant plnts[][DIMEN]);
int InitNabes(int *nabes, int dist);
void CollectPollen(Plant plnts[DIMEN][DIMEN]);
void DispersePollen(Plant plnts[DIMEN][DIMEN], int idim, int jdim);
int NewInd(int ind);
void SeedDisperse(Plant plnts[DIMEN][DIMEN]);
void PutSeed(Plant plnts[DIMEN][DIMEN], int *nabes, int num_nabes, 
		int idim, int jdim, int type);
void KillPlants(Plant *plnt);
void Germinate(Plant *plnt);
void MixPlants(Plant *plnt1);
void CountEm(Plant *plnt);
void EPS_Header(FILE *, int, int);
void EPS_Image(FILE *fileid, Plant plnts[][DIMEN]);
void EPS_Trailer(FILE *);

int main(int argc, char **argv)
{
	int	ttt, j;
	char	fname[20];
	FILE	*fopen(const char *, const char *), *fileid;

	deathrate = DEATHRATE;
	polldist = POLLDIST; hsdist = HSEEDDIST; fsdist = FSEEDDIST;
	if (argc>1)  ArgumentControl(argc,argv);
	InitStuff();

	for(ttt=0;ttt<MAXGENS;ttt++) /* measurement period */
	{
		if(ttt%PRNTTIME==0)
		{
			CountEm(&plants[0][0]);
			printf(" %f  %f\n",(double)numh/(double)NPLANTS,
				(double)numf/(double)NPLANTS);
		}

#if(PSOUT==1)
		if(ttt%PSOUTTIME==0)
		{
			sprintf(fname,"xy_%-1d.eps",ttt);
			fileid = fopen(fname,"w");
			EPS_Header(fileid,DIMEN,DIMEN);
			EPS_Image(fileid,plants);
			EPS_Trailer(fileid);
			fclose(fileid);
		}
#endif
		for(j=0;j<STEPSPERT;j++)
		{
			CollectPollen(plants);
			SeedDisperse(plants);
			KillPlants(&plants[0][0]);
			Germinate(&plants[0][0]);
#if(MIX==1)
			MixPlants(plants);
#endif
		}
	}
#if(PSOUT==1)
	fileid = fopen("xy_end.eps","w");
	EPS_Header(fileid,DIMEN,DIMEN);
	EPS_Image(fileid,plants);
	EPS_Trailer(fileid);
#endif
	return (0);
}

void ArgumentControl(int argc, char **argv)
{
	int argcount, skip;

	for(argcount=1; argcount<argc; argcount+=skip)
	{
		if(argv[argcount][0] == '-')
		{
			skip = 1;
			switch(argv[argcount][1])
			{
				case 'm':
					sscanf(argv[argcount+1],"%lf",&deathrate);
					skip++;
					break;
				case 'd':
					if(argv[argcount][2]=='p')
						sscanf(argv[argcount+1],"%d",&polldist);
					else if(argv[argcount][2]=='h')
						sscanf(argv[argcount+1],"%d",&hsdist);
					else if(argv[argcount][2]=='f')
						sscanf(argv[argcount+1],"%d",&fsdist);
					else exit(0);
					skip++;
					break;
				default:
					printf("**invalid command line argument >%c< \n",
						    argv[argcount][1]);
				case 'h':
					printf("command line arguments (and default values):\n");
					printf(" -h >> help\n");
					printf(" -m %f >> death rate\n", deathrate);
					printf(" -dp %d >> pollen dispersal\n", polldist);
					printf(" -dh %d >> herm seed dispersal\n", hsdist);
					printf(" -df %d >> female seed dispersal\n", fsdist);
					exit(0);
			}
		}
		else
		{
			printf("**invalid command line flag >%c<\n",argv[argcount][0]);
			printf("try -h for help.\n");
			exit(0);
		}
	}
}


void InitStuff(void)
{
	InitSeed();
	InitPlants(plants);
	num_pnabes = InitNabes(&pnabes[0][0],polldist);
	num_hnabes = InitNabes(&hnabes[0][0],hsdist);
	num_fnabes = InitNabes(&fnabes[0][0],fsdist);
}


void InitSeed(void)
{
	int	seed;
	time_t	nowtime;
	struct	tm *preztime;

	time(&nowtime);
	preztime = localtime(&nowtime);
	seed = (int)((preztime->tm_sec+1)*(preztime->tm_min+1)*
		(preztime->tm_hour+1)*(preztime->tm_year)*(preztime->tm_year));
	if(seed%2==0) seed++;

	srand48(seed);
}

void InitPlants(Plant plnts[DIMEN][DIMEN])
{
	int i,j;
	double xxx;

	for(i=0;i<DIMEN;i++)
	{
		for(j=0;j<DIMEN;j++)
		{
			plnts[i][j].alive = plnts[i][j].herm = 0;
			if(i<DIMEN/2 && j<DIMEN/2)  /* populate only one quadrat */
			{
				xxx = drand48();
				if(xxx<(INITHERMS+INITFEMS))
				{
					plnts[i][j].alive = 1;
					if(xxx<INITHERMS) plnts[i][j].herm = 1;
				}
			}
		}
	}
}

int InitNabes(int *nabes, int dist)
{
	int i, j, tot=0, tester, ii;

	tester = dist*dist;
	for(i=-dist;i<=dist;i++) /* finds nabes within dist */
	{
		ii = i*i;
		for(j=-dist;j<=dist;j++)
		{
			if(ii+j*j<=tester && !(i==0 && j==0))
			{
				nabes[tot++] = i;
				nabes[tot++] = j;
			}
		}
	}
	return(tot/2);
}

void CollectPollen(Plant plnts[DIMEN][DIMEN])
{
	int iplant, i, j;
	Plant *plnt;

	plnt = &plnts[0][0];
	for(iplant=0;iplant<NPLANTS;iplant++,plnt++)
	{
		plnt->inpollen = 0.0;
		plnt->seeds = plnt->inseeds[0] = plnt->inseeds[1] = 0;
	}

	plnt = &plnts[0][0];
	for(i=0;i<DIMEN;i++)
		for(j=0;j<DIMEN;j++,plnt++)
			if(plnt->alive && plnt->herm)
				DispersePollen(plnts,i,j);
}

void DispersePollen(Plant plnts[DIMEN][DIMEN], int idim, int jdim)
{
	int nn, i, j;
	double amnt;

	amnt = POLLPROD/(double)num_pnabes;
	for(nn=0;nn<num_pnabes;nn++)
	{
		i = NewInd(idim + pnabes[nn][0]);
		j = NewInd(jdim + pnabes[nn][1]);
		plnts[i][j].inpollen += amnt;
	}
}

int NewInd(int ind)	/* use periodic boundary conditions */
{
	if(ind<0) ind += DIMEN;
	if(ind>=DIMEN) ind -= DIMEN;
	return(ind);
}

void SeedDisperse(Plant plnts[DIMEN][DIMEN])
{
	int i, j, k;
	double fertprob;
	Plant *plnt;

	plnt = &plnts[0][0];
	for(i=0;i<DIMEN;i++)
	{
		for(j=0;j<DIMEN;j++,plnt++)
		{
			if(plnt->alive==1)
			{
				fertprob = DELTAT*plnt->inpollen/(POLLZERO+plnt->inpollen);
				if(plnt->herm==0)
				{
					for(k=0;k<FSEEDS;k++)
						if(drand48()<fertprob)
							PutSeed(plnts,&fnabes[0][0],num_fnabes,i,j,0);
				}
				else
				{
					for(k=0;k<HSEEDS;k++)
						if(drand48()<fertprob)
							PutSeed(plnts,&hnabes[0][0],num_hnabes,i,j,1);
				}
			}
		}
	}
}

void PutSeed(Plant plnts[][DIMEN], int *nabes, int num_nabes, int idim, int jdim, int type)
{
	int newspot, i, j;
						/* uniform dispersal */
	newspot = 2 * (int)(num_nabes*drand48());
	i = NewInd(idim + nabes[newspot]);
	j = NewInd(jdim + nabes[newspot+1]);
	plnts[i][j].inseeds[type]++;
}

void KillPlants(Plant *plnt)
{
	int i;

	for(i=0;i<NPLANTS;i++,plnt++)
		if(plnt->alive)
			if(drand48()<deathrate*DELTAT)
				plnt->alive = -1;	/* flag for death */
}

void Germinate(Plant *plnt)
{
	int iplant;
	double xxx;

	for(iplant=0;iplant<NPLANTS;iplant++,plnt++)
	{
		if(plnt->alive==0)
		{
			xxx = (double)(plnt->inseeds[0]+plnt->inseeds[1]);
			if(xxx>0.0)
			{
				plnt->alive = plnt->herm = 1;
				if(xxx*drand48()<(double)plnt->inseeds[0]) plnt->herm = 0;
			}
		}
		else if(plnt->alive==-1) plnt->alive = 0;
	}

}

void MixPlants(Plant *plnt1)
{
	int i, j, k;
	Plant *plnt;
	plnt = plnt1;

	for(i=0;i<NPLANTS;i++,plnt++)
	{
		j = NPLANTS*drand48();
		if(plnt1[j].alive)
		{
			plnt1[j].alive = plnt->alive;
			plnt->alive = 1;

			k = plnt1[j].herm;
			plnt1[j].herm = plnt->herm;
			plnt->herm = k;
		}
	}
}

void CountEm(Plant *plnt)
{
	int iplant;

	numh = numf = 0;
	for(iplant=0;iplant<NPLANTS;iplant++,plnt++)
	{
		if(plnt->alive)
		{
			if(plnt->herm)	numh++;
			else		numf++;
		}
	}

}



#define	PSWIDTHMAX	(300.0) /* 8.5x11 page with no margins */
#define	PSLENGTHMAX	(500.0) /*     is 612x792 points */

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

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

	rat = PSWIDTHMAX/neww;
	neww = PSWIDTHMAX;
	newl = rat*newl;

/*	if(neww>PSWIDTHMAX)
	{
		rat = PSWIDTHMAX/neww;
		neww = PSWIDTHMAX;
		newl = rat*newl;
	}
*/
	if(newl>PSLENGTHMAX)
	{
		rat = PSLENGTHMAX/newl;
		newl = PSLENGTHMAX;
		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 EPS_Image(FILE *fileid, Plant (*plnts)[DIMEN])
{
	int	i, j, k=0;

	for(i=0;i<DIMEN;i++)
	{
		for(j=0;j<DIMEN;j++)
		{
			if(plnts[i][j].alive)
			{
				if(plnts[i][j].herm) 	fprintf(fileid,"0"); /* dark */
				else 			fprintf(fileid,"9");
			}
			else 				fprintf(fileid,"F"); /* blank */
			if((++k)%60==0) fprintf(fileid,"\n");
		}
		fprintf(fileid,"\n");
	}
}


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