code: def compute_nutation(t):
      
      oblm, oblt, eqeq, psi, eps = t._earth_tilt
  
      cobm = cos(oblm * DEG2RAD)
      sobm = sin(oblm * DEG2RAD)
      cobt = cos(oblt * DEG2RAD)
      sobt = sin(oblt * DEG2RAD)
      cpsi = cos(psi * ASEC2RAD)
      spsi = sin(psi * ASEC2RAD)
  
      return array(((cpsi,
                    -spsi * cobm,
                    -spsi * sobm),
                    (spsi * cobt,
                     cpsi * cobm * cobt + sobm * sobt,
                     cpsi * sobm * cobt - cobm * sobt),
                    (spsi * sobt,
                     cpsi * cobm * sobt - sobm * cobt,
                     cpsi * sobm * sobt + cobm * cobt)))
docs: Generate the nutation rotations for Time `t`.
  
      If the Julian date is scalar, a simple ``(3, 3)`` matrix is
      returned; if the date is an array of length ``n``, then an array of
      matrices is returned with dimensions ``(3, 3, n)``.