code: def shore_matrix(radial_order, zeta, gtab, tau=1 / (4 * np.pi**2)):
      r"""Compute the SHORE matrix for modified Merlet's 3D-SHORE [1]_
  
      ..math::
              :nowrap:
                  \begin{equation}
                      \textbf{E}(q\textbf{u})=\sum_{l=0, even}^{N_{max}}
                                              \sum_{n=l}^{(N_{max}+l)/2}
                                              \sum_{m=-l}^l c_{nlm}
                                              \phi_{nlm}(q\textbf{u})
                  \end{equation}
  
      where $\phi_{nlm}$ is
      ..math::
              :nowrap:
                  \begin{equation}
                      \phi_{nlm}^{SHORE}(q\textbf{u})=\Biggl[\dfrac{2(n-l)!}
                          {\zeta^{3/2} \Gamma(n+3/2)} \Biggr]^{1/2}
                          \Biggl(\dfrac{q^2}{\zeta}\Biggr)^{l/2}
                          exp\Biggl(\dfrac{-q^2}{2\zeta}\Biggr)
                          L^{l+1/2}_{n-l} \Biggl(\dfrac{q^2}{\zeta}\Biggr)
                          Y_l^m(\textbf{u}).
                  \end{equation}
  
      Parameters
      ----------
      radial_order : unsigned int,
          an even integer that represent the order of the basis
      zeta : unsigned int,
          scale factor
      gtab : GradientTable,
          gradient directions and bvalues container class
      tau : float,
          diffusion time. By default the value that makes q=sqrt(b).
  
      References
      ----------
      .. [1] Merlet S. et al., "Continuous diffusion signal, EAP and
      ODF estimation via Compressive Sensing in diffusion MRI", Medical
      Image Analysis, 2013.
  
      """
  
      qvals = np.sqrt(gtab.bvals / (4 * np.pi**2 * tau))
      qvals[gtab.b0s_mask] = 0
      bvecs = gtab.bvecs
  
      qgradients = qvals[:, None] * bvecs
  
      r, theta, phi = cart2sphere(
          qgradients[:, 0], qgradients[:, 1], qgradients[:, 2]
      )
      theta[np.isnan(theta)] = 0
      F = radial_order / 2
      n_c = int(np.round(1 / 6.0 * (F + 1) * (F + 2) * (4 * F + 3)))
      M = np.zeros((r.shape[0], n_c))
  
      counter = 0
      for l in range(0, radial_order + 1, 2):
          for n in range(l, int((radial_order + l) / 2) + 1):
              for m in range(-l, l + 1):
                  M[:, counter] = (
                      real_sh_descoteaux_from_index(m, l, theta, phi)
                      * genlaguerre(n - l, l + 0.5)(r**2 / zeta)
                      * np.exp(-(r**2) / (2.0 * zeta))
                      * _kappa(zeta, n, l)
                      * (r**2 / zeta) ** (l / 2)
                  )
                  counter += 1
      return M
  
imports: from warnings import warn
  from math import factorial
  import numpy
  from scipy.special import genlaguerre, gamma, hyp2f1
  from dipy.reconst.cache import Cache
  from dipy.reconst.multi_voxel import multi_voxel_fit
  from dipy.reconst.shm import real_sh_descoteaux_from_index
  from dipy.core.geometry import cart2sphere
  from dipy.utils.optpkg import optional_package
label_desc: Use image node_id 1 for calling a global function with example usage: shore_matrix(radial_order, zeta, gtab, tau) and returns: M