orbit

Implementation of formulae for calculating orbits

earth = PlanetParameters.earth() module-attribute

An instance of PlanetParameters representing Earth

OrbitDateAndAngle dataclass

Pairing of a date and the corresponding orbit angle

Source code in src/analemma/orbit.py

@dataclass
class OrbitDateAndAngle:
    """
    Pairing of a date and the corresponding orbit angle
    """

    date: datetime.date
    sigma: float

PlanetParameters dataclass

Parameters defining a planet for sundial calculation purposes

Parameters:

Name	Type	Description	Default
N	float	Number of mean days in a year	required
T_d	int	Number of seconds in a mean day	required
rho	float	Angle between axes of the ellipse and the equinoxes / solstices	required
alpha	float	Inclination of the earths axis of rotation	required
a	float	Length of the planet's orbit's semi-major axis	required
e	float	Eccentricity of the planet's orbit	required

The following attributes are calculated given the above parameters.

Attributes:

Name	Type	Description
T_y	int	Number of seconds in a mean year
om_y	float	Mean angular speed of the earth's centre of mass in its orbit
om_d	float	Angular speed of a point on the earth about the earth's centre of mass
om_sd	float	Angular speed of a point on an earth that revolves once per siderial day
T_sd	float	Number of seconds in a siderial day
Om	float	Angular speed parameter used in the spinor orbit formalism ( = om_y / 2 * a )

Source code in src/analemma/orbit.py

@dataclass
class PlanetParameters:
    """
    Parameters defining a planet for sundial calculation purposes

    Parameters:
        N: Number of mean days in a year
        T_d: Number of seconds in a mean day
        rho: Angle between axes of the ellipse and the equinoxes / solstices
        alpha: Inclination of the earths axis of rotation
        a: Length of the planet's orbit's semi-major axis
        e: Eccentricity of the planet's orbit

    The following attributes are calculated given the above parameters.

    Attributes:
        T_y: Number of seconds in a mean year
        om_y: Mean angular speed of the earth's centre of mass in its orbit
        om_d: Angular speed of a point on the earth about the earth's centre of mass
        om_sd: Angular speed of a point on an earth that revolves once per siderial day
        T_sd: Number of seconds in a siderial day
        Om: Angular speed parameter used in the spinor orbit formalism ( = om_y / 2 * a )
    """

    N: float
    T_d: int
    rho: float
    alpha: float
    a: float
    e: float
    T_y: int = field(init=False)
    om_y: float = field(init=False)
    om_d: float = field(init=False)
    om_sd: float = field(init=False)
    T_sd: float = field(init=False)
    Om: float = field(init=False)

    def __post_init__(self):
        self.T_y = self.N * self.T_d
        self.om_y = 2 * pi / self.T_y
        self.om_d = 2 * pi / self.T_d
        self.om_sd = (self.N + 1) / self.N * self.om_d
        self.T_sd = self.N / (self.N + 1) * self.T_d
        self.Om = pi / self.T_y * self.a

    def clone_with_eccentricity(self, e: float):
        """
        Return a new instance with the same parameters but a different eccentricity
        """
        return PlanetParameters(
            N=self.N,
            T_d=self.T_d,
            rho=self.rho,
            alpha=self.alpha,
            a=self.a,
            e=e,
        )

    def daily_noons(self) -> np.array:
        """
        Daily time samples in seconds from noon at perihelion
        """
        return self.T_d * np.arange(int(self.N))

    def rotation_angle(self, t: npt.ArrayLike) -> npt.ArrayLike:
        r"""
        Angle of planetary rotation increases linearly with time, one complete revolution
        after one siderial day, and an offset of rho

        Parameters:
            t: Collection of time in seconds starting at perihelion at which to calculate the planet's angle of rotation

        Returns:
            The angle $\psi$ measuring planetary rotation at the input times
        """
        return np.mod(self.rho + self.om_sd * t, 2 * pi)

    def orbit_angle(self, t: npt.ArrayLike) -> npt.ArrayLike:
        """

        Parameters:
            t: Collection of time in seconds starting at perihelion at which to calculate the planet's angle of rotation
        """
        phi = orbital_angle(spinor_time(t))
        return np.mod(
            pi + self.rho + phi, 2 * pi
        )  # phi starts at perihelion, sigma starts at winter solstice

    def __hash__(self):
        # this is not as efficient as it could be but covers the
        # common case of reuse of the same object for a given planet
        return id(self)

    @classmethod
    def earth(cls):
        """
        Return a PlanetParameters instance representing Earth
        """
        return PlanetParameters(
            N=365.2422,
            T_d=24 * 3600,
            rho=12.25 / 180 * pi,
            alpha=23.5 / 180 * pi,
            a=149598000000,
            e=0.017,
        )

clone_with_eccentricity(e)

Return a new instance with the same parameters but a different eccentricity

Source code in src/analemma/orbit.py

def clone_with_eccentricity(self, e: float):
    """
    Return a new instance with the same parameters but a different eccentricity
    """
    return PlanetParameters(
        N=self.N,
        T_d=self.T_d,
        rho=self.rho,
        alpha=self.alpha,
        a=self.a,
        e=e,
    )

daily_noons()

Daily time samples in seconds from noon at perihelion

Source code in src/analemma/orbit.py

def daily_noons(self) -> np.array:
    """
    Daily time samples in seconds from noon at perihelion
    """
    return self.T_d * np.arange(int(self.N))

earth() classmethod

Return a PlanetParameters instance representing Earth

Source code in src/analemma/orbit.py

@classmethod
def earth(cls):
    """
    Return a PlanetParameters instance representing Earth
    """
    return PlanetParameters(
        N=365.2422,
        T_d=24 * 3600,
        rho=12.25 / 180 * pi,
        alpha=23.5 / 180 * pi,
        a=149598000000,
        e=0.017,
    )

orbit_angle(t)

Parameters:

Name	Type	Description	Default
t	ArrayLike	Collection of time in seconds starting at perihelion at which to calculate the planet's angle of rotation	required

Source code in src/analemma/orbit.py

def orbit_angle(self, t: npt.ArrayLike) -> npt.ArrayLike:
    """

    Parameters:
        t: Collection of time in seconds starting at perihelion at which to calculate the planet's angle of rotation
    """
    phi = orbital_angle(spinor_time(t))
    return np.mod(
        pi + self.rho + phi, 2 * pi
    )  # phi starts at perihelion, sigma starts at winter solstice

rotation_angle(t)

Angle of planetary rotation increases linearly with time, one complete revolution after one siderial day, and an offset of rho

Parameters:

Name	Type	Description	Default
t	ArrayLike	Collection of time in seconds starting at perihelion at which to calculate the planet's angle of rotation	required

Returns:

Type	Description
ArrayLike	The angle \(\psi\) measuring planetary rotation at the input times

Source code in src/analemma/orbit.py

def rotation_angle(self, t: npt.ArrayLike) -> npt.ArrayLike:
    r"""
    Angle of planetary rotation increases linearly with time, one complete revolution
    after one siderial day, and an offset of rho

    Parameters:
        t: Collection of time in seconds starting at perihelion at which to calculate the planet's angle of rotation

    Returns:
        The angle $\psi$ measuring planetary rotation at the input times
    """
    return np.mod(self.rho + self.om_sd * t, 2 * pi)

orbit_date_to_day(the_date, year=None)

Convert from the date to the number of days since perihelion

Source code in src/analemma/orbit.py

def orbit_date_to_day(the_date: datetime.date, year: int = None) -> int:
    """
    Convert from the date to the number of days since perihelion
    """
    if not year:
        year = datetime.date.today().year
    return (the_date - _perihelion_date(year)).days

orbit_day_to_date(orbit_day, year=None)

Convert from the number of days since perihelion to the date

Source code in src/analemma/orbit.py

def orbit_day_to_date(orbit_day: int, year: int = None) -> datetime.date:
    """
    Convert from the number of days since perihelion to the date
    """
    if not year:
        year = datetime.date.today().year
    return _perihelion_date(year) + datetime.timedelta(days=int(orbit_day))

orbital_angle(s, planet=earth)

Calculate orbital angular coordinate given time parameter, phi(s)

Source code in src/analemma/orbit.py

def orbital_angle(s: npt.ArrayLike, planet: PlanetParameters = earth):
    """
    Calculate orbital angular coordinate given time parameter, phi(s)
    """
    A, B, Om, _ = _kepler_params(planet)
    tanSigY = (A**2 - B**2) * sin(2 * Om * s)
    tanSigX = (A**2 + B**2) * cos(2 * Om * s) + 2 * A * B
    return np.arctan2(tanSigY, tanSigX) + pi

orbital_radius(s, planet=earth)

Calculate orbital radial coordinate given spinor time parameter, r(s)

Source code in src/analemma/orbit.py

def orbital_radius(s: npt.ArrayLike, planet: PlanetParameters = earth):
    """
    Calculate orbital radial coordinate given spinor time parameter, r(s)
    """
    A, B, Om, _ = _kepler_params(planet)
    return A**2 + B**2 + 2 * A * B * cos(2 * Om * s)

orbital_time(s, planet=earth)

Calculate orbital time given time parameter, t(s)

Parameters:

Name	Type	Description	Default
s	ArrayLike	Spinor time parameter (note, called $ au$ in the paper)	required
planet	PlanetParameters	Planet whose orbit is being analyzed	earth
e		Optional override for the orbit's eccentricity	required

Source code in src/analemma/orbit.py

def orbital_time(s: npt.ArrayLike, planet: PlanetParameters = earth):
    """
    Calculate orbital time given time parameter, t(s)

    Parameters:
        s: Spinor time parameter (note, called $\tau$ in the paper)
        planet: Planet whose orbit is being analyzed
        e: Optional override for the orbit's eccentricity
    """
    A, B, Om, T_y = _kepler_params(planet)
    return (A**2 + B**2) * s + A * B / Om * sin(2 * Om * s) + T_y / 2

season_event_info(season_value, year)

Return the date and orbit angle for an equinox or solstice in a given year, identified by the season's value as per analemma.geometry.Season.

Source code in src/analemma/orbit.py

def season_event_info(season_value: int, year: int) -> OrbitDateAndAngle:
    """
    Return the date and orbit angle for an equinox or solstice in a given year, identified by
    the season's value as per [analemma.geometry.Season][].
    """
    season_events = _skyfield_season_events(year)
    # S S A W (seasons)
    # 1 2 3 0 (Season enum)
    # 0 1 2 3 (Skyfield)
    skyfield_season_value = (season_value + 3) % 4
    season_event_angles = [pi / 2, 0, 3 * pi / 2, pi]
    return OrbitDateAndAngle(
        season_events[skyfield_season_value].utc_datetime().date(),
        season_event_angles[skyfield_season_value],
    )

spinor_time(t, planet=earth)

Invert t(s), the relationship of orbital time t with the parameter in the spinor treatment of the Kepler problem, s, to give s(t).

Source code in src/analemma/orbit.py

def spinor_time(t: npt.ArrayLike, planet: PlanetParameters = earth):
    """
    Invert t(s), the relationship of orbital time t with the parameter in the spinor
    treatment of the Kepler problem, s, to give s(t).
    """
    s_points, t_points = _finegrained_interp_points(planet)
    return np.interp(t, t_points, s_points)