/* PHYS 307:  Solution to Homework #4, problems 1(a) and 1(c).
 *
 * Default is to compute the values of epsilon in the problem set,
 * comparing the analytic solutions with perturbation theory.
 *
 * For a plot of analytic versus perturbation theory results, give
 * a command-line argument and pipe into gpl:
 *
 *	hw4_1ac 1 | gpl -k -c 1 2 3 -x epsilon -y "dE/E0"
 */

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

typedef double	real;
#define TOL	1.e-8

real E0 = (M_PI*M_PI/4);
real epsilon;

real g(real z)
{
    real eta = z*z - epsilon*E0;
    if (eta <= 0) eta = 0;
    else eta = sqrt(eta);

    return eta*sin(0.5*eta)*sin(0.5*z) - z*cos(0.5*eta)*cos(0.5*z);
}

real solve(real z1, real z2, real tolerance)
{
    /* Use bisection until the range is less than the specified tolerance. */

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

    /* Perform one final false-position iteration. */

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

main(int argc, char *argv[])
{
    real z1 = 0.01*sqrt(E0), z2 = 1.99*sqrt(E0), z;
    if (argc == 1) {
	real eps[7] = {0.0, 0.01, 0.05, 0.1, 0.2, 0.5, 1.0};
	int i;
	for (i = 0; i < 7; i++) {
	    epsilon = eps[i];
	    z = solve(z1, z2, TOL);
	    printf("epsilon = %f, g(z) = %f, dE/E0 = %f, dE/E0(pert) = %f\n",
		   epsilon, g(z), z*z/E0 - 1, (0.5+1/M_PI)*epsilon);
	}
    } else {
	for (epsilon = 0; epsilon <= 1; epsilon += 1/256.) {
	    z = solve(z1, z2, TOL);
	    printf("%f %f %f\n",
		   epsilon, z*z/E0-1, (0.5+1/M_PI)*epsilon);
	}
    }
}










//