DeltaT and sunrise/sunset times

[sgbmzm]

Skyfield includes DeltaT in the time calculation. but there is a problem. DeltaT is very important when trying to calculate ancient astronomical events such as Halley's Comet, etc., which are trying to calculate them backwards according to the known rate at which they occur. And also regarding solar eclipses and the like. But when trying to calculate daily times, such as sunrise, sunset, and meridian transit, it creates a problem, and I will give an example. Let's take Jerusalem as an example. Every year on the date of Vernal Equinox the sun sets around 18:00 according to standard local time which is: UTC+2.

But when I try to do the same calculation for 23 AD, that is 2000 years back, the sun sets at about 21:00 at the Vernal Equinox in Jerusalem, according to UTC+2.

The reason why the sun supposedly sets at 21:00 is because deltaT comes into the calculation. And in 23 AD deltaT was about 10,417 seconds which is 173 minutes which is close to three hours.

But no person in Jerusalem thought in 23 that the sun sets at 9:00 p.m. It is clear that they thought that the sun sets at 18:00 and of course the middle of the day was at 12:00 and not at 15:00

[brandon-rhodes]

Could you include sample code, and show which time scale you are using? If you are using a timescale that includes ∆T, you should find a much better match between solar events and the time.

[sgbmzm]

Here is a simple example of the sunset.

from pytz import timezone

israel = timezone('Israel')

eph = load('de441s.bsp')
bluffton = wgs84.latlon(31.940826 * N, 35.037057 * E)
t0 = ts.utc(23, 3, 21, 4)
t1 = ts.utc(23, 3, 22, 4)
t, y = almanac.find_discrete(t0, t1, almanac.sunrise_sunset(eph, bluffton))

print(t.utc_iso())
print(y)
print(t.astimezone(israel))

['0023-03-21T06:33:07Z', '0023-03-21T18:41:05Z'] [1 0] [datetime.datetime(23, 3, 21, 8, 54, 7, 70091, tzinfo=<DstTzInfo 'Israel' LMT+2:21:00 STD>) datetime.datetime(23, 3, 21, 21, 2, 4, 900269, tzinfo=<DstTzInfo 'Israel' LMT+2:21:00 STD>)]

from pytz import timezone
israel = timezone('Israel')
eph = load('de441s.bsp')
bluffton = wgs84.latlon(31.940826 * N, 35.037057 * E)
t0 = ts.utc(2023, 3, 21, 4)
t1 = ts.utc(2023, 3, 22, 4)
t, y = almanac.find_discrete(t0, t1, almanac.sunrise_sunset(eph, bluffton))

print(t.utc_iso())
print(y)
print(t.astimezone(israel))

['2023-03-21T15:51:42Z', '2023-03-22T03:41:43Z'] [0 1] [datetime.datetime(2023, 3, 21, 17, 51, 42, 320825, tzinfo=<DstTzInfo 'Israel' IST+2:00:00 STD>) datetime.datetime(2023, 3, 22, 5, 41, 43, 473581, tzinfo=<DstTzInfo 'Israel' IST+2:00:00 STD>)]

[brandon-rhodes]

Good! That output makes sense: because the UTC time scale ignores ∆T, it fails to follow the change in day length, both back into the past, and also forward into the future. Only over the narrow timespan of the past few decades, when leap seconds are active, does it follow sunrise and sunset and solar noon.

So what you really want to do is compute sunrise and sunset, and then use those two times to define what ancient peoples would have called ‘the first hour’ and ‘the twelfth hour’.

Failing that, you could get an approximate match using UT1, which does follow ∆T:

https://rhodesmill.org/skyfield/time.html#ut1-and-downloading-iers-data

But in my opinion it's best to ignore modern time scales and name ancient times using ancient sunrise and sunset, since that's what people at the time would have used.

Note that even in modern times when UTC follows ∆T, sunrise on the equinox and sunset aren't at 6:00 and 18:00 because of the analemma:

https://en.wikipedia.org/wiki/Analemma

I'm going to close this issue as I don't think there's a bug with Skyfield involved here, but please comment further on this issue if you have any more questions!

[sgbmzm]

Thanks! If UTC time scale ignores ∆T how can you use another time scale that does not ignore ∆T?

[brandon-rhodes]

You could get an approximate match using UT1, which does follow ∆T:

https://rhodesmill.org/skyfield/time.html#ut1-and-downloading-iers-data

[sgbmzm]

There are two matters:

1. Give Skyfield an input of date ut1. I understand that I do this using ts.ut1(20, 6, 1, 16, 5) instead of ts.utc(20, 6, 1, 16, 5) but I usually get dates in datetime format so how do I load UT1 using ts.from_datetime()?
2. Get from Skyfield an output of date ut1. I understand that this is done using t.ut1_strftime(), and it does give an output that doesn't ignore delta-t but the thing is that I usually need to get a datetime output. So how can you get an output of ut1 as a datetime?
3. And just like that, referring to the example I gave above, I get a result that is really close to UT1 even when I use the following option: t.astimezone(israel) - timedelta(seconds=int(t[0].delta_t))

t.ut1_strftime(): ['0023-03-21 03:43:55 UT1', '0023-03-21 15:51:53 UT1']

t.astimezone(israel) - timedelta(seconds=int(t[0].delta_t)): [datetime.datetime(23, 3, 21, 6, 4, 13, 70131, tzinfo=<DstTzInfo 'Israel' LMT+2:21:00 STD>) datetime.datetime(23, 3, 21, 18, 12, 10, 900249, tzinfo=<DstTzInfo 'Israel' LMT+2:21:00 STD>)]

Thanks again

[sgbmzm]

Another serious problem I encountered because of delta-t is the calculation of the equation of time. I defined it as the difference between midnight local average time and 12:00. So in contemporary years it works, but if UTC does not consider DELTA-T and makes it so that midnight is for example at 14:40 then of course the time equation comes out wrong. If I try to omit DELTA-T from the result I still get times that are slightly different than what I get here: https://gml.noaa.gov/grad/solcalc/ So how can I calculate the equation of time in Skyfield another way? The way I used it is:

Calculation of meridian transit:

time_6 = israel.localize(datetime(2023, 9, 30, 6))
meridian_transit = almanac.meridian_transits(eph, eph["sun"], location)
transit, y = almanac.find_discrete(ts.from_datetime(time_6), ts.from_datetime(time_6 + dt.timedelta(hours=24)), meridian_transit)
day_transit = transit[0].astimezone(location_timezone)

convert time to LMT

def localmeantime(utc_time,longitude):
    lmt = utc_time + timedelta(seconds=round(4*60*longitude))
    lmt = lmt.replace(tzinfo=None)
    return lmt

convert day_transit to LMT

transit_lmt = localmeantime(day_transit.astimezone(timezone('utc')),location.longitude.degrees)
transit_12 = transit_lmt.replace(hour=12, minute=0, second=0, microsecond=0)
if transit_lmt < transit_12:
        equation_of_time = transit_12-transit_lmt
else:
        equation_of_time = transit_lmt-transit_12

[brandon-rhodes]

I usually get dates in datetime format so how do I load UT1 using ts.from_datetime()?

Each Python datetime object has attributes offering the month, year, and so forth. I think you can simply pass those to the Skyfield ts.ut1() method?

So how can you get an output of ut1 as a datetime?

Take the output of ut1_calendar() and pass the components to the datetime() constructor.

So how can I calculate the equation of time in Skyfield another way?

My approach would be to compute afresh the timing of the “mean Sun” for the ancient year I am interested in. So, I would compute the times of all the solar transits for the year, then take an average across all of them, and call that the average solar noon for that year. Then, the equation of time would be computed compared to that.

[sgbmzm]

Hi. Thank you! For a long time I could not make progress on this matter, so I am only now responding.

1. Can you define exactly from which year to which year the UTC time scale does not ignore delta-T?
2. I did as you suggested regarding the calculation of "mean Sun" and to derive the "equation of time". The problem is that it is a lot of calculations, and takes a lot of time, much more than I can afford for my software. Is there a way to calculate this faster? I am attaching an example of what I did

from skyfield.api import N, S, E, W, wgs84, load, load_file
from skyfield.framelib import ecliptic_frame
from skyfield.earthlib import refraction
from skyfield import almanac
from skyfield.units import Angle
from skyfield.nutationlib import iau2000b_radians

from skyfield.searchlib import find_discrete

find_discrete

import datetime as dt
from datetime import datetime, timedelta, date
from pytz import timezone

from numpy import cos, zeros_like, arccos

eph = load('de441.bsp')

israel = timezone('Israel')
utc = timezone('utc')

location_timezone = israel

location = wgs84.latlon(31.940826 * N, 35.037057 * E)

ts = load.timescale()

target_time = datetime(3,6,15) 

time = israel.localize(target_time)

# Setting the time to 6 am (local)
time_6 = time.replace(hour=6, minute=0, second=0, microsecond=0)

# Calculation of transits in the search between six in the morning and six in the morning the next day
meridian_transit = almanac.meridian_transits(eph, eph["sun"], location)
transit, y = almanac.find_discrete(ts.from_datetime(time_6), ts.from_datetime(time_6 + dt.timedelta(hours=24)), meridian_transit)
day_transit = transit[0].astimezone(location_timezone)
night_transit = transit[1].astimezone(location_timezone)

# Setting the start and end date of the requested Gregorian year, in local time
datetime_greg_year_start = location_timezone.localize(datetime(target_time.year, 1, 1, 0, 0 ,0))
datetime_greg_year_end = location_timezone.localize(datetime(target_time.year, 12, 31, 23, 59 ,59))

# Calculation of midnight of the day and night for each day of the requested year
transits, y = almanac.find_discrete(ts.from_datetime(datetime_greg_year_start), ts.from_datetime(datetime_greg_year_end), meridian_transit)
# Collection every day-transits only according to the time in Greenwich time for each of the days of the year into a new array
transits = [transits[i].astimezone(timezone('utc')) for i in range(len(transits)) if y[i] == 1]

# Conversion of the day-transits hour of each day of the year from hours minutes and seconds - to seconds
total_seconds = sum(time.hour * 3600 + time.minute * 60 + time.second for time in transits)

# Calculate average seconds
# Average calculation of the number of seconds which is basically a calculation of average day-transits for the whole year
average_seconds = total_seconds // len(transits)

# Convert average seconds back to datetime object
mean_time = datetime(1900, 1, 1) + timedelta(seconds=average_seconds)

# Convert mean time back to string format
mean_time_str = mean_time.strftime('%H:%M:%S')

# Real day-transit calculation of the requested day in Greenwich Mean Time
transit_utc = day_transit.astimezone(timezone('utc'))
# Average day-transit calculation in Greenwich Mean Time
transit_mean = transit_utc.replace(hour=mean_time.hour, minute=mean_time.minute, second=mean_time.second, microsecond=mean_time.microsecond)

print("transit_utc",transit_utc)
print("transit_mean",transit_mean)

# Checking the time gap, i.e. the delta between the local day-transit time and the average day-transit time, which is the time equation, and checking whether the time equation adds minutes or subtracts minutes
if transit_utc < transit_mean:
    equation_of_time = transit_mean-transit_utc
    seconds_equation_of_time = equation_of_time.total_seconds()
    str_equation_of_time = f'+{str(equation_of_time)[2:7]}'
else:
    equation_of_time = transit_utc-transit_mean
    seconds_equation_of_time = - equation_of_time.total_seconds()
    str_equation_of_time = f'-{str(equation_of_time)[2:7]}'

print("str_equation_of_time",str_equation_of_time)

transit_utc 0003-06-15 12:28:14.912553+00:00 transit_mean 0003-06-15 12:32:53+00:00 str_equation_of_time +04:38

[brandon-rhodes]

Apparently the holidays and travel led to my losing track of your question — the tab I opened it in must have gotten closed accidentally! In case answers are still useful, more than a month later:

Can you define exactly from which year to which year the UTC time scale does not ignore delta-T?

Yes! Leap seconds have been maintained from 1972 to the present:

https://data.iana.org/time-zones/tzdb-2018a/leap-seconds.list

For all previous eras, UTC is simply a count of SI seconds backwards into history that ignore UT1.

The problem is that it is a lot of calculations, and takes a lot of time, much more than I can afford for my software.

Then let's look for something quicker. I just added new rising, setting, and transit routines to Skyfield in the most recent release, that are much faster, so there's a good chance that we can speed up your calculation.

But can we speed it up enough? We should start by asking: how much time can you afford for your software? How quickly does the computation need to finish?

Also, could you outline what problem you are trying to solve which these calculations — what question you are trying to answer? Re-reading this thread, I'm not sure what your goal is. It might be that you could accomplish it with fewer calculations, and thus have a faster script.

[sgbmzm]

Thanks for your response. Well, here's the story: my main problem started when I noticed that the times I was getting for early years (for example AD 10) didn't fit a time scale that seemed appropriate to me. I especially noticed this when I tried to derive the equation of time from the difference between 12:00 local time, and the day transit time of the local. I then understood from you that this is happening because I am using UTC timescale instead of using UT1 timescale. But it was awkward for me. So I tried to solve the problem by reducing DELTA-T but it didn't help enough. So, per your advice, I built the function I brought above to find the annual average local day-transit time. And the difference between it and real Transit will be the equation of time. The calculation took about 6 seconds per year, while I can give a maximum of 1 second to this calculation. But I have already moved on. Instead of reducing DELTA-T I used the following definition: ts = load.timescale(delta_t=0.0) So I took the code above, and changed a few things:

1. ts = load.timescale(delta_t=0.0)
2. I changed the geographic location to longitude 0: location = wgs84.latlon(0 * N, 0 * E) Correct me if I'm wrong, but I think the result of the day-transit time averaged over a whole year, must be 12:00. The result I got was different. So I don't know what the problem is, but I looped through all the years from year 1 to year 3000 which is the relevant year range for me (the loop ran for 3 hours), and checked how much the difference is between 12:00 UTC at longitude 0 and the average yearly result. The difference is between 50 and 72 seconds. So I saved the second difference for each year into a file and added it on DELTA-T 0.0. For me it is a corrected time scale, suitable for UT1.

So the topic title right now is, in my opinion, An annual mean day-transit at longitude 0 does not point to 12:00 even though delta-T is set to zero. What do you think?

I am attaching the corrected code that I used, and sorry that the explanations are in Hebrew

# סקייפילד
from skyfield.api import N, S, E, W, wgs84, load, load_file
from skyfield.framelib import ecliptic_frame
from skyfield.earthlib import refraction
from skyfield import almanac
from skyfield.units import Angle
from skyfield.nutationlib import iau2000b_radians

from skyfield.searchlib import find_discrete

# השורה הזו חשובה מאוד והיא מועתקת מתוך ספריית האלמנך
find_discrete

# תאריכים
import datetime as dt
from datetime import datetime, timedelta, date
from pytz import timezone

# חישובים
from numpy import cos, zeros_like, arccos

eph = load('de441.bsp')

utc = timezone('utc')

location_timezone = utc

location = wgs84.latlon(0 * N, 0 * E)

ts = load.timescale(delta_t=0.0)

# יצירת מילון לאחסון תוצאות
results_dict = {}

# לולאה על כל שנה מ-2 עד 2999
for year in range(1, 3000):
    # יצירת תאריך ל-17 ליוני בשנה הנוכחית
    target_time = datetime(year, 6, 17)

    time = utc.localize(target_time)

    # הגדרת הזמן לשעה שש בבוקר
    time_6 = time.replace(hour=6, minute=0, second=0, microsecond=0)

    # חישוב חצות היום והלילה בחיפוש בין שש בבוקר לשש בבוקר למחורת: נזכיר כי הזמן הוא כבר בשעון מקומי
    meridian_transit = almanac.meridian_transits(eph, eph["sun"], location)
    transit, y = almanac.find_discrete(ts.from_datetime(time_6), ts.from_datetime(time_6 + dt.timedelta(hours=24)), meridian_transit)
    day_transit = transit[0].astimezone(location_timezone)
    night_transit = transit[1].astimezone(location_timezone)

    # הגדרת תאריך ההתחלה והסיום של השנה הגרגוריאנית המבוקשת, בשעון מקומי
    datetime_greg_year_start = location_timezone.localize(datetime(target_time.year, 1, 1, 0, 0 ,0))
    datetime_greg_year_end = location_timezone.localize(datetime(target_time.year, 12, 31, 23, 59 ,59))

    # חישוב חצות היום והלילה עבור כל יום מימי השנה המבוקשת
    transits, y = almanac.find_discrete(ts.from_datetime(datetime_greg_year_start), ts.from_datetime(datetime_greg_year_end), meridian_transit)
    # איסוף כל חצות היום בלבד לפי השעה בשעון גריניץ עבור כל אחד מימי השנה לתוך מערך חדש
    transits = [transits[i].utc for i in range(len(transits)) if y[i] == 1]

    # המרה של שעת חצות היום של כל יום מימי השנה מ-שעות דקות ושניות - ל-שניות
    total_seconds = sum(time.hour * 3600 + time.minute * 60 + time.second for time in transits)

    # חישוב ממוצע של מספר השניות שזה בעצם חישוב של חצות ממוצע עבור כל השנה
    average_seconds = total_seconds // len(transits)

    # הפיכת הזמן שהתקבל לזמן פייתון (תחילת יממת שנת 1900 זה לא בדווקא)
    mean_time = datetime(1900, 1, 1) + timedelta(seconds=average_seconds)

    # Convert mean time back to string format
    mean_time_str = mean_time.strftime('%H:%M:%S')

    # חישוב חצות מקומי אמיתי של היום המבוקש בשעון גריניץ
    transit_utc = day_transit.astimezone(timezone('utc'))

    # שעון גריניץ של חצות הממוצע
    transit_mean = transit_utc.replace(hour=mean_time.hour, minute=mean_time.minute, second=mean_time.second, microsecond=mean_time.microsecond)

    delta_t = ts.from_datetime(day_transit).delta_t

    transit_12 = transit_utc.replace(hour=12, minute=0, second=0, microsecond=0)

    # בדיקת פער הזמן כלומר הדלתא בין שעה מקומית ממוצעת של חצות לבין השעה 12 בצהריים ובדיקה האם יש להוסיף דקות או להפחית דקות
    if transit_mean < transit_12:
        error_delta_t = transit_12-transit_mean
        seconds_error_delta_t = error_delta_t.total_seconds()
        str_error_delta_t = f'+{str(error_delta_t)[2:7]}'
    else:
        error_delta_t = transit_mean-transit_12
        seconds_error_delta_t = - error_delta_t.total_seconds()
        str_error_delta_t = f'-{str(error_delta_t)[2:7]}'

    repaired_delta_t = delta_t + seconds_error_delta_t # חיבור כולל בתוכו חיסור כשיש מינוס

    # עדכון המילון בערך המתוקן של דלטא_טי עבור השנה המתאימה
    results_dict[year] = repaired_delta_t

    #print("original_delta_t",delta_t)
    #print("seconds_eror_delta_t", seconds_error_delta_t)
    #print("repaired_delta_t",repaired_delta_t)
    #print("transit_utc",transit_utc)
    #print("transit_mean",transit_mean)

# מכאן והלאה שמירת הנתונים לקובץ סי.אס.וי
# נתיב לקובץ CSV
csv_file_path = 'ts_seconds.csv'

# כותב CSV
with open(csv_file_path, 'w', encoding=cu_encod, newline='') as csvfile:
    # הגדרת עמודות ה-CSV
    fieldnames = ['Year', 'Seconds']

    # יצירת אובייקט מכתב עבור CSV
    writer = csv.DictWriter(csvfile, fieldnames=fieldnames)

    # כתיבת כותרות עמודות
    writer.writeheader()

    # כתיבת נתונים לקובץ CSV
    for year, seconds in results_dict.items():
        writer.writerow({'Year': year, 'Seconds': seconds})

[brandon-rhodes]

In case any other contributors want to try out your script — I needed to remove the encoding= parameter because I got an error:

Traceback (most recent call last):
  File "tmp904.py", line 204, in <module>
    with open(csv_file_path, 'w', encoding=cu_encod, newline='') as csvfile:
NameError: name 'cu_encod' is not defined

And then there was another error, fixed with import csv:

Traceback (most recent call last):
  File "tmp904.py", line 209, in <module>
    writer = csv.DictWriter(csvfile, fieldnames=fieldnames)
NameError: name 'csv' is not defined
python tmp904.py  14.33s user 9.74s system 317% cpu 7.587 total

The script then ran successfully! Running it over three years, I got this output — is this similar to what you're seeing?

Year,Seconds
1999,65.0
2000,65.0
2001,66.0

By the way, you can switch to using Skyfield's new transit function, if you want your script to run a bit faster. Here is what the central part of your loop would look like:

    # חישוב חצות היום והלילה בחיפוש בין שש בבוקר לשש בבוקר למחורת: נזכיר כי הזמן הוא כבר בשעון מקומי
    observer = eph['earth'] + location
    start = ts.from_datetime(time_6)
    end = ts.from_datetime(time_6 + dt.timedelta(hours=24))
    transits = almanac.find_transits(observer, eph['sun'], start, end)
    day_transit = transits[0].astimezone(location_timezone)

    # הגדרת תאריך ההתחלה והסיום של השנה הגרגוריאנית המבוקשת, בשעון מקומי
    datetime_greg_year_start = location_timezone.localize(datetime(target_time.year, 1, 1, 0, 0 ,0))
    datetime_greg_year_end = location_timezone.localize(datetime(target_time.year, 12, 31, 23, 59 ,59))

    # חישוב חצות היום והלילה עבור כל יום מימי השנה המבוקשת
    start = ts.from_datetime(datetime_greg_year_start)
    end = ts.from_datetime(datetime_greg_year_end)
    transits = almanac.find_transits(observer, eph['sun'], start, end)

    # איסוף כל חצות היום בלבד לפי השעה בשעון גריניץ עבור כל אחד מימי השנה לתוך מערך חדש
    transits = [t.utc for t in transits]

So, to address your question:

An annual mean day-transit at longitude 0 does not point to 12:00 even though delta-T is set to zero.

With ∆T set to zero, we would expect solar noon to occur somewhere around noon Terrestrial Time. And since the difference between TT and UTC in 1999 was around 64 or 65 seconds (~32.1s for the difference between TT and TAI, and then 32 leap seconds between TAI and UTC), the value 65.0 prints out.

Does that explain the result you are seeing from the script?

[sgbmzm]

Yes. Thanks! I think that explains it. The difference between TT and TAI is what caused the problem for me. Now, if we go to the year 500 the difference between TT and UTC will only be the difference between TT and TAI because in the year 500 there were no leap seconds. If so, without whatever I've done, is there another way in Skyfield to get the number of seconds between TT and UTC to add to delta t zeros? And basically I want to get three data:

1. Difference between TT and TAI
2. Difference between TAI and UTC
3. All together: the difference between TT and UTC

[brandon-rhodes]

I'll look and see how you can most easily compute those offsets!

But — would things be simpler if you just used TT, and stop using UTC? There's a method tt_calendar() on times that would give you a tuple whose year/month/day/hour/minute/second you could use.

[sgbmzm]

You wrote:

"With ∆T set to zero, we would expect solar noon to occur somewhere around noon Terrestrial Time. And since the difference between TT and UTC in 1999 was around 64 or 65 seconds (~32.1s for the difference between TT and TAI, and then 32 leap seconds between TAI and UTC), the value 65.0 prints out."

Try to do this calculation for the year 500. Is the data explained even then? In my csv file I got a difference of 47 seconds for year 500

day_transit by utc_datetime() comes out: datetime.datetime(500, 6, 17, 11, 56, 21, 893619, tzinfo=datetime.timezone.utc)

ut1_calendar() and tt_calendar() give the same result: (500, 6, 17, 11, 57, 4.07761886715889)

If you take utc_datetime() and add 47 seconds to it you will get: 11:57:08 (not 11:57:04)

[brandon-rhodes]

The very first leap second was in 1972, and moved the TAI - UTC offset from 10 to 11 seconds. Here's how you can print the leap second table:

from skyfield.api import load
ts = load.timescale(builtin=True)
for jd, offset in zip(ts.leap_dates, ts.leap_offsets):
    ymd = ts.tt_jd(jd).tt_strftime('%Y-%m-%d')
    print(ymd, '{:+}'.format(int(offset)))

The result:

1972-07-01 +11
1973-01-01 +12
1974-01-01 +13
1975-01-01 +14
1976-01-01 +15
1977-01-01 +16
1978-01-01 +17
1979-01-01 +18
1980-01-01 +19
1981-07-01 +20
1982-07-01 +21
1983-07-01 +22
1985-07-01 +23
1988-01-01 +24
1990-01-01 +25
1991-01-01 +26
1992-07-01 +27
1993-07-01 +28
1994-07-01 +29
1996-01-01 +30
1997-07-01 +31
1999-01-01 +32
2006-01-01 +33
2009-01-01 +34
2012-07-01 +35
2015-07-01 +36
2017-01-01 +37

Does that explain the offset you're seeing?

[sgbmzm]

Well, even if we use tt_calendar() the average transit at longitude 0 does not come out at 12:00 if we go to a year like year 500 (difference of about 5 seconds) or year 2 (difference of 8 seconds) then your first explanation might have been correct For 1999 but for very early years I still have a second gap that I don't know how to explain. And as I said: the right time scale for me is a time scale where the average transit will be at 12:00.

For this purpose, I changed the following lines from the code I provided above as follows:

transits = [transits[i].tt_calendar() for i in range(len(transits)) if y[i] == 1]
total_seconds = sum(time[3] * 3600 + time[4] * 60 + time[5] for time in transits)

[brandon-rhodes]

Let me know when you have the chance to try out the new find_transits() routine. I think it's important, before you spend more time with your attempt to come up with a "year-specific mean time", for you to be able to test your ideas more quickly than with the old transmit routine.

[sgbmzm]

wow It's really wonderful. I did a test on 50 years using the old method including exporting a csv file, and it took 3 minutes and 10 seconds (from the moment I confirmed the restart kernel). The same using the new method took only about 15 seconds (from the moment I confirmed restart kernel). So if you now run the following code for 3000 years, you get in 15 minutes a csv file where you can see the difference between the average transit according to tt_calendar() and the time 12:00. It becomes slightly more significant before 800 AD. Do you not agree with me that day-transit, on an annual average, in national time, at zero longitude, should be exactly at - 12:00? For example: Regarding the year 500 I think, maybe the reason for the second difference is because of a slightly different delta-t calculation from reality. Because I understand that when I set Skyfield's timescale to be with delta-t zero, then Skyfield subtracts the delta-t from the result. So if delta t is minus 5 seconds for year 500, then mean transit according to tt_calendar() will be exactly at 12:00. What do you think?

Here is the code I used. I put the old code inside ''' ''' so it won't be active.

import csv

# סקייפילד
from skyfield.api import N, S, E, W, wgs84, load, load_file
from skyfield.framelib import ecliptic_frame
from skyfield.earthlib import refraction
from skyfield import almanac
from skyfield.units import Angle
from skyfield.nutationlib import iau2000b_radians

from skyfield.searchlib import find_discrete

# השורה הזו חשובה מאוד והיא מועתקת מתוך ספריית האלמנך
find_discrete

# תאריכים
import datetime as dt
from datetime import datetime, timedelta, date
from pytz import timezone

# חישובים
from numpy import cos, zeros_like, arccos

eph = load('de441s.bsp')

utc = timezone('utc')

location_timezone = utc

location = wgs84.latlon(0 * N, 0 * E)

ts = load.timescale(delta_t=0.0)

# יצירת מילון לאחסון תוצאות
results_dict = {}

# לולאה על כל שנה מ-2 עד 2999
for year in range(1, 3000):
    # יצירת תאריך ל-17 ליוני בשנה הנוכחית
    target_time = datetime(year, 6, 17)

    time = utc.localize(target_time)

    # הגדרת הזמן לשעה שש בבוקר
    time_6 = time.replace(hour=6, minute=0, second=0, microsecond=0)

    '''
    # חישוב חצות היום והלילה בחיפוש בין שש בבוקר לשש בבוקר למחורת: נזכיר כי הזמן הוא כבר בשעון מקומי
    meridian_transit = almanac.meridian_transits(eph, eph["sun"], location)
    transit, y = almanac.find_discrete(ts.from_datetime(time_6), ts.from_datetime(time_6 + dt.timedelta(hours=24)), meridian_transit)
    day_transit = transit[0].astimezone(location_timezone)
    night_transit = transit[1].astimezone(location_timezone)

    # הגדרת תאריך ההתחלה והסיום של השנה הגרגוריאנית המבוקשת, בשעון מקומי
    datetime_greg_year_start = location_timezone.localize(datetime(target_time.year, 1, 1, 0, 0 ,0))
    datetime_greg_year_end = location_timezone.localize(datetime(target_time.year, 12, 31, 23, 59 ,59))

    # חישוב חצות היום והלילה עבור כל יום מימי השנה המבוקשת
    transits, y = almanac.find_discrete(ts.from_datetime(datetime_greg_year_start), ts.from_datetime(datetime_greg_year_end), meridian_transit)
    # איסוף כל חצות היום בלבד לפי השעה בשעון גריניץ עבור כל אחד מימי השנה לתוך מערך חדש
    transits = [transits[i].tt_calendar() for i in range(len(transits)) if y[i] == 1]
    '''

     # חישוב חצות היום והלילה בחיפוש בין שש בבוקר לשש בבוקר למחורת: נזכיר כי הזמן הוא כבר בשעון מקומי
    observer = eph['earth'] + location
    start = ts.from_datetime(time_6)
    end = ts.from_datetime(time_6 + dt.timedelta(hours=24))
    transits = almanac.find_transits(observer, eph['sun'], start, end)
    day_transit = transits[0].astimezone(location_timezone)

    # הגדרת תאריך ההתחלה והסיום של השנה הגרגוריאנית המבוקשת, בשעון מקומי
    datetime_greg_year_start = location_timezone.localize(datetime(target_time.year, 1, 1, 0, 0 ,0))
    datetime_greg_year_end = location_timezone.localize(datetime(target_time.year, 12, 31, 23, 59 ,59))

    # חישוב חצות היום והלילה עבור כל יום מימי השנה המבוקשת
    start = ts.from_datetime(datetime_greg_year_start)
    end = ts.from_datetime(datetime_greg_year_end)
    transits = almanac.find_transits(observer, eph['sun'], start, end)

    # איסוף כל חצות היום בלבד לפי השעה בשעון גריניץ עבור כל אחד מימי השנה לתוך מערך חדש
    transits = [t.tt_calendar() for t in transits]

    # המרה של שעת חצות היום של כל יום מימי השנה מ-שעות דקות ושניות - ל-שניות
    total_seconds = sum(time[3] * 3600 + time[4] * 60 + time[5] for time in transits)

    # חישוב ממוצע של מספר השניות שזה בעצם חישוב של חצות ממוצע עבור כל השנה
    average_seconds = total_seconds // len(transits)

    # הפיכת הזמן שהתקבל לזמן פייתון (תחילת יממת שנת 1900 זה לא בדווקא)
    mean_time = datetime(1900, 1, 1) + timedelta(seconds=average_seconds)

    # Convert mean time back to string format
    mean_time_str = mean_time.strftime('%H:%M:%S')

    # חישוב חצות מקומי אמיתי של היום המבוקש בשעון גריניץ
    transit_utc = day_transit.astimezone(timezone('utc'))

    # שעון גריניץ של חצות הממוצע
    transit_mean = transit_utc.replace(hour=mean_time.hour, minute=mean_time.minute, second=mean_time.second, microsecond=mean_time.microsecond)

    delta_t = ts.from_datetime(day_transit).delta_t

    transit_12 = transit_utc.replace(hour=12, minute=0, second=0, microsecond=0)

    # בדיקת פער הזמן כלומר הדלתא בין שעה מקומית ממוצעת של חצות לבין השעה 12 בצהריים ובדיקה האם יש להוסיף דקות או להפחית דקות
    if transit_mean < transit_12:
        error_delta_t = transit_12-transit_mean
        seconds_error_delta_t = error_delta_t.total_seconds()
        str_error_delta_t = f'+{str(error_delta_t)[2:7]}'
    else:
        error_delta_t = transit_mean-transit_12
        seconds_error_delta_t = - error_delta_t.total_seconds()
        str_error_delta_t = f'-{str(error_delta_t)[2:7]}'

    repaired_delta_t = delta_t + seconds_error_delta_t # חיבור כולל בתוכו חיסור כשיש מינוס

    # עדכון המילון בערך המתוקן של דלטא_טי עבור השנה המתאימה
    results_dict[year] = repaired_delta_t

    #print("original_delta_t",delta_t)
    #print("seconds_eror_delta_t", seconds_error_delta_t)
    #print("repaired_delta_t",repaired_delta_t)
    #print("transit_utc",transit_utc)
    #print("transit_mean",transit_mean)

#results_dict

# מכאן והלאה שמירת הנתונים לקובץ סי.אס.וי
# נתיב לקובץ CSV
csv_file_path = 'tttttttccccc_seconds.csv'

# כותב CSV
with open(csv_file_path, 'w', newline='') as csvfile:
    # הגדרת עמודות ה-CSV
    fieldnames = ['Year', 'Seconds']

    # יצירת אובייקט מכתב עבור CSV
    writer = csv.DictWriter(csvfile, fieldnames=fieldnames)

    # כתיבת כותרות עמודות
    writer.writeheader()

    # כתיבת נתונים לקובץ CSV
    for year, seconds in results_dict.items():
        writer.writerow({'Year': year, 'Seconds': seconds})

[sgbmzm]

Here is the new code in clean form, with an English translation. No need to create a csv. The data is displayed as a dictionary. ("almanac.find_transits" needs the latest version of Skyfield)

from skyfield.api import wgs84, load
from skyfield import almanac
from datetime import datetime, timedelta, date

#Ephemeris
eph = load('de441.bsp')

# יצירת מילון לאחסון תוצאות
#Creating a dictionary to store results
results_dict = {}

# לולאה על כל שנה מ-1 עד 2999
# Loop over each year from 1 to 2999
for year in range(1, 2999):

    # סולם זמן של סקייפילד בלי דלטא-טי
    # Timescale of Skyfield without delta-T
    ts = load.timescale(delta_t = 0.0)

    # חישוב חצות היום בלבד בקו אורך אפס בשיטה החדשה של סקייפילד עבור כל יום מימי השנה המבוקשת
    # Calculation of transit at zero longitude using Skyfield's new method for each day of the requested year
    transits = almanac.find_transits(eph['earth'] + wgs84.latlon(0, 0), eph['sun'], ts.utc(year, 1, 1, 0, 0 ,0), ts.utc(year, 12, 31, 23, 59 ,59))

    # המרת כל זמני החצות לזמנים שלהם ב-יו.טי.סי. לתוך מערך חדש
    # Convert all transit times to their UTC times. into a new array
    transits = [t.utc for t in transits] #or t.ut1_calendar() or t.tt_calendar()

    # המרה של שעת חצות היום של כל יום מימי השנה מ-שעות דקות ושניות - ל-שניות ואיסוף סכום סופי של כל השניות
    # Conversion of the transit time of each day of the year from hours, minutes and seconds - to seconds and collecting a final sum of all seconds
    total_seconds = sum(time[3] * 3600 + time[4] * 60 + time[5] for time in transits)

    # חישוב ממוצע של מספר השניות שזה בעצם חישוב של חצות ממוצע עבור כל השנה
    # Average calculation of the number of seconds which is basically a calculation of an average transit for the whole year
    average_seconds = total_seconds // len(transits)

    # הפיכת הזמן שהתקבל לזמן פייתון (תחילת יממה ראשונה של השנה זה לא בדווקא אלא יכול להיות כל שנה)
    # Converting the resulting time to Python time
    transit_mean = datetime(year, 1, 1) + timedelta(seconds=average_seconds)

    # הפיכת הזמן הממוצע לשעה 12:00
    # Converting the average time to 12:00
    transit_12 = transit_mean.replace(hour=12, minute=0, second=0, microsecond=0)

    # בדיקת פער הזמן כלומר הדלתא בין שעה מקומית ממוצעת של חצות לבין השעה 12 בצהריים ובדיקה האם יש להוסיף דקות או להפחית דקות
    # Checking the time gap, i.e. the delta between the average local time of the transit and 12 noon and checking whether minutes should be added or minutes reduced
    if transit_mean < transit_12:
        error_delta_t = transit_12-transit_mean
        seconds_error_delta_t = error_delta_t.total_seconds()
    else:
        error_delta_t = transit_mean-transit_12
        seconds_error_delta_t = - error_delta_t.total_seconds()

    # עדכון המילון בערך המתוקן של דלטא_טי עבור השנה המתאימה
    # Updating the dictionary with the corrected value of delta_t for the corresponding year
    results_dict[year] = seconds_error_delta_t

# הדפסת המילון
# Print the dictionary
results_dict

[brandon-rhodes]

I spent a couple of hours this weekend looking at your sun-averaging idea, and I would indeed have thought that averaging the Sun's position would have resulted in a match for the ∆T value. But, like you, I see a several-second offset between them.

One adjustment I made was to not use transits, because they are not even samples — they are clustered closer together at aphelion, and are farther apart at perihelion, so they don't evenly sample the Sun's position.

But that only helped a little.

If someday I have time for a deep dive into all the math, I'll let you know what I find; but until then, the few-second different at year 500 (and the growing difference before that) will, alas, remain a mystery!

For your application, I would leave ∆T at its default value, and for each day, record its anti transit, sunrise, solar noon, and sunset, because those are the four things that the ancient observer would be able to measure; they wouldn't, for example, be able to track the difference between apparent time and our modern concept of mean time. Then use those four reference points to define the hours of the day and watches of the night.

[sgbmzm]

I am very grateful for all the helpful discussion. The new transit function was very helpful. For myself, I think the method I suggested is the best method for me right now, and I use it in my software to fix the delta-t zero time scale. By the way, all the astronomical calculations in my software are according to skyfield, so I thank you for everything. My software was written in Hebrew, but I translated parts of it into English. My software is free and can be downloaded from my website here: https://sites.google.com/view/cochavim-uzmanim/ Here is a picture of the front of my software. I took the general inspiration from an Android app called "sol et umbra" but it calculates data only for the sun, and my software also for the moon, the planets, and some bright stars. I recommend also looking at the "now updating" option in the "times and events" menu. It takes a new calculation every second, and then you can see the movement of the sun and other bodies.