code: 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