#include "OdeFunc.h"
#include "ode-steppers.h"

using namespace blitz;

OdeFunc::OdeFunc(Range r)
    : _range(r)
{
    _stepper = euler_stepper;
}

ComplexVector OdeFunc::Step(double x, double& step_size, 
                            ComplexVector y, double prec)
{
    assert (_stepper);
    return _stepper(*this,x,step_size,y,prec);
}


// Integrate in one shot.
ComplexVector OdeFunc::Solve(double x0, double xf, double step_size,
                             ComplexVector y0, double prec)
{
    assert(_stepper);
    ComplexVector y(_range);
    y = y0;
    for (double x = x0; x <= xf; x += step_size) {
        y = this->Step(x,step_size, y, prec);
//        cerr << "x=" << x << ", y= " << y << endl;
        double mag =  vector_magnitude(y);
        if (fabs(1-mag) > 0.001)
            cerr << "|y|=" << mag << endl;
    }
    return y;
}

// Integrate over discrete regions.  This will integrate from x0[0] to
// xf[0] with step_size[0], then restart at x0[1], etc.  This can be
// used to give more steps to certain regions, or to break the
// integration at a discontinuity.
ComplexVector OdeFunc::Solve(double x0[], double xf[], double step_size[],
                             int nregions, ComplexVector y0, double prec)
{
    assert(_stepper);
    ComplexVector y(_range);
    y = y0;
    for (int region = 0; region < nregions; ++region) {
        int steps=0;
        for (double x = x0[region]; x <= xf[region]; x += step_size[region]) {
            y = this->Step(x,step_size[region], y, prec);
            ++steps;
        }
        cerr << "Region " << region << ": " << steps << " steps\n";
    }
    return y;
}