> In examples/pse.py:
> 
> > +        """
> > +        # Set up the scales of the sky and pixels
> > +        self.N = N
> > +        self.sky_size = np.radians(sky_size_deg)
> > +        self.dx = self.sky_size / N
> > +
> > +        # Define the possible ell range, the bin edges and effective ell values.
> > +        # This is necessary for binning the power spectrum in ell.
> > +        lmin = 2*pi / self.sky_size
> > +        lmax = np.sqrt(2.)*pi / self.dx # in 2 dimensions
> > +        self.bin_edges = np.logspace(np.log10(lmin), np.log10(lmax), nbin+1)
> > +        # By default, report an area-averaged value of ell, which should be fine if there is
> > +        # no weighting (in which case it's recomputed) and if there are many ell modes in
> > +        # each bin.  The latter assumption is most likely to break down at low ell.
> > +        self.ell = (2./3.)*(self.bin_edges[1:]**3-self.bin_edges[:-1]**3) \
> > +                                   / (self.bin_edges[1:]**2-self.bin_edges[:-1]**2)
> It has just occurred to me that this area averaged bin weighting is not correct for the higher frequency bins in which the lower ell edge is greater than pi / self.dx, due to the annular ell region being cut off by the square grid edges. This seems like something that is a pain to calculate and correct for, but maybe we should mention this in the dox? I'm happy to add a few words if people agree.
> 

Yes, please do add some warning about that, if you don't mind.