Calculation of Julian day is off for negative dates

Question

I bought the book Astronomical Algorithms by Jean Meuss and I am reading it in order to build a small library I will use for another project. I made a program in Rust to calculate the JD given a date see here, this program uses the following formula from the book:

Given Y as year, M as month and D as day;

• if M > 2 Y and M are unchanged
• if M == 1 or 2 replace Y with Y-1 and M with M+12

If the date is from a Gregorian calendar then:

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

else

• B = 0

So JD = INT(365.25*(Y+4716)) + INT(30.6001*(M+1)) + D + B - 1524.5

The program seems to work correctly for positive dates (see tests) but whenever the date is a negative one the result is wrong; for example JD for -4712 January 1.5 return 37 instead of 0 or JD for -1000 July 12.5 return 1356010 instead of 1356001

I made the calculation also on paper and the program confirms what I get on paper, what is wrong with the formula?

[EDIT] I reworked my gist using the source code of libnova, the code now works all the time, but I am still confused on why

new gist -> https://gist.github.com/MattBlack85/c2030dc56c97f5f7fd06dd27a97413bb rust playground -> https://play.rust-lang.org/?version=stable&mode=debug&edition=2021&gist=98e513a4fa196d9e6eccbf6b83509cc2

Answer 1

This answer offers in two parts, (1) a first thought that might possibly help you to get debugged (but no promises), and (2) some ways to set up reliable date-conversion without having to 're-invent the wheel'.

(1) Your code adopts Meeus' floating-point-based algorithm, but I see (at least on your linked web-page) that you have included the constant 30.6 . Meeus explains on page 61 that this value, though it would be technically correct in a 'perfect' arithmetic, sometimes gives incorrect rounding/truncation to integer in real-life computer floating-point arithmetic. This can cause 'one-day-off' errors. So he instructs to use 30.601 or 30.6001 instead, to be always on the correct side of each integer boundary for rounding/truncation purposes.

So a first suggestion is that you try the recommended 30.601 or 30.6001 instead of 30.6, and see if that solves any problem. If not, the debugging will take a little longer.

(2) The programming problem for calendar-to-Julian-day conversions and the reverse was effectively addressed in the 1960s at the JPL (Jet Propulsion Laboratory) -- (and also elsewhere, as in PM2Ring's reference). I suppose that a reason why this wheel has been so often reinvented is the never-ending stream of new programming languages and interface requirements.

The JPL codes to convert between Julian Day Numbers and calendar dates (in both the Julian and Gregorian calendars) rely solely on integer arithmetic. The JPL codes assume separation of date+time combinations into two parts, one a date at Greenwich noon, to be converted with integer arithmetic only, and the other part is a floating-point difference or day-fraction from mean noon, in whatever units may be chosen suitably for the application in hand.

Much of the 'point' of using exclusively integer arithmetic for calendar conversions is that this avoids operations of rounding or truncation, which can give rise to obscure problems and 'one-day-off' errors.

The reliable JPL codes were copied in the 1992 'Explanatory Supplement to the Astronomical Almanac', and as they are quite short and cover both Julian and Gregorian calendars, I include them here below (in the form of page-extract-images, to avoid that I commit any typo errors here).

An alternative worth considering (for Gregorian-form calendar dates only) is to use the calendar subroutines included in the IAU 'Standards of Fundamental Astronomy' program-suite available in C or Fortran at (http://iausofa.org/). (Look for SOFA functions Cal2jd and Jd2cal).

If you continue to suffer from any obscure problem, it's possible that the shortest way in the end could be to install* (or simply port) the JPL or other integer routines into your environment, and test them out on suitable test-data (se below). (* Ensure that you use an integer type big enough to contain the numbers that will be processed, e.g. 'long integer'.)

Otherwise, if you still need to debug your code and can't find a shorter way than the following, I suggest you go through the calculation for two dates, one on the 'right' side and the other on the 'wrong' side of the boundary where your problem date-range begins. For each date work through both the algorithm that you installed, and the corresponding JPL algorithm, and isolate the point of divergence to localize the error in your code.

First, some sample test data from 'Explanatory Supplement to the Astronomical Ephemeris' (1961). Years before AD 1 are 'astronomical', counted with inclusion of a year '0'. For each year, the given JD is for Greenwich noon on 'January 0.5' (i.e. noon on Dec 31 of the year preceding).

>   Year      JD at Jan 0.5
>             for Julian Jan 0.5:
>   -2000       990557
>   -1900      1027082
> [...]
>    -200      1648007
>    -100      1684532
>       0      1721057
>     100      1757582
> [...]
>              Julian Jan 0.5   Gregorian Jan 0.5
>    1500      2268932          2268923
>    1600      2305457          2305447
>    1900      2415032          2415020

The following algorithms convert from and to Julian calendar dates:--

The following algorithms convert from and to Gregorian calendar dates:--

Answer 2

Sorry, I don't know Rust, so I won't try to debug your code, but it looks like there's some issue with the Julian calendar vs the Gregorian calendar. However, there's an off-by-one error too.

I see that you're using the arithmetic year numbering convention (sometimes known as the astronomical convention), eg 1 AD = year 1, 1 BC = year 0, 2 BC = year -1, etc.

Now -1000 Jul 12.5 Gregorian is JD 1356010, -1000 Jul 12.5 Julian is JD 1356001. But -4712 Jan 1.5 Gregorian is JD 38, -4712 Jan 1.5 Julian is JD 0. Perhaps you're getting 37 instead of 38 due to floating-point issues. I don't know how Rust converts negative floats to ints; some languages round towards zero, some round towards negative infinity.

Personally, I prefer to avoid using floats when the calculations can be done using integer arithmetic.

I see from your profile that you're familiar with Python. I have Python code that does Julian day number conversion for the Gregorian and Julian calendars on Github. It was derived from RG Tantzen (1963), ACM. Algorithm 199: conversions between calendar date and Julian day number. doi: 10.1145/366707.390020.

This code only uses integer arithmetic, and it should work correctly for any date. It has been exhaustively tested against all the dates that the standard Python datetime module accepts.

Here's a live version running on the SageMathCell server. When inputting a calendar date, you must give the month in numeric form, but you can use anything non-numeric to separate the year, month & day fields.

---

JPL Horizons provides a calendar date / Julian day number conversion tool, but it doesn't let you choose a calendar, it just silently switches to the Julian calendar for dates prior to 1582 Oct 15.