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

#define	DELTAT		0.1
#define	NPERT		((int)(1.0/DELTAT))
#define	MAXTIME 	(1000*NPERT)
#define	SIZEDEL		(0.025)
#define	NUMPNTS		((int)(KKKK/SIZEDEL)+1)

#define	NPATCH		20.0
#define ALPHA		0.05
#define	KKKK		1.0
#define	BETA		0.1
#define	MOVECOST	0.1
#define	TIMECOST	0.01

double	sizes1[NUMPNTS], sizes2[NUMPNTS], growth[NUMPNTS],
	graze[NUMPNTS], dubsize[NUMPNTS];
double	grazerate;

void InitStuff(void);
void InitSizeDist(double *sizedist);
void UpdateSizeDist(double *, double *);
double	Eaten(double, double *);
double	grazed(double, double);

int main(void)
{
	int	ttt;
	double	staytime, aveeaten, *s1, *s2, *stmp;

	InitStuff();

/* Determine the steady-state size distribution */

	s1 = sizes1; s2 = sizes2;
	InitSizeDist(s1);

	for(ttt=0;ttt<MAXTIME;ttt++)
	{
		UpdateSizeDist(s1,s2);
		stmp = s1; s1 = s2; s2 = stmp;
	}

/* Calculate intake rates */

	for(staytime=2.0;staytime<40.0;staytime+=1.0)
	{
		aveeaten = Eaten(staytime,s1);
		printf(" %4.1f %f\n", staytime, aveeaten);
	}
	return (0);
}

void InitStuff(void)
{
	int	size;
	double	val;

	grazerate = 1.0/NPATCH;			/* one patch has a grazer */
	for(size=1;size<NUMPNTS;size++)	/* tabulate growth/grazing functions */
	{
		val = ((double)size)*(double)SIZEDEL;
		dubsize[size] = val;		/* relates index to resource size */
		growth[size] = ALPHA*DELTAT*val*(1.0 - val/KKKK)/SIZEDEL;
		graze[size] = -(growth[size] - BETA*DELTAT*val/SIZEDEL);
	}
}

void InitSizeDist(double *sizedist)
{
	int	size;

	for(size=1;size<NUMPNTS-1;size++)
		sizedist[size] = 0.0;
	sizedist[NUMPNTS-1] = 1.0;	/* start with patches at largest resource value */
}

void UpdateSizeDist(double *old, double *new)
{
	int	size;

	new[1] = old[1] * (1.0 - (1.0-grazerate)*growth[1])
			+ grazerate*graze[2]*old[2];

	for(size=2;size<NUMPNTS-1;size++)
	{
		new[size] = old[size] *
			(1.0 - (1.0-grazerate)*growth[size] - grazerate*graze[size])
			+ (1.0-grazerate)*growth[size-1]*old[size-1]
			+ grazerate*graze[size+1]*old[size+1];
	}
	new[NUMPNTS-1] = old[NUMPNTS-1] * (1.0 - grazerate*graze[NUMPNTS-1])
		+ (1.0-grazerate)*growth[NUMPNTS-2]*old[NUMPNTS-2];
}

double Eaten(double stime, double *sizedist)
{
	int	size;
	double	xxx;

	xxx = -MOVECOST-TIMECOST*stime;
	for(size=0;size<NUMPNTS;size++)
		xxx += grazed(dubsize[size],stime)*sizedist[size];
	return (xxx/stime);
}

double	grazed(double n0, double stime)
{
	return( (BETA*KKKK/ALPHA) *
		( log( ((ALPHA*n0+KKKK*(BETA-ALPHA))*exp((BETA-ALPHA)*stime)
			- ALPHA*n0) / (KKKK*(BETA-ALPHA)) ) )
		- BETA*(BETA-ALPHA)*stime*KKKK/ALPHA );
}