/* integrate.c
 * Integrate using the FORTRAN routine QRDMETH78.
 * (This is just a quad precision version of rkd78/rdmeth78)
 */

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

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

QRKD78 *Qrkd78Alloc(int n, long double time, long double *y, QF_EVAL f, int stiff_flag,
    long double tol, long double zero)
{
    QRKD78 *r = Calloc(1, sizeof(QRKD78));

    if(stiff_flag != 0 || zero != 0.0)
	Fatal("Qrkd78Alloc: 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(long double));
    r->W = Malloc(r->NW * 23 * sizeof(long double));

    return r;
}

void Qrkd78Free(QRKD78 *r)
{
    Free(r->C);
    Free(r->W);
    Free(r);
}


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


static QRKD78 *_gr;

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


long double Qrkd78Integrate(QRKD78 *r, long double TOUT) /* returns actual tout */
{
    int direction = TOUT - r->T > 0 ? 1 : -1;
    extern void qrdmeth_();

     _gr = r;

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