#include <math.h>
#include <assert.h>
#include <malloc.h>
#include "misc.h"
#include "bsode.h"
#include "matvec.h"

extern int odeint(double ystart[], int nvar, double *x1, double x2,
    double eps, double h1, double hmin, int *nok, int *nbad,
    void(*derivs)(int,double,double[],double[]),
    int(*step)(double[],double[],int,double*,double,double,double[],
	double*,double*,void(*)(int,double,double[],double[])));

extern int bsstep(double y[], double dydx[], int nv, double *xx, double htry,
    double eps, double yscal[], double *hdid, double *hnext,
    void(*derivs)(int,double,double*,double*));


BSODE *BsodeAlloc(int n, double t, double *y, F_EVAL f, int stiff_flag,
    double eps, double zero)
{
    BSODE *b = Calloc(1, sizeof(BSODE));
    b->N = n;
    b->T = t;
    b->Y = y;
    b->F = f;
    if(stiff_flag != 0)
	Fatal("BsodeAlloc: stiff_flag(%d) must be 0", stiff_flag);
    b->ISTATE = 0;
    b->EPS = eps;
    if(zero == 0.0)
	b->HMIN = 1e-16*VecLength(n,y);
    else
	b->HMIN = zero;
    b->H0 = 1;
    return b;
}

/* Do the actual integration.  Return the actual TOUT */
double BsodeIntegrate(BSODE *b, double TOUT)
{
    int direction = TOUT - b->T > 0 ? 1 : -1;
    int nok, nbad;

    while(direction*(b->T - TOUT) < 0)
    {
	int status;
	status = odeint(b->Y, b->N, &b->T, TOUT, b->EPS, b->H0, b->HMIN,
	    &nok, &nbad, b->F, bsstep);

	switch(status)
	{
	case 0:		/* OK */
	    assert(b->T == TOUT);
	    break;
	case 1:		/* stepsize underflow */
	    assert(b->H0 > 0.0);
	    b->H0 *= 100;
	    break;
	case 2:		/* too many steps */
	    break;
	}
    }
    return b->T;
}


void BsodeFree(BSODE *b)
{
    free(b);
}
