// PHYS 307:  Solution to Homework #1, problem 2.

#include <iostream>
#include <cmath>
using namespace std;

typedef double real;		// convenient alias

real g(real z)
{
    return (1+z)*pow(sin(z),2) - 1;
}

real solve(real z1, real z2, real epsilon)
{
    // Use bisection until the range is less than the specified tolerance.
    // Global n_iter is used to count iterations in the bisection loop.

    real zm;
    while (z2 - z1 > epsilon) {
	zm = (z1 + z2)/2;
	if (g(zm)*g(z1) > 0)
	    z1 = zm;
	else			// note no check for equality
	    z2 = zm;
    }

    // Perform a final false-position iteration.

    return z1 + (z2 - z1) * (-g(z1)) / (g(z2) - g(z1));
}

main()
{
    // Natural scale is ~pi/2.

    real z1, z2, dz = 0.1*M_PI, epsilon = 1.e-4, z;
    int count = 0;

    z1 = 0;
    z2 = dz;

    while (count < 5) {

	// Find the next range (z1, z2) bracketing a root.

	while (g(z2)*g(z1) > 0) {
	    z1 = z2;
	    z2 += dz;
	}

	// The range (z1, z2) brackets a root.  Solve for the root.

	z = solve(z1, z2, epsilon);

	count++;
	cout << "root #" << count << ": z = " << z
	     << ", g(z) = " << g(z) << endl;

	z1 = z2;
	z2 += dz;

    }
}










//