def _derZ(self, w, x, y, z):
    '\n        Returns the derivative with respect to z of the interpolated function\n        at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeZ.\n        '
    if _isscalar(w):
        x_pos = max(min(np.searchsorted(self.x_list, x), (self.x_n - 1)), 1)
        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 = ((x - self.x_list[(x_pos - 1)]) / (self.y_list[x_pos] - self.x_list[(x_pos - 1)]))
        beta = ((y - self.y_list[(y_pos - 1)]) / (self.y_list[y_pos] - self.y_list[(y_pos - 1)]))
        dfdz = ((((((((1 - alpha) * (1 - beta)) * self.wInterpolators[(x_pos - 1)][(y_pos - 1)][z_pos](w)) + (((1 - alpha) * beta) * self.wInterpolators[(x_pos - 1)][y_pos][z_pos](w))) + ((alpha * (1 - beta)) * self.wInterpolators[x_pos][(y_pos - 1)][z_pos](w))) + ((alpha * beta) * self.wInterpolators[x_pos][y_pos][z_pos](w))) - ((((((1 - alpha) * (1 - beta)) * self.wInterpolators[(x_pos - 1)][(y_pos - 1)][(z_pos - 1)](w)) + (((1 - alpha) * beta) * self.wInterpolators[(x_pos - 1)][y_pos][(z_pos - 1)](w))) + ((alpha * (1 - beta)) * self.wInterpolators[x_pos][(y_pos - 1)][(z_pos - 1)](w))) + ((alpha * beta) * self.wInterpolators[x_pos][y_pos][(z_pos - 1)](w)))) / (self.z_list[z_pos] - self.z_list[(z_pos - 1)]))
    else:
        m = len(x)
        x_pos = np.searchsorted(self.x_list, x)
        x_pos[(x_pos > (self.x_n - 1))] = (self.x_n - 1)
        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
        dfdz = (np.zeros(m) + np.nan)
        for i in range(1, self.x_n):
            for j in range(1, self.y_n):
                for k in range(1, self.z_n):
                    c = np.logical_and(np.logical_and((i == x_pos), (j == y_pos)), (k == z_pos))
                    if np.any(c):
                        alpha = ((x[c] - self.x_list[(i - 1)]) / (self.x_list[i] - self.x_list[(i - 1)]))
                        beta = ((y[c] - self.y_list[(j - 1)]) / (self.y_list[j] - self.y_list[(j - 1)]))
                        dfdz[c] = ((((((((1 - alpha) * (1 - beta)) * self.wInterpolators[(i - 1)][(j - 1)][k](w[c])) + (((1 - alpha) * beta) * self.wInterpolators[(i - 1)][j][k](w[c]))) + ((alpha * (1 - beta)) * self.wInterpolators[i][(j - 1)][k](w[c]))) + ((alpha * beta) * self.wInterpolators[i][j][k](w[c]))) - ((((((1 - alpha) * (1 - beta)) * self.wInterpolators[(i - 1)][(j - 1)][(k - 1)](w[c])) + (((1 - alpha) * beta) * self.wInterpolators[(i - 1)][j][(k - 1)](w[c]))) + ((alpha * (1 - beta)) * self.wInterpolators[i][(j - 1)][(k - 1)](w[c]))) + ((alpha * beta) * self.wInterpolators[i][j][(k - 1)](w[c])))) / (self.z_list[k] - self.z_list[(k - 1)]))
    return dfdz
Returns the derivative with respect to z of the interpolated function
at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeZ.
def _derZ(self, w, x, y, z):
    '\n        Returns the derivative with respect to z of the interpolated function\n        at each value in w,x,y,z. Only called internally by HARKinterpolator4D.derivativeZ.\n        '
    if _isscalar(w):
        x_pos = max(min(np.searchsorted(self.x_list, x), (self.x_n - 1)), 1)
        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 = ((x - self.x_list[(x_pos - 1)]) / (self.y_list[x_pos] - self.x_list[(x_pos - 1)]))
        beta = ((y - self.y_list[(y_pos - 1)]) / (self.y_list[y_pos] - self.y_list[(y_pos - 1)]))
        dfdz = ((((((((1 - alpha) * (1 - beta)) * self.wInterpolators[(x_pos - 1)][(y_pos - 1)][z_pos](w)) + (((1 - alpha) * beta) * self.wInterpolators[(x_pos - 1)][y_pos][z_pos](w))) + ((alpha * (1 - beta)) * self.wInterpolators[x_pos][(y_pos - 1)][z_pos](w))) + ((alpha * beta) * self.wInterpolators[x_pos][y_pos][z_pos](w))) - ((((((1 - alpha) * (1 - beta)) * self.wInterpolators[(x_pos - 1)][(y_pos - 1)][(z_pos - 1)](w)) + (((1 - alpha) * beta) * self.wInterpolators[(x_pos - 1)][y_pos][(z_pos - 1)](w))) + ((alpha * (1 - beta)) * self.wInterpolators[x_pos][(y_pos - 1)][(z_pos - 1)](w))) + ((alpha * beta) * self.wInterpolators[x_pos][y_pos][(z_pos - 1)](w)))) / (self.z_list[z_pos] - self.z_list[(z_pos - 1)]))
    else:
        m = len(x)
        x_pos = np.searchsorted(self.x_list, x)
        x_pos[(x_pos > (self.x_n - 1))] = (self.x_n - 1)
        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
        dfdz = (np.zeros(m) + np.nan)
        for i in range(1, self.x_n):
            for j in range(1, self.y_n):
                for k in range(1, self.z_n):
                    c = np.logical_and(np.logical_and((i == x_pos), (j == y_pos)), (k == z_pos))
                    if np.any(c):
                        alpha = ((x[c] - self.x_list[(i - 1)]) / (self.x_list[i] - self.x_list[(i - 1)]))
                        beta = ((y[c] - self.y_list[(j - 1)]) / (self.y_list[j] - self.y_list[(j - 1)]))
                        dfdz[c] = ((((((((1 - alpha) * (1 - beta)) * self.wInterpolators[(i - 1)][(j - 1)][k](w[c])) + (((1 - alpha) * beta) * self.wInterpolators[(i - 1)][j][k](w[c]))) + ((alpha * (1 - beta)) * self.wInterpolators[i][(j - 1)][k](w[c]))) + ((alpha * beta) * self.wInterpolators[i][j][k](w[c]))) - ((((((1 - alpha) * (1 - beta)) * self.wInterpolators[(i - 1)][(j - 1)][(k - 1)](w[c])) + (((1 - alpha) * beta) * self.wInterpolators[(i - 1)][j][(k - 1)](w[c]))) + ((alpha * (1 - beta)) * self.wInterpolators[i][(j - 1)][(k - 1)](w[c]))) + ((alpha * beta) * self.wInterpolators[i][j][(k - 1)](w[c])))) / (self.z_list[k] - self.z_list[(k - 1)]))
    return dfdz