/* integrate.c
 * Integrate using the FORTRAN routine DVERK.
 */

#include <math.h>
#include <assert.h>
#include <malloc.h>
#include "misc.h"
#include "dverk.h"

#define NED 0

/* init_integrator()
 *
 * We must pass a whole bunch of stuff to DVERK, all of which must be
 * variables (no constants), because FORTRAN passes by reference.  These
 * must be globals.  I'm using all CAPs (even though it's ugly) to ensure
 * the names don't conflict with Isaac variable names.
 */


DVERK *DverkAlloc(int n, double time, double *y, F_EVAL f, int stiff_flag,
    double tol, double zero)
{
    DVERK *r = Calloc(1, sizeof(DVERK));

    if(stiff_flag != 0 || zero != 0.0)
	Fatal("DverkAlloc: stiff_flag must be 0 && zero must be 0.0");

    r->NW = r->N = n;
    r->T = time;
    r->Y = y;
    r->F = f;

    r->TOL = tol;
    if(r->C) Free(r->C);
    r->C = Malloc(24 * sizeof(double));
#if NED
    r->IND = 2;
    r->C[0] = 1.0;
#else
    r->IND = 1;	/* no options */
#endif
    if(r->W) Free(r->W);
    r->W = Malloc(r->NW * 9 * sizeof(double));

    return r;
}

void DverkFree(DVERK *d)
{
#if NED
    printf("# steps = %g\n", d->C[23]);
#endif
    Free(d->C);
    Free(d->W);
    Free(d);
}


/* The function called by DVERK to compute derivatives.
 *
 * If something is horribly wrong, set (*NN) to zero, and
 * DVERK will notify the calling program.
 */


static DVERK *gr;

static void FORT_F(int *NN, double *TT, double *YY, double *YP)
{
    (*(gr->F))( *NN, *TT, YY, YP);
}


double DverkIntegrate(DVERK *r, double TOUT) /* returns actual tout */
{
    int direction = TOUT - r->T > 0 ? 1 : -1;
    extern dverk_();

    gr = r;

    while(direction*(r->T - TOUT) < 0)
    {
	dverk_(&r->N, FORT_F, &r->T, r->Y, &TOUT, &r->TOL, &r->IND, r->C,
	    &r->NW, r->W);

	switch(r->IND)
	{
	case 3:           /* successful integration */
	    assert(r->T == TOUT);
	    break;
	case 4:
	case 5:
	case 6:
	    Fatal("DVERK thinks it's handling interrupts");
	    break;
	case -3:          /* EPS too small */
	    Warning("DVERK step size too small: increasing TOL to %g",
		r->TOL*=10);
	    r->IND = 1;
	    break;
	case -2:
	    Warning("DverkIntegrate: stepsize too big? increasing TOL to %g",
		r->TOL*=10);
	    r->IND = 1;
	    break;
	case -1:          /* too many steps */
	    Warning("DVERK working hard, perhaps wrong mf");
	    r->IND = 3;    /* Just continue */
	    break;
	default:
	    Fatal("AAAAHHGHHGHRHGHH!  DVERK returned bogus IND!");
	    break;
	}
    }
    return r->T;
}
