docs: def sample_radius(self, n):
          """
          Sample the radial distribution (deg) from the 2D stellar density.
          Output is elliptical radius in true projected coordinates.
          """
          edge = self.edge if self.edge<20*self.extension else 20*self.extension
          radius = np.linspace(0, edge, 1.e5)
          pdf = self._pdf(radius) * np.sin(np.radians(radius))
          cdf = np.cumsum(pdf)
          cdf /= cdf[-1]
          fn = scipy.interpolate.interp1d(cdf, list(range(0, len(cdf))))
          index = np.floor(fn(np.random.uniform(size=n))).astype(int)
          return radius[index]