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

from __future__ import division

import sys
import math
from numpy import(abs, amax, amin, arange, arccos, arctan, array, cos, cosh,
                  cross, exp, log, ndarray, newaxis, ones_like, pi, power,
                  repeat, sin, sinh, sqrt, sum, tan, tile, zeros_like)

from skyfield.constants import AU_KM, DAY_S, DEG2RAD
from skyfield.functions import dots, length_of, mxv
from skyfield.descriptorlib import reify
from skyfield.elementslib import OsculatingElements, normpi
from skyfield.units import Distance, Velocity
from skyfield.vectorlib import VectorFunction

_CONVERT_GM = DAY_S * DAY_S / AU_KM / AU_KM / AU_KM

class _KeplerOrbit(VectorFunction):
    def __init__(self,
                 position,
                 velocity,
                 epoch,
                 mu_au3_d2,
                 center=None,
                 target=None,
        ):
        """ Calculates the position of an object using 2 body propagation

        Parameters
        ----------
        position : Distance
            Position vector at epoch with shape (3,)

# TODO: This is expensive, and should be extensively trimmed to only
    # include the most important terms underlying GAST.  But it improves
    # 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 km_per_s(self):
        return self.au_per_d * AU_KM / DAY_S

def ITRF_position_velocity_error(self, t):
        """Return the ITRF position, velocity, and error at time `t`.

        The position is an x,y,z vector measured in au, the velocity is
        an x,y,z vector measured in au/day, and the error is a vector of
        possible error messages for the time or vector of times `t`.

        """
        rTEME, vTEME, error = self._position_and_velocity_TEME_km(t)
        rTEME /= AU_KM
        vTEME /= AU_KM
        vTEME *= DAY_S
        rITRF, vITRF = TEME_to_ITRF(t.whole, rTEME, vTEME, 0.0, 0.0,
                                    t.ut1_fraction)
        return rITRF, vITRF, error

def find_maxima(start_time, end_time, f, epsilon=1.0 / DAY_S, num=12):
    """Find the local maxima in the values returned by a function of time.

    This routine is used to find events like highest altitude and
    maximum elongation.  See :doc:`searches` for how to use it yourself.

    """
    #    @@       @@_@@       @@_@@_@@_@@
    #   /  \     /     \     /           \
    # @@    @@ @@       @@ @@             @@
    # +1 -1    +1  0 -1    +1  0  0  0 -1    sd = sign(diff(y))
    # -2       -1 -1       -1  0  0 -1       diff(sign(diff(y))

    ts = start_time.ts
    jd0 = start_time.tt
    jd1 = end_time.tt

Supply the Terrestrial Time (TT) as a proleptic Gregorian
        calendar date:

        >>> t = ts.tt(2014, 1, 18, 1, 35, 37.5)
        >>> t.tt
        2456675.56640625
        >>> t.tt_calendar()
        (2014, 1, 18, 1, 35, 37.5)

        """
        if jd is not None:
            return self.tt_jd(jd)  # deprecate someday
        a = _to_array
        whole = julian_day(a(year), a(month), a(day)) - 0.5
        fraction = (a(second) + a(minute) * 60.0 + a(hour) * 3600.0) / DAY_S
        return Time(self, whole, fraction)

dist * cdc * cra,
            dist * cdc * sra,
            dist * sdc,
            ))

        # Compute Doppler factor, which accounts for change in light
        # travel time to star.

        k = 1.0 / (1.0 - self.radial_km_per_s / C * 1000.0)

        # Convert proper motion and radial velocity to orthogonal
        # components of motion with units of au/day.

        pmr = self.ra_mas_per_year / (parallax * 365.25) * k
        pmd = self.dec_mas_per_year / (parallax * 365.25) * k
        rvl = self.radial_km_per_s * DAY_S / self.au_km * k

        # Transform motion vector to equatorial system.

        self._velocity_au_per_d = array((
            - pmr * sra - pmd * sdc * cra + rvl * cdc * cra,
              pmr * cra - pmd * sdc * sra + rvl * cdc * sra,
              pmd * cdc + rvl * sdc,
              ))

# Compute local sidereal time factors at the observer's longitude.

    stlocl = 15.0 * DEG2RAD * gast + longitude
    sinst = sin(stlocl)
    cosst = cos(stlocl)

    # Compute position vector components in kilometers.

    ac = ach * cosphi
    acsst = ac * sinst
    accst = ac * cosst
    pos = array((accst, acsst, zero + ash * sinphi))

    # Compute velocity vector components in kilometers/sec.

    vel = ANGVEL * DAY_S * array((-acsst, accst, zero))

    return pos, vel

from skyfield.keplerlib import KeplerOrbit

ts = load.timescale(builtin=True)

from numpy import sqrt

df['semimajor_axis_au'] = (
    df['perihelion_distance_au'] / (1.0 - df['eccentricity'])
)

from skyfield.data.gravitational_parameters import GM_dict
from skyfield.constants import AU_KM, DAY_S

mu_km3_s2 = GM_dict[10]
mu_au3_d2 = mu_km3_s2 / (AU_KM**3.0) * (DAY_S**2.0)

row = df.ix[0]

t_perihelion = ts.tt(
    row.perihelion_year, row.perihelion_month, row.perihelion_day
)

df['mean_anomaly_degrees'] = (
    sqrt(mu_au3_d2 / (row.semimajor_axis_au ** 3.0))
    *
    (ts.J2000.tt - t_perihelion.tt)
    * 360.0 / tau
)

comet = df.iloc[0:1]

def periapsis_time(self):
        M = mean_anomaly(self.eccentric_anomaly.radians, self.eccentricity, shift=False)
        tp = time_since_periapsis(M,
                                  self.mean_motion_per_day.radians/DAY_S,
                                  self.true_anomaly.radians,
                                  self.semi_latus_rectum.km,
                                  self._mu)
        ts = self.time.ts
        times = self.time.tdb - tp/DAY_S
        return ts.tdb(jd=times)