How to use the skyfield.functions.mxm function in skyfield

def rotation_and_rate_at(self, t):
        """Return rotation and rate matrices for this frame at time ``t``."""
        components, rates = self._segment.compute(t.tdb, 0.0, True)
        ra, dec, w = components
        R = mxm(rot_z(-w), mxm(rot_x(-dec), rot_z(-ra)))

        zero = w * 0.0
        one = 1.0 + zero
        ca = cos(w)
        sa = sin(w)
        u = cos(dec)
        v = -sin(dec)

        solutn = array((
            (one, zero, u),
            (zero, ca, -sa * v),
            (zero, sa, ca * v),
        ))

        domega = mxv(solutn, rates[::-1])

def build_ecliptic_matrix(t):
    """Build the matrix to rotate an ICRF vector into ecliptic coordinates."""
    _, d_eps = t._nutation_angles_radians
    true_obliquity = t._mean_obliquity_radians + d_eps
    return mxm(rot_x(- true_obliquity), t.M)

def rotation_at(self, t):
        """Return the rotation matrix for this frame at time ``t``."""
        ra, dec, w = self._segment.compute(t.tdb, 0.0, False)
        R = mxm(rot_z(-w), mxm(rot_x(-dec), rot_z(-ra)))
        if self._matrix is not None:
            R = mxm(self._matrix, R)
        return R

def C(self):
        # Calculate the Equation of Origins in cycles
        eq_origins = (earth_rotation_angle(self.ut1) - self.gast / 24.0)
        R = rot_z(2 * pi * eq_origins)
        return mxm(R, self.M)

solutn = array((
            (one, zero, u),
            (zero, ca, -sa * v),
            (zero, sa, ca * v),
        ))

        domega = mxv(solutn, rates[::-1])

        drdtrt = array((
            (zero, domega[0], domega[2]),
            (-domega[0], zero, domega[1]),
            (-domega[2], -domega[1], zero),
        ))

        dRdt = mxm(drdtrt, R)

        if self._matrix is not None:
            R = mxm(self._matrix, R)
            dRdt = mxm(self._matrix, dRdt)

        return R, dRdt

(zero, ca, -sa * v),
            (zero, sa, ca * v),
        ))

        domega = mxv(solutn, rates[::-1])

        drdtrt = array((
            (zero, domega[0], domega[2]),
            (-domega[0], zero, domega[1]),
            (-domega[2], -domega[1], zero),
        ))

        dRdt = mxm(drdtrt, R)

        if self._matrix is not None:
            R = mxm(self._matrix, R)
            dRdt = mxm(self._matrix, dRdt)

        return R, dRdt

center = self._get_assignment('FRAME_{0}_CENTER'.format(integer))
        spec = self.assignments.get('TKFRAME_{0}_SPEC'.format(integer))
        if spec is None:
            matrix = None
        else:
            if spec == 'ANGLES':
                angles = self.assignments['TKFRAME_{0}_ANGLES'.format(integer)]
                axes = self.assignments['TKFRAME_{0}_AXES'.format(integer)]
                units = self.assignments['TKFRAME_{0}_UNITS'.format(integer)]
                scale = _unit_scales[units]
                matrix = 1,0,0, 0,1,0, 0,0,1
                matrix = array(matrix)
                matrix.shape = 3, 3
                for angle, axis in list(zip(angles, axes)):
                    rot = _rotations[axis]
                    matrix = mxm(rot(angle * scale), matrix)
            elif spec == 'MATRIX':
                matrix = self.assignments['TKFRAME_{0}_MATRIX'.format(integer)]
                matrix = array(matrix)
                matrix.shape = 3, 3
            else:
                raise NotImplementedError('spec %r not yet implemented' % spec)
            relative = self.assignments['TKFRAME_{0}_RELATIVE'.format(integer)]
            integer = self.assignments['FRAME_{0}'.format(relative)]

        if _segment is None:
            segment = self._segment_map.get(integer)
        else:
            segment = _segment

        if segment is None:
            raise LookupError('you have not yet loaded a binary PCK file that'

(zero, sa, ca * v),
        ))

        domega = mxv(solutn, rates[::-1])

        drdtrt = array((
            (zero, domega[0], domega[2]),
            (-domega[0], zero, domega[1]),
            (-domega[2], -domega[1], zero),
        ))

        dRdt = mxm(drdtrt, R)

        if self._matrix is not None:
            R = mxm(self._matrix, R)
            dRdt = mxm(self._matrix, dRdt)

        return R, dRdt