def _evaluate(self,w,x,y,z):
        '''
        Returns the level of the interpolated function at each value in w,x,y,z.
        Only called internally by HARKinterpolator4D.__call__ (etc).
        '''
        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.x_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])
            gamma = (z - self.z_list[z_pos-1])/(self.z_list[z_pos] - self.z_list[z_pos-1])
            f = (
                (1-alpha)*(1-beta)*(1-gamma)*self.wInterpolators[x_pos-1][y_pos-1][z_pos-1](w)
              + (1-alpha)*(1-beta)*gamma*self.wInterpolators[x_pos-1][y_pos-1][z_pos](w)
              + (1-alpha)*beta*(1-gamma)*self.wInterpolators[x_pos-1][y_pos][z_pos-1](w)
              + (1-alpha)*beta*gamma*self.wInterpolators[x_pos-1][y_pos][z_pos](w)
              + alpha*(1-beta)*(1-gamma)*self.wInterpolators[x_pos][y_pos-1][z_pos-1](w)
              + alpha*(1-beta)*gamma*self.wInterpolators[x_pos][y_pos-1][z_pos](w)
              + alpha*beta*(1-gamma)*self.wInterpolators[x_pos][y_pos][z_pos-1](w)
              + alpha*beta*gamma*self.wInterpolators[x_pos][y_pos][z_pos](w))
        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
            f = 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])
                            gamma = (z[c] - self.z_list[k-1])/(self.z_list[k] - self.z_list[k-1])
                            f[c] = (
                                (1-alpha)*(1-beta)*(1-gamma)*self.wInterpolators[i-1][j-1][k-1](w[c])
                              + (1-alpha)*(1-beta)*gamma*self.wInterpolators[i-1][j-1][k](w[c])
                              + (1-alpha)*beta*(1-gamma)*self.wInterpolators[i-1][j][k-1](w[c])
                              + (1-alpha)*beta*gamma*self.wInterpolators[i-1][j][k](w[c])
                              + alpha*(1-beta)*(1-gamma)*self.wInterpolators[i][j-1][k-1](w[c])
                              + alpha*(1-beta)*gamma*self.wInterpolators[i][j-1][k](w[c])
                              + alpha*beta*(1-gamma)*self.wInterpolators[i][j][k-1](w[c])
                              + alpha*beta*gamma*self.wInterpolators[i][j][k](w[c]))
        return f
Returns the level of the interpolated function at each value in w,x,y,z.
        Only called internally by HARKinterpolator4D.__call__ (etc).