/* root-finding algorithm */

#include	<stdlib.h>
#include 	<stdio.h>
#include	<math.h>
	
#define ALPHA		0.05
#define	KK		1.0
#define BETA		0.1
#define	MOVECOST	0.1
#define	TIMECOST	0.01

double	RootFn(double);
double	tgprime(double);
double	ggg(double);
double	alpha, beta, k, cost1, cost2, rstar;

int main(void)
{
	int	step;
	double	lowt, midt, hight;
	double	lowf, midf, highf, n, nfac;

	alpha = ALPHA;
	beta = BETA;
	k = KK;
	cost1 = MOVECOST;
	cost2 = TIMECOST;
	
	for(n=3.0,nfac=1.0;n<1000.0;n+=nfac*0.1,nfac*=2.0)
	{
		rstar = (alpha-beta/n)*k/alpha;
		lowt = 0.0001;
		hight = 1000.0;
		lowf = RootFn(lowt);
		highf = RootFn(hight);

		for(step=0;step<50;step++)
		{
			midt = (lowt + hight)/2.0;
			midf = RootFn(midt);

			if(lowf*midf<0.0)	hight = midt;
			else if(highf*midf<0.0)	lowt = midt;
			else break;
		}
		printf("%g %g \n",midt,ggg(midt)/midt); /* intake rates vs residence time */
	}
	return (0);
}

double	RootFn(double x)
{
	return( tgprime(x) - ggg(x)/x );
}

double	tgprime(double x)
{
	return( (ggg(1.001*x)-ggg(x))/(0.001*x) );
}

double	ggg(double x)
{
	double mess;

	mess = ( (alpha*rstar+k*(beta-alpha))*exp((beta-alpha)*x)
				-alpha*rstar) / k / (beta-alpha);
	return( (beta*k/alpha)*(log( mess ) - (beta-alpha)*x) 
			- (cost1+cost2*x) );
}