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

#include <math.h>
#include <assert.h>
#include <malloc.h>
#include "misc.h"
#include "iqrkd78.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

IQRKD78 *Iqrkd78Alloc(int n, double t0, double *y, F_EVAL f, int stiff_flag,
    double tol, double zero)
{
    IQRKD78 *r = Calloc(1, sizeof(IQRKD78));
    int i;

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

    r->NW = r->N = n;
    r->QT = r->T = t0;
    r->Y = y;
    r->QY = Calloc(n, sizeof(long double));
    for(i=0; i<n; i++)
	r->QY[i] = y[i];
    r->F = f;

    r->QTOL = 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 Iqrkd78Free(IQRKD78 *r)
{
    Free(r->C);
    Free(r->W);
    Free(r->QY);
    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 IQRKD78 *_gr;

static void QFORT_F(int *NN, long double *TT, long double *YY, long double *YP)
{
    int i;
    double tt, yy[*NN], yp[*NN];

    /* convert internal QUADs to doubles before calling user's F */
    tt = *TT;
    for(i=0; i<*NN; i++)
	yy[i] = YY[i];

    (*(_gr->F))(*NN, tt, yy, yp);

    /* convert answer back to QUAD */
    for(i=0; i<*NN; i++)
	YP[i] = yp[i];
}


double Iqrkd78Integrate(IQRKD78 *r, double TOUT) /* returns actual tout */
{
    int direction = TOUT - r->T > 0 ? 1 : -1, i;
    long double QTOUT = TOUT;
    extern void qrdmeth78_();

     _gr = r;

    while(direction*(r->QT - QTOUT) < 0)
    {
	qrdmeth78_(&r->N, QFORT_F, &r->QT, r->QY, &QTOUT, &r->QTOL, &r->IND, r->C,
	    &r->NW, r->W);

	switch(r->IND)
	{
	case 3:           /* successful integration */
	    assert(r->QT == QTOUT);
	    assert((long double)TOUT == QTOUT);
	    r->T = r->QT;
	    break;
	case 4:
	case 5:
	case 6:
	    Fatal("QRDMETH78 thinks it's handling interrupts");
	    break;
	case -3:          /* EPS too small */
	    r->TOL*=10;
	    r->QTOL = r->TOL;
	    fprintf(stderr, "QRDMETH78 step size too small: increasing TOL to %g\n", r->TOL);
	    r->IND = IND_ENTRY;
	    break;
	case -2:
	    r->TOL*=10;
	    r->QTOL = r->TOL;
	    fprintf(stderr, "Qrkd78Integrate: stepsize too big? increasing TOL to %g\n", r->TOL);
	    r->IND = IND_ENTRY;
	    break;
	case -1:          /* too many steps */
	    fprintf(stderr, "QRDMETH78 working hard, perhaps wrong mf");
	    r->IND = 3;    /* Just continue */
	    break;
	default:
	    Fatal("AAAAHHGHHGHRHGHH!  QRDMETH78 returned bogus IND!");
	    break;
	}
    }
    for(i=0; i<r->N; i++)
	r->Y[i] = r->QY[i];
    return r->T;
}
