docs: def sphere_points_from_angles_and_tilt(angles, tilted_axis):
      """
      For a given tilt of the rotational axis `tilted_axis`, compute
      the points on a unit sphere that correspond to the distribution
      `angles` along the great circle about this axis.
  
      Parameters
      ----------
      angles: 1d ndarray
          The angles that will be distributed on the great circle.
      tilted_axis: list of length 3
          The tilted axis of rotation that determines the great
          circle.
  
      Notes
      -----
      The reference axis is always [0,1,0].
      `theta` is the azimuthal angle measured down from the y-axis.
      `phi` is the polar angle in the x-z plane measured from z towards x.
  
      """
      assert len(angles.shape) == 1
      # Normalize tilted axis.
      tilted_axis = norm_vec(tilted_axis)
      [u, v, w] = tilted_axis
      # Initial distribution of points about great circle (x-z).
      newang = np.zeros((angles.shape[0], 3), dtype=float)
      # We subtract angles[0], because in step (a) we want that
      # newang[0]==[0,0,1]. This only works if we actually start
      # at that point.
      newang[:, 0] = np.sin(angles-angles[0])
      newang[:, 2] = np.cos(angles-angles[0])
  
      # Compute rotational angles w.r.t. [0,1,0].
      # - Draw a unit sphere with the y-axis pointing up and the
      #   z-axis pointing right
      # - The rotation of `tilted_axis` can be described by two
      #   separate rotations. We will use these two angles:
      #   (a) Rotation from y=1 within the y-z plane: theta
      #       This is the rotation that is critical for data
      #       reconstruction. If this angle is zero, then we
      #       have a rotational axis in the imaging plane. If
      #       this angle is PI/2, then our sinogram consists
      #       of a rotating image and 3D reconstruction is
      #       impossible. This angle is counted from the y-axis
      #       onto the x-z plane.
      #   (b) Rotation in the x-z plane: phi
      #       This angle is responsible for matching up the angles
      #       with the correct sinogram images. If this angle is zero,
      #       then the projection of the rotational axis onto the
      #       x-y plane is aligned with the y-axis. If this angle is
      #       PI/2, then the axis and its projection onto the x-y
      #       plane are identical. This angle is counted from the
      #       positive z-axis towards the positive x-axis. By default,
      #       angles[0] is the point that touches the great circle
      #       that lies in the x-z plane. angles[1] is the next point
      #       towards the x-axis if phi==0.
  
      # (a) This angle is the azimuthal angle theta measured from the
      #     y-axis.
      theta = np.arccos(v)
  
      # (b) This is the polar angle measured in the x-z plane starting
      #     at the x-axis and measured towards the positive z-axis.
      if np.allclose(u, 0) and np.allclose(w, 0):
          # Avoid flipping the axis of rotation due to numerical
          # errors during its computation.
          phi = 0
      else:
          phi = np.arctan2(u, w)
  
      # Determine the projection points on the unit sphere.
      # The resulting circle meets the x-z-plane at phi, and
      # is tilted by theta w.r.t. the y-axis.
  
      # (a) Create a tilted data set. This is achieved in 3 steps.
  
      # a1) Determine radius of tilted circle and get the centered
      #     circle with a smaller radius.
      rtilt = np.cos(theta)
      newang *= rtilt
  
      # a2) Rotate this circle about the x-axis by theta
      #     (right-handed/counter-clockwise/basic/elemental rotation)
      Rx = np.array([
          [1,             0,              0],
          [0, np.cos(theta), -np.sin(theta)],
          [0, np.sin(theta),  np.cos(theta)]
      ])
      for ii in range(newang.shape[0]):
          newang[ii] = np.dot(Rx, newang[ii])
  
      # a3) Shift newang such that newang[0] is located at (0,0,1)
      newang = newang - (newang[0] - np.array([0, 0, 1])).reshape(1, 3)
  
      # (b) Rotate the entire thing with phi about the y-axis
      #     (right-handed/counter-clockwise/basic/elemental rotation)
      Ry = np.array([
          [+np.cos(phi), 0, np.sin(phi)],
          [0,            1,           0],
          [-np.sin(phi), 0, np.cos(phi)]
      ])
  
      for jj in range(newang.shape[0]):
          newang[jj] = np.dot(Ry, newang[jj])
  
      # For visualiztaion:
      # import matplotlib.pylab as plt
      # from mpl_toolkits.mplot3d import Axes3D
      # from matplotlib.patches import FancyArrowPatch
      # from mpl_toolkits.mplot3d import proj3d
      #
      # class Arrow3D(FancyArrowPatch):
      #     def __init__(self, xs, ys, zs, *args, **kwargs):
      #         FancyArrowPatch.__init__(self, (0,0), (0,0), *args, **kwargs)
      #         self._verts3d = xs, ys, zs
      #
      #     def draw(self, renderer):
      #         xs3d, ys3d, zs3d = self._verts3d
      #         xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
      #         self.set_positions((xs[0],ys[0]),(xs[1],ys[1]))
      #         FancyArrowPatch.draw(self, renderer)
      #
      # fig = plt.figure(figsize=(10,10))
      # ax = fig.add_subplot(111, projection='3d')
      # for vec in newang:
      #     u,v,w = vec
      #     a = Arrow3D([0,u],[0,v],[0,w],
      #                 mutation_scale=20, lw=1, arrowstyle="-|>")
      #     ax.add_artist(a)
      #
      # radius=1
      # ax.set_xlabel('X')
      # ax.set_ylabel('Y')
      # ax.set_zlabel('Z')
      # ax.set_xlim(-radius*1.5, radius*1.5)
      # ax.set_ylim(-radius*1.5, radius*1.5)
      # ax.set_zlim(-radius*1.5, radius*1.5)
      # plt.tight_layout()
      # plt.show()
  
      return newang