<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def make_radial_kernel(psf, fn, sigma, npix, cdelt, xpix, ypix, psf_scale_fn=None, normalize=False, klims=None, sparse=False):
 """Make a kernel for a general radially symmetric 2D function. Parameters psf : `~fermipy.irfs.PSFModel` fn : callable Function that evaluates the kernel at a radial coordinate r. sigma : float 68% containment radius in degrees. """


    if klims is None:
        egy = psf.energies
    else:
        egy = psf.energies[klims[0]:klims[1] + 1]
    ang_dist = make_pixel_distance(npix, xpix, ypix) * cdelt
    max_ang_dist = np.max(ang_dist) + cdelt
    #dtheta = np.linspace(0.0, (np.max(ang_dist) * 1.05)**0.5, 200)**2.0
    # z = create_kernel_function_lookup(psf, fn, sigma, egy,
    #                                  dtheta, psf_scale_fn)

    shape = (len(egy), npix, npix)
    k = np.zeros(shape)

    r99 = psf.containment_angle(energies=egy, fraction=0.997)
    r34 = psf.containment_angle(energies=egy, fraction=0.34)

    rmin = np.maximum(r34 / 4., 0.01)
    rmax = np.maximum(r99, 0.1)
    if sigma is not None:
        rmin = np.maximum(rmin, 0.5 * sigma)
        rmax = np.maximum(rmax, 2.0 * r34 + 3.0 * sigma)
    rmax = np.minimum(rmax, max_ang_dist)

    for i in range(len(egy)):

        rebin = min(int(np.ceil(cdelt / rmin[i])), 8)
        if sparse:
            dtheta = np.linspace(0.0, rmax[i]**0.5, 100)**2.0
        else:
            dtheta = np.linspace(0.0, max_ang_dist**0.5, 200)**2.0

        z = eval_radial_kernel(psf, fn, sigma, i, dtheta, psf_scale_fn)
        xdist = make_pixel_distance(npix * rebin,
                                    xpix * rebin + (rebin - 1.0) / 2.,
                                    ypix * rebin + (rebin - 1.0) / 2.)
        xdist *= cdelt / float(rebin)
        #x = val_to_pix(dtheta, np.ravel(xdist))

        if sparse:
            m = np.ravel(xdist) < rmax[i]
            kk = np.zeros(xdist.size)
            #kk[m] = map_coordinates(z, [x[m]], order=2, prefilter=False)
            kk[m] = np.interp(np.ravel(xdist)[m], dtheta, z)
            kk = kk.reshape(xdist.shape)
        else:
            kk = np.interp(np.ravel(xdist), dtheta, z).reshape(xdist.shape)
            # kk = map_coordinates(z, [x], order=2,
            #                     prefilter=False).reshape(xdist.shape)

        if rebin > 1:
            kk = sum_bins(kk, 0, rebin)
            kk = sum_bins(kk, 1, rebin)

        k[i] = kk / float(rebin)**2

    k = k.reshape((len(egy),) + ang_dist.shape)
    if normalize:
        k /= (np.sum(k, axis=0)[np.newaxis, ...] * np.radians(cdelt) ** 2)

    return k