def lininterp_3d(x_grid, y_grid, z_grid, vals, s):
    """
    Fast 3D interpolation.  Uses linear extrapolation for points outside the
    grid.  Note that the grid must be regular (i.e., evenly spaced). 

    Based on linear interpolation code written by @albop.

    Parameters
    ----------

    x_grid: np.ndarray
        grid points for x, one dimensional regular grid

    y_grid: np.ndarray
        grid points for y, one dimensional regular grid

    z_grid: np.ndarray
        grid points for z, one dimensional regular grid

    vals: np.ndarray
        vals[i, j, k] = f(x[i], y[j], z[k])
     
    s: np.ndarray
        3D point at which to evaluate function

    """


    d = 3
    smin = (x_grid[0], y_grid[0], z_grid[0])
    smax = (x_grid[-1], y_grid[-1], z_grid[-1])

    order_0 = len(x_grid)
    order_1 = len(y_grid)
    order_2 = len(z_grid)

    # (s_1, ..., s_d) : evaluation point
    s_0 = s[0]
    s_1 = s[1]
    s_2 = s[2]

    # normalized evaluation point (in [0,1] inside the grid)
    s_0 = (s_0-smin[0])/(smax[0]-smin[0])
    s_1 = (s_1-smin[1])/(smax[1]-smin[1])
    s_2 = (s_2-smin[2])/(smax[2]-smin[2])

    # q_k : index of the interval "containing" s_k
    q_0 = max( min( int(s_0 *(order_0-1)), (order_0-2) ), 0 )
    q_1 = max( min( int(s_1 *(order_1-1)), (order_1-2) ), 0 )
    q_2 = max( min( int(s_2 *(order_2-1)), (order_2-2) ), 0 )

    # lam_k : barycentric coordinate in interval k
    lam_0 = s_0*(order_0-1) - q_0
    lam_1 = s_1*(order_1-1) - q_1
    lam_2 = s_2*(order_2-1) - q_2

    # v_ij: values on vertices of hypercube "containing" the point
    v_000 = vals[(q_0), (q_1), (q_2)]
    v_001 = vals[(q_0), (q_1), (q_2+1)]
    v_010 = vals[(q_0), (q_1+1), (q_2)]
    v_011 = vals[(q_0), (q_1+1), (q_2+1)]
    v_100 = vals[(q_0+1), (q_1), (q_2)]
    v_101 = vals[(q_0+1), (q_1), (q_2+1)]
    v_110 = vals[(q_0+1), (q_1+1), (q_2)]
    v_111 = vals[(q_0+1), (q_1+1), (q_2+1)]

    # interpolated/extrapolated value
    output = (1-lam_0)*((1-lam_1)*((1-lam_2)*(v_000) + (lam_2)*(v_001)) + (lam_1)*((1-lam_2)*(v_010) + (lam_2)*(v_011))) + (lam_0)*((1-lam_1)*((1-lam_2)*(v_100) + (lam_2)*(v_101)) + (lam_1)*((1-lam_2)*(v_110) + (lam_2)*(v_111)))
    return output