code: def interpolate(p, axis_values, pixelgrid, order=1, mode='constant', cval=0.0):
      
      # convert requested parameter combination into a coordinate
      p_ = np.array([np.searchsorted(av_,val) for av_, val in zip(axis_values,p)])
      lowervals_stepsize = np.array([[av_[p__-1], av_[p__]-av_[p__-1]] \
                           for av_, p__ in zip(axis_values,p_)])
      p_coord = (p-lowervals_stepsize[:,0])/lowervals_stepsize[:,1] + p_-1
      
      # interpolate
      if order > 1:
          prefilter = False
      else:
          prefilter = False
  
      return [ndimage.map_coordinates(pixelgrid[...,i],p_coord, order=order,
                                      prefilter=prefilter, mode=mode, cval=cval) \
                  for i in range(np.shape(pixelgrid)[-1])]
docs: Interpolates in a grid prepared by create_pixeltypegrid().
      
      p is an array of parameter arrays
      
      @param p: Npar x Ninterpolate array
      @type p: array
      @return: Ndata x Ninterpolate array
      @rtype: array