code: def principal_rotation_axis(gyro_data):
      
      N = np.zeros((3,3))
      for x in gyro_data.T: # Transpose because samples are stored as columns
          y = x.reshape(3,1)
          N += y.dot(y.T)
          
      (eig_val, eig_vec) = np.linalg.eig(N)
      i = np.argmax(eig_val)
      v = eig_vec[:,i]
      
      # Make sure v has correct sign
      s = 0
      for x in gyro_data.T: # Transpose because samples are stored as columns
          s += v.T.dot(x.reshape(3,1))
          
      v *= np.sign(s)
      
      return v
docs: Get the principal rotation axis of angular velocity measurements.    
      
      Parameters
      -------------
      gyro_data : (3, N) ndarray
              Angular velocity measurements
             
      Returns
      -------------
      v : (3,1) ndarray
              The principal rotation axis for the chosen sequence