code: def compute_coeffs(density_func, nmax, lmax, M, r_s, args=(),
                     skip_odd=False, skip_even=False, skip_m=False,
                     S_only=False, progress=False, **nquad_opts):
      from gala._cconfig import GSL_ENABLED
      if not GSL_ENABLED:
          raise ValueError("Gala was compiled without GSL and so this function "
                           "will not work.  See the gala documentation for more "
                           "information about installing and using GSL with "
                           "gala: http://gala.adrian.pw/en/latest/install.html")
      lmin = 0
      lstride = 1
      if skip_odd or skip_even:
          lstride = 2
      if skip_even:
          lmin = 1
      Snlm = np.zeros((nmax+1, lmax+1, lmax+1))
      Snlm_e = np.zeros((nmax+1, lmax+1, lmax+1))
      Tnlm = np.zeros((nmax+1, lmax+1, lmax+1))
      Tnlm_e = np.zeros((nmax+1, lmax+1, lmax+1))
      nquad_opts.setdefault('limit', 256)
      nquad_opts.setdefault('epsrel', 1E-10)
      limits = [[0, 2*np.pi],
                [-1, 1.],
                [-1, 1.]]
      nlms = []
      for n in range(nmax+1):
          for l in range(lmin, lmax+1, lstride):
              for m in range(l+1):
                  if skip_m and m > 0:
                      continue
                  nlms.append((n, l, m))
      if progress:
          try:
              from tqdm import tqdm
          except ImportError as e:
              raise ImportError('tqdm is not installed - you can install it '
                                'with `pip install tqdm`.\n' + str(e))
          iterfunc = tqdm
      else:
          iterfunc = lambda x: x
      for n, l, m in iterfunc(nlms):
          Snlm[n, l, m], Snlm_e[n, l, m] = si.nquad(
              Snlm_integrand, ranges=limits,
              args=(density_func, n, l, m, M, r_s, args),
              opts=nquad_opts)
          if not S_only:
              Tnlm[n, l, m], Tnlm_e[n, l, m] = si.nquad(
                  Tnlm_integrand, ranges=limits,
                  args=(density_func, n, l, m, M, r_s, args),
                  opts=nquad_opts)
      return (Snlm, Snlm_e), (Tnlm, Tnlm_e)