def _evaluate(self, x, y, z):
    '\n        Returns the level of the interpolated function at each value in x,y,z.\n        Only called internally by HARKinterpolator3D.__call__ (etc).\n        '
    if _isscalar(x):
        y_pos = max(min(np.searchsorted(self.y_list, y), (self.y_n - 1)), 1)
        z_pos = max(min(np.searchsorted(self.z_list, z), (self.z_n - 1)), 1)
        alpha = ((y - self.y_list[(y_pos - 1)]) / (self.y_list[y_pos] - self.y_list[(y_pos - 1)]))
        beta = ((z - self.z_list[(z_pos - 1)]) / (self.z_list[z_pos] - self.z_list[(z_pos - 1)]))
        f = ((((((1 - alpha) * (1 - beta)) * self.xInterpolators[(y_pos - 1)][(z_pos - 1)](x)) + (((1 - alpha) * beta) * self.xInterpolators[(y_pos - 1)][z_pos](x))) + ((alpha * (1 - beta)) * self.xInterpolators[y_pos][(z_pos - 1)](x))) + ((alpha * beta) * self.xInterpolators[y_pos][z_pos](x)))
    else:
        m = len(x)
        y_pos = np.searchsorted(self.y_list, y)
        y_pos[(y_pos > (self.y_n - 1))] = (self.y_n - 1)
        y_pos[(y_pos < 1)] = 1
        z_pos = np.searchsorted(self.z_list, z)
        z_pos[(z_pos > (self.z_n - 1))] = (self.z_n - 1)
        z_pos[(z_pos < 1)] = 1
        f = (np.zeros(m) + np.nan)
        for i in range(1, self.y_n):
            for j in range(1, self.z_n):
                c = np.logical_and((i == y_pos), (j == z_pos))
                if np.any(c):
                    alpha = ((y[c] - self.y_list[(i - 1)]) / (self.y_list[i] - self.y_list[(i - 1)]))
                    beta = ((z[c] - self.z_list[(j - 1)]) / (self.z_list[j] - self.z_list[(j - 1)]))
                    f[c] = ((((((1 - alpha) * (1 - beta)) * self.xInterpolators[(i - 1)][(j - 1)](x[c])) + (((1 - alpha) * beta) * self.xInterpolators[(i - 1)][j](x[c]))) + ((alpha * (1 - beta)) * self.xInterpolators[i][(j - 1)](x[c]))) + ((alpha * beta) * self.xInterpolators[i][j](x[c])))
    return f
Returns the level of the interpolated function at each value in x,y,z.
Only called internally by HARKinterpolator3D.__call__ (etc).
def _evaluate(self, x, y, z):
    '\n        Returns the level of the interpolated function at each value in x,y,z.\n        Only called internally by HARKinterpolator3D.__call__ (etc).\n        '
    if _isscalar(x):
        y_pos = max(min(np.searchsorted(self.y_list, y), (self.y_n - 1)), 1)
        z_pos = max(min(np.searchsorted(self.z_list, z), (self.z_n - 1)), 1)
        alpha = ((y - self.y_list[(y_pos - 1)]) / (self.y_list[y_pos] - self.y_list[(y_pos - 1)]))
        beta = ((z - self.z_list[(z_pos - 1)]) / (self.z_list[z_pos] - self.z_list[(z_pos - 1)]))
        f = ((((((1 - alpha) * (1 - beta)) * self.xInterpolators[(y_pos - 1)][(z_pos - 1)](x)) + (((1 - alpha) * beta) * self.xInterpolators[(y_pos - 1)][z_pos](x))) + ((alpha * (1 - beta)) * self.xInterpolators[y_pos][(z_pos - 1)](x))) + ((alpha * beta) * self.xInterpolators[y_pos][z_pos](x)))
    else:
        m = len(x)
        y_pos = np.searchsorted(self.y_list, y)
        y_pos[(y_pos > (self.y_n - 1))] = (self.y_n - 1)
        y_pos[(y_pos < 1)] = 1
        z_pos = np.searchsorted(self.z_list, z)
        z_pos[(z_pos > (self.z_n - 1))] = (self.z_n - 1)
        z_pos[(z_pos < 1)] = 1
        f = (np.zeros(m) + np.nan)
        for i in range(1, self.y_n):
            for j in range(1, self.z_n):
                c = np.logical_and((i == y_pos), (j == z_pos))
                if np.any(c):
                    alpha = ((y[c] - self.y_list[(i - 1)]) / (self.y_list[i] - self.y_list[(i - 1)]))
                    beta = ((z[c] - self.z_list[(j - 1)]) / (self.z_list[j] - self.z_list[(j - 1)]))
                    f[c] = ((((((1 - alpha) * (1 - beta)) * self.xInterpolators[(i - 1)][(j - 1)](x[c])) + (((1 - alpha) * beta) * self.xInterpolators[(i - 1)][j](x[c]))) + ((alpha * (1 - beta)) * self.xInterpolators[i][(j - 1)](x[c]))) + ((alpha * beta) * self.xInterpolators[i][j](x[c])))
    return f