/*
** This is the generic integrator tester.  To test a new integrator, the
** following global substitutions must be made with an editor:

:%s/bsode/inc_name/g
:%s/Bsode/fcn_prefix/g
:%s/BSODE/object_name/g

*/
#include "bsode.h"

/* Make this "long double" for QUAD routines, otherwise "double", and
 * revise TOL appropriately (1e-15 for double, smaller for long double).
 */
#define REAL double
#define TOL 1e-16

#include "misc.h"
#include <assert.h>
#include <stdio.h>
#include <math.h>

#define ONLY_1 0	/* only do one ODE? */
#define TRIES 1
#define GO_BACK 0
#define FIND_ROOTS 0
#define COMPUTE_OBSERVED_ORDER 0

/*
** A test to ensure the BSODE routines can do multiple integrations at once.
*/
void f_y(int n, REAL t, REAL y[n], REAL ydot[n])
{
    assert(n==1);
    ydot[0] = y[0];
}

void f_y2(int n, REAL t, REAL y[n], REAL ydot[n])
{
    assert(n==1);
    ydot[0] = y[0]*y[0];
}


#if 0
/* Kepler problem */
void f_k(int n, REAL t, REAL y[n], REAL ydot[n])
{
    double ecc,r;
    assert(n==4);
    data ecc / 0.5d0 /
    integer n
    if(n .ne. 4) then
    print*,'subroutine D3: n must be 4, not ',n
    stop
    endif
    yp(1)=y(3)
    yp(2)=y(4)
    r=dsqrt(y(1)**2+y(2)**2)
    yp(3)=-y(1)/r**3
    yp(4)=-y(2)/r**3
}
#endif

void f_y4(int n, REAL t, REAL y[n], REAL ydot[n])
{
    assert(n==1);
    ydot[0] = y[0]/4 * (1-y[0]/20);
}

void f_y3_2(int n, REAL t, REAL y[n], REAL ydot[n])
{
    assert(n==1);
    ydot[0] = -y[0]*y[0]*y[0]/2;
}

void f_1_y(int n, REAL t, REAL y[n], REAL ydot[n])
{
    assert(n==1);
    ydot[0] = 1/y[0];
}

void f_one(int n, REAL t, REAL y[n], REAL ydot[n])
{
    assert(n==1);
    ydot[0] = 1;
}
REAL g_one(int n, REAL t, REAL y[n], int r)
{
    assert(n==1);
    assert(r==0);
    return y[0];
}

void f_two(int n, REAL t, REAL y[n], REAL ydot[n])
{
    assert(n==1);
    ydot[0] = 2;
}
REAL g_two(int n, REAL t, REAL y[n], int r)
{
    assert(n==1);
    assert(r==0);
    return y[0];
}

void f_t(int n, REAL t, REAL y[n], REAL ydot[n])
{
    assert(n==1);
    ydot[0] = t;
}
REAL g_t(int n, REAL t, REAL y[n], int r)
{
    assert(n==1);
    assert(r==0);
    return y[0];
}

void f_t2(int n, REAL t, REAL y[n], REAL ydot[n])
{
    assert(n==1);
    ydot[0] = t*t;
}
REAL g_t2(int n, REAL t, REAL y[n], int r)
{
    assert(n==1);
    assert(r==0);
    return y[0];
}

int main(void)
{
    int i, try;
#if COMPUTE_OBSERVED_ORDER
    REAL error[512];
#endif
    REAL approx, exact;
for(try=0; try<TRIES; try++)
{
    REAL eps = TOL, t;
    REAL y_y=1;
#if !ONLY_1
    REAL y_y2_0=0.1, y_y2=y_y2_0;
    REAL y_y4_0=1, y_y4=y_y4_0;
    REAL y_y3_2_0=1, y_y3_2=y_y3_2_0;
    REAL y_1_y_0=1, y_1_y=y_1_y_0;
    REAL y_one=0;
    REAL y_two=0;
    REAL y_t=0;
    REAL y_t2=0;
#endif
    BSODE *l_y
#if !ONLY_1
    , *l_y2, *l_y4, *l_y3_2, *l_1_y, *l_one, *l_two, *l_t, *l_t2
#endif
    ;

    l_y   = BsodeAlloc(1, 0., &y_y,   f_y,   0, eps, 0.);
#if !ONLY_1
    l_y2  = BsodeAlloc(1, 0., &y_y2,  f_y2,  0, eps, 0.);
    l_y4  = BsodeAlloc(1, 0., &y_y4,  f_y4,  0, eps, 0.);
    l_y3_2= BsodeAlloc(1, 0., &y_y3_2,f_y3_2,0, eps, 0.);
    l_1_y = BsodeAlloc(1, 0., &y_1_y, f_1_y, 0, eps, 0.);
    l_one = BsodeAlloc(1, 0., &y_one, f_one, 0, eps, 0.);
    l_two = BsodeAlloc(1, 0., &y_two, f_two, 0, eps, 0.);
    l_t   = BsodeAlloc(1, 0., &y_t,   f_t,   0, eps, 0.);
    l_t2  = BsodeAlloc(1, 0., &y_t2,  f_t2,  0, eps, 0.);
#endif

    for(i=1; i<=10; i+=1)
    {
	t = 0.1*i;
	assert(t == BsodeIntegrate(l_y, t));
#if !ONLY_1
	assert(t == BsodeIntegrate(l_y2, t));
	assert(t == BsodeIntegrate(l_y4, t));
	assert(t == BsodeIntegrate(l_y3_2, t));
	assert(t == BsodeIntegrate(l_1_y, t));
	assert(t == BsodeIntegrate(l_one, t));
	assert(t == BsodeIntegrate(l_two, t));
	assert(t == BsodeIntegrate(l_t, t));
	assert(t == BsodeIntegrate(l_t2, t));
#endif
    }
#if GO_BACK
    for(i=1; i>=0; i-=1)
    {
	t = 1.0*i;
	assert(t == BsodeIntegrate(l_y, t));
#if !ONLY_1
	assert(t == BsodeIntegrate(l_y2, t));
	assert(t == BsodeIntegrate(l_y4, t));
	assert(t == BsodeIntegrate(l_y3_2, t));
	assert(t == BsodeIntegrate(l_1_y, t));
	assert(t == BsodeIntegrate(l_one, t));
	assert(t == BsodeIntegrate(l_two, t));
	assert(t == BsodeIntegrate(l_t, t));
	assert(t == BsodeIntegrate(l_t2, t));
#endif /*GO_BACK*/
    }
#endif

    BsodeFree(l_y);
#if !ONLY_1
    BsodeFree(l_y2);
    BsodeFree(l_y4);
    BsodeFree(l_y3_2);
    BsodeFree(l_1_y);
    BsodeFree(l_one);
    BsodeFree(l_two);
    BsodeFree(l_t);
    BsodeFree(l_t2);
#endif

    exact=exp(t); approx=y_y;
    printf("y: %.16g = %.16g (%.8g)\n", (double)exact, (double)approx, (double)(approx/exact-1));

#if !ONLY_1
    exact=1/(1/y_y2_0 - t); approx=y_y2;
    printf("y2:  %.16g = %.16g (%.8g)\n", (double)exact, (double)approx, (double)(approx/exact-1));

    exact=(20/(1+19*exp(-t/4))); approx=y_y4;
    printf("y4:  %.16g = %.16g (%.8g)\n", (double)exact, (double)approx, (double)(approx/exact-1));

    exact=(1/sqrt(1+t)); approx=y_y3_2;
    printf("y3_2:%.16g = %.16g (%.8g)\n", (double)exact, (double)approx, (double)(approx/exact-1));

    exact=sqrt(2*t + y_1_y_0); approx=y_1_y;
    printf("1/y: %.16g = %.16g (%.8g)\n", (double)exact, (double)approx, (double)(approx/exact-1));

    exact=t; approx=y_one;
    printf("one: %.16g = %.16g (%.8g)\n", (double)exact, (double)approx, (double)(approx/exact-1));

    exact=(2*t); approx=y_two;
    printf("two: %.16g = %.16g (%.8g)\n", (double)exact, (double)approx, (double)(approx/exact-1));

    exact=(.5*t*t); approx=y_t;
    printf("t: %.16g = %.16g (%.8g)\n", (double)exact, (double)approx, (double)(approx/exact-1));

    exact=(t*t*t/3); approx=y_t2;
    printf("t^2: %.16g = %.16g (%.8g)\n", (double)exact, (double)approx, (double)(approx/exact-1));
#endif

#if FIND_ROOTS
    y_one=0;
    y_two=0;
    y_t2=0;
    l_one = BsodeAlloc(1, 0., &y_one, f_one, 0, eps, 0.);
    l_two = BsodeAlloc(1, 0., &y_two, f_two, 0, eps, 0.);
    l_t2  = BsodeAlloc(1, 0., &y_t2,  f_t2,  0, eps, 0.);

    printf("Out to 1, then back looking for y=0\n");

    assert(1.0 == BsodeIntegrate(l_one, 1.0));
    assert(1.0 == BsodeIntegrate(l_two, 1.0));
    assert(1.0 == BsodeIntegrate(l_t2,  1.0));

    /* Now go back looking for 0 */
    BsodePrepRoot(l_one, 1, g_one);
    BsodePrepRoot(l_two, 1, g_two);
    BsodePrepRoot(l_t2,  1, g_t2);

    t = BsodeIntegrate(l_one, -1.);
    printf("one: 0 = %.16g (%.16g)\n", (double)y_one, (double)t);
    t = BsodeIntegrate(l_two, -1.);
    printf("two: 0 = %.16g (%.16g)\n", (double)y_two, (double)t);
    t = BsodeIntegrate(l_t2,  -1.);
    printf("t^2: 0 = %.16g (%.16g)\n", (double)y_t2,  (double)t);

    assert(BsodeWhichRoot(l_one) == 0);
    assert(BsodeWhichRoot(l_two) == 0);
    assert(BsodeWhichRoot(l_t2) == 0);

    BsodeFree(l_one);
    BsodeFree(l_two);
    BsodeFree(l_t2);
#endif	/* FIND_ROOTS */

#if COMPUTE_OBSERVED_ORDER
    eps = 1;
    for(i=0; i<8; i++)
    {
	double order=0.0, base=100;
	y_y=1;
	l_y = BsodeAlloc(1, 0., &y_y, f_y, 0, eps, 0.);
	t = 1;
	assert(t == BsodeIntegrate(l_y, t));
	BsodeFree(l_y);
	error[i] = fabs(y_y - exp(t));
	if(i > 0)
	    order = log(error[i-1]/error[i])/log(base);
	printf("eps=%g y=%.16g err=%g err/eps=%g order=%g\n", (double)eps, y_y, (double)error[i], (double)(error[i]/eps), order);
	eps /= base;
    }
    
#endif	/* COMPUTE_OBSERVED_ORDER */
}
    return 0;
}
