def _evaluate(self, w, x, y, z):
    '\n        Returns the level of the interpolated function at each value in x,y,z.\n        Only called internally by HARKinterpolator4D.__call__ (etc).\n        '
    if _isscalar(w):
        w_pos = max(min(self.wSearchFunc(self.w_list, w), (self.w_n - 1)), 1)
        x_pos = max(min(self.xSearchFunc(self.x_list, x), (self.x_n - 1)), 1)
        y_pos = max(min(self.ySearchFunc(self.y_list, y), (self.y_n - 1)), 1)
        z_pos = max(min(self.zSearchFunc(self.z_list, z), (self.z_n - 1)), 1)
    else:
        w_pos = self.wSearchFunc(self.w_list, w)
        w_pos[(w_pos < 1)] = 1
        w_pos[(w_pos > (self.w_n - 1))] = (self.w_n - 1)
        x_pos = self.xSearchFunc(self.x_list, x)
        x_pos[(x_pos < 1)] = 1
        x_pos[(x_pos > (self.x_n - 1))] = (self.x_n - 1)
        y_pos = self.ySearchFunc(self.y_list, y)
        y_pos[(y_pos < 1)] = 1
        y_pos[(y_pos > (self.y_n - 1))] = (self.y_n - 1)
        z_pos = self.zSearchFunc(self.z_list, z)
        z_pos[(z_pos < 1)] = 1
        z_pos[(z_pos > (self.z_n - 1))] = (self.z_n - 1)
    i = w_pos
    j = x_pos
    k = y_pos
    l = z_pos
    alpha = ((w - self.w_list[(i - 1)]) / (self.w_list[i] - self.w_list[(i - 1)]))
    beta = ((x - self.x_list[(j - 1)]) / (self.x_list[j] - self.x_list[(j - 1)]))
    gamma = ((y - self.y_list[(k - 1)]) / (self.y_list[k] - self.y_list[(k - 1)]))
    delta = ((z - self.z_list[(l - 1)]) / (self.z_list[l] - self.z_list[(l - 1)]))
    f = (((1 - alpha) * (((1 - beta) * ((((((1 - gamma) * (1 - delta)) * self.f_values[((i - 1), (j - 1), (k - 1), (l - 1))]) + (((1 - gamma) * delta) * self.f_values[((i - 1), (j - 1), (k - 1), l)])) + ((gamma * (1 - delta)) * self.f_values[((i - 1), (j - 1), k, (l - 1))])) + ((gamma * delta) * self.f_values[((i - 1), (j - 1), k, l)]))) + (beta * ((((((1 - gamma) * (1 - delta)) * self.f_values[((i - 1), j, (k - 1), (l - 1))]) + (((1 - gamma) * delta) * self.f_values[((i - 1), j, (k - 1), l)])) + ((gamma * (1 - delta)) * self.f_values[((i - 1), j, k, (l - 1))])) + ((gamma * delta) * self.f_values[((i - 1), j, k, l)]))))) + (alpha * (((1 - beta) * ((((((1 - gamma) * (1 - delta)) * self.f_values[(i, (j - 1), (k - 1), (l - 1))]) + (((1 - gamma) * delta) * self.f_values[(i, (j - 1), (k - 1), l)])) + ((gamma * (1 - delta)) * self.f_values[(i, (j - 1), k, (l - 1))])) + ((gamma * delta) * self.f_values[(i, (j - 1), k, l)]))) + (beta * ((((((1 - gamma) * (1 - delta)) * self.f_values[(i, j, (k - 1), (l - 1))]) + (((1 - gamma) * delta) * self.f_values[(i, j, (k - 1), l)])) + ((gamma * (1 - delta)) * self.f_values[(i, j, k, (l - 1))])) + ((gamma * delta) * self.f_values[(i, j, k, l)]))))))
    return f
Returns the level of the interpolated function at each value in x,y,z.
Only called internally by HARKinterpolator4D.__call__ (etc).
def _evaluate(self, w, x, y, z):
    '\n        Returns the level of the interpolated function at each value in x,y,z.\n        Only called internally by HARKinterpolator4D.__call__ (etc).\n        '
    if _isscalar(w):
        w_pos = max(min(self.wSearchFunc(self.w_list, w), (self.w_n - 1)), 1)
        x_pos = max(min(self.xSearchFunc(self.x_list, x), (self.x_n - 1)), 1)
        y_pos = max(min(self.ySearchFunc(self.y_list, y), (self.y_n - 1)), 1)
        z_pos = max(min(self.zSearchFunc(self.z_list, z), (self.z_n - 1)), 1)
    else:
        w_pos = self.wSearchFunc(self.w_list, w)
        w_pos[(w_pos < 1)] = 1
        w_pos[(w_pos > (self.w_n - 1))] = (self.w_n - 1)
        x_pos = self.xSearchFunc(self.x_list, x)
        x_pos[(x_pos < 1)] = 1
        x_pos[(x_pos > (self.x_n - 1))] = (self.x_n - 1)
        y_pos = self.ySearchFunc(self.y_list, y)
        y_pos[(y_pos < 1)] = 1
        y_pos[(y_pos > (self.y_n - 1))] = (self.y_n - 1)
        z_pos = self.zSearchFunc(self.z_list, z)
        z_pos[(z_pos < 1)] = 1
        z_pos[(z_pos > (self.z_n - 1))] = (self.z_n - 1)
    i = w_pos
    j = x_pos
    k = y_pos
    l = z_pos
    alpha = ((w - self.w_list[(i - 1)]) / (self.w_list[i] - self.w_list[(i - 1)]))
    beta = ((x - self.x_list[(j - 1)]) / (self.x_list[j] - self.x_list[(j - 1)]))
    gamma = ((y - self.y_list[(k - 1)]) / (self.y_list[k] - self.y_list[(k - 1)]))
    delta = ((z - self.z_list[(l - 1)]) / (self.z_list[l] - self.z_list[(l - 1)]))
    f = (((1 - alpha) * (((1 - beta) * ((((((1 - gamma) * (1 - delta)) * self.f_values[((i - 1), (j - 1), (k - 1), (l - 1))]) + (((1 - gamma) * delta) * self.f_values[((i - 1), (j - 1), (k - 1), l)])) + ((gamma * (1 - delta)) * self.f_values[((i - 1), (j - 1), k, (l - 1))])) + ((gamma * delta) * self.f_values[((i - 1), (j - 1), k, l)]))) + (beta * ((((((1 - gamma) * (1 - delta)) * self.f_values[((i - 1), j, (k - 1), (l - 1))]) + (((1 - gamma) * delta) * self.f_values[((i - 1), j, (k - 1), l)])) + ((gamma * (1 - delta)) * self.f_values[((i - 1), j, k, (l - 1))])) + ((gamma * delta) * self.f_values[((i - 1), j, k, l)]))))) + (alpha * (((1 - beta) * ((((((1 - gamma) * (1 - delta)) * self.f_values[(i, (j - 1), (k - 1), (l - 1))]) + (((1 - gamma) * delta) * self.f_values[(i, (j - 1), (k - 1), l)])) + ((gamma * (1 - delta)) * self.f_values[(i, (j - 1), k, (l - 1))])) + ((gamma * delta) * self.f_values[(i, (j - 1), k, l)]))) + (beta * ((((((1 - gamma) * (1 - delta)) * self.f_values[(i, j, (k - 1), (l - 1))]) + (((1 - gamma) * delta) * self.f_values[(i, j, (k - 1), l)])) + ((gamma * (1 - delta)) * self.f_values[(i, j, k, (l - 1))])) + ((gamma * delta) * self.f_values[(i, j, k, l)]))))))
    return f