code: def minimal_rotation(R, t, iterations=2):
      
      from scipy.interpolate import InterpolatedUnivariateSpline as spline
      if iterations == 0:
          return R
      R = quaternion.as_float_array(R)
      Rdot = np.empty_like(R)
      for i in range(4):
          Rdot[:, i] = spline(t, R[:, i]).derivative()(t)
      R = quaternion.from_float_array(R)
      Rdot = quaternion.from_float_array(Rdot)
      halfgammadot = quaternion.as_float_array(Rdot * quaternion.z * R.conjugate())[:, 0]
      halfgamma = spline(t, halfgammadot).antiderivative()(t)
      Rgamma = np.exp(quaternion.z * halfgamma)
      return minimal_rotation(R * Rgamma, t, iterations=iterations-1)
docs: Adjust frame so that there is no rotation about z' axis
  
      The output of this function is a frame that rotates the z axis onto the same z' axis as the
      input frame, but with minimal rotation about that axis.  This is done by pre-composing the input
      rotation with a rotation about the z axis through an angle gamma, where
  
          dgamma/dt = 2*(dR/dt * z * R.conjugate()).w
  
      This ensures that the angular velocity has no component along the z' axis.
  
      Note that this condition becomes easier to impose the closer the input rotation is to a
      minimally rotating frame, which means that repeated application of this function improves its
      accuracy.  By default, this function is iterated twice, though a few more iterations may be
      called for.
  
      Parameters
      ==========
      R: quaternion array
          Time series describing rotation
      t: float array
          Corresponding times at which R is measured
      iterations: int [defaults to 2]
          Repeat the minimization to refine the result