def interpolate_lagrange_3D(x, y, z, x_nodes, y_nodes, z_nodes, field_nodes):
    """Interpolate values using 3D Lagrangian polynomials

    The order of the polynomial interpolation in each direction depends on the
    number of field nodes passed to this function

    Parameters
    ----------
    x, y, z : (N,) array_like
        Coordinates used to interpolate
    x_nodes, y_nodes, z_nodes, field_nodes : array_like
        Nodal coordinates and field to be interpolated at corresponding nodes

    Returns
    -------
    field : (N,) array_like
        Interpolated field

    """
    field = 0
    for m, (xm, ym, zm, theta) in enumerate(zip(x_nodes, y_nodes, z_nodes,
        field_nodes)):
        tmp = 1
        for k, xk in enumerate(x_nodes):
            if k != m and xm - xk != 0:
                tmp *= (x - xk) / (xm - xk)
        for p, yp in enumerate(y_nodes):
            if p != m and ym - yp != 0:
                tmp *= (y - yp) / (ym - yp)
        for p, zp in enumerate(z_nodes):
            if p != m and zm - zp != 0:
                tmp *= (z - zp) / (zm - zp)
        field += tmp * theta
    return field
