Skip to content

Instantly share code, notes, and snippets.

@arodland
Last active August 18, 2023 02:12
Show Gist options
  • Save arodland/9c2341dcda8214e6c0badc4a6abb7749 to your computer and use it in GitHub Desktop.
Save arodland/9c2341dcda8214e6c0badc4a6abb7749 to your computer and use it in GitHub Desktop.
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()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment