linear_1d_interpolation.hpp

Go to the documentation of this file.

#ifndef LINEAR_1D_INTERPOLATION_HPP
#define LINEAR_1D_INTERPOLATION_HPP

#include "space_settings.hpp"


/* Simple 1D linear interpolation
 * It is assumed that x is sorted in ascending order
 * @param(in) x  View of input x values
 * @param(in) y  View of input y(x) values
 * @param(in) x1  Single x value
 * @return T2 Single y value
 */
template<typename T, typename T2>
KOKKOS_INLINE_FUNCTION T2 linear_1d_interpolation(const T& x, const T& y, const T2 x1){
    const int n = x.size();
    if (x1 >= x(0) && x1 <= x(n-1)){
        // use bisection to find points in x closest to x1
        int is = 0;
        int ie = n-1;

        while(ie-is > 1){
            int ip = (is+ie)/2;
            if (x1 > x(ip)){
                is=ip;
            }else{
                ie=ip;
            }
        }

        double slope = (y(ie)-y(is))/(x(ie)-x(is));
        return (y(is) + slope * (x1-x(is)));
    }else{
        // Extrapolate if outside range
        if (x1 < x(0)){
            double slope = (y(1)-y(0)) / (x(1)-x(0));
            return (y(0) + slope * (x1-x(0)));
        }else if (x1 > x(n-1)){
            double slope = (y(n-1)-y(n-2)) / (x(n-1)-x(n-2));
            return (y(n-1) + slope * (x1-x(n-1)));
        }else{
            return 0.0;
            // NaN?
        }
    }
}


/* Simple 1D linear interpolation
 * It is assumed that x is sorted in ascending order
 * @param(in) x  View of input x values
 * @param(in) y  View of input y(x) values
 * @param(in) x1  View of output x values
 * @param(out) y1  View of output y(x1) values
 */
template<typename T>
void linear_1d_interpolation(const T& x, const T& y, const T& x1, const T& y1){
    Kokkos::parallel_for("interpol_1d", Kokkos::RangePolicy<typename T::execution_space>( 0,x1.size() ), KOKKOS_LAMBDA( const int i){
        y1(i) = linear_1d_interpolation(x, y, x1(i));
    });
    Kokkos::fence();
}

#endif