//=========================================================================
#define PROGRAMMER "SOLUTIONS, NoFrills"
#define PROG_CODE  "soln"
#define COURSE     "ECE-1021"
#define YEAR       (2003)
#define TERM       "Fall"
#define SECTION    (0)
#define ASSIGNMENT "HW #2A"
#define REVISION   (0)
#define TITLE      "Beyond Numeracy"
#define SUBTITLE   "Study in Chaos"
#define EMAIL      "wbahn@eas.uccs.edu"
#define FILENAME   "02a0soln.c"
//=========================================================================

// PROBLEM:
//
// Simulate the population of some item (egrets) according to the
// Feigenbaum model presented in the text.
//
// The model is supposed to demonstrate that, for some combinations of
// inputs, the system stabilizes to a single value, for others it
// oscillates between two values after it settles down, and for others
// it appears to generate random fluctuations.
//
// PSEUDOCODE:
//
// 1) TASK: Get data from the user
// 1.1) GET: Rate (rate: floating point value between 0 and 4)
// 1.2) GET: Population (pop: integer between 0 and 1,000,000)
// 1.3) GET: Years (pop: positive integer)
// 2) TASK: Compute and Display Year-by-Year Population
// 2.1) SET: year = 0
// 2.2) PUT: year and pop
// 2.3) WHILE: (year < years)
// 2.3.1) SET: year++
// 2.3.2) SET: pop *= (rate)*(1- pop/N)  (N = 1,000,000)
// 2.3.3) PUT: year and pop

//== QUESTIONS TO BE TURNED IN WITH SOURCE CODE ===========================
/*
3) QUESTION
	As the problem states, for certain values of R this system oscillates
	between two values each year (after some number of iterations). One
	such value, 3.14, was given. Find another value and determine what the
	two population values are that the system oscillates between. Discuss
	these results in view of the answer to Question #2 above.

	ANSWER
	Following the basic approach used to answer Qeustion #2 (that was
	turned in with the pseudocode), in order to the system to oscillate
	between two values requires that, after some number of years, we end
	up with a situation where:

	B = rate * A * ( 1 - (A/N) )  (for one year)
	A = rate * B * ( 1 - (B/N) )  (for the next year)

	Since we have two equations and three unknowns, what we can get are
	the required relationships between the three unknowns - and then
	by picking one of them we can find the other two.

	B = A*rate - (A^2)*(rate/N)
	A = B*rate - (B^2)*(rate/N)

	Eliminating the right hand terms, we get:

	(B^3 - A^3) = A*(B^2)*rate - (A^2)*B*rate = (rate)(AB)(B-A)

	We know that (B-A) is a factor of (B^3 - A^3), so dividing both
	side by (B-A) gives us:

	B^2 + AB + A^2 = (rate)(AB)

	B^2 + A^2 = (rate)(AB) - (AB) = (AB)(rate - 1)

	solving for rate, we get:

	rate = 1 + (B/A) + (A/B)

	This is one constraint that must be satisfied. Notice the second
	two terms are recipricols of each other. Let's define the parameter
	x = B/A and substitute it in:

	rate = 1 + (x + 1/x)

	Either x is greater than one or 1/x is greater than one. This
	fact all by itself means that rate must be greater than two. The
	minimum value of rate can be found by taking the derivative and
	setting it equal to zero. Doing this we find:

	d(rate)/dx = 1 - 1/x^2 = 0 which occurs at x = 1

	So the minimum value of rate is 3.

	Given a rate, we need to find the value of x such that:

	(x + 1/x) = (rate - 1)

	Or, in other words, we need to find the zero of the function:

	(x + 1/x) - (rate - 1) = 0

	This is simply a quadratic equation (multiply through by x) or
	we could use Newton's method from HW#1.

	Keep in mind that (A, B, and rate) are a matched set of values (for
	a given oscillating condition). At this point, if we are given
	a rate, we can find what the ratio is between the two oscillating
	values of the population, but we don't know that those two values
	actually are. We can get this by simply plugging the results thus
	far back into one of the original equations:

	B = rate * A * ( 1 - (A/N) )
	x = B/A  => B = A*x

	A*x = rate*A*(1 - A/N)
	x = rate*(1 - A/N)

	Solving for A:

	A = N*(1 - x/rate)

	As a check, let's see if these results agree with the example value
	given in the text:

	rate = 3.14
	A = 538007
	B = 780464

	x = B/A = 1.45066
	rate = 1 + x + 1/x = 3.14000

	x/rate = 0.461992
	A = N * (1 - x/rate) = 538007

	Life is good.

	So lets find another oscillating value:

	let's pick a value of x

	x = 1.75

	===============================================
	SECOND STABLE OSCILLATING POINT:

	rate = 1 + x + 1/x = 3.32142
	A = N * (1 - x/rate) = 473118
	B = A*x = 827957

	Running the prpgram with this rate yielded an
	oscillation between 473121 and 827955

	Note that the answer to Question #2 says that
	for this same rate, there should be a stable
	population of pop = N * ( 1 - (1/rate) )
	which is 698924. If we start the system out
	at this point, is stays close to this point for
	several iterations and then slowly (but at
	an accelerating pace) diverges and finally
	settles in the oscillating point.

	So there IS a static value that is "stable"
	but that point is very sensitive. It is what
	is called an "unstable equilibrium" like when
   a pencil is balanced on one end. 
	===============================================
*/
//== INCLUDE FILES ========================================================
#include <stdio.h> // printf()

//== MACROS ===============================================================
#define MAX_POP (1000000L)

//== FUNCTION PROTOTYPES ==================================================
void PrintHeader(void);

//== MAIN FUNCTION ========================================================
int main(void)
{
	long int pop, year, years;
	double rate;

	PrintHeader();
	printf("\n"); // Print a blank line

	// 1) TASK: Get data from the user
	printf("Enter the Growth Rate (0 <= rate <= 4): ");
	scanf("%lf", &rate);

	printf("Enter the initial population: ");
	scanf("%li", &pop);

	printf("Enter the number of years to simulate: ");
	scanf("%li", &years);

	// 2) TASK: Compute and Display Year-by-Year Population
	year = 0;
	printf("\n");
	printf("Growth Rate: %7.4f\n", rate);
	printf("\n");
	printf("       Year   Population\n");
	printf(" %10li    %7li\n", year, pop);

	while(year < years)
	{
   	year++;
		pop *= rate * (1 - ((double) pop)/((double) MAX_POP) );
		printf(" %10li    %7li\n", year, pop);
	}

	return(0);
}

//== SUPPORT FUNCTIONS ====================================================
void PrintHeader(void)
{
	printf("========================================"
			 "=======================================\n");
	printf("Course....... %s-%i (%s %i)\n", COURSE, SECTION, TERM, YEAR);
	printf("Programmer... %s (%s)\n", PROGRAMMER, PROG_CODE);
	printf("Assignment... %s (Rev %i) (Source Code in %s)\n",
		ASSIGNMENT, REVISION, FILENAME);
	printf("Description.. %s\n", TITLE);
	printf("              %s\n", SUBTITLE);
	printf("========================================"
			 "=======================================\n");

	return;
}