interp1Dexpl.c

Go to the documentation of this file.

#include <stdlib.h>
#include <math.h>
#include "../ascot5.h"
#include "../math.h"
#include "interp.h"
#include "spline.h"

int interp1Dexpl_init_coeff(real* c, real* f, int n_x, int bc_x,
                            real x_min, real x_max) {
 
    if(c == NULL) {
        return 1;
    }
 
    /* Calculate cubic spline coefficients. For each grid cell i_x, there are
       four coefficients. Note how we account for normalized grid. */
 
    /* Cubic spline along x, using f values to get a total of four
       coefficients */
    splineexpl(f, n_x, bc_x, c);
 
    return 0;
}

void interp1Dexpl_init_spline(interp1D_data* str, real* c,
                              int n_x, int bc_x, real x_min, real x_max) {
 
    /* Calculate grid interval. For periodic boundary condition, grid maximum
       value and the last data point are not the same. Take this into account
       in grid interval. */
    real x_grid = (x_max - x_min) / ( n_x - 1 * (bc_x == NATURALBC) );
 
    /* Initialize the interp1D_data struct */
    str->n_x    = n_x;
    str->bc_x   = bc_x;
    str->x_min  = x_min;
    str->x_max  = x_max;
    str->x_grid = x_grid;
    str->c      = c;
}

a5err interp1Dexpl_eval_f(real* f, interp1D_data* str, real x) {
 
    /* Make sure periodic coordinates are within [min, max] region. */
    if(str->bc_x == PERIODICBC) {
        x = fmod(x - str->x_min, str->x_max - str->x_min) + str->x_min;
        x = x + (x < str->x_min) * (str->x_max - str->x_min);
    }
 
    /* Index for x variable. The -1 needed at exactly grid end. */
    int i_x = (x-str->x_min)/str->x_grid - 1*(x==str->x_max);
    /* Normalized x coordinate in current cell */
    real dx = (x-(str->x_min+i_x*str->x_grid))/str->x_grid;
    /* Helper variables */
    real dx2 = dx*dx;
    real dx3 = dx2*dx;
 
    int n = i_x*4; /* Index jump to cell */
 
    int err = 0;
 
    /* Check that the coordinate is within the grid. */
    if( str->bc_x == NATURALBC && !(x >= str->x_min && x <= str->x_max) ) {
        err = 1;
    }
 
    if(!err) {
        *f = str->c[n+0]+dx*str->c[n+1]+dx2*str->c[n+2]+dx3*str->c[n+3];
    }
 
    return err;
}

a5err interp1Dexpl_eval_df(real* f_df, interp1D_data* str, real x) {
 
    /* Make sure periodic coordinates are within [min, max] region. */
    if(str->bc_x == PERIODICBC) {
        x = fmod(x - str->x_min, str->x_max - str->x_min) + str->x_min;
        x = x + (x < str->x_min) * (str->x_max - str->x_min);
    }
 
    /* Index for x variable. The -1 needed at exactly grid end. */
    int i_x = (x-str->x_min)/str->x_grid - 1*(x==str->x_max);
    /* Normalized x coordinate in current cell */
    real dx = (x-(str->x_min+i_x*str->x_grid))/str->x_grid;
    /* Helper variables */
    real dx2 = dx*dx;
    real dx3 = dx2*dx;
    real xgi = 1.0/str->x_grid;
 
    int n = i_x*4; /* Index jump to cell */
 
    int err = 0;
 
    /* Check that the coordinate is within the grid. */
    if( str->bc_x == NATURALBC && !(x >= str->x_min && x <= str->x_max) ) {
        err = 1;
    }
 
    if(!err) {
        /* f */
        f_df[0] = str->c[n+0]+dx*str->c[n+1]+dx2*str->c[n+2]+dx3*str->c[n+3];
 
        /* df/dx */
        f_df[1] = xgi*(str->c[n+1]+2*dx*str->c[n+2]+3*dx2*str->c[n+3]);
 
        /* d2f/dx2 */
        f_df[2] = xgi*xgi*(2*str->c[n+2]+6*dx*str->c[n+3]);
    }
 
    return err;
}