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  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){
    int n = x.size();
    int n1 = x1.size();

    Kokkos::parallel_for("interpol_1d", Kokkos::RangePolicy<typename T::execution_space>( 0,n1 ), KOKKOS_LAMBDA( const int i){
        if (x1(i) >= x(0) && x1(i) <= x(n-1)){
            // use bisection to find points in x closest to x1(i)
            int is = 0;
            int ie = n-1;

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

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

#endif