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