#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include "misc.h"
#include "lsode.h"

/*
LSODE *LsodeAlloc(int n, double t, double *y, F_EVAL f, int stiff_flag,
    double eps, double zero);
double LsodeIntegrate(LSODE *l, double TOUT);
void LsodeFree(LSODE *l);
*/

#define k 1
#define mu 0.1
#define mass 1

void fric(int n, double t, double *y, double *Ydot)
{
    assert(n==3);
    Ydot[0] = y[1];
    Ydot[1] = -k*y[0]/mass - mu*y[1];
    Ydot[2] = mu*SQR(y[1]);
}

int main(void)
{
    double y[3], t = 0;
    LSODE *l;
    y[0] = -1;	/* x */
    y[1] = 0;	/* v */
    y[2] = 0;	/* energy lost due to friction */

    l = LsodeAlloc(3, 0.0, y, fric, 0, 1e-10, 0.0);

    while(1)
    {
	double KE, PE, FE, TE;
	t = LsodeIntegrate(l, t+0.1);
	KE = 0.5*mass*SQR(y[1]);
	PE = 0.5*k*SQR(y[0]);
	FE = y[2];
	TE = KE+PE+FE;
	printf("t=%g x=%g v=%g KE=%g PE=%g FE=%g TE=%g\n",
	    t, y[0], y[1], KE, PE, FE, TE);
    }
    return 0;
}
