template <typename T>
class Grid {
public:
    int m, n;         // number of grid nodes
    vec2 x0;          // world-space position of node (0, 0)
    float dx;         // grid node spacing
    vector<T> values; // grid node values
    Grid() {}
    Grid(int m, int n, vec2 x0, float dx):
        m(m), n(n), x0(x0), dx(dx), values(m*n) {}
    // value at node (i, j)
    const T& get(int i, int j) const;
    T& operator()(int i, int j) {
        return (T&)get(i, j);
    }
    const T& operator()(int i, int j) const {
        return get(i, j);
    }
    // assign all values to a constant
    void assign(T value);
    // interpolated value at world-space position x
    T interpolate(vec2 x) const;
};

template <typename T>
const T& Grid<T>::get(int i, int j) const {
    return values[i+m*j];
}

template <typename T>
void Grid<T>::assign(T value) {
    values.assign(values.size(), value);
}

template <typename T>
T lerp(float x, T y0, T y1) { // linear interpolation between y0 and y1
    return y0 + (y1 - y0)*x;
}

template <typename T>
T Grid<T>::interpolate(vec2 x) const {
    vec2 index = (x - x0)/dx;
    int i = floor(index(0)), j = floor(index(1)); // integer part of index
    float fi = index(0) - i, fj = index(1) - j;   // fractional part
    // clamp to interior of grid if at or outside boundary
    if (i < 0)
        {i = 0; fi = 0;}
    else if (i >= m-1)
        {i = m-2; fi = 1;}
    if (j < 0)
        {j = 0; fj = 0;}
    else if (j >= n-1)
        {j = n-2; fj = 1;}
    return lerp(fj, lerp(fi, get(i, j),   get(i+1, j)),
                    lerp(fi, get(i, j+1), get(i+1, j+1)));
}