/* If the timestep changes by a factor "base" and the
 * errors are e1 and e2, respectively, then the order
 * is log(e1/e2)/log(base)
 */

#include "misc.h"
#include "variable_leapfrog.h"
#include "lsode.h"
#include "matvec.h"
#include <assert.h>
#include <math.h>

#define NUM_H 20

/* Set T_FINAL to DT to measure local error, or a constant to measure global */
#define T_FINAL 1.0

void pendulum(int n, double t, double *y, double *yp)
{
    assert(n==2);
    yp[0] = y[1];
    yp[1] = -y[0];
}

void pendulum2(int n, double t, double *r, double *rpp)
{
    assert(n==1);
    rpp[0] = -r[0];
}

static double DT;

double ComputeDT(int n, double T, double current_DT, double *R, double *V)
{
#if 0
    static int half = 0;
    extern double drand48(void);
    if(drand48() < 0.00)
	half = !half;
    if(half)
	return DT/2;
    else
#endif
	return DT;
}

int main(void)
{
    double base=2.0, ytrue[2], yleap[2], error[NUM_H];
    int i;
    LSODE *lsode;
    
    srand48(time(NULL)+getpid());

    for(i=0; i<NUM_H; i++)
    {
	double diff[2];

	DT = IntPow(base, -i);

	ytrue[0] = 1.0;
	ytrue[1] = 0.0;
	lsode = LsodeAlloc(2, 0.0, ytrue, pendulum, 0, 4e-16, 0.0);
	assert(T_FINAL ==  LsodeIntegrate(lsode, T_FINAL));
	LsodeFree(lsode);

	yleap[0] = 1.0;
	yleap[1] = 0.0;
	init_variable_leapfrog(1, 0.0, ComputeDT, yleap, yleap+1, pendulum2);
	assert(T_FINAL ==  integrate_variable_leapfrog(T_FINAL));

	VecDiff(2, diff, ytrue, yleap);
	error[i] = VecNorm1(2, diff);
	if(i > 0)
	{
	    double order = log(error[i-1]/error[i])/log(base);
	    printf("%2d %g %g %g\n", i, DT, error[i], order);
	}
    }

    return 0;
}
