#!/usr/bin/env python

import sys
from math import *

# Solution to Homework 2.

TOL = 1.e-4

# Want to solve xi sin xi = eta cos xi (even), or -xi cos xi = eta
# sin xi (odd), subject to xi^2 + etai^2 = u0.

def g(z):				# equation to solve, z is xi

    eta = u0 - z*z
    if eta < 0: eta = 0.0		# avoid problems with z^2 ~ u0
    eta = sqrt(eta)

    if even == 1:
	return z*sin(z) - eta*cos(z)
    else:
	return z*cos(z) + eta*sin(z)

def solve(z1, z2, epsilon):

    # Use bisection until the range is less than the specified
    # tolerance.

    while z2 - z1 > epsilon:
	zm = 0.5*(z1 + z2)
	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))


if __name__ == '__main__':

    oddeven = ("odd", "even")
    start = 0
    finish = 1

    u_max = 100.0		       	# maximum scaled potential
    du = 0.1

    argc = len(sys.argv)
    if argc > 1:			# specify u_max on command line
	u_max = float(sys.argv[1])

    if argc > 2:			# select even/odd with 2nd argument
	which = int(sys.argv[2])
	if which <= 0: finish = 0
	elif which == 1: start = 1

    u0 = du
    while u0 <= u_max + TOL*du:

	dz = 0.025*sqrt(u0)		# should be small enough...

        even = start
        while even <= finish:

	    # Systematically seek new roots between 0 and u0.

	    zl = 0.0
	    while zl < sqrt(u0)-TOL*dz:
		zr = zl + dz
		if zr > sqrt(u0): zr = sqrt(u0)

		if g(zl)*g(zr) <= 0 and (even or zl > 0):
		    xi = solve(zl, zr, TOL)
		    print u0, xi, xi*xi

		zl = zr

            even += 1

        u0 += du