skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.39k stars 209 forks source link

Ancient hebrew calendar investigations #728

Closed hanetzer closed 2 years ago

hanetzer commented 2 years ago

Hello. As you may be aware, prior to the construction of the modern mathematical metonic hebrew calendar, months were began via the sighting of a new moon crescent, and the year based around the spring equinox and the celebration of passover.

To this end, my intent was to create an almanac of sorts. the idea was to have seven 'days' held in variables, like so:

ts = load.timescale()
t1 = ts.utc(a sunday date)
t2 = ts.utc(a monday date)

and so on, incrementing the variables t1-t7 by seven days each time, writing the data to a csv file for additional manipulation in say libreoffice or the like. Can one 'add' to a skyfield.timelib.Time object in this manner?

brandon-rhodes commented 2 years ago

Yes! Take a look at the documentation where it talks about "Date arithmetic":

https://rhodesmill.org/skyfield/time.html#date-arithmetic

Does that help you write your Python script?

hanetzer commented 2 years ago

This looks workable, thank you. Another question: as I'm starting each 'week' on sunday, and not every year starts on a sunday, is there a mechanism for getting, for example, the sunday before the start of year AD1, unless AD1 itself starts on a sunday, in which case it would return Jan 1, AD1 (all dates gregorian, so no need for julian trickery)

brandon-rhodes commented 2 years ago

If you ask Skyfield which day of the week it is, you can subtract that. So Sunday = 0 has no effect; Monday = 1 backs the date up one day to Sunday; Tuesday = 2 backs the date up 2 days to Sunday; and so forth. Hopefully that will work!

hanetzer commented 2 years ago

ahhhh, I getcha. so if 1/1/1 is a sunday, you subtract 0, and get 0, so still a sunday, but monday would subtract 1 and get the sunday prior. Nice. How does one ask skyfield the day of the week?

brandon-rhodes commented 2 years ago

Well — drat. None of the time tuples listed in the API docs for a Time return a day-of-the-week, do they?

https://rhodesmill.org/skyfield/api.html#time-objects

And you can't use a Python datetime object, because they don't work for ancient dates. Hmm. Can you ask for its plain .ut1 float? Something like int(t.ut1) should give you a day, and if you knew which day of the week day 0 was, you could produce a day of the week from that.

hanetzer commented 2 years ago

well, for instance, ts.utc(1,1,1) as an int gives 1721425. I could manually account for centuries and such, but I'm hoping to avoid as much manual adjustment as possible. And no, I don't see much of a way to get a weekday out cleanly, unless I parse something like utc_strftime() manually.

hanetzer commented 2 years ago

However, I can do t0.utc_strftime(format="%a") which returns Mon for 1/1/1. Suppose I could use some indexed object like Sun:0, Mon:1, and so on and subtract based on that.

hanetzer commented 2 years ago

Ran into a possible snafu. There is no year 0 (so the day after Dec 31, 1BC is Jan 1, 1AD). In principle the math works, however.

t0 = ts.utc(1,1,1)
weekday = { 'Sun': 0, 'Mon': 1, 'Tue': 2, 'Wed': 3, 'Thu': 4, 'Fri': 5, 'Sat': 6 }
t1 = t0 - weekday[t0.utc_strftime(format="%a")]
print(t1.utc_strftime())
#=> 0-12-31 00:00:00 UTC
brandon-rhodes commented 2 years ago

I could manually account for centuries and such, but I'm hoping to avoid as much manual adjustment as possible.

Happily, the "Julian date" time scale used by .ut1 and .tt ignores centuries, and so requires no manual adjustment at all. It's simply a count of days forward into history and back into the past. My idea was that you could apply an offset to account for the fact that maybe day 0.0 isn't the first day of the week; in a quick experiment I just did, the offset 5.0 works. So maybe:

day_of_week = (t.ut1 // 1.0 - 5.0) % 7.0
print(dow)

Try that out with some different dates t and see what you get. (If you don't want UT1, then change to .tt or whatever you want.)

hanetzer commented 2 years ago

Oh, I just realized, astronomical year numbering does include a year 0. Interesting.

brandon-rhodes commented 2 years ago

Oh, I just realized, astronomical year numbering does include a year 0. Interesting.

True. But, note that the year does not affect the progress of the week — the week repeats in a seven-day cycle that ignores months and years.

hanetzer commented 2 years ago

Yep, aware of that. Even the Julian->Gregorian shift didn't impact it.

hanetzer commented 2 years ago

Is there a defined mechanism for not using a zero year? Makes converting between astronomical years and 'historical years' a bit annoying, if nothing else.

brandon-rhodes commented 2 years ago

No, without a year zero, math would not work between years. Maybe you could write a small function of your own that takes a Skyfield Time and gives you the year in the format you want?

hanetzer commented 2 years ago

Ended up doing something along those lines. Another question: I know I can use pytz and such to set the timezone to israel, but that makes my 'days' begin and end two hours early/late, forget which. Is there some mechanism I'm missing that can give me days that start at 00:00 and end at 24:00 in a given time zone?

brandon-rhodes commented 2 years ago

I would need to see your example code to figure out what went wrong — I always just cut and paste the pytz examples from Skyfield's Time documentation and they have worked for me.

hanetzer commented 2 years ago

If you'll ignore the bad python (not a pythonist), this is what I got. Currently it prints the time of the events in UTC, but I'd like them in Jerusalem time.

#! /usr/bin/env python3
# -*- coding: utf-8 -*-
# vim:fenc=utf-8:spell:spelllang=en
#
# Copyright © 2022 hanetzer <hanetzer@proprietary-killer>
#
# Distributed under terms of the MIT license.

"""
Generates a CSV calendar for the given years. Includes Y/M/D by the Proleptic Gregorian Calendar,
the date and time of the Solstices and Equinoxes, and
"""
from datetime import datetime, timedelta
from pytz import timezone
from skyfield import almanac, eclipselib
from skyfield.api import load, wgs84, Topos, Star
import numpy as np

ts = load.timescale()
israel_tz = timezone("Israel")
a_day = timedelta(hours=24)

eph = load('de431t.bsp')
# print(eph)
sun, moon, earth = eph['sun'], eph['moon'], eph['earth']

# jerusalem = earth + Topos(latitude_degrees=(31, 46, 19), longitude_degrees=(35, 13, 1))
jerusalem = wgs84.latlon(31.7683, 35.2137)

WEEKDAYS = {
    'Sun': 0,
    'Mon': 1,
    'Tue': 2,
    'Wed': 3,
    'Thu': 4,
    'Fri': 5,
    'Sat': 6
}

RISESET = {
    0: "Sunset",
    1: "Sunrise"
}

def gen_first_week(time):
    """
    Finds the first week of a year, backdating to the previous Sunday.
    """
    t1 = time - WEEKDAYS[time.utc_strftime(format="%a")]
    t2 = t1 + 1
    t3 = t2 + 1
    t4 = t3 + 1
    t5 = t4 + 1
    t6 = t5 + 1
    t7 = t6 + 1

    return [t1, t2, t3, t4, t5, t6, t7]

def check_year_heads(time):
    """
    Finds the Equinoxes or Solstices that may occur in a given week.
    """
    t, y = almanac.find_discrete(time[WEEKDAYS['Sun']], time[WEEKDAYS['Sat']] + a_day, almanac.seasons(eph))
    if np.size(y):
        for yi, ti in zip(y, t):
            return f"{almanac.SEASON_EVENTS[yi]} {ti.utc_strftime('%m-%d %H:%M')}"
    else:
        return ""

def check_moon(time):
    """
    Finds Conjunctions or Full Moons which may occur in a given week.
    """
    returnval = None
    for t in time:
        ti, yi = almanac.find_discrete(t, t + a_day, almanac.moon_phases(eph))
        for tip, yip in zip(ti, yi):
            if yip == 0:
                returnval = f"Conjunction {tip.utc_strftime('%m-%d %H:%M')}"
            elif yip == 2:
                returnval = f"{almanac.MOON_PHASES[yip]} {tip.utc_strftime('%m-%d %H:%M')}"

    if returnval is None:
        returnval = ""

    return returnval

def gen_lunar_eclipse(time):
    t, y, details = eclipselib.lunar_eclipses(time[WEEKDAYS['Sun']], time[WEEKDAYS['Sat']] + a_day, eph)
    if np.size(y):
        for yi, ti, de in zip(y, t, details):
            return f"{eclipselib.LUNAR_ECLIPSES[yi]} {ti.utc_strftime('%m-%d %H:%M')}"
    else:
        return ""

def gen_sunrise_sunset(time):
    returnval = []
    for t in time:
        ti, yi = almanac.find_discrete(t, t + a_day, almanac.sunrise_sunset(eph, jerusalem))
        returnval.append(f'{ti[0].utc_strftime("%H:%M")}/{ti[1].utc_strftime("%H:%M")}')

    return ",".join(returnval)

def gen_week_string(time):
    year = int(time[WEEKDAYS['Sun']].utc_strftime('%Y')) # 0
    if year < 1:
        yearstr = "{}BC".format((year * -1) + 1)
    else:
        yearstr = "{}AD".format(year)

    return "{0},{1},{2},{3},{4},{5},{6},{7},{8},{9}".format(
            yearstr, # 0
            time[WEEKDAYS['Sun']].utc_strftime('%Y'), # 1
            time[WEEKDAYS['Sun']].utc_strftime('%m'), # 2
            time[WEEKDAYS['Sun']].utc_strftime('%d'), # 3
            time[WEEKDAYS['Mon']].utc_strftime('%d'), # 4
            time[WEEKDAYS['Tue']].utc_strftime('%d'), # 5
            time[WEEKDAYS['Wed']].utc_strftime('%d'), # 6
            time[WEEKDAYS['Thu']].utc_strftime('%d'), # 7
            time[WEEKDAYS['Fri']].utc_strftime('%d'), # 8
            time[WEEKDAYS['Sat']].utc_strftime('%d') # 9
    )

def increase_week(time):
    """
    """
    newdates = []
    for day in time:
        newdates.append(day + 7)

    return newdates

print('Common Year,Astronomical Year,Month,Sun,Mon,Tue,Wed,Thu,Fri,Sat,"Sun Rise/Set","Mon Rise/Set","Tue Rise/Set","Wed Rise/Set","Thu Rise/Set","Fri Rise/Set","Sat Rise/Set","Equinox and Solstice","Conjunction and Full Moon","Lunar Eclipses"')

time = ts.utc(-3,1,1,0,0)
dates = gen_first_week(time)
for i in range(13): # approximately 5218 weeks in 100 years
    # print(dates[WEEKDAYS['Sun']].utc_strftime())
    weekstring = gen_week_string(dates)
    riseset = gen_sunrise_sunset(dates)
    heads = check_year_heads(dates)
    moons = check_moon(dates)
    lunar_eclipses = gen_lunar_eclipse(dates)
    print(f"{weekstring},{riseset},{heads},{moons},{lunar_eclipses}")

    dates = increase_week(dates)
brandon-rhodes commented 2 years ago

My guess is that instead of:

time = ts.utc(-3,1,1,0,0)

— you will want to ask for midnight in the time zone israel_tz, using the first big block of example code in the section:

https://rhodesmill.org/skyfield/time.html#utc-and-your-timezone

Try it out, and let me know if you get a usable result.

hanetzer commented 2 years ago

Well, starting at Jan 1, 1AD, with the following code snippet:

d = datetime(1,1,1,0,0)
jt = israel_tz.localize(d)
time = ts.from_datetime(jt)

I get the following error on execution:

Traceback (most recent call last):
  File "calends/gencal.py", line 155, in <module>
    jt = israel_tz.localize(d)
  File "calends/venv/lib/python3.9/site-packages/pytz/tzinfo.py", line 323, in localize
    loc_dt = dt + delta
OverflowError: date value out of range
brandon-rhodes commented 2 years ago

Oh. I suppose that makes sense: Python datetime cannot go back to years before the year 1, and trying to go back 2 hours from the date you gave it exceeds that limit.

I wish that Python datetimes weren't so limited.

Well, drat. I guess you might have to subtract 2.0 / 24.0 from your first datetime, then, if you want two hours earlier, since it looks like Python time zone support won't help. Which I guess they could defend on the grounds that the time zone was not really in force that many years ago!

If you are doing serious work, you might want to search for the true times of sunrise and sunset in Jerusalem on those ancient dates, and use those times to begin and end each day, since the ancient authorities would have used their view of the sun, and not a continuous perfect clock, as their reference.

hanetzer commented 2 years ago

Unfortunately I don't think I'm going to find a hebrew almanac that old, lmao. I'm just using the skyfield calculated times as a 'general reference', which is good enough for my purposes.

hanetzer commented 2 years ago

that being said, I'd like to be able to specify the time on the command line instead of manually editing the file. However, I'm not seeing any mechanism to pass a string-like object to ts.something() to construct it. Tips?

brandon-rhodes commented 2 years ago

Unfortunately I don't think I'm going to find a hebrew almanac that old, lmao. I'm just using the skyfield calculated times as a 'general reference', which is good enough for my purposes.

Oh — I meant, you could ask Skyfield for the sunrise and sunset times — I wasn't suggesting you find an external almanac source.

However, I'm not seeing any mechanism to pass a string-like object…

I would suggest a quick re-read of the Python Tutorial on the language's website, it includes all the basic operations you need to split a string into pieces then turn those pieces into numbers.

hanetzer commented 2 years ago

Yeah. figured it out. small amount of argparse and I just pass it a year (since I'm doing 100y chunks from Jan 1 of each chunk) and str() it.

brandon-rhodes commented 2 years ago

Great, I'm glad you worked out a solution! Given that the original issue now sounds like it's resolved, I'm going to go ahead and close the issue, but please feel free to add further comments to the closed issue.

hanetzer commented 2 years ago

Spoke too soon. adjusting for israel time manually ended up giving me two dec 30's in a row one year and not adjusting gives me borked sunrise/sunset dates. If datetime worked with BC dates I think this would be no issue.

brandon-rhodes commented 2 years ago

If datetime worked with BC dates I think this would be no issue.

Without an example, I can't speak with much confidence about the results you're seeing, but I don't understand your idea that duplicate results, or sunrise and sunset times you don't like, could result from as minor an issue as whether a year is called "-3" or "4 BCE"? The name of the year would not normally have the power to change sunrise and sunset results.

hanetzer commented 2 years ago

To account for Israel time I was started from 22:00 UTC(which is midnight israel time) and adjusting on the fly while printing. Have you considered rolling a tz aware variation of Time? It just strikes me as odd I can calc sunrise/sunset based on location but can't easily produce a time for it relative to that location. And I do -3, not 4BCE, but I do provide conversions .

brandon-rhodes commented 2 years ago

Have you considered rolling a tz aware variation of Time?

Alas, no.

It's possible that we could explore a routine that would take a Time and determine what local hour-of-the-day an ancient person would have computed based on sunrise, transit, and sunset. Would that help with the historical question you are trying to answer?

hanetzer commented 2 years ago

Well yes, I'm intimately familiar with the idea of 'at the xth hour'. Its just that at certain points in time, using utc 0 thru utc 24 as a 'day' will fail to find local sunrise/sunset for a given location. Plus, even if python timezones don't work that far back, the fact of the matter is israel's sunrise/set is offset by two hours compared to UTC/GMT.

It makes sense to have local times for local phenomenons.

brandon-rhodes commented 2 years ago

Spoke too soon. adjusting for israel time manually ended up giving me two dec 30's in a row one year and not adjusting gives me borked sunrise/sunset dates.

Backing up, to remind myself of the actual problem: could you share this result you got, with two dec 30's? That should not happen because of a manual adjustment; and if it did, then support for time zones or support for BC wouldn't fix it; the underlying math that Skyfield did would be the same. So I think we should focus on the duplicate dec 30, which must be a bug from somewhere else; it can't be related to subtracting two hours from your start and stop times (unless, of course, something very subtle is taking place that I haven't yet grasped).