/*
cc seir.c -o bseir -lm
*/

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

#define	HEALTHY		0
#define	EXPOSED		1
#define	INFECTED	2
#define	IMMUNE		3

#define	EXP_TIME	3
#define	INF_TIME	3
#define	IMM_TIME	3

#define	INF_RATE	(0.45)
#define	INIT_EXP	(0.1)

#define	MAXTIME		500
#define	CELLS		1000

int	state[CELLS], timer[CELLS];
double	inf_rate, scnt, ecnt, icnt, rcnt;

void InitStuff();
void SpreadDisease(void);
void StoE(int, int, int);
void ChangeStates(void);
void Measurements(void);
void Spit_EPS_Header(FILE *, int, int);
void Spit_EPS_Row(FILE *);
void Spit_EPS_Trailer(FILE *);

int main(void)
{
	int i;
	FILE	*fopen(const char *, const char *), *fileid;

	fileid = fopen("seirim.eps","w");

	InitStuff();

	Measurements();

	Spit_EPS_Header(fileid,CELLS,MAXTIME);

	for(i=0;i<MAXTIME;i++)
	{
		Spit_EPS_Row(fileid);

		SpreadDisease();
		ChangeStates();

		Measurements();
		printf("%d %f %f %f %f\n",
			i, scnt, ecnt, icnt, rcnt);
	}
	Spit_EPS_Trailer(fileid);

	return (0);
}

void InitStuff(void)
{
	int i;
	unsigned long seed=145739853;
	srand48(seed);

	inf_rate = INF_RATE;
	for(i=0;i<CELLS;i++)
	{
		state[i] = HEALTHY;
		timer[i] = 0;
		if(drand48()<INIT_EXP)
		{
			state[i] = EXPOSED;
			timer[i] = EXP_TIME * drand48() + 1;
		}
	}
}

void SpreadDisease(void)
{
	int i;

	if(state[0]==HEALTHY) StoE(0,CELLS-1,1);
	for(i=1;i<CELLS-1;i++)
		if(state[i]==HEALTHY) StoE(i,i-1,i+1);
	if(state[CELLS-1]==HEALTHY) StoE(CELLS-1,CELLS-2,0);
}

void StoE(int cell, int cellL, int cellR)
{
	int cnt=0;

	if(state[cellL]==INFECTED) cnt = 1;
	if(state[cellR]==INFECTED) cnt++;
	if(cnt)
	{
		if(drand48()<(double)cnt*inf_rate)
		{
			state[cell] = EXPOSED;
			timer[cell] = EXP_TIME;
		}
	}
}

void ChangeStates(void)
{
	int i;

	for(i=0;i<CELLS;i++)
	{
		switch(state[i])
		{
			case IMMUNE:
				timer[i]--;
				if(timer[i]==0)
					state[i] = HEALTHY;
				break;

			case INFECTED:
				timer[i]--;
				if(timer[i]==0)
				{
					state[i] = IMMUNE;
					timer[i] = IMM_TIME;
				}
				break;

			case EXPOSED:
				timer[i]--;
				if(timer[i]==0)
				{
					state[i] = INFECTED;
					timer[i] = INF_TIME;
				}
				break;

			default:
				break;
		}
	}
}

void Measurements(void)
{
	int i;

	scnt = ecnt = icnt = rcnt = 0.0;

	for(i=0;i<CELLS;i++)
	{
		if(state[i]==HEALTHY) 		scnt += 1.0;
		else if(state[i]==EXPOSED) 	ecnt += 1.0;
		else if(state[i]==INFECTED) 	icnt += 1.0;
		else if(state[i]==IMMUNE) 	rcnt += 1.0;
	}
	scnt /= (double)(CELLS);
	ecnt /= (double)(CELLS);
	icnt /= (double)(CELLS);
	rcnt /= (double)(CELLS);
}


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

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

	if(neww>600.0)
	{
		rat = 600.0/neww;
		neww = 600.0;
		newl = rat*newl;
	}
	if(newl>780.0)
	{
		rat = 780.0/newl;
		newl = 780.0;
		neww = rat*neww;
	}
/* 8.5x11 page with no margins is 612x792 points */

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

	fprintf(fileid,"%%!PS-Adobe-2.0 EPSF-2.0\n");
	fprintf(fileid,"%%%%BoundingBox: %f %f %f %f\n",
			xorigin, yorigin,xorigin+neww,yorigin+newl);
	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);

	fprintf(fileid,"/bufstr %d string def\n\n",width);
	fprintf(fileid,"%f %f translate\n",xorigin,yorigin);
	fprintf(fileid,"%f %f scale\n\n",neww,newl);
	fprintf(fileid,"%d %d 4\n",width,length);
	fprintf(fileid,"[%d 0 0 %d 0 %d]\n",width,-length,length);
	fprintf(fileid,"{currentfile bufstr readhexstring pop} bind image\n");
}


void Spit_EPS_Row(FILE *fileid)
{
	int	i, j=0;

	for(i=0;i<CELLS;i++,j++)
	{
		if(state[i]==HEALTHY) 		fprintf(fileid,"F");
		else if(state[i]==EXPOSED) 	fprintf(fileid,"5");
		else if(state[i]==INFECTED) 	fprintf(fileid,"0");
		else if(state[i]==IMMUNE) 	fprintf(fileid,"B");
		if((j+1)%60==0) fprintf(fileid,"\n");
	}
	fprintf(fileid,"\n");
}


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