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

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

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


LDBSODE *LDBsodeAlloc(int n, long double t, long double *y, LD_EVAL f, int stiff_flag,
    long double eps, long double zero)
{
    LDBSODE *b = Calloc(1, sizeof(LDBSODE));
    double temp_y[n];
    int i;
    b->N = n;
    b->T = t;
    b->Y = y;
    b->F = f;
    if(stiff_flag != 0)
	Fatal("LDBsodeAlloc: stiff_flag(%d) must be 0", stiff_flag);
    b->ISTATE = 0;
    b->EPS = eps;
    for(i=0; i<n; i++)temp_y[i]=y[i];
    if(zero == 0.0)
	b->HMIN = 1e-16*VecLength(n,temp_y);
    else
	b->HMIN = zero;
    b->H0 = 1;
    return b;
}

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

    while(direction*(b->T - TOUT) < 0)
    {
	int status;
	status = ld_odeint(b->Y, b->N, &b->T, TOUT, b->EPS, b->H0, b->HMIN,
	    &nok, &nbad, b->F, ld_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 LDBsodeFree(LDBSODE *b)
{
    free(b);
}
