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

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

/* init_integrator()
 *
 * We must pass a whole bunch of stuff to RKD78, 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.
 */

#define IND_ENTRY 1

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

    if(stiff_flag != 0 || zero != 0.0)
	Fatal("Rkd78Alloc: 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; /*MAX(tol, 1e-10);*/
    r->IND = IND_ENTRY;
    r->C = Malloc(24 * sizeof(double));
    r->W = Malloc(r->NW * 23 * sizeof(double));

    return r;
}

void Rkd78Free(RKD78 *r)
{
    Free(r->C);
    Free(r->W);
    Free(r);
}


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


static RKD78 *_gr;

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


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

     _gr = r;

    while(direction*(r->T - TOUT) < 0)
    {
	rdmeth78_(&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("RDMETH78 thinks it's handling interrupts");
	    break;
	case -3:          /* EPS too small */
	    fprintf(stderr, "RDMETH78 step size too small: increasing TOL to %g\n",
		r->TOL*=10);
	    r->IND = IND_ENTRY;
	    break;
	case -2:
	    fprintf(stderr, "Rkd78Integrate: HMIN too big, increasing TOL to %g\n",
		r->TOL*=10);
	    r->IND = IND_ENTRY;
	    break;
	case -1:          /* too many steps */
	    fprintf(stderr, "RDMETH78 working hard, perhaps wrong mf");
	    r->IND = 3;    /* Just continue */
	    break;
	default:
	    Fatal("AAAAHHGHHGHRHGHH!  RDMETH78 returned bogus IND!");
	    break;
	}
    }
    return r->T;
}
