cosinekitty / astronomy

Astronomy Engine: multi-language calculation of Sun, Moon, and planet positions. Predicts lunar phases, eclipses, transits, oppositions, conjunctions, equinoxes, solstices, rise/set times, and other events. Provides vector and angular coordinate transforms among equatorial, ecliptic, horizontal, and galactic orientations.
MIT License
496 stars 64 forks source link

Date calculations #325

Open sbooth opened 1 year ago

sbooth commented 1 year ago

The documentation for astro_time_t says that ut is the UT1/UTC number of days since noon on January 1, 2000. My understanding is that means that the value ut is more or less the Julian date with 2451545 subtracted to form a year 2000 variant. This seems to be the case for dates after the Gregorian calendar changeover:

Creation JD JD - 2451545 ut
Astronomy_MakeTime(2000, 1, 1, 0, 0, 0) 2451544.5 -0.5 -0.5
Astronomy_MakeTime(2000, 1, 31, 0, 0, 0) 2451574.5 29.5 29.5
Astronomy_MakeTime(1582, 10, 15, 0, 0, 0) 2299160.5 -152384.5 -152384.5

The value diverges for earlier dates (at least using the NASA JD calculator, which automatically uses the Julian calendar dates before 1582-10-15):

Creation JD JD - 2451545 ut Δ
Astronomy_MakeTime(1582, 10, 4, 0, 0, 0) 2299159.5 -152385.5 -152395.5 10
Astronomy_MakeTime(1, 1, 1, 0, 0, 0) 1721423.5 -730121.5 -730119.5 -2

I think all dates passed to Astronomy_MakeTime are being interpreted as Gregorian dates?

Using the proleptic Gregorian calendar, 0001-01-01 has a JD of 1721425.5, which with the Y2K adjustment applied matches the value of ut (1721425.5 - 2451545 = -730119.5).

This may be the meaning of the documentation for the Astronomy_MakeTime function which states that it uses the UTC calendar.

I think it would be helpful if the documentation stated that the dates are interpreted as Gregorian dates.

If you're interested in handling both Julian and Gregorian dates (which I think would simplify things for historical calculations) I have written code that could be easily ported to Astronomy Engine.

cosinekitty commented 1 year ago

Yes, you are absolutely right. The intention was to have a continuous time scale without any discontinuities, so I chose to use the modern Gregorian calendar extrapolated backwards indefinitely. I can see why people might want to be able to convert old Julian calendar dates into astro_time_t. So if you want to add a new function that can convert Julian dates when appropriate, that would be great as a separate function. I don't want to change the meaning of Astronomy_MakeTime because there is code that uses it over spans of thousands of years and expects it to be a continuous function.

There are other functions that also assume Gregorian calendars over all time:

I also like to keep all programming languages in sync (C/C++, C#, JavaScript, Kotlin, Python, and soon Go). It would be nice to add support to the other languages as well. Currently all languages use the Gregorian calendar consistently throughout.

So you can see we have at least 4 functions, multiplied by 6 programming languages, which could end up being a lot of work. I'm not sure how eager you are to do all this. Things in Astronomy Engine often end up being far more hassle than new contributors realize!

sbooth commented 1 year ago

I can definitely write the functions in C/C++, Javascript, and Python. I haven't used the other languages before but it's mostly integer math for the algorithms so maybe it would not be terribly difficult. Or maybe it would be, who knows!

I'll code something in C as an experiment (Astronomy_MakeTime_Julian?) and post it here for you to decide if it's something worth pursuing.

sbooth commented 1 year ago

Here's a C version that seems to work but is not extensively tested:

const int64_t y = 4716;
const int64_t j = 1401;
const int64_t m = 2;
const int64_t n = 12;
const int64_t r = 4;
const int64_t p = 1461;
const int64_t q = 0;
const int64_t v = 3;
const int64_t u = 5;
const int64_t s = 153;
const int64_t t = 2;
const int64_t w = 2;

astro_time_t Astronomy_MakeTime_Julian(int year, int month, int day, int hour, int minute, double second)
{
    astro_time_t time;
    int64_t y = (int64_t)year;
    int64_t m = (int64_t)month;
    int64_t d = (int64_t)day;
    int64_t cycles = 0;

    if (y < -4712) {
        cycles = (-4713 - y) / 4 + 1;
        y += cycles * 4;
    }

    const int64_t h = m - 2;
    const int64_t g = y + 4716 - (12 - h) / 12;
    const int64_t f = (h - 1 + 12) % 12;
    const int64_t e = (1461 * g + 0) / 4 + d - 1 - 1401;
    int64_t J = e + (153 * f + 2) / 5;

    if (cycles > 0)
        J -= cycles * 1461;

    int64_t y2000 = J - 2451545;

    time.ut = (y2000 - 0.5) + (hour / 24.0) + (minute / 1440.0) + (second / 86400.0);
    time.tt = TerrestrialTime(time.ut);
    time.psi = time.eps = time.st = NAN;

    return time;
}

It is based on E.G. Richards' algorithm with modifications by me to support years before -4712.

cosinekitty commented 1 year ago

That's great! We just need a few test dates/times to know whether the numbers are correct. Then I can help write a unit test to validate the C version. I can help implement some of the versions in the other languages too.

I just created a new branch julian that you can feel free to submit a PR into. Then we can work offline without disrupting the master branch.

Thanks again for taking an interest in this!