#include "shooting.h"

#include <stdexcept>

ivp<vec<2>> bvp( vec<2> f( double x, vec<2> u ),
     std::pair<double, double> x_range,
     std::pair<double, double> u_bndry,
     double eps_abs,
     unsigned int max_iterations,
     double h,
     std::pair<double, double> h_range
) {
    // We must find the appropriate initial slope 's0'
    //   - the ivp solver associated with this initial slope will be 'soln0'
    //   - the error of this solution at u(b) will be 'ub0err'

    double s0{ (u_bndry.second - u_bndry.first)/(x_range.second - x_range.first) };

    ivp<vec<2>> soln0{ f, x_range.first, vec<2>{u_bndry.first, s0},
                       h, h_range, eps_abs, vec<2>::norm };

    double ub0err{ soln0( x_range.second )[0] - u_bndry.second };

    std::clog << "u(" << x_range.second << ") = "
              << soln0( x_range.second )[0]
              << " when u'(" << x_range.first << ") = "
              << s0 << std::endl;

    // We must find the appropriate next slope 's1'
    //   - the ivp solver associated with this initial slope will be 'soln1'
    //   - the error of this solution at u(b) will be 'ub1err'

    double s1{ (u_bndry.second - ub0err - u_bndry.first)/
                        (x_range.second - x_range.first) };

    ivp<vec<2>> soln1{ f, x_range.first, vec<2>{u_bndry.first, s1},
                       h, h_range, eps_abs, vec<2>::norm };

    double ub1err{ soln1( x_range.second )[0] - u_bndry.second };
 
    std::clog << "u(" << x_range.second << ") = "
              << soln1( x_range.second )[0]
              << " when u'(" << x_range.first << ") = "
              << s1 << std::endl;

    // Now we have two approximations: (s0, ub0err) and (s1, ub1err )
    //  - iterate using the secant method!

    for ( unsigned int k{0}; k < max_iterations; ++k ) {
        // Use the secant method to find the next initial slope
        double s2{ (s1*ub0err - s0*ub1err)/(ub0err - ub1err) };
        
        ivp<vec<2>> soln2{ f, x_range.first, vec<2>{u_bndry.first, s2},
                           h, h_range, eps_abs, vec<2>::norm };

        double ub2err{ soln2( x_range.second )[0] - u_bndry.second };

        std::clog << "u(" << x_range.second << ") = "
                  << soln2( x_range.second )[0]
                  << " when u'(" << x_range.first << ") = "
                  << s2 << std::endl;

        if ( std::abs( ub2err ) < eps_abs ) {
            // Return the solver that has the required
            // initial slope to closely approximate the
            // solution to the boundary-value problem
            return soln2;
        }

        s0 = s1;
        s1 = s2;
        ub0err = ub1err;
        ub1err = ub2err;
    }

    throw std::runtime_error( "The shooting method did not converge after "
                            + std::to_string( max_iterations )
                            + " steps" );
}