docs: def normalize(u):
      '''
      normalize(u) yields a vetor with the same direction as u but unit length, or, if u has zero
      length, yields u.
      '''
      u = np.asarray(u)
      unorm = np.sqrt(np.sum(u**2, axis=0))
      z = np.isclose(unorm, 0)
      c = np.logical_not(z) / (unorm + z)
      return u * c