////////////////////////////////////////////////////////////////////////
//
// $Id: OdeFunc.h,v 1.1 2001/10/26 19:19:13 bviren Exp $
//
// OdeFunc
//
// Package: elbo
//
// Solve dy/dx = F(x,y) where y and F are vectors.  OdeFunc is the
// base class.  Subclasses should define the function via operator().
//
// Contact: bv@bnl.gov
//
// Created on: Mon Sep 17 12:52:45 2001
//
////////////////////////////////////////////////////////////////////////


#ifndef ODEFUNC_H
#define ODEFUNC_H

#include "matrix.h"

class OdeFunc {
public:
    typedef ComplexVector (*OdeStepper)(OdeFunc& ode_func,
                                        double x, double& step_size, 
                                        ComplexVector y, double prec);

    OdeFunc(blitz::Range r = blitz::Range(1,3));

    // user overrides this operator.
    virtual ComplexVector operator()(double x, ComplexVector y) = 0;

    // Get the range used to create any related matrices/vectors :
    // defaults to (1,3)
    virtual blitz::Range GetRange(void) const { return _range; }

    // Returns vector y0 transported from x0 to xf by integrating in
    // one shot.  Not all steppers use prec.
    ComplexVector Solve(double x0, double xf, double step_size, 
                        ComplexVector y0, double prec);
    // Like above, but 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.  Not
    // all steppers use prec.
    ComplexVector Solve(double x0[], double xf[], double step_size[],
                        int nregions, ComplexVector y0, double prec);

    void SetStepper(OdeStepper s) { _stepper = s; }
    ComplexVector Step(double x, double& step_size, 
                       ComplexVector y, double prec);

private:
    blitz::Range _range;
    OdeStepper _stepper;
};

#endif  // ODE_H