planets.py

import datetime as dt
import pytz

from skyfield import api, jpllib, vectorlib, timelib
from skyfield.toposlib import wgs84


def determine_object_visibilities(
        object: vectorlib.VectorSum |jpllib.ChebyshevPosition,
        observer: vectorlib.VectorSum,
        observing_times: timelib.Time
):
    observations = observer.at(observing_times).observe(object)

    alt_angles, *_ = observations.apparent().altaz()

    positive_altitudes = [alt >= 0 for alt in alt_angles.degrees]
    return positive_altitudes


def main():
    load = api.Loader(".")
    ephemeris = load("de422.bsp")
    timescale = load.timescale()

    earth = ephemeris['earth']
    other_planets = [
        ephemeris['mercury'],
        ephemeris['venus'],
        ephemeris['mars'],
        ephemeris['jupiter_barycenter'],
        ephemeris['saturn_barycenter'],
        ephemeris['uranus_barycenter'],
        ephemeris['neptune_barycenter']
    ]
    sun = ephemeris['sun']

    tz = pytz.timezone("UTC")

    start_time = tz.localize(dt.datetime(1001, 1, 1))
    end_time = tz.localize(dt.datetime(2001, 1, 1))

    t0, t1 = timescale.from_datetimes([start_time, end_time])

    observation_count = 1_000_000
    observation_times = timescale.linspace(t0, t1, observation_count)

    longitude, latitude, elevation = 0, 0, 0
    observer = earth + wgs84.latlon(longitude, latitude, elevation)

    sun_visible = determine_object_visibilities(sun, observer, observation_times)

    planetary_visibilities = tuple(
        determine_object_visibilities(planet, observer, observation_times)
        for planet in other_planets
    )

    all_other_planets_visible = tuple(
        all(visible_planets) for visible_planets in zip(*planetary_visibilities, sun_visible)
    )

    visibility_count = sum(all_other_planets_visible)
    sun_count = sum(sun_visible)

    print(f"Observing latitude: {abs(latitude)}°{'N' if latitude >= 0 else 'S'}\n"
          f"Observing longitude {abs(longitude)}°{'E' if longitude >= 0 else 'W'}\n"
          f"Starting time: {start_time:%Y-%m-%d %H:%M %z}\n"
          f"Ending time: {end_time:%Y-%m-%d %H:%M %z}\n"
          f"Number of observations with sun visible: {sun_count:,}\n"
          f"Observations with sun and all planets visible: {visibility_count:,}\n"
          f"Percentage of total: {visibility_count / sun_count:5.3%}")


if __name__ == '__main__':
    main()