DEV Community: Rémy Hannequin

Introducing Caelus

Rémy Hannequin — Fri, 31 Oct 2025 10:05:32 +0000

It's not often you finally build something you've been dreaming about for more than a decade.
For me, that moment came today.

While what I'm presenting today is far from complete and still only offers a few features, the way it is built is an achievement I have been dreaming of for many years.

How to compute astronomical data

You open your favourite weather app. It tells you the Sun will rise at 07:34 today, and you trust it. There's no reason not to. But how does it know? Then, you change your location to see how the weather is in another part of the country, and you realise the sunset time is different. Finally, you also realise the app tells you what is the current Moon phase, and when the next full Moon will happen. How does it know what will happen tonight, tomorrow, a year from now?

I wasn't satisfied with simply coding and manipulating astronomical data. I wanted a deeper understanding of the entire process, from how the data is gathered to how the sky fundamentally operates.

I believe astronomical data is universal, pun intended. If scientists managed to understand how the sky works, then this knowledge should be accessible to everyone.

Astronoby

Four years ago, I started working on Astronoby, initially basing my work on many astronomy and astrometry books. More recently, Astronoby was refactored to use raw position data of Solar System objects directly from the best sources available: IMCCE and NASA. The former provides astronomical data for EU countries, while the latter is the US space agency.

Astronoby is open source, so all the code and logic to transform the data, all the algorithms to compute event times, everything is accessible for anyone to understand how it works and spot any mistakes I might have made. It is precise enough to enable any amateur telescope to target planets for the next 50 years, accurate to the second.

Finally, with access to accurate astronomical data based on logic that is open to anyone, I can finally tackle my decade-long goal: build my own astronomical tool and make it accessible to the public. This is Caelus.

Caelus

Caelus is heavily inspired by the amazing Heavens Above. Aside from the fact that Heavens Above currently provides many more features, the main difference is that everything in Caelus is traceable.

If Caelus shows a sunrise time of 07:34, you can trace back to the source data and understand the entire chain of logic that produced it. This includes the algorithm calculating the Sun's position, the code handling the observer's location, and even adjustments for atmospheric effects. Every single piece of data can be checked by the user, making it truly universal.

What's now, what's next

Currently Caelus offers three main webpages that provide most of the general information about tonight for the observer.

You will see information about the Moon

Dedicated information about the Sun

Information about the visibility of planets

Accessible deep-sky objects visible from your location

The two main things I want to focus on now are the documentation and an observation planner.

The first is creating educational content that explains the astronomy behind the data. This will include definitions for common terms, explanations of physical effects (like atmospheric refraction), and a breakdown of the calculations. If the user doesn't know how to read Ruby code but still wants to learn more about how we compute this data and what it means, I want Caelus to help.

The observation planner will be a way to specify what you're interested in, when you're going to observe the sky and where you will be observing from, and Caelus will provide all the right data so that the user can plan their observation night.

I'm also keen to hear from the community and try to produce whatever data people are interested in. I want Caelus to become a useful and educational open source tool.

The entire project is open source and hosted on GitHub: https://github.com/rhannequin/caelus. I welcome feedback, ideas, issues, and pull requests from anyone interested in collaborating.

Cheers! 🌌

Cover picture credits: altargods.com

How to calculate the best days for planetary observation

Rémy Hannequin — Mon, 01 Sep 2025 08:14:56 +0000

A few months ago, I started a tradition: each time I release a new version of Astronoby, I write a blog post showing how Ruby and the amateur astronomy communities can use it to compute astronomical events.

Today, to celebrate v0.8.0, we'll compute the best days for planetary observation. This illustrates Astronoby's new caching and performance improvements, and hints at features coming in future versions.

But first, a bit of context.

What is planetary observation?

In amateur astronomy, when we say planetary, we usually mean three types of bodies:

the planets
the Moon
the Sun

If you jumped because the Moon and the Sun are not planets, you're right, but please hear me out.

Planetary observation differs from deep sky observation. Deep sky objects are very different from planets, usually much fainter. But also they always appear at the same place in the sky at the same time of the year. They represent nebulae, galaxies, stars clusters from our own galaxy or surrounding. On a human timescale, they appear fixed in the sky. Only Earth's revolution around the Sun changes their visibility throughout the year.

Planets on the other hand are very close and very bright, and they have their own motion which is noticeable on a daily basis. The Moon changes phases, while some planets like Mars is best to be observed every two years.

Credits: BBC Sky at Night Magazine

Planetary observation depends on many constantly changing factors.

Sun observation

Let's start with the easiest one on paper, but also the most dangerous. The "astronomer code of conduct" requires me to remind you, even if you might already know it, to never look at the Sun with the naked eye or any equipment, unless you are using a very secure and safe one.

You might think, why would we observe something that can burn the retina? Because, with a dedicated equipment, the Sun is full of great things to see like sun spots, which are cooler areas on the surface, or flares, which are ejections of coronal mass, and many more.

The Sun is easy because it doesn't really move, our perspective of movement is mainly due to the Earth's rotation and revolution.

The Earth's tilt makes the Sun appear more or less high in the sky across the year. If you want to know more about the why and how, check out my previous blog post on the sub-solar point.

Therefore, the best time to observe the Sun is around the summer solstice, June in the northern hemisphere, December in the southern, when it is highest in the sky and visible for the longest time.

Credits: Sky & Telescope

You don't need an astronomy library to know when are the few weeks around the solstice. But you might be interested to know what is the highest altitude the Sun could go as seen from your location. Let's use Astronoby to compute the day of solstice, and then compute the highest altitude of the Sun that day, when it transits. As someone living in France, I'm going to use the June solstice.

# Download an ephemeris if you don't have one yet.
# Check out https://github.com/rhannequin/astronoby/wiki/Ephem

# Load an ephemeris
ephem = Astronoby::Ephem.load("inpop19a.bsp")

# Set up the observer with your GPS coordinates
# For the example: Paris, France
observer = Astronoby::Observer.new(
  latitude: Astronoby::Angle.from_degrees(48.8575),
  longitude: Astronoby::Angle.from_degrees(2.3514)
)

# Compute solstice day
solstice_day = Astronoby::EquinoxSolstice
  .june_solstice(2025, ephem)
  .to_date
# => #<Date: 2025-06-21 ((2460848j,0s,0n),+0s,2299161j)>

calculator = Astronoby::RiseTransitSetCalculator.new(
  body: Astronoby::Sun,
  observer: observer,
  ephem: ephem
)

# Compute time of Sun's transit
transit_time = calculator
  .event_on(solstice_day)
  .transit_time

# Compute Sun' state at this instant
sun = Astronoby::Sun.new(
  instant : Astronoby::Instant.from_time(transit_time),
  ephem: ephem
)

# Get the altitude in degrees
altitude = sun
  .observed_by(observer)
  .horizontal
  .altitude
  .degrees
  .round(2)
# => 64.58

Paris is decently North, so the Sun gets pretty high during the solstice but not as high as observed from the Tropic of Cancer.

Planets observation

Things get a bit more complicated from here. First, because the planets are constantly moving, and so does Earth, so they cannot be expected to be at the same location in the sky every year.

But also, some planets have orbits closer to the Sun than Earth, and some are further, which changes the orientation of looking at them.

Inferior planets

Let's start with the only two planets whose orbits lie closer to the Sun than Earth's: Mercury and Venus.

From our perspective, they are always quite close to the Sun. Sometimes, when they are between the Sun and the Earth, or when the Sun is between them and the Earth, they cannot be seen at all. That's why some cultures call Venus the Morning Star: it's often visible just before sunrise (or just after sunset).

Credits: University of Nebraska-Lincoln

In conclusion, the best moment to observe these planets is a range of a few days or weeks around when the planet is the further from the Sun, from our perspective. This moment is called the greatest elongation. It is the maximal angular distance between the planet and the Sun, as observed from Earth.

We are going to use Astronoby to compute the top 30 days when Venus is the further from the Sun.

Astronoby doesn't have a built-in function (yet) for this, so we're going to do it ourselves. However, it means we are going to compute Venus' position for every day of the year, this can be quite expensive.

That's where I can advertise on Astronoby's latest feature: internal caching. You just have to enable it and Astronoby will internally save some calculations that are redundant, making the total computation faster. Then, we are going to use some trigonometry to compute the daily elongation, and save the 30 largest.

# Assuming you have already downloaded
# and loaded the ephemeris, cf above.

# Enable caching and better performance
Astronoby.configuration.cache_enabled = true

# Iterate over every day of the year
results = (Date.new(2025, 1, 1)..Date.new(2025, 12, 31))
  .each_with_object({}) do |date, results|
    # Let's check elongation at midnight UTC
    time = Time.utc(date.year, date.month, date.day)
    instant = Astronoby::Instant.from_time(time)

    # Compute Sun and Venus positions
    sun = Astronoby::Sun.new(instant: instant, ephem: ephem)
    venus = Astronoby::Venus.new(instant: instant, ephem: ephem)

    # Extract their astrometric unit vectors
    # We normalize the position vectors of Venus and the Sun
    # to compute the angle between them.
    venus_unit_vector = venus
      .astrometric
      .position
      .map { |x| x.m / venus.astrometric.position.norm.m }
    sun_unit_vector = sun
      .astrometric
      .position
      .map { |x| x.m / sun.astrometric.position.norm.m }

    # Compute the elongation angle and store it
    results[date] = Astronoby::Angle.acos(
      Astronoby::Util::Maths.dot_product(
        venus_unit_vector,
        sun_unit_vector
      )
    )
  end

top_days = results
  .sort_by { |_date, angle| -angle }
  .first(30)
  .map { |date, _angle| date }
# => [
#   #<Date: 2025-01-10 ((2460686j,0s,0n),+0s,2299161j)>,
#   ...

In 2025, the greatest elongation of Venus happened on January 10 and the best days to observe Venus were around this day.

Note that if you look at the whole array, you will notice the presence of some days from June, it's because Venus reached its greatest elongation twice in 2025 (East in January, West in June).

Superior planets

Now, superior planets are the ones which orbit is further from the Sun than the Earth's, which include Jupiter, Saturn, Uranus and Neptune.

For these planets, the best time to observe them is at opposition, when Earth is exactly between the Sun and the planet. At that point, the planet is visible all night and reaches transit near midnight, giving the best conditions for observation.

Credits: University of Nebraska-Lincoln

The geometrical and physical definition of a planet's opposition is to have the difference between the Sun and the planet's geocentric apparent ecliptic longitude to equal 180°.

Good news, we also have this information available straight from Astronoby. Let's calculate the day of opposition by finding the day when this angle is the closest to 180°.

# Assuming you have already downloaded
# and loaded the ephemeris, cf above.

# Enable caching and better performance
Astronoby.configuration.cache_enabled = true

# Iterate over every day of the year
results = (Date.new(2025, 1, 1)..Date.new(2025, 12, 31))
  .each_with_object({}) do |date, results|
    # Let's check opposition at midnight UTC
    time = Time.utc(date.year, date.month, date.day)
    instant = Astronoby::Instant.from_time(time)

    # Compute Sun and Saturn positions
    sun = Astronoby::Sun.new(instant: instant, ephem: ephem)
    saturn = Astronoby::Saturn.new(instant: instant, ephem: ephem)

    # Compute the opposition angle
    # This calculates the difference in ecliptic longitude,
    # opposition is when it's 180°.
    angle = sun.apparent.ecliptic.longitude - saturn.apparent.ecliptic.longitude

    # Store the absolute angle reduced by 180° for ease of use
    angle = Astronoby::Angle.from_degrees((angle.degrees.abs - 180).abs)
    results[date] = angle
  end

top_days = results
  .sort_by { |_date, angle| angle }
  .first(30)
  .map { |date, _angle| date }
# => [
#   #<Date: 2025-09-21 ((2460940j,0s,0n),+0s,2299161j)>,
#   ...

In 2025, the opposition of Saturn will happen on September 10 and the best days to observe Saturn were around this day.

And you don't have to trust my math, you can also use Astronoby to get the transit time on the opposition day to ensure Saturn will reach its maximum altitude in the middle of the night:

opposition_date = Date.new(2025, 9, 21)

# Set up the observer with your GPS coordinates
# For the example: Paris, France
observer = Astronoby::Observer.new(
  latitude: Astronoby::Angle.from_degrees(48.8575),
  longitude: Astronoby::Angle.from_degrees(2.3514)
)

calculator = Astronoby::RiseTransitSetCalculator.new(
  body: Astronoby::Saturn,
  observer: observer,
  ephem: ephem
)

# Compute time of Saturn's transit, Paris local time
transit_time = calculator
  .event_on(opposition_date)
  .transit_time
  .localtime("+01:00")
# => 2025-09-22 00:44:30 +0100

Moon observation

Last stop of our journey: The Moon.

You might think this is easy: the Moon is visible almost every day, so observing it should be straightforward. But it's not that simple.

It depends on what you want to observe, when, and in what conditions.

First, the Moon rotates around the Earth in around 29 days, its phase changes every day. Whether you want to observe a full Moon or a thin crescent, the desired phase isn't available every day.

Also, you probably want to observe it when it's high in the sky, visible for a long time. Observing it close to the horizon might make the observation hard due to the atmospheric refraction bending the light.

You might have noticed it, but sometimes the Moon can be seen very high in the sky, and sometimes very close to the horizon. This is because the Moon's orbit is inclined by 5.14° to the ecliptic, in addition to the Earth's axial tilt to orbit of 23.44°. The combination of these tilts constantly changes how the Moon appears in the sky. There are also many small and long-term effects that have an influence but let's not get too complex.

The main thing to understand is that the perfect configuration for you might happen very rarely, that's why computing when it happens with Astronoby is very helpful.

For the example, I'm going to be very selective to show how precise you can get with Astronoby. Let's say I want to know the days when I can observe the Moon in the evening, not too close to the Full Moon because I want to observe craters with grazing light, when it will be high enough in the sky at transit.

To be more technical, this translates into the following requirements:

Altitude at transit at least 45°
Phase angle below 110° so that the Moon doesn't appear too close to the Sun
Illuminated fraction below 0.85 to avoid the Full Moon and the days around it
Transit happens after noon to allow for a nice observation in the evening

Here's how the Moon would look like when meeting all the requirements. Some craters are beautiful during this phase.

Last note before we jump into the code. I'm going to use the tzinfo gem to have an accurate local noon and dates across the year.

require "tzinfo"

# Assuming you have already downloaded
# and loaded the ephemeris, cf above.

# Enable caching and better performance
Astronoby.configuration.cache_enabled = true

# Set up the observer with your GPS coordinates
# For the example: Paris, France
observer = Astronoby::Observer.new(
  latitude: Astronoby::Angle.from_degrees(48.8575),
  longitude: Astronoby::Angle.from_degrees(2.3514)
)

# Use Paris' time zone for accurate time and date conversion
tz = TZInfo::Timezone.get("Europe/Paris")

# Define requirements
minimum_altitude = Astronoby::Angle.from_degrees(45)
maximum_phase_angle = Astronoby::Angle.from_degrees(110)
maximum_illuminated_fraction = 0.85

calculator = Astronoby::RiseTransitSetCalculator.new(
  body: Astronoby::Moon,
  observer: observer,
  ephem: ephem
)

# Compute all the transit times of year 2025
# Thanks to the caching feature, this is quite fast
transit_times = calculator
  .events_between(
    Time.utc(2025, 1, 1),
    Time.utc(2026, 1, 1)
  )
  .transit_times

results = transit_times
  .each_with_object([]) do |transit_time, results|
    date = tz.to_local(transit_time).to_date
    transit_instant = Astronoby::Instant.from_time(transit_time)

    # Compute the Moon's position at transit time
    moon = Astronoby::Moon.new(
      instant: transit_instant,
      ephem: ephem
    )

    # Get the Moon's topocentric altitude
    # Altitude as seen by an observer on Earth
    altitude = moon
      .observed_by(observer)
      .horizontal
      .altitude

    # Define when noon happens in Paris that day
    local_noon = tz
      .local_time(date.year, date.month, date.day, 12, 0, 0)

    # Decide whether or not the requirements are met
    requirements_pass = altitude > minimum_altitude &&
      moon.phase_angle < maximum_phase_angle &&
      moon.illuminated_fraction < maximum_illuminated_fraction &&
      transit_time > local_noon

    results << date if requirements_pass
  end

results.size
# => 28

There are only 28 days in 2025 where I can observe the Moon from Paris in the particular configuration I require.

Of course I was very strict for the sake of the example and the Moon is nice to observe any time it is visible.

In this post, we saw how Astronoby can compute the best days for observing the Sun, planets, and Moon, from solstices and elongations to oppositions and lunar phases. These tools help astronomers plan sessions months in advance. In future releases, Astronoby will support even more events. Until then, I'd love to hear what you'd like to compute with Astronoby!

Cover picture credits: NASA

Where on Earth do shadows disappear? Finding the sub-solar point with Ruby

Rémy Hannequin — Fri, 20 Jun 2025 14:58:49 +0000

There is constantly a place on Earth where you would cast no shadow. Or more exactly, if you were a straight object like a traffic cone, you would not cast any visible shadow.

This happens when you right at the sub-solar point, meaning "right under the Sun". At this instant, the Sun is at the zenith, right above your head, 90° above the horizon.

This place on Earth constantly moves, act fast if you want to snap a strangely fascinating photo. But what exactly is happening?

First, let's have a quick astronomical recap.

Obliquity of the ecliptic

The Earth does its revolution around the Sun on a plane called the ecliptic. But as we know, the Earth also rotates on its axis in about 24 hours, along what we call the equatorial plane.

Our planet is tilted, meaning the ecliptic and equatorial planes are not the same. The angle between them is called the obliquity of the ecliptic. Its current value is 23.44° and it slightly changes over time (millennia), but this part will be for another day.

This tilt is what causes the seasons. As we go through the year, one hemisphere will have longer or shorter days.

But what does this have to do with having the Sun at zenith?

Sub-solar point

Because the Earth is round (if you don't agree, the rest of this article might be a bit less enjoyable) and tilted on its axis (that's the obliquity), we don't all receive the light rays with the same angle.

This also means there's always a place on Earth that receives the sun light perpendicularly to the surface. That's the sub-solar point!

Obviously this point doesn't stay there, let's not forget the Earth's rotation - it's always noon somewhere.

Now let's combine the two pieces of information we have: There is a place where the sunlight hits the surface perpendicularly, and this place changes constantly because of the Earth rotation and its revolution around the Sun.

Finding the exact spot: the theory

Now comes the fun part: finding that elusive spot.

To compute it, we need 2 pieces of information:

the Sun's apparent declination
the longitude of the meridian where the real solar noon happens

Lots of unfamiliar terms here, let's dig a little bit deeper and then we should be good with the theoretical stuff.

Declination

It's very similar as the geographic latitude, which is how North or South the location is compared to the equator. If you extended Earth's equator into space, it would trace the celestial equator. The declination is then how up or down the object is compared to the celestial equator.

Solar noon

Yes, there's a real noon... and a fake one. 😱

Actually, it's not fake, it's mean. Noon is commonly known as mid-day, when the Sun is at its highest in the sky. Mean noon is a standardised time that everyone in the same time zone agrees to use. Solar noon is when the Sun is at it's highest in the sky for a given observer, regardless of its legal time.

Actually finding the exact spot

Finally we can move to the practical part and actually compute the sub-solar point's coordinates.

We're using Astronoby, a library that enables to compute astronomical event and data in pure Ruby. It relies on ephemerides published by official administrations and institutes like NASA/JPL or the IMCCE.

The first action is to download one of these ephemerides, they are the fuel of the library. This needs to be done only once and then the file is stored locally, if you already used Astronoby it might not be necessary.

Astronoby::Ephem.download(
  name: "inpop19a.bsp",
  target: "tmp/inpop19a.bsp"
)

From there, the ephemeris can be loaded and used by Astronoby in any computation.

ephem = Astronoby::Ephem.load("tmp/inpop19a.bsp")

Then, we need to specify what time we want to know where the sub-solar point is. Astronoby abstracts time through a Instant object.

time = Time.utc(2025, 6, 20)
instant = Astronoby::Instant.from_time(time)

Up to this point, we have everything we need to compute the Sun's apparent coordinates.

sun = Astronoby::Sun.new(ephem: ephem, instant: instant)

The sub-solar point's latitude is right there, described by th sun's apparent declination.

sub_solar_latitude = sun.apparent.equatorial.declination

The longitude is slightly more complicated as it is the difference between the apparent right ascension and the Greenwich sidereal time. We didn't explain these two new notions for simplicity's sake, and it might be confusing to compare and angle with a time, but astronomy allows that. Let's just say right ascension is the celestial equivalent of longitude, paired with declination it forms a full coordinate system for the sky. Sidereal time is a time scale that is based on the Earth's rate of rotation measured relative to the fixed stars.

The longitude is the difference between right ascension and sidereal time, normalised in a range between -180° and 180°.

right_ascension = sun.apparent.equatorial.right_ascension
sidereal_time = Astronoby::GreenwichSiderealTime
  .from_utc(time.utc)
  .time
longitude = right_ascension - Astronoby::Angle.from_hours(sidereal_time)

# Normalise to range [-180°, 180°]
sub_solar_longitude = Astronoby::Angle.from_degrees(
  ((longitude.degrees + 180) % 360 - 180)
)

That's it, we have our coordinates!

puts "At #{time}, the sub-solar point is at:"
puts "  #{sub_solar_latitude.degrees.round(2)}°,#{sub_solar_longitude.degrees.round(2)}°"

# Output:
# At 2025-06-20 00:00:00 UTC, the sub-solar point is at:
#   23.43°,-179.62°

But wait, there's more!

You might be aware, but this post is supposedly published the day before June 21, 2025's June solstice. What's a solstice you may say? This is the moment in time where the Sun will reach it's highest or lowest (depending on the observer's hemisphere) altitude in the sky for the year.

The first thing, is that you can compute this instant easily with Astronoby.

Astronoby::EquinoxSolstice.june_solstice(2025, ephem)
# => 2025-06-21 02:42:19 UTC

Remember when we said at the beginning of this post that the Earth's tilt was 23.44°. Well, the second thing I want to show you is that with the beauty of geometry, we can compute this angle back by looking at the sub-solar point's coordinates during solstice!

time = Astronoby::EquinoxSolstice.june_solstice(2025, ephem)
instant = Astronoby::Instant.from_time(time)

# ... same instructions as before

# Output:
# At 2025-06-21 02:42:19 UTC, the sub-solar point is at:
#   23.44°,139.86°

This makes sense as solstices are moments where the Earth's tilt as the most effect on how we receive sunlight.

And guess what, we've actually named the extreme sub-solar latitudes 23.44° and -23.44° and you might have already heard them: Tropic of Cancer and Tropic of Capricorn.

Who would've guessed we could play with astronomy and geography, in Ruby? Now you can!

Check out on Astronoby for even more celestial wonders seen from Earth!

Credits:

[Boost]

Rémy Hannequin — Thu, 19 Jun 2025 21:08:38 +0000

George Kosmopoulos

Jun 17 '25

Check the code health of your rails app

#rails #ruby #gem #health

Comments 1

5 min read

Astronoby v0.7.0: Planets and ephemerides

Rémy Hannequin — Mon, 12 May 2025 10:19:43 +0000

tl;dr at the end

After months of refactoring, Astronoby v0.7.0 is released!

This new version introduces two major changes: Planets and ephemerides.

Quick glossary:

An ephemeris is a software or file that helps generating positions and velocities of celestial bodies
The plural form is ephemerides

Astronoby ❤️ Ephem

5 months ago, Ephem was released, a Ruby library to compute position and velocity vectors of the major planets of the Solar System using SPICE Kernel files (SPK).

Ephem is now v0.3.0 and supports two sources of SPK:

The new version of Astronoby adds Ephem as a dependency to enable the extreme precision of calculating planet coordinates. The Development Ephemeris files from the JPL are, for instance, used for actual spacecraft navigation.

Planets!

While only the Sun and the Moon were supported as celestial bodies in previous Astronoby versions, now all the major Solar System planets can be manipulated, from Mercury to Neptune.

As an observer on Earth, it is possible to compute with great precision where a planet is in the sky, like the Sun and the Moon.

Usage

All the changes are documented in the release notes and the whole library is documented in the repository's Wiki, but let me give you a quick look at what is now possible with Astronoby.

Download and load an ephemeris

The ephemeris is the main source of data for Astronoby, which uses it to compute different kind of positions for different needs.

Ephemerides can be downloaded directly using Astronoby.

Astronoby::Ephem.download(name: "de440s.bsp", target: "data/ephem.bsp")

This will download the file and store it at the specified location. You only need to do it once, you can then load it any time you want to compute data and events with Astronoby.

ephem = Astronoby::Ephem.load("data/ephem.bsp")

Deal with time with `Instant`

Because ephemerides require extra precision and manipulate a different time scale, the new version of Astronoby introduced Astronoby::Instant, a value object to deal with time scales.

You can still use Time and Date, but it is necessary to go through Instant in most calculations.

time = Time.utc(2025, 11, 5)

instant = Astronoby::Instant.from_time(time)

instant.to_time
# => 2025-11-05 00:00:00 UTC
instant.to_date
# => #<Date: 2025-11-05 ((2460985j,0s,0n),+0s,2299161j)>

Instantiate a planet

Once you have an ephem and an instant, you have everything you need to instantiate a Solar System body.

saturn = Astronoby::Saturn.new(ephem: ephem, instant: instant)

From there, you can compute 5 different positions: geometric, astrometric, mean_of_date, apparent and topocentric.

Most users will be interested in apparent and topocentric. The apparent position is corrected and geocentric, it describes how the object would be seen from the centre of Earth, while the topocentric position is centred as an observer on Earth.

A topocentric position can be computed using an observer object that already existed in the previous versions of Astronoby:

observer = Astronoby::Observer.new(
  latitude: Astronoby::Angle.from_degrees(55.8652),
  longitude: Astronoby::Angle.from_degrees(-4.3062)
)

topocentric = saturn.observed_by(observer)

From a topocentric position, you can access horizontal coordinates, which are basically where the object is up/down and left/right in the sky:

horizontal = topocentric.horizontal

horizontal.altitude.degrees.round(4)
# => 20.5112

horizontal.azimuth.degrees.round(4)
# => 226.8246

Rising, transit and setting times

This new version also introduces a new calculator for rising, transit and setting times, performant and supporting ranges of dates.

calculator = Astronoby::RiseTransitSetCalculator.new(
  body: Astronoby::Saturn,
  observer: observer,
  ephem: ephem
)

events = calculator.events_between(
  Time.utc(2025, 11, 1),
  Time.utc(2025, 11, 3)
)

events.rising_times
# => [2025-11-01 15:40:09 UTC, 2025-11-02 15:36:09 UTC]

events.transit_times
# => [2025-11-01 21:19:55 UTC, 2025-11-02 21:15:49 UTC]

events.setting_times
# => [2025-11-01 03:03:49 UTC, 2025-11-02 02:59:38 UTC]

It is still possible to compute Moon phases, twilight times and other data, check out the Wiki to discover all the library's features.

Precision

Not only does this new version enable the user to manipulate Solar System planets, but it also does it with extreme precision.

Topocentric positions from Astronoby have an error margin of less than 10 arc seconds, which corresponds to half the size of Saturn in the sky when it is the closest to Earth.

This precision is by far enough to guide automatically any amateur telescope.

Conclusion

This new version enables Astronoby to enter a new phase: computing any phenomenon with accuracy. Thanks to the help of ephemerides, Astronoby has access to extremely precise data and can use it in algorithms and provide many more features in the future.

While I was limited to the data provided by the few books the library was based on, with ephemerides Astronoby can be autonomous and write its own story.

Don't hesitate to reach out to me here or on GitHub if you have any question or for support to integrate this new version, I will be more than happy to help.

tl;dr

Astronoby v0.7.0 is out
It uses Ephem for extreme accuracy
All planets of the Solar System can be observed from Earth

Credit cover picture: https://gostargazing.co.uk

Introducing Ephem

Rémy Hannequin — Mon, 13 Jan 2025 15:07:36 +0000

tl;dr at the end

I am glad to announce a new Ruby gem on astronomy has been published: Ephem, a library that produces ephemerides from NASA's JPL Development Ephemeris series.

JPL DE4XX

The Jet Propulsion Laboratory (JPL) produces a series of mathematical models called Development Ephemeris (DE), which are used for spacecraft navigation and astronomy.

Such an ephemeris is available through a SPK file, or SPICE Kernel file, a complex binary file containing mathematical parameters to compute the position of Solar System bodies based on time.

The precision delivered by these models is so high that space agencies use it to send spacecraft into space. As an order of magnitude, while the planet Mars is usually located within a range of 50 to 250 million kilometres from the Earth, they are expected to predict its location with a margin error of one kilometre.

These kernels are large and complex and follow an arbitrary format of headers, segments, and coefficients that is unfriendly to beginners. However, once the right coefficients are found and used in the right mathematical formula, barycentric Cartesian coordinates are produced.

Solar System Cartesian coordinates

While we are used to coordinates relative to one's town, country or even the Earth, a DE kernel produces a more absolute series of coordinates, relative to barycentres in the Solar System.

A [barycentre] is the centre of mass of two or more bodies. In the same way, the point of balance for two different objects on a plank depends on their mass, Solar System bodies such as planets or even the Sun are related to their barycentre.

Also, Cartesian coordinates are a 3-D space coordinate system of three axes going through a single origin. Often named x, y, z, these are coordinates where an object is located in space relative to the origin.

Here is an example in more practical words: JPL DE405 enables us to compute the Earth-Moon system's (barycentre) location from the Solar System's barycentre. It can also calculate the Earth's location from the Earth-Moon barycentre. Therefore, through a simple formula, it is possible to compute the Earth's location from the Solar System's barycentre.

How to use these coordinates

Now the question is: What can we do with such raw data?

As we can see these coordinates as absolute, we can then convert them into coordinates that are useful to an observer.

Through a list of trigonometry functions, time and referential conversions and many corrections due to how complex the motion of celestial bodies is, the barycentric coordinates can be converted into topocentric coordinates, also known as geodetic coordinates or even horizontal coordinates. From Ephem's data, it is possible to calculate the altitude and azimuth of a celestial object in the sky for a given time and place on Earth.

The coordinates produced by Ephem can be used to observe the sky with extreme precision.

While Ephem remains its own gem so that the community can access raw data, it will become a dependency of Astronoby to exploit its full potential. Having access to such precise and accurate data will open the door to all possible computations of celestial events with the greatest accuracy available to the public. This will enable Astronoby to truly become a scientific library.

Usage

Here is a common use case of the library:

# Download and store the SPICE binary kernel file
Ephem::IO::Download.call(name: "de421.bsp", target: "tmp/de421.bsp")

# Load the kernel
spk = Ephem::SPK.open("tmp/de421.bsp")

# Define the center and target bodies
center = 0 # Solar system barycenter
target = 5 # Jupiter

# Get the right segment
segment = spk[center, target]

# Get the position and velocity of the target body at a given time
# The time is expressed in Julian Date TDB
time = 2460676.5 # 2025-01-01 00:00 UTC
state = segment.compute_and_differentiate(time)

# Display the position and velocity vectors
puts "Position: #{state.position}"
# => Position: Vector[157123190.6507038, 684298787.8592143, 289489366.7262833]
# The position is expressed in km

puts "Velocity: #{state.velocity}"
# => Velocity: Vector[-1117315.1437825128, 254177.26336095092, 136149.03901534996]
# The velocity is expressed in km/day

What if I want to use it right now?

If you want to use this data immediately and not wait for Astronoby to use it, or if you don't want to rely on Astronoby at all, I got you covered. You will find here a simple Ruby script that computes the position of Saturn on a specified date and time and displays the horizontal coordinates (altitude and azimuth) for an observer on Earth.

Please be aware the results will be approximate. Many tiny corrections should be added to the calculations to compute accurate coordinates. However, for amateur purposes, you will probably find the results accurate enough.

Conclusion

Access to astronomical data with this level of precision is a first time in the Ruby ecosystem. While Ephem is only v0.1, the main usage is covered and accurate data can already be used by the community.

It will help Astronoby reach its full potential and provide more and more features while ensuring accuracy and precision.

Check out the GitHub repository for more usage information and please spread the word so that we can improve the library from feedback.

tl;dr

Ephem v0.1 is out: https://rubygems.org/gems/ephem
Supports computing Solar System's planets Cartesian coordinates
Based on NASA JPL Development Ephemeris files
Accuracy from our point of view is estimated to be below 0.1 arcseconds
Will be used in Astronoby to compute geocentric coordinates
Such accuracy allows for computing any astronomical event and data with confidence

Credit cover picture: https://www.highpointscientific.com/astronomy-hub/post/astronomy-101/8-things-you-need-to-know-about-the-8-planets-in-our-solar-system
Credit JPL center: https://www.jpl.nasa.gov/annual-report/today-tomorrow/
Credit barycenter animation: https://en.wikipedia.org/wiki/Barycenter_(astronomy)
Credit Cartesian coordinates: https://fr.m.wikiversity.org/wiki/Fichier:3D_coordinate_system.svg

Vision for Astronoby - Call for contributors and maintainers

Rémy Hannequin — Wed, 18 Dec 2024 10:13:16 +0000

Genesis

I started building Astronoby in April 2022. My initial ambition was to provide astronomical data in pure Ruby without the need for another programming language or for integrating a third-party API. Astronomical data are by definition universal, they should be accessible and understandable by anyone. They don't belong to a copy-righted website and they should not be accessible only to a niche of programmers.

My goal was to build a Ruby library that could provide data and be self-documented so that the calculations, geometrical and physical notions could be understood by just reading the code.

I was also aiming for decent accuracy, knowing I could never compare with official agencies and institutes or libraries with thousands of contributors for several decades, but good enough to cover most non-professional needs.

Current state

One week ago, version 0.6.0 was released and with it several bug fixes, but also a benchmark to start having a proper look at the accuracy of the library.

The current benchmark could greatly benefit from a more scientific approach and compare more data, but it turns out Astronoby is currently decently accurate. It is a few minutes away from the Sun and the Moon setting and rising times in most cases. The Sun and Moon's coordinates are usually less than half an arc-minute of margin error which represents approximately an angular distance 60 times smaller than the full Moon in the sky (this is not documented at the moment).

However, I have to admit not everything is perfect.

On the self-documented aspect, I have to be honest and say the current implementation does not enable anyone to understand celestial mechanics by just reading the code. There are also some notions I still don't understand myself that I managed to implement anyway.

This certainly didn't help encourage external help from the community. If it requires having access to the same books as me, reading the algorithms and trying to understand how I implemented them to contribute to Astronoby, I know why no one so far submitted a pull request on GitHub for a fix or a new feature. Contributing currently requires more than just Ruby knowledge; it is not newcomers-friendly.

On the accuracy aspect, the library is going to face some limits sooner or later. Because the algorithms are based on books with simplified calculations and definitions, there will always be some imprecision compared to what the current state of knowledge allows us. Even the different books that serve as sources don't always share the same level of precision, or the same methods, which doesn't help find the best and yet still most accurate implementation.

Finally, I regularly discover new sources of knowledge and I understand there are tools and models available for programmers to reach high accuracy but to the cost of having formulas with dozens of terms and absolutely no documentation -readable by a non-PhD person- to help understand celestial mechanics.

Vision for the future of Astronoby

I don't want to lose the battle on both sides. I want Astronoby to be the most accurate and performant astronomy library for Ruby. Still, I also want Astronoby to be a source of knowledge and inspiration for anyone who would like to learn more about how it works or even how to implement it themselves. The Ruby community deserves a Ruby gem that fulfils both of these goals.

I have the feeling Astronoby should be a library made of two distinct yet similar projects.

One could be a project for accuracy. By integrating physical models and with the inspiration of existing important projects like Skyfield or Astropy, this project could focus on providing the most accurate and performant results possible in Ruby. Contributors could help optimise the code, running benchmarks, and covering as many use cases as possible.

Another project could be education. With a focus on the quality of implementation, maintainability, and documentation of both the code and the astronomical concepts, this project would provide "good enough" results with code and documentation that would enable any contributor to help make astronomy more accessible to the community.

To give a practical example, let's take the location of a planet in the sky. To have a decent estimation of coordinates, it is necessary to calculate the mean and true anomaly, to include perturbations of the orbit or perturbations from the observer, and also to define all these concepts so that it's easier to understand what's going on and apply a detailed algorithm. To have a high-quality prediction, there are some models like VSOP or short-period ephemeris like DE 421 to apply formulas and generate extremely high-accuracy results. Both aspects seem relevant and useful to me.

Call for contributors

For too long I have kept the project unwelcoming. It is already difficult to dig into an unknown library and try to help, if it also requires understanding algorithms that I don't always understand myself, there is no way anyone would like to join.

With two projects and two focuses, I want to make Astronoby welcoming to newcomers, either interested in having a high-accuracy astronomy library in Ruby or having an educational astronomy Ruby gem for the community.

To the courageous folks who read until this point, join me and help me build Astronoby for the community. Join a group effort to make Astronoby not only the best astronomy Ruby library but also an efficient, accurate and educational astronomy library for the whole programming community. You can be either a curious developer, a scientist looking for a way to share your knowledge or someone who wants to help the community in any way.

I'm happy to share more ideas I have in mind and to discover ideas from others in the comments section or on GitHub issues. We could also consider having periodic pair-programming sessions work directly on the codebase or Q&A events to learn more about the project and its current state.

Thank you.

Cover picture credit: Pixabay

I counted the number of visible sunsets with Ruby

Rémy Hannequin — Tue, 02 Apr 2024 12:30:24 +0000

I live in a flat in the city. There are buildings around me, so I'm not exposed to direct Sunlight all day long.

My living room is West-facing, but when I moved there I quickly discovered that I wasn't able to enjoy and beautiful sunset every day (if you exclude the fact that the weather in Paris is 💩). This could look counter-intuitive. We learn in school that the Sun rises in the East and sets in the West. So first, let's have a quick recap about the apparent motion of the Sun in the sky.

Axis tilt and seasons

Credit: timeanddate.com

The Earth orbits around the Sun. But the Earth's rotation axis is not perpendicular to the orbits' plane, we call this angle the axial tilt, or obliquity.

The axial tilt is approximately 23.44°. It slightly changes throughout millennials, but this has no impact on the rest of this article.

This angle is responsible for the seasons we encounter on Earth, especially in temperate regions. One hemisphere is always more exposed to the rays of the Sun than the other one, except for the two equinoxes.

The other consequence is the apparent motion of the Sun in the sky. The duration of sunlight exposure is higher in the summer than in the winter, but the Earth always rotates at the same speed: this means the Sun appears higher in the sky during the summer than in the winter.

Evolution of the Sun's azimuth

Credit: Between Solstices by György Soponyaim

The Sun rises almost exactly in the East and almost exactly in the West only during the days of equinoxes. The rest of the year, the azimuth of these phenomena changes over time.

The azimuth is the angle with the North direction on the horizontal plane, or more simply the horizon for an observer. So, the Sun's azimuth is its location on the horizontal axis.

During the summer, the Sun rises and sets more in the North direction. For example, from Paris, during the June solstice, the Sun rises at 51° (or North-East) and sets at 308° (or North-West). During the December solstice, it rises at 126° (or South-East) and sets at 233° (or South-West).

This difference in time is highly influenced by the observer's latitude. At the equator, the azimuth difference between the two solstices is less than 50°, while it's more than 100° in northern locations like Oslo.

Observing the sunset from my living room

After all this explanation, let's get back to business.

From my living room, the South direction is obstructed by a large building. For a significant part of the year, the Sun sets behind it. To be more precise, as long as the sunset azimuth is too much in the South, the Sun is hidden.

As you might know, I am building a Ruby gem called astronoby to compute astronomical phenomena. I recently released a new version that includes information about the Sun, and the sunset azimuth is part of it.

I approximately estimated the border of the building to be 273° on the horizon, when I sit on my couch. In theory, as long as the sunset azimuth is greater than 273°, it should not be hidden by the building.

With astronoby, I can compute the sunset azimuth of each day of the year and select only those within the limit of visibility for me.

require "astronoby"

# Actual coordinates hidden for privacy reasons
observer = Astronoby::Observer.new(
  latitude: Astronoby::Angle.from_degrees(2),
  longitude: Astronoby::Angle.from_degrees(48),
  elevation: 40
)

azimuth_limit = Astronoby::Angle.from_degrees(273)
range = (Date.new(2024, 1, 1)..Date.new(2024, 12, 31))

visible_dates = range.to_a.select do |date|
  noon = Time.utc(date.year, date.month, date.day, 12, 0, 0)
  epoch = Astronoby::Epoch.from_time(noon)
  sun = Astronoby::Sun.new(epoch: epoch)
  setting_azimuth = sun.setting_azimuth(observer: observer)

  setting_azimuth > azimuth_limit
end

dates.first.to_s
# => "2024-03-26"

dates.last.to_s
# => "2024-09-16"

Final results

In theory, if the weather is fine, I should be able to enjoy a nice sunset from my living room couch from the 26th of March to the 16th of September.

The Ruby part of this article might feel a bit disappointing, we spent two sections talking about astronomy while the dates can be computed in barely 10 lines of code. This is however an opportunity for me to brag and share that the astronoby library is powerful.

Knowing the location of the Sun in the sky at any time as an observer anywhere in the world has many applications. From the Sun's exposure to the quantity of energy received from the sky, the astronoby gem opens new doors to Ruby developers.

Cover picture credits: Dramatic City Sunset Sky Pro Vector by John Alberton

astronoby gem v0.1.0 released

Rémy Hannequin — Wed, 28 Feb 2024 10:30:40 +0000

I am proud to announce the first official non-major version of astronoby.

What is Astronoby?

Astronoby is a Ruby gem I created almost two years ago. It is a project meant to enable Ruby developers to compute astronomical data and events.

Why?

Astronomy is universal. Astronomical knowledge doesn't belong to anyone. But it is also hard to grasp, calculate and comprehend.

Almost 10 years ago, I tried to make web apps with astronomical data by scrapping websites. It felt wrong. Why should we rely on proprietary content, often protected with techniques against scrapping, to access such universal data?

I then decided to go deeper and compute the data myself, with existing Python libraries such as PyEphem and Skyfield. I have nothing against Python, and to be honest the Python ecosystem is incredibly more advanced on this kind of scientific topic than Ruby will probably ever be. But I wasn't satisfied not to be able to the the programming languages I wanted to.

This couldn't be the end of it. Is astronomical computing in Ruby really something impossible? When getting started with astronomy, astrometry and programming, do we always have to use tools so powerful and complicated that even professional astronomers use them?

I created Astronoby for these reasons:

Make the formulas and calculations accessible to developers of any level
Bring a new tool to the list of scientific projects to the Ruby community
Have fun and understand how these giant things up there move

What now?

This first official version is only the first step. For the first time, this gem enables developers to compute one data: the Sun's location in the sky based on an observer's location and a time. This only real feature is the reason this first release is not a major one.

time = Time.utc(2023, 2, 17, 11, 0, 0)
epoch = Astronoby::Epoch.from_time(time)

latitude = Astronoby::Angle.as_degrees(48.8566)
longitude = Astronoby::Angle.as_degrees(2.3522)

sun = Astronoby::Sun.new(epoch: epoch)

horizontal_coordinates = sun.horizontal_coordinates(
  latitude: latitude,
  longitude: longitude
)

horizontal_coordinates.altitude.str(:dms)

It is not ground breaking, but it is a first step. Later will come more features associated with the Sun, like sunrise and sunset, equinoxes and solstices. Then, some data related with the Moon will come up like phases and eclipses. After that, the planets will be introduced with multiple events to predict.

The list is long and the remaining work is tremendous. This will take time, especially as one of the main goal is not only to make it work -libraries in other languages already do this- but to make it understandable, maintainable and pleasant to use for anyone enjoying astronomy or Ruby, or both.

It already took me two years of inconsistent spikes to reach this first step. But I believe the new ones will be easier to get to now that some foundations are built.

Contributing

I am aware this kind of project is not easy to contribute to, as most of the new features will be a Ruby translation from books I read by authors like Jean Meeus, J. L. Lawrence, Peter Duffet-Smith and Jonathan Zwart.

The Ruby code though, all the good practices, performance improvements, bux fixes, documentation, are accessible to anyone. This is my first gem, I already made many mistakes I tried to slowly fix and I am aware the gem could already benefit from many improvements unrelated to adding features.

I am also interested in learning more from the community what this library could be. It doesn't have to stop to features covered by a few books. A lot of things can happen, and this is the reason I am releasing this first official version: to make the project finally open to the community.

Please spread the word, share your thoughts, join the project and let's enjoy our excitement for astronomy and Ruby.

Cheers.

Cover picture credits: Between Solstices by György Soponyai

Astronomical computing – Ep 6 - Solving Kepler's equation in Ruby

Rémy Hannequin — Tue, 09 May 2023 09:46:44 +0000

One of the main goals of this series is to understand the math behind astronomical events calculations, and translate them into simple and maintainable code.

In order to forecast where a planet will be visible in the night sky at a given moment we need to rely on topics that we haven’t yet covered such as orbital elements.

Unlike far away stars which can be considered motionless^[1], planets change position on the celestial sphere from one night to the other. Even the Sun, from our point of view, seems to be slowly moving alongside the zodiac belt throughout the year. Our first step towards finding out where a planet will be visible in the sky is to calculate its equatorial coordinates.

Kepler's second law states that the orbital speed of an object varies throughout its orbit. With an elliptic orbit and a changing speed, calculating the planet's location on its orbit for a known date is not a linear process.

We won't cover the calculations needed to compute the coordinates of a planet in this article. We will make sure to cover those in a later installment though. Still, this is an essential step to understand a planet’s motion on its orbit.

Some definitions

Almost all of the following content can be traced to one single source: Celestial Calculations, A Gentle Introduction to Computational Astronomy, by J. L. Lawrence, editions MIT Press.

This is not the first time I mention this book as it is one of my main sources for this series, and I cannot be more grateful of the incredible work of communication and explanation the author has produced to make all these phenomenons accessible.

Section of ellipse showing eccentric and true anomaly^[2]

The illustration above will be useful to understand the different definitions:

The red ellipse is the true anomaly
The blue circle is the mean anomaly
F is the main focus
C is the geometrical center of both the red and blue circles
f is the true anomaly
E is the eccentric anomaly

True anomaly

We can consider it as being the angle between the body's current position and its periapsis, as seen from the main focus of the ellipse. I can imagine how this definition is not helping, so let's go a bit further.

The main focus of the ellipse is basically the point around which the body rotates. When talking about planets of the Solar System, this point can be considered as the Sun^[3].

The periapsis is the nearest point on the orbit to the main focus. For planets, it's called perihelion and is the point where it's the closest to the Sun.

In other words, the true anomaly is the angular distance between the current position and a fixed point, as seen from the Sun.

Mean anomaly

As a planet's orbit is an ellipse and its speed changes constantly depending on how close it is from the Sun, it's pretty difficult to calculate the true anomaly.

To tackle this non trivial task, scholars imagined a fictive orbit, circular, on which the body would move with a constant speed, called the mean orbit. The real orbit would share its geometrical center and its periapsis with the mean orbit's center, which means the body would have the same orbital period on both orbits.

The mean orbit is the angle, from the fictive circle's center, between the body's current position and the periapsis.

The mean orbit is considerably easier to use. A body's mean position on its mean orbit with an orbital period of n sidereal days, after t days from a starting point (the periapsis), can be described this way:

mean_anomaly = 360 / n * t

Eccentric anomaly

The eccentric anomaly is another angular parameter that helps defining the position of a body on its orbit.

It is the projection of the body's current position on the mean orbit, parallel to the line that goes through the mean orbit's center and the periapsis.
Technically we don't use the mean orbit but a circle called the auxiliary circle, but they're the same circle.

Knowing the orbital eccentricity, and the eccentric anomaly, the true anomaly can be defined as:

true_anomaly = arc_tangent(
  square_root(
    (1 + orbital_eccentricity) / (1 - orbital_eccentricity)
  ) * tangent(eccentric_anomaly / 2)
) * 2

Okay, yes, this is getting weirder and weirder but please bear with me, we're getting closer.

Kepler's equation

Kepler discovered that the eccentric and the mean anomalies are related by the following equation:

mean_anomaly = eccentric_anomaly -
  orbital_eccentricity * sine(eccentric_anomaly)

where the mean and eccentric anomalies are expressed in radians.

Knowing this relation, all we have to do is to use the orbit’s orbital eccentricity along with the mean and eccentric anomalies to finally calculate the true anomaly. The calculation for the orbital eccentricity of an orbit is not covered in this article, so we can consider that we already know its value.

Iterations

Unfortunately, Kepler's equation is a transcendental equation, there is no known closed-form solution for it. To compute the eccentric anomaly, we must resolve the equation iteratively and get the most accurate approximation as possible.

If two consecutive iterations have a difference inferior to a pre-defined accuracy parameter, then we can use the last iteration' solution as the most accurate approximation.

Instead of resolving the equation defined above, we will resolve another one, defined as the Newton/Raphson method, which is known to be more performant and more accurate.

Iteratively, we will resolve the following equation:

eccentric_anomaly(i) = eccentric_anomaly(i-1) - (
  (
    eccentric_anomaly(i-1) -
      orbital_eccentricity * sin(eccentric_anomaly(i-1)) -
      mean_anomaly
  ) / (
    1 -
    orbital_eccentricity * cosine(eccentric_anomaly(i-1))
  )
)

where i represents the following the current iteration.

For the first iteration, where i is 0, we will have the following fixed solution:

mean_anomaly if the orbital eccentricity is 0.75 or less,
π if the orbital eccentricity is more than 0.75.

Finally, we will also introduce a maximum number of iterations, as some configuration may lead to values that never converge, and to prevent from a too small accuracy parameter.

Implementation in Ruby

Finally the part that actually matches the article's title.

The implementation is actually quite simple and uses recursion. This is one solution, other interesting implementations might be available out there.

def kepler_equation_newton_raphson_iteration_scheme(
  mean_anomaly_in_radians,
  orbital_eccentricity,
  precision,
  maximum_iteration_count,
  current_iteration = 0,
  solution_on_previous_interation = nil
)
  solution = if current_iteration == 0
    if orbital_eccentricity <= 0.75
      mean_anomaly_in_radians
    else
      Math::PI
    end
  else
    solution_on_previous_interation - (
      (
        solution_on_previous_interation -
        orbital_eccentricity *
        Math.sin(solution_on_previous_interation) -
        mean_anomaly_in_radians
      ) / (
        1 -
        orbital_eccentricity *
        Math.cos(solution_on_previous_interation)
      )
    )
  end

  if solution_on_previous_interation
    delta = solution - solution_on_previous_interation
    puts "E#{current_iteration} ≈ #{solution}, " \
    "Δ = #{delta.ceil(6).abs}"
  else
    puts "E#{current_iteration} ≈ #{solution}"
  end

  if current_iteration >= maximum_iteration_count ||
      (
        solution_on_previous_interation &&
        delta.abs <= precision
      )
    return solution
  end

  kepler_equation_newton_raphson_iteration_scheme(
    mean_anomaly_in_radians,
    orbital_eccentricity,
    precision,
    maximum_iteration_count,
    current_iteration + 1,
    solution
  )
end

Nothing fancy here, we compute the current iteration' solution from the previous iteration' solution, unless it is the first iteration. Recursion helps computing each iteration and passing the previous iteration' solution.
As long as we have two solutions to compare, we check if their difference is significant, and if we have some iterations left before having to stop manually.
I also added some logs to have a better view of what's happening.

To illustrate this implementation, let's solve Kepler's equation with the Newton/Raphson method for an orbital eccentricity of 0.5 and a mean anomaly of 24.742896°.

First we have to convert the mean anomaly from degrees to radians:

radians = degrees * π / 180
        = 0.431845

Then, we set the accuracy parameter to 0.000002 radians, which is around 0.4 arc seconds.
We set the maximum number of iterations to 10.

Finally, we convert the solution from radians to degrees.

kepler_equation_newton_raphson_iteration_scheme(
  0.431845,
  0.5,
  2e-06,
  10
) * 180.0 / Math::PI

# E0 ≈ 0.431845
# E1 ≈ 0.8151984038210852, Δ = 0.383354
# E2 ≈ 0.7856420972684516, Δ = 0.029556
# E3 ≈ 0.7853985310762203, Δ = 0.000243
# E4 ≈ 0.7853985148507631, Δ = 0.0

# => 45.00002013679163

The expected answer is 45°, which differs of only 0.073" from the solution we computed.

Conclusion

The true anomaly is one step towards our ability to compute where a planet is supposed to appear on the celestial sphere for a given time from a given place on Earth. We are one step closer to this goal and we already have the Ruby implementation ready.

This also helped us understand why we need this value and in what way it is not straightforward to compute it.

^[1] This is absolutely false as everything in space moves, all the time, depending on the referential. Far away stars also rotate around the galaxy's center. The only difference with planets is that their motion appear extremely slow from our point of view, so slow that we often looks like they are fixed compared with short periods of time such as a Human life.
^[2] Source: Wikimedia Commons. Author: Brews ohare.
^[3] In reality the planets orbit around the Solar System barycenter, which moves through time due to the different positions of planets. It is always located between the Sun's center and its corona, so we usually consider the Sun itself as the barycenter.

Astronomical computing – Ep 5 - Locate an object in the sky

Rémy Hannequin — Tue, 02 Aug 2022 11:40:00 +0000

I am still reading multiple books at the same time, as they all provide different knowledge, methods and levels of detail.
I recently finished Chapter 5 from Celestial Calculations, "Stars in the Nighttime sky".

This chapter describes a very interesting method to locate an object from its equatorial coordinates, to its horizontal coordinates from an observer's point of view, and I would like to detail it here and explain all the steps. This methods is particularly interesting as it mixes many important notions, and it is one primary step to compute, in the future, planets locations in the night sky.
Objects from the solar system move, their position in space varies from one day to another from our point of view and in real life. But, for a fixed time, they have fixed equatorial coordinates, which means we are going to be able to locate them in the sky. This is what we are going to do.

Now let's forget about solar system objects and use ones that have fixed equatorial coordinates at (almost) all time: distant stars. Actually, all stars except the Sun, as they are all very distant from us, so distant that their motion cannot be perceived to us in such a short period of time as a human lifetime.

With this very long introduction, here are the steps to convert equatorial to horizontal coordinates:

Convert local civil time (LCT) to universal time (UT)
Convert universal time to Greenwich sidereal time (GST)
Convert Greenwich sidereal time to local sidereal time (LST)
Convert right ascension to hour angle
Convert equatorial to horizontal coordinates

We are now going to explain all the new terms and why these steps are necessary. Also, it will be the opportunity to see how to do so in Ruby, using the standard library or not.

LCT to UT

Local civil time, or LCT, is the time from an observer's watch. It is the legal time from their current time zone.

Universal time, or UT, used to be also knows as the Greenwich mean time, or GMT. It is the solar mean time at The Royal Observatory in Greenwich, London. The solar mean time is just a way to define time with days having the same amount of time. Before time zones, it used to be the way to define time in geographical areas such as cities or regions.

Converting LCT to UT is actually just a step towards having the local sidereal time, which will be defined below.

# It is 10:10:05 pm on June 26th, 2022, in Paris, France

local_civil_time = Time.new(2022, 6, 26, 3, 10, 5, "+02:00")
# => 2022-06-26 03:10:05 +0200

universal_time = local_civil_time.utc
# => 2022-06-26 01:10:05 UTC

UT to GST

Sidereal time is a time scale that is based on Earth's rate of rotation measured relative to fixed stars. Because of Earth's revolution around the Sun, it takes a few minutes longer to have the Sun at culmination than a distant star.

As equatorial coordinates are a way to define the position of an object regardless of the observer's position or time, it is linked to the apparent motion of stars, on the one of the Sun. Our usual time is based of the motion of the Sun, therefore we need to have a conversion of units.

Converting UT to GST is a multiple steps process with constants defined in Celestial Calculations by J. L. Lawrence in chapter 3.10 (see Ruby example below). It also requires calculating the universal time's associated Julian day number. In the previous posts of this series if discovered that a custom implementation of the Julian day calculation fast more performant than the one provided by the Ruby standard library. However, for the sake of simplicity in this post, we will use Date#ajd.

Let's calculate gst, the Greenwich sidereal time, from a defined universal time:

require "bigdecimal"
require "date"

universal_time = Time.new(2022, 6, 26, 3, 10, 5, "+02:00").utc
jd = universal_time.to_date.ajd.to_f
jd0 = Time.utc(universal_time.year, 1, 1, 0, 0, 0).to_date.ajd - 1
days = jd - jd0
t = (jd0 - BigDecimal("2415020")) / BigDecimal("36525")
r = BigDecimal("6.6460656") +
      BigDecimal("2400.051262") * t +
      BigDecimal("0.00002581") * t * t
b = 24 - r + 24 * (universal_time.year - 1900)
t0 = BigDecimal("0.0657098") * days - b
ut = universal_time.hour +
       universal_time.min / BigDecimal("60") +
       universal_time.sec / BigDecimal("3600")
gst = t0 + BigDecimal("1.002738") * ut

gst += 24 if gst.negative?
gst -= 24 if gst > 24

gst.to_f
# => 19.444822123178245

GST to LST

As GST defines the sidereal time at Greenwich, we need to convert if to the observer's local sidereal time, or LST, in order to take into account the observer's longitude.

To convert GST to LST, we only need to convert GST to a decimal format (already the case in the Ruby example above) and add a time zone adjustment defined by:

adjustment = (observer's longitude in decimal degrees) / 15

Let's take the example of an observer at the beautiful beach of Étretat, France.

gst = BigDecimal("14.496844123178246")
longitude = BigDecimal("0.20271537957527094")
adjustment = longitude / 15
lst = (gst + adjustment).to_f
# => 19.4583364818166

Right ascension to hour angle

Converting the right ascension to hour angle is the final step to link the observer's longitude and the right ascension.

The hour angle is a measure of how long it has been since the object crossed an observer's meridian. It varies both with time of day (relative to a sidereal location, explaining why we need the LST) and an observer's location.

It is a time measurement, while the right ascension is an angular measurement, but both can be expressed in hour-minute-second (HMS) format. Therefore, the formula that links them is the following:

hour_angle  = lst - right_ascension

On June 26th, 2022, at 3:10:05 AM, Saturn's right ascension was Right 21h 49m 08.6s. We can compute in Ruby the hour angle with the following algorithm:

lst = BigDecimal("19.4583364818166")
hour = 21
minute = 49
second = 8.6

hour_angle = lst - (
  hour +
    Rational(minute, 60) +
    Rational(second, 3600)
  )

hour_angle += 24 if hour_angle.negative?
hour_angle -= 24 if hour_angle > 24

hour_angle.to_f
# => 21.63928092626104

Equatorial to horizontal coordinates

The final step, having horizontal coordinates to be able, as an observer, to detect the position of the object, on the celestial sphere.

Not much to say here as converting these two coordinates systems is actually essentially trigonometry.
We have all the data we need:

declination
hour angle
latitude

One thing to keep in mind is that trigonometry functions work with angles in radians. Here, we are mostly manipulating angles in HMS format or degrees, so we will need to convert degrees into radians each time we'll use cosine or sine functions.

Let's keep our previous data. Note that Saturn's declination associated to our data was -14° 26' 57.4".

declination_degree = -14
declination_minute = 26
declination_second = 57.4
hour_angle = BigDecimal("21.63928092626104")
latitude = BigDecimal("49.70911954641343")

to_radians = Math::PI / 180
to_degrees = 180 / Math::PI
hour_angle_degree = hour_angle * 15
hour_angle_radians = hour_angle_degree * to_radians
declination_decimal = (
  declination_degree +
    declination_minute / 60.0 +
    declination_second / 3600.0
).to_f
declination_radians = declination_decimal * to_radians
latitude_radians = latitude * to_radians
t0 = Math.sin(declination_radians) *
       Math.sin(latitude_radians) +
       Math.cos(declination_radians) *
       Math.cos(latitude_radians) *
       Math.cos(hour_angle_radians)
h = Math.asin(t0) * to_degrees
h_radians = h * to_radians
t1 = Math.sin(declination_radians) -
       Math.sin(latitude_radians) *
       Math.sin(h_radians)
t2 = t1 / (Math.cos(latitude_radians) * Math.cos(h_radians))
a = Math.acos(t2) * to_degrees
sin_hour_angle = Math.sin(hour_angle_radians)

a = 360 - a if sin_hour_angle.positive?

a.to_f
# => 19.490674347775652

h.to_f
# => 143.30560194994217

Conclusion

On June 26th, 2022, at 03:10:05 AM, at Étretat, France, planet Saturn was visible in the night sky by looking 19° above the horizon and 143° from the North cardinal point.

I made a public gist to sum up all our algorithm. Note that this is only an example gist, it is not optimized nor considered as good code. But it works.

Of course this algorithm only works when we have fixed equatorial coordinates. Stars, for example, have "fixed" equatorial coordinates as they are so far from us that their motion is barely perceivable. Planets however are very much closer to us, their position in the Solar System changes all the time. From one day to another, Saturn's position would have changed a little bit from our point of view, and its equatorial coordinates would have been slightly different, so would have its horizontal coordinates.

What next

This algorithm works and the results are really close to what I compared with several simulation software. However it doesn't include multiple corrections that we'll take time to explain in future posts, such as:

atmospheric refraction, present especially with objects close to the horizon, like planets
precession, which slowly changes the celestial poles positions
nutation, which is another phenomenon that changes the celestial poles positions as well

Also, the final goal is still to build a nice and maintainable Ruby library. The current algorithm will need to be improved and developed with Ruby objects, ready to welcome future features.

Astronomical computing – Ep 4 - Julian day

Rémy Hannequin — Mon, 23 May 2022 12:11:53 +0000

Now that we are managing angles in degrees or radians, we have one of the elementary blocks of the library.

One other important notion is time. Time is complicated, we are going to talk about it a lot in this series. But first, we can tackle a time notion in astronomy that is going to follow us all the way.

Julian day number

Most of us now use the Gregorian calendar^[1] that was introduced on October 1582 and starts on the 15th of October, 1582. Before this day, the Julian Calendar was used. They are both very similar calendars, but they space leap years differently, which makes the Gregorian calendar more accurate regarding the real motion of the Earth around the Sun.

Other calendars still exist, like the Hebrew calendar, the Hijri calendar. Some existed but are not used anymore, like the French Revolutionary calendar or the Maya calendar.

We didn't talk about leap years, time zones or leap seconds but we can understand from here that time is a complicated, if not ambiguous, notion to grasp. While ambiguity is not something that we want to deal with in astronomy and science in general.

Therefore, the Julian day^[2] was introduced by Joseph Scaliger in 1583 to manipulate time as a universal entity. Its definition is the number of days that have elapsed since noon at Greenwich, England, on January 1, 4713 BC. Usually, the Julian day is written as a decimal number, with the decimal part being the fraction of the day that passed. The Julian day is the Julian day number plus the fractional day.

Using the Julian day is very useful in astronomical calculations, especially when comparing two different dates, in order to know the time that separates them.

Converting a date and time to a Julian day

Algorithm

In Astronomical Algorithms, Jean Meeus details a simple formula^[3] to calculate the Julian day associated to a date. The following method works for any positive Julian day, which means any date starting from January 1, 4713 BC.

Given a date, sliced in year (Y), number of the month (M), day of the month (D), hour (h), minute (m) and second (s).

If M > 2, then Y and M stays the same.
If M equals 1 or 2, then Y equals Y - 1 and M equals M + 12.

If the date is on the Gregorian calendar, meaning equal of greater than 15th October 1582, then:

A = INT(Y / 100)
B = 2 - A + INT(A/4)

With INT() being a truncation function so that it represents the whole part of the number.

If the date is on the Julian calendar, then A is not necessary and B equals 0.

Then, the Julian day number is defined by the following formula:

INT(365.25 (Y + 4716)) + INT(30.6001 (M + 1)) + D + B - 1524.5

The Julian day, which is the Julian day number plus the fraction of the day that passed, is the same formula, replacing D as such:

D = D + h / 24 + m / 1440 + s / 86400

One notion to keep in mind is that the Julian day dating started at noon, which means a time at noon will give a whole Julian day, while any other time will give a Julian day with decimals.

Also, the dating started at noon at Greenwich, so any time defined on another time zone than UTC will need to be converted into UTC before being converted into a Julian day. But we will talk more about time zones in future posts.

Implementation in Ruby

`DateTime#jd` and `DateTime#ajd`

The Ruby standard library^[4] provides two methods to compute the Julian day number and the Julian day.

DateTime#jd returns the corresponding Julian day number of a date_time. It is the "chronological" Julian day number, which means if suffers from half a day advance from the original Julian day number. Also, it doesn't include the fraction of time passed in the day. To have the original Julian day from DateTime#jd, it is required to add a few instructions to it:

DateTime.new(2022, 4, 1, 13, 30, 0).jd +
  13 / 24.0 + 30 / 1440.0 - 0.5
# => 2459671.0625

DateTime#ajd returns the "astronomical" Julian day, which is the original one. However, it returns it as a fractional number, so it needs to be converted to be displayed in a way that is understandable to the human eye:

ajd = DateTime.new(2022, 4, 1, 13, 30, 0).ajd
# => (39354737/16)

ajd.to_f
# => 2459671.0625

Using the most performant method

Even if it is tempting to use something ready out of the box, what the Ruby standard library provides us is not perfect. It needs complementary work, or conversion, and this may sometimes reflect reduced performance.

I developed a simple benchmark to compare DateTime#jd with adjustments, DateTime#ajd with conversion, and the custom algorithm provided by Jean Meeus.

require "benchmark"
require "bigdecimal"
require "date"

year = 1957
month = 10
day = 4
hour = 19
minute = 28
second = 34

iterations = 10_000_000

Benchmark.bm do |benchmark|
  benchmark.report("DateTime#jd") do
    iterations.times do
      DateTime.new(year, month, day, hour, minute, second).jd +
        hour / 24.0 + minute / 1440.0 + second / 86400.0 - 0.5
    end
  end

  benchmark.report("DateTime#ajd") do
    iterations.times do
      DateTime
        .new(year, month, day, hour, minute, second)
        .ajd
        .to_f
    end
  end

  benchmark.report("Meeus") do
    iterations.times do
      if month > 2
        m = month
        y = year
      else
        m = month + 12
        y = year - 1
      end

      if (
        year > 1582 ||
          (year == 1582 && month > 10) ||
          (year == 1582 && month == 10 && day >= 15)
      )
        a = y / 100
        b = 2 - a + a / 4
      else
        b = 0
      end

      julian_day_number = (365.25 * (y + 4716)).floor +
        (30.6001 * (m + 1)).floor + day + b - 1524.5
      time_of_day = hour / 24.0 + minute / 1440.0 + second / 86400.0

      julian_day_number + time_of_day
    end
  end
end

#                user       system     total      real
# DateTime#jd    2.766334   0.003965   2.770299   (2.770522)
# DateTime#ajd   5.054064   0.000000   5.054064   (5.054338)
# Meeus          1.978400   0.000000   1.978400   (1.978430)

I also ran this benchmark on multiple computers and Ruby versions, with very similar results. Thanks to my teammates at Getaround for their help.

The custom algorithm from Jean Meeus is the most performant, it is the one I am going to implement in astronoby.

Converting a Julian day to a date and time

Algorithm

Still from Jean Meeus' book^[5], given J the Julian day augmented by 0.5.

Z = INT(J) and F = DEC(J)

With DEC() being a truncation function so that it represents the decimal part of the number.

If Z < 2299161, then we define A as the same value as Z. If Z is greater or equal to 2299161, then we define A as such:

α = INT((Z - 1867216.25) / 36524.25)
A = Z + 1 + α - INT(α / 4)

We can now define B, C, D and E as such:

B = A + 1524
C = INT((B - 122.1) / 365.25)
D = INT(265.25 * C)
E = INT((B - D) / 30.6001)

Finally, we have the day of the month, as a decimal number, equal to B - D - INT(30.6001 * E) + F.

The number of the month (M) is equal to E - 1 if E < 14, and equal to E - 13 if E is equal to 14 or 15.

The year number is equal to C - 4716 if M > 2, and equal to C - 4715 if M is equal to 1 or 2.

Implementation in Ruby

`DateTime::jd`

Similarly, DateTime provides a class method jd to instantiate Datetime from a chronological Julian day. Adding 0.5 makes it astronomical, and will result into the accurate date and time associated.

DateTime.jd(2436116.31 + 0.5)
# => #<DateTime: 1957-10-04T19:26:24+00:00 ((2436116j,69984s,4828n),+0s,2299161j)>

Using the most performant method

I developed a second benchmark to compare DateTime#jd with a custom implementation.

require "benchmark"
require "bigdecimal"
require "date"

j = BigDecimal("2436116.31")

iterations = 1_000_000

Benchmark.bm do |benchmark|
  benchmark.report("DateTime::jd") do
    iterations.times do
      DateTime.jd(j + 0.5)
    end
  end

  benchmark.report("Meeus") do
    iterations.times do
      julian_day = j + 0.5
      z = julian_day.floor
      f = julian_day.abs - julian_day.floor

      if z < 2299161
        a = z
      else
        aa = ((z - 1867216.25) / 36524.25 ).floor
        a = z + 1 + aa - aa / 4
      end

      b = a + 1524
      c = ((b - 122.1) / 365.25).floor
      d = (365.25 * c).floor
      e = ((b - d) / 30.6001).floor

      day = b - d - (30.6001 * e).floor
      month = e < 14 ? e - 1 : e - 13
      year = month > 2 ? c - 4716 : c - 4715

      hour = f * 24
      minute = (hour - hour.floor) * 60
      second = ((minute - minute.floor) * 60).floor

      DateTime.new(
        year,
        month,
        day,
        hour.floor,
        minute.floor,
        second.floor,
      )
    end
  end
end

#               user       system     total      real
# DateTime::jd  6.562874   0.000000   6.562874   (6.563703)
# Meeus         6.316613   0.000000   6.316613   (6.317376)

This time, the custom implementation is only slightly more performant than the standard one. I would normally use the standard implementation for such a small performance difference, however the custom one will let me have access to each part of the date time (year, month, ...) without having to extract it afterwards, and it will probably be useful in the future.

Conclusion

I am really satisfied with this research as one of my goal with this series is to both implement and explain astronomical calculations. By having custom algorithms being as performant or more performant than those implemented in the standard library, it is a logical choice to implement the custom ones and therefore to have an explanatory source code.

However, it takes a lot of time to translate my research into code, and then write about the theory and the Ruby implementation. It already takes me quite some time to understand these notions, try them and add them to astronoby. Duplicating the Ruby code into articles is going to take too much time, and will probably be a bit difficult to read.
In future articles, I will try to find something in the middle, with Ruby code only when it's appropriate and necessary.

Happy hacking

Sources

^[1] Gregorian calendar
^[2] Julian day
^[3] Astronomical Algorithms by Jean Meeus, Chapter 7, page 59
^[4] Ruby DateTime
^[5] Astronomical Algorithms by Jean Meeus, Chapter 7, page 63

DEV Community: Rémy Hannequin

Introducing Caelus

How to compute astronomical data

Astronoby

Caelus

What's now, what's next

How to calculate the best days for planetary observation

What is planetary observation?

Sun observation

Planets observation

Inferior planets

Superior planets

Moon observation

Where on Earth do shadows disappear? Finding the sub-solar point with Ruby

Obliquity of the ecliptic

Sub-solar point

Finding the exact spot: the theory

Declination

Solar noon

Actually finding the exact spot

But wait, there's more!

[Boost]

Check the code health of your rails app

Astronoby v0.7.0: Planets and ephemerides

Astronoby ❤️ Ephem

Planets!

Usage

Download and load an ephemeris

Deal with time with Instant

Instantiate a planet

Rising, transit and setting times

Precision

Conclusion

tl;dr

Introducing Ephem

JPL DE4XX

Solar System Cartesian coordinates

How to use these coordinates

Usage

What if I want to use it right now?

Conclusion

tl;dr

Vision for Astronoby - Call for contributors and maintainers

Genesis

Current state

Vision for the future of Astronoby

Call for contributors

I counted the number of visible sunsets with Ruby

Axis tilt and seasons

Evolution of the Sun's azimuth

Observing the sunset from my living room

Final results

astronoby gem v0.1.0 released

What is Astronoby?

Why?

What now?

Contributing

Astronomical computing – Ep 6 - Solving Kepler's equation in Ruby

Some definitions

True anomaly

Mean anomaly

Eccentric anomaly

Kepler's equation

Iterations

Implementation in Ruby

Conclusion

Astronomical computing – Ep 5 - Locate an object in the sky

LCT to UT

UT to GST

GST to LST

Right ascension to hour angle

Equatorial to horizontal coordinates

Conclusion

What next

Astronomical computing – Ep 4 - Julian day

Julian day number

Converting a date and time to a Julian day

Algorithm

Implementation in Ruby

DateTime#jd and DateTime#ajd

Deal with time with `Instant`

`DateTime#jd` and `DateTime#ajd`

`DateTime::jd`