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

if n_vec.ndim == 1:
        if length_of(e_vec) < 1e-15: # circular
            return 0

        elif length_of(n_vec) < 1e-15: # equatorial and not circular
            angle = arctan2(e_vec[1], e_vec[0])%tau
            return angle if cross(pos_vec, vel_vec, 0, 0).T[2] >= 0 else -angle%tau

        else: # not circular and not equatorial
            angle = angle_between(n_vec, e_vec)
            return angle if e_vec[2] > 0 else -angle%tau
    else:
        w = zeros_like(pos_vec[0]) # defaults to 0 for circular orbits

        equatorial = length_of(n_vec) < 1e-15
        circular = length_of(e_vec) < 1e-15

        inds = ~circular*equatorial
        angle = arctan2(e_vec[1][inds], e_vec[0][inds])%tau
        condition = cross(pos_vec[:, inds], vel_vec[:, inds], 0, 0).T[2] >= 0
        w[inds] = where(condition, angle, -angle%tau)

        inds = ~circular*~equatorial
        angle = angle_between(n_vec[:, inds], e_vec[:, inds])
        condition = e_vec[2][inds] > 0
        w[inds] = where(condition, angle, -angle%tau)
        return w

distance = length_of(tposition - cposition)
    light_time0 = 0.0
    for i in range(10):
        light_time = distance / C_AUDAY
        delta = light_time - light_time0
        if abs(max(delta)) < 1e-12:
            break

        # We assume a light travel time of at most a couple of days.  A
        # longer light travel time would best be split into a whole and
        # fraction, for adding to the whole and fraction of TDB.
        t2 = ts.tdb_jd(whole, tdb_fraction - light_time)

        tposition, tvelocity, gcrs_position, message = target._at(t2)
        distance = length_of(tposition - cposition)
        light_time0 = light_time
    else:
        raise ValueError('light-travel time failed to converge')
    return tposition - cposition, tvelocity - cvelocity, t, light_time

Parameters
    ----------
    position : ndarray
       Position vector with shape (3,)
    velocity : ndarray
        Velocity vector with shape (3,)
    t0 : float
        Time corresponding to `position` and `velocity`
    t1 : float or ndarray
        Time or times to propagate to
    gm : float
        Gravitational parameter in units that match the other arguments
    """
    if gm <= 0:
        raise ValueError("'gm' should be positive")
    if length_of(velocity) == 0:
        raise ValueError('Velocity vector has zero magnitude')
    if length_of(position) == 0:
        raise ValueError('Position vector has zero magnitude')

    r0 = length_of(position)
    rv = dots(position, velocity)

    hvec = cross(position, velocity)
    h2 = dots(hvec, hvec)

    if h2 == 0:
        raise ValueError('Motion is not conical')

    eqvec = cross(velocity, hvec)/gm + -position/r0
    e = length_of(eqvec)
    q = h2 / (gm * (1+e))

"""Compute distance to intersections of a line and a sphere.

    Given a line through the origin (0,0,0) and an (x,y,z) ``endpoint``,
    and a sphere with the (x,y,z) ``center`` and scalar ``radius``,
    return the distance from the origin to their two intersections.

    If the line is tangent to the sphere, the two intersections will be
    at the same distance.  If the line does not intersect the sphere,
    two ``nan`` values will be returned.

    """
    # See http://paulbourke.net/geometry/circlesphere/index.html#linesphere
    # Names "b" and "c" designate the familiar values from the quadratic
    # formula; happily, a = 1 because we use a unit vector for the line.

    minus_b = 2.0 * (endpoint / length_of(endpoint) * center).sum(axis=0)
    c = (center * center).sum(axis=0) - radius * radius
    discriminant = minus_b * minus_b - 4 * c
    dsqrt = discriminant ** where(discriminant < 0, nan, 0.5)  # avoid sqrt(<0)
    return (minus_b - dsqrt) / 2.0, (minus_b + dsqrt) / 2.0

def light_time_difference(position, observer_position):
    """Returns the difference in light-time, for a star,
      between the barycenter of the solar system and the observer (or
      the geocenter).

    """
    # From 'pos1', form unit vector 'u1' in direction of star or light
    # source.

    dis = length_of(position)
    u1 = position / dis

    # Light-time returned is the projection of vector 'pos_obs' onto the
    # unit vector 'u1' (formed from 'pos1'), divided by the speed of light.

    diflt = einsum('a...,a...', u1, observer_position) / C_AUDAY
    return diflt

def speed(self):
        """Compute the magnitude of the velocity vector.

        >>> v = ICRF([0, 0, 0], [1, 2, 3])
        >>> print(v.speed())
        3.74166 au/day

        """
        return Velocity(length_of(self.velocity.au_per_d))

the positions of an `observer` and a `deflector` of reciprocal mass
    `rmass`, this function updates `position` in-place to show how much
    the presence of the deflector will deflect the image of the object.

    """
    # Construct vector 'pq' from gravitating body to observed object and
    # construct vector 'pe' from gravitating body to observer.

    pq = observer + position - deflector
    pe = observer - deflector

    # Compute vector magnitudes and unit vectors.

    pmag = length_of(position)
    qmag = length_of(pq)
    emag = length_of(pe)

    phat = position / where(pmag, pmag, 1.0)  # where() avoids divide-by-zero
    qhat = pq / where(qmag, qmag, 1.0)
    ehat = pe / where(emag, emag, 1.0)

    # Compute dot products of vectors.

    pdotq = dots(phat, qhat)
    qdote = dots(qhat, ehat)
    edotp = dots(ehat, phat)

    # If gravitating body is observed object, or is on a straight line
    # toward or away from observed object to within 1 arcsec, deflection
    # is set to zero set 'pos2' equal to 'pos1'.

    make_no_correction = abs(edotp) > 0.99999999999

def semi_latus_rectum(h_vec, mu):
    return length_of(h_vec)**2/mu

def speed(self):
        """Compute the magnitude of the velocity vector.

        >>> v = ICRF([0, 0, 0], [1, 2, 3])
        >>> print(v.speed())
        3.74166 au/day

        """
        return Velocity(length_of(self.velocity.au_per_d))