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)