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

def frame_latlon(self, frame):
        """Return as longitude, latitude, and distance in the given frame."""
        R = frame.rotation_at(self.t)
        vector = mxv(R, self.position.au)
        d, lat, lon = to_spherical(vector)
        return (Angle(radians=lat, signed=True),
                Angle(radians=lon),
                Distance(au=d))

def _to_altaz(position_au, observer_data, temperature_C, pressure_mbar):
    """Compute (alt, az, distance) relative to the observer's horizon.

    """
    elevation_m = observer_data.elevation_m
    R = observer_data.altaz_rotation

    if (elevation_m is None) or (R is None):
        raise ValueError('to compute an altazimuth position, you must'
                         ' observe from a specific Earth location or from'
                         ' a position on another body loaded from a set'
                         ' of planetary constants')

    # TODO: wobble

    position_au = mxv(R, position_au)
    r_au, alt, az = to_spherical(position_au)

    if temperature_C is None:
        alt = Angle(radians=alt)
    else:
        if temperature_C == 'standard':
            temperature_C = 10.0
        if pressure_mbar == 'standard':
            pressure_mbar = 1010.0 * exp(-elevation_m / 9.1e3)
        alt = refract(alt * RAD2DEG, temperature_C, pressure_mbar)
        alt = Angle(degrees=alt)

    return alt, Angle(radians=az), Distance(r_au)

def _at(self, t):
        R, dRdt = self._frame.rotation_and_rate_at(t)
        r = mxv(_T(R), self._position_au)
        v = mxv(_T(dRdt), self._position_au) * DAY_S
        return r, v, None, None

# the precision by something like 1e5 times when compared to using
    # the round number skyfield.constants.ANGVEL!
    #
    # See the test `test_velocity_in_ITRF_to_GCRS2()`.
    #
    if _high_accuracy:
        _one_second = 1.0 / DAY_S
        t_later = t.ts.tt_jd(t.whole, t.tt_fraction + _one_second)
        angvel = (t_later.gast - t.gast) / 24.0 * tau
    else:
        angvel = ANGVEL

    velocity[0] += DAY_S * angvel * - position[1]
    velocity[1] += DAY_S * angvel * position[0]

    position = mxv(t.MT, position)
    velocity = mxv(t.MT, velocity)

    return position, velocity

def ITRF_to_GCRS(t, rITRF):

    # Todo: wobble

    spin = rot_z(t.gast / 24.0 * tau)
    position = mxv(spin, array(rITRF))
    return mxv(t.MT, position)

def _at(self, t):
        R, dRdt = self._frame.rotation_and_rate_at(t)
        r = mxv(_T(R), self._position_au)
        v = mxv(_T(dRdt), self._position_au) * DAY_S
        return r, v, None, None

def frame_xyz(self, frame):
        """Express this position as an (x,y,z) vector in a particular frame."""
        R = frame.rotation_at(self.t)
        return Distance(au=mxv(R, self.position.au))

(CIRS), a dynamical coordinate system referenced to the Celestial
        Intermediate Origin (CIO). As this is a dynamical system it must be
        calculated at a specific epoch.
        """
        if isinstance(epoch, Time):
            pass
        elif isinstance(epoch, float):
            epoch = Time(None, tt=epoch)
        elif epoch == 'date':
            epoch = self.t
        else:
            raise ValueError('the epoch= must be a Time object,'
                             ' a floating point Terrestrial Time (TT),'
                             ' or the string "date" for epoch-of-date')

        vector = mxv(epoch.C, self.position.au)
        return Distance(vector)

"""
    # TODO: are xp and yp the values from the IERS?  Or from general
    # nutation theory?

    theta, theta_dot = theta_GMST1982(jd_ut1, fraction_ut1)
    angular_velocity = multiply.outer(_zero_zero_minus_one, theta_dot)

    R = rot_z(-theta)

    if len(rTEME.shape) == 1:
        rPEF = (R).dot(rTEME)
        vPEF = (R).dot(vTEME) + _cross(angular_velocity, rPEF)
    else:
        rPEF = mxv(R, rTEME)
        vPEF = mxv(R, vTEME) + _cross(angular_velocity, rPEF)

    if xp == 0.0 and yp == 0.0:
        rITRF = rPEF
        vITRF = vPEF
    else:
        W = (rot_x(yp)).dot(rot_y(xp))
        rITRF = (W).dot(rPEF)
        vITRF = (W).dot(vPEF)
    return rITRF, vITRF