/* allocation simulation 

cc simpall.c -o bsimp -lm
*/

#include	<stdlib.h>
#include 	<stdio.h>
#include 	<math.h>
#include 	<time.h>
	
#define	MAXGENS	 	100
#define	NPLANTS		1000
#define	SEASONLEN	100.0
#define GROWTH		0.02
#define	MUTATERATE	0.05
#define LZERO		1.0
#define SOLNTYPE	1 /* 1-constant, 2-switch */

typedef struct plant_def
{
	double	trait;
	double	size;
	double	fitness;
} Plant;
Plant plants1[NPLANTS], plants2[NPLANTS];

void InitSeed(void);
void InitPlants(Plant *plnt);
void ZeroPlants(Plant *plnt);
void GrowPlants(Plant *plnt);
void ReproducePlants(Plant *parents, Plant *babies);
void PrintPlants(Plant *plnt);

int main(void)
{
	int	ttt;
	Plant	*babes, *adults, *tmpy;

	adults = plants2;
	babes = plants1;
	InitSeed();
	InitPlants(babes);

	for(ttt=0;ttt<MAXGENS;ttt++)
	{
		tmpy = adults;
		adults = babes;
		babes = tmpy;

		ZeroPlants(adults);
		GrowPlants(adults);
		ReproducePlants(adults,babes);
	}
	PrintPlants(adults);
	return(0);
}

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 *plnt)
{
	int i;

	for(i=0;i<NPLANTS;i++,plnt++)
	{
#if(SOLNTYPE==1) 
		plnt->trait = drand48();
#endif
#if(SOLNTYPE==2)
		plnt->trait = drand48()*SEASONLEN;
#endif
	}
}

void ZeroPlants(Plant *plnt)
{
	int i;
	for(i=0;i<NPLANTS;i++,plnt++)
	{
		plnt->size = LZERO;
		plnt->fitness = 0.0;
	}
}

void GrowPlants(Plant *plnt)
{
	int i;
	double efac;

	for(i=0;i<NPLANTS;i++,plnt++)
	{
#if(SOLNTYPE==1)
		efac = exp(GROWTH*(1.0-plnt->trait)*SEASONLEN);
		plnt->size = LZERO*efac;
		plnt->fitness = (GROWTH*plnt->trait*LZERO)*(efac-1.0) 
					/(GROWTH*(1.0-plnt->trait));
#endif
#if(SOLNTYPE==2)
		plnt->size = LZERO*exp(GROWTH*plnt->trait);
		plnt->fitness = GROWTH*plnt->size*(SEASONLEN-plnt->trait);
#endif
	}
}

void ReproducePlants(Plant *parents, Plant *babies)
{
	int	i, j;
	double	totrep, rnum;
	double	brakes[NPLANTS];
	Plant *plnt;

	plnt = parents; /* measure the total fitness in the habitat */
	totrep = 0.0;
	for(i=0;i<NPLANTS;i++,plnt++)
	{
		totrep += plnt->fitness;
		brakes[i] = totrep;
	}

	for(i=0;i<NPLANTS;i++,babies++) /* Clonal reproduction */
	{
		rnum =  drand48() * totrep; /* choose adults proportional to fitness */
		j = (int)((rnum*NPLANTS)/totrep);

		if(rnum<brakes[0]) j = 0;
		else if(brakes[j]<rnum)
			while(rnum > brakes[j]) j++;
		else 
		{
			while(rnum < brakes[j]) j--;
			j++;
		}
		babies->trait = parents[j].trait;

		if(drand48()<MUTATERATE)
#if(SOLNTYPE==1)
			babies->trait = drand48();
#endif
#if(SOLNTYPE==2)
			babies->trait = drand48()*SEASONLEN;
#endif
	}
}

/* dist[] holds freq distribution of trait vals
   flag[] marks counted individuals (0-uncounted,-1 marks halotype(?),
					1-counted) */
void PrintPlants(Plant *plnt)
{
	int i, j, flag[NPLANTS], dist[NPLANTS], biggie=0;

	for(i=0;i<NPLANTS;i++)
		flag[i] = dist[i] = 0;

	for(i=0;i<NPLANTS;i++)
	{
		if(flag[i]==0) /* uncounted */
		{
			dist[i] = 1;
			flag[i] = -1;
			for(j=i+1;j<NPLANTS;j++) /* seek out clones */
			{
				if(plnt[i].trait==plnt[j].trait)
				{
					dist[i]++;
					flag[j] = 1;
				}
			}
			if(dist[i]>biggie) biggie = dist[i];
		}
	}

	for(i=0;i<NPLANTS;i++,plnt++) /* print unique clones */
	{
		if(flag[i]==-1)
			printf("%f  %f  %f %f\n",
				plnt->trait,plnt->size,plnt->fitness,
				5.0*dist[i]/(double)biggie);
	}
}