How to use the skyfield.constants.tau function in skyfield

def season_at(t):
        """Return season 0 (Spring) through 3 (Winter) at time `t`."""
        t._nutation_angles_radians = iau2000b_radians(t)
        e = earth.at(t)
        _, slon, _ = e.observe(sun).apparent().ecliptic_latlon('date')
        return (slon.radians // (tau / 4) % 4).astype(int)

def longitude_of_ascending_node(i, h_vec):
    if i.ndim == 0:
        return arctan2(h_vec[0], -h_vec[1])%tau if i != 0 else float64(0)
    else:
        return arctan2(h_vec[0], -h_vec[1], out=zeros_like(i), where=i!=0)%tau

def convergeEccentricAnomaly(mean_anomaly, eccentricity, precision):
    # calculate the delta
    delta = 10 ** -precision

    # normalize the mean anomaly
    m = mean_anomaly % constants.tau

    # set up the first guess
    eccentric_anomaly = constants.tau
    if eccentricity < 0.8:
        eccentric_anomaly = m

    # do the initial test
    test = eccentric_anomaly - eccentricity * sin(m) - m

    # while we're not happy with the result, and we haven't been dawdling too long
    max_iterations = 30
    count = 0
    while ((math.fabs(test) > delta) and (count < max_iterations)):
        # calculate the next guess for an eccentric anomaly
        eccentric_anomaly = (
            eccentric_anomaly - test /

* 1 — Satellite culminated and started to descend again.
        * 2 — Satellite fell below ``altitude_degrees``.

        Note that multiple culminations in a row are possible when,
        without setting, the satellite reaches a second peak altitude
        after descending partway down the sky from the first one.

        """
        # First, we find the moments of maximum altitude over the time
        # period.  Some of these maxima will be negative, meaning the
        # satellite failed to crest the horizon.

        ts = t0.ts
        at = (self - topos).at
        half_second = 0.5 / DAY_S
        orbits_per_minute = self.model.no_kozai / tau
        orbits_per_day = 24 * 60 * orbits_per_minute

        # Note the protection against zero orbits_per_day.  I have not
        # yet learned with which satellite caused a user to run into a
        # ZeroDivisionError here, and without a concrete example I have
        # little sense of whether 1.0 is a good fallback value or not.
        rough_period = 1.0 / orbits_per_day if orbits_per_day else 1.0

        # Long-period satellites might rise each day not because of
        # their own motion, but because the Earth rotates under them, so
        # check position at least each quarter-day.  We might need to
        # tighten this even further if experience someday shows it
        # missing a pass of a particular satellite.
        if rough_period > 0.25:
            rough_period = 0.25

e = self.eccentricity
        I = self.inclination
        Om = self.longitude_ascending
        #n = 0.230605479
        n = 0.230652907

        w = self.argument_perihelion

        M = self.mean_anomaly
        d = date.tdb - self.epoch.tdb

        Om = Om / 360.0 * constants.tau
        w = w / 360.0 * constants.tau
        I = I / 360.0 * constants.tau
        M = M / 360.0 * constants.tau
        n = n / 360.0 * constants.tau

        M += d * n

        # calculate the mean anomaly in rads
        E = convergeEccentricAnomaly(M, e, 30)

        # calculate the initial primes
        x_prime = a * (cos(E) - e)
        y_prime = a * (1 - e ** 2.0) ** (0.5) * sin(E)

        """
        http://ssd.jpl.nasa.gov/txt/aprx_pos_planets.pdf
        x_ecl = cos(w)cos(Om)-sin(w)sin(Om)cos(I) * x_prime +
            (-sin(w)cos(Om) - cos(w)sin(Om)cos(I)) * y_prime
        y_ecl = cos(w)sin(Om)-sin(w)cos(Om)cos(I) * x_prime +
            (-sin(w)cos(Om) - cos(w)sin(Om)cos(I)) * y_prime

def reverse_terra(xyz_au, gast, iterations=3):
    """Convert a geocentric (x,y,z) at time `t` to latitude and longitude.

    Returns a tuple of latitude, longitude, and elevation whose units
    are radians and meters.  Based on Dr. T.S. Kelso's quite helpful
    article "Orbital Coordinate Systems, Part III":
    https://www.celestrak.com/columns/v02n03/

    """
    x, y, z = xyz_au
    R = sqrt(x*x + y*y)

    lon = (arctan2(y, x) - 15 * DEG2RAD * gast - pi) % tau - pi
    lat = arctan2(z, R)

    a = ERAD / AU_M
    f = 1.0 / IERS_2010_INVERSE_EARTH_FLATTENING
    e2 = 2.0*f - f*f
    i = 0
    C = 1.0
    while i < iterations:
        i += 1
        C = 1.0 / sqrt(1.0 - e2 * (sin(lat) ** 2.0))
        lat = arctan2(z + a * C * e2 * sin(lat), R)
    elevation_m = ((R / cos(lat)) - a * C) * AU_M
    return lat, lon, elevation_m

"""Open a BPC file, read its angles, and produce rotation matrices."""

import re
from numpy import array, cos, nan, sin
from jplephem.pck import DAF, PCK
from .constants import ASEC2RAD, AU_KM, DAY_S, tau
from .functions import _T, mxv, mxm, mxmxm, rot_x, rot_y, rot_z
from .units import Angle, Distance
from .vectorlib import VectorFunction

_TEXT_MAGIC_NUMBERS = b'KPL/FK', b'KPL/PCK'
_NAN3 = array((nan, nan, nan))
_halftau = tau / 2.0
_quartertau = tau / 4.0

class PlanetaryConstants(object):
    """Planetary constants manager.

    You can use this class to build working models of Solar System
    bodies by loading text planetary constant files and binary
    orientation kernels.  For a full description of how to use this, see
    :doc:`planetary`.

    """
    def __init__(self):
        self.assignments = {}
        self._binary_files = []
        self._segment_map = {}

(715923.2178 +
              (    31.8792 +
              (     0.051635 +
              (    -0.00024470)
              * t) * t) * t) * t) * ASEC2RAD
              + (1325.0*t % 1.0) * tau)

    # Mean Anomaly of the Sun.

    fa[1] = ((1287104.793048 +
              (1292581.0481 +
              (     -0.5532 +
              (     +0.000136 +
              (     -0.00001149)
              * t) * t) * t) * t) * ASEC2RAD
              + (99.0*t % 1.0) * tau)

    # Mean Longitude of the Moon minus Mean Longitude of the Ascending
    # Node of the Moon.

    fa[2] = (( 335779.526232 +
              ( 295262.8478 +
              (    -12.7512 +
              (     -0.001037 +
              (      0.00000417)
              * t) * t) * t) * t) * ASEC2RAD
              + (1342.0*t % 1.0) * tau)

    # Mean Elongation of the Moon from the Sun.

    fa[3] = ((1072260.703692 +
              (1105601.2090 +