/* 

finds energy eigenvalues for a 1-d square well via the bisection
method for root finding.  the user needs
to supply a (closed) energy interval [a,b].
the depth of the well is V0 (V0 > 0).
*/

#include <iostream>
#include <cstdlib>
#include <cmath>
#include <gsl/gsl_math.h>

using std::endl;
using std::cout;

// define function:


double func_f(double x){
  const double V0 = 83.0;   // well depth on MeV
  return 2.0*pow(0.0483*(x+V0),0.5)*tan(2.0*pow(0.0483*(x+V0),0.5))
    - 2.0*pow((-0.0483*x), 0.5);
}


int main() 
{    
   int N= 0;
   int N_limit = 1000;           // max number of trys

   double a = -76;   // lower bound for [a,b]
   double b = -73.0;      // upper bound for [a,b]

   double tol = 1.0e-4;          // tolerance for finding zero

   if (func_f(a) == 0.0) {
     cout << a << " is a root." << endl;
     return 0;
   }

   if (func_f(b) == 0.0) {
     cout << b << " is a root." << endl;
     return 0;
   }

   cout << "f(a) = " << func_f(a) << endl;
   cout << "f(b) = " << func_f(b) << endl;


   if (GSL_SIGN(func_f(a)) * GSL_SIGN(func_f(b)) > 0.0){
     cout << "[" << a << "," << b << "] does not bracket a root." << endl;
     return 0;
   }
   
   double x_low = a;
   double x_high = b;
   double x_mid;

   while ((fabs(x_high - x_low) >= tol) && (N < N_limit)){
     x_mid = 0.5*(x_low + x_high);
     if ( GSL_SIGN(func_f(x_mid)) * GSL_SIGN(func_f(x_high)) < 0){
       x_low = x_mid;
     }
     else {
       x_high = x_mid;
     }
     ++N;
   }

   if(N > N_limit)
     cout << "iteration limit exceeded" << endl;
   else
     cout << "estimated root = " << x_mid << " MeV with " << N << " iterations." << endl;
   
   return 0;
}