How to use the pixell.utils.interpol function in pixell

def __call__(self, x):
		ix = ((x.T-self.box[0])/(self.box[1]-self.box[0])*(np.array(self.y.shape[1:])-1)).T
		return utils.interpol(self.y, ix, *self.args, **self.kwargs)

def at(map, pos, order=3, mode="constant", cval=0.0, unit="coord", prefilter=True, mask_nan=False, safe=True):
	if unit != "pix": pos = sky2pix(map.shape, map.wcs, pos, safe=safe)
	return utils.interpol(map, pos, order=order, mode=mode, cval=cval, prefilter=prefilter, mask_nan=mask_nan)

"""Given an array d[...,nt] of data that has been sampled at times t[nt],
	return an array that has been resampled to have a constant sampling rate."""
	# Find the typical sampling rate of the input. We will lose information if
	# we don't use a sampling rate that's higher than the highest rate in the
	# input. But we also don't want to exaggerate the number of samples. Use a
	# configurable quantile as a compromise.
	dt   = np.percentile(np.abs(t[1:]-t[:-1]), quantile*100)
	# Modify so we get a whole number of samples
	nout = utils.nint(np.abs(t[-1]-t[0])/dt)+1
	dt   = (t[-1]-t[0])/(nout-1)
	# Construct our output time steps
	tout = np.arange(nout)*dt + t[0]
	# To interpolate, we need the input sample number as a function of time
	samples = np.interp(tout, t, np.arange(len(t)))
	# Now that we have the samples we can finally evaluate the function
	dout = utils.interpol(d, samples[None], mode="nearest", order=order, mask_nan=mask_nan)
	return dout, tout

model = enmap.zeros(amps.shape[-1:]+posmap.shape[-2:], posmap.wcs, dtype)
	for cy in range(ncy):
		for cx in range(ncx):
			nsrc = nhit[cy,cx]
			if verbose: print("map cell %5d/%d with %5d srcs" % (cy*ncx+cx+1, ncy*ncx, nsrc))
			if nsrc == 0: continue
			srcs  = cell_srcs[cy,cx,:nsrc]
			y1,y2 = (cy+0)*cres[0], (cy+1)*cres[0]
			x1,x2 = (cx+0)*cres[1], (cx+1)*cres[1]
			pixpos = posmap[:,y1:y2,x1:x2]
			srcpos = poss[srcs].T # [2,nsrc]
			srcamp = amps[srcs].T # [ncomp,nsrc]
			r      = utils.angdist(pixpos[::-1,None,:,:],srcpos[::-1,:,None,None])
			bpix   = (r - beam[0,0])/(beam[0,1]-beam[0,0])
			# Evaluate the beam at these locations
			bval   = utils.interpol(beam[1], bpix[None], mode="constant", order=1, mask_nan=False) # [nsrc,ry,rx]
			cmodel = srcamp[:,:,None,None]*bval
			cmodel = op.reduce(cmodel,-3)
			op(model[:,y1:y2,x1:x2], cmodel, model[:,y1:y2,x1:x2])
	return model

def aberrate(imap, dir, beta, mode="wrap", order=3, recenter=False, modulation=True):
	pol = imap.ndim > 2
	pos = imap.posmap()
	# The ::-1 stuff switches between dec,ra and ra,dec ordering.
	# It is a bit confusing to have different conventions in enmap
	# and coordinates.
	pos = remap(pos[::-1], dir, beta, pol=pol, recenter=recenter, modulation=modulation)
	pos[:2] = pos[1::-1]
	pix = imap.sky2pix(pos[:2], corner=True) # interpol needs corners
	omap= enmap.ndmap(utils.interpol(imap, pix, mode=mode, order=order), imap.wcs)
	if pol:
		c,s = np.cos(2*pos[2]), np.sin(2*pos[2])
		omap[1] = c*omap[1] + s*omap[2]
		omap[2] =-s*omap[1] + c*omap[2]
	if modulation:
		omap *= pos[2+pol]
	return omap

ls = np.sum(lmap(oshape, wcs, oversample=oversample)**2,0)**0.5
	if smooth == "auto":
		# Determine appropriate fourier-scale smoothing based on 2d fourer
		# space resolution. We wish to smooth by about this width to approximate
		# averaging over sub-grid modes
		smooth = 0.5*(ls[1,0]+ls[0,1])
		smooth /= 3.41 # 3.41 is an empirical factor
	if smooth > 0:
		cov = smooth_spectrum(cov, kernel="gauss", weight="mode", width=smooth)
	# Translate from steradians to pixels
	cov = cov * np.prod(shape[-2:])/area(shape,wcs)
	if exp != 1.0: cov = multi_pow(cov, exp)
	cov[~np.isfinite(cov)] = 0
	# Use order 1 because we will perform very short interpolation, and to avoid negative
	# values in spectra that must be positive (and it's faster)
	res = ndmap(utils.interpol(cov, np.reshape(ls,(1,)+ls.shape),mode=mode, mask_nan=False, order=1),wcs)
	res = downgrade(res, oversample)
	res = res.reshape(oshape[:-2]+res.shape[-2:])
	return res