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

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:

1. Convert local civil time (LCT) to universal time (UT)
2. Convert universal time to Greenwich sidereal time (GST)
3. Convert Greenwich sidereal time to local sidereal time (LST)
4. Convert right ascension to hour angle
5. 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.