code: def generate_momentum(d, massmatrix, alpha, beta):
      """draw momentum from a circle, with amplitude following the beta distribution"""
      momentum = np.random.multivariate_normal(
          np.zeros(d), np.dot(massmatrix, np.eye(d))
      )
      # generate normalisation from beta distribution
      # add a bit of noise in the step size
      # norm *= np.uniform(0.2, 2)
      betainc = scipy.special.betainc
      auxnorm = (
          -betainc(alpha + 1, beta, 1)
          + betainc(alpha + 1, beta, 0)
          + betainc(alpha, beta, 1)
      )
      u = np.random.uniform()
      if u > 0.9:
          norm = 1.0
      else:
          u /= 0.9
          norm = betainc(alpha, beta, u)
      momnorm = -np.log((norm + 1e-10) / auxnorm)
      assert momnorm >= 0, (momnorm, norm, auxnorm)
  
      momentum *= momnorm / (momentum**2).sum() ** 0.5
  
      return momentum
  
imports: import numpy
  import matplotlib.pyplot
  import scipy.special
  import scipy.stats
label_desc: Use image node_id 8 for calling a global function with example usage: generate_momentum(d, massmatrix, alpha, beta) and returns: momentum