/* Solution to Homework 4, question 1(b), based on Homework 3 solution.
 *
 * Default is to compute the values of epsilon in the problem set,
 * comparing the numerical solutions with perturbation theory.
 *
 * For a plot of numerical versus perturbation theory results, give
 * a command-line argument and pipe into gpl:
 *
 *	hw4_1b 1 | gpl -k -c 1 2 3 -x epsilon -y "dE/E0"
 */

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

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

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

int  even = 1;

/*------------------------------------------------------------------------*/

/* Integrator taken from Exercise 3.1. */

real xi;

real acc(real y, real yp, real x)	/* determine d2y/dx2, given x, y, yp */
{
    real pert = 0;
    if (fabs(x) < 0.5)
	pert = epsilon*E0;

    return (-xi*xi + pert) * y;
}

void take_a_step(real *y, real *yp, real *x, real dx)
{
    real a = acc(*y, *yp, *x);

    /* Predict: */

    (*y) += dx * ((*yp) + 0.5*dx*a);
    (*yp) += dx * a;

    (*x) += dx;

    /* Correct: */

    (*yp) += 0.5*dx*(acc(*y, *yp, *x) - a);
}

real integrate(real x0, real y0, real yp0, real dx, real x_max)
{
    real x = x0, y = y0, yp = yp0;

    while (x < x_max) {

	take_a_step(&y, &yp, &x, dx);	/* note: passing addresses here
					     -- function expects pointers */
    }

    return y;
}

/*----------------------------------------------------------------------*/

/* Function g ties the pieces together. */

real g(real z)				/* equation to solve, z is xi */
{
    xi = z;
    if (even)
	return integrate(0, 1, 0, 0.001, 1);
    else
	return integrate(0, 0, 1, 0.001, 1);
}

/*----------------------------------------------------------------------*/

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);
	}
    }
}











//