/* implementation of DE-filling */

#include <math.h>
#include <float.h>
#include "diffeq.h"

/* exact function we are approximating */
static double fexact(double x, int K)
{
    return exp(-(K+1)*x) + exp(-x);
}

/*
 * approximation of f' using parameter K and
 * approximation of f(x) (fx)
 */
static double fprime(double x, double fx, int K)
{
    return -fx - K*(fx - exp(-x));
}

/*
 * approximation of function f, using the approximation of of f
 * at the previous step (oldf), and the approximation of f'
 * at the previous step (oldfprime).
 */
static double nextf(double oldf, double oldfprime)
{
    return oldf + STEPSIZE * oldfprime;
}


void delist(double initx, double init_f_of_x, int K, struct deline *detable)
{
    int i;

    /* initial conditions */
    detable[0].x = initx;
    detable[0].f_x = init_f_of_x;
    detable[0].f_of_x = fexact(detable[0].x, K);
    detable[0].error = (detable[0].f_x - detable[0].f_of_x)
                       / detable[0].f_of_x;

    /* remaining entries */
    for (i = 1; i < TABLELINES; i++) {
        detable[i].x = detable[i-1].x + STEPSIZE;
        detable[i].f_x = nextf(detable[i-1].f_x,
                               fprime(detable[i-1].x, detable[i-1].f_x, K));
        detable[i].f_of_x = fexact(detable[i].x, K);
        if (detable[i].f_of_x == 0 && detable[i].f_x == 0)
            detable[i].error = 0;
        else if (detable[i].f_of_x == 0)
            detable[i].error = FLT_MAX;
        else
            detable[i].error = (detable[i].f_x - detable[i].f_of_x)
                               / detable[i].f_of_x;
    }
}