def compute_coeffs_discrete(xyz, mass, nmax, lmax, r_s,
                            skip_odd=False, skip_even=False, skip_m=False,
                            compute_var=False):
    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))
    Tnlm = np.zeros((nmax+1, lmax+1, lmax+1))
    if compute_var:
        Snlm_var = np.zeros((nmax+1, lmax+1, lmax+1))
        Tnlm_var = np.zeros((nmax+1, lmax+1, lmax+1))
    xyz = np.ascontiguousarray(np.atleast_2d(xyz))
    mass = np.ascontiguousarray(np.atleast_1d(mass))
    r = np.sqrt(np.sum(xyz**2, axis=-1))
    s = r / r_s
    phi = np.arctan2(xyz[:,1], xyz[:,0])
    X = xyz[:,2] / r
    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
                Snlm[n,l,m], Tnlm[n,l,m] = STnlm_discrete(s, phi, X, mass, n, l, m)
                if compute_var:
                    Snlm_var[n,l,m], Tnlm_var[n,l,m] = STnlm_var_discrete(s, phi, X, mass, n, l, m)
    if compute_var:
        return (Snlm,Snlm_var), (Tnlm,Tnlm_var)
    else:
        return Snlm, Tnlm