skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.38k stars 208 forks source link

Incorrect values reported by almanac.dark_twilight_day() #879

Closed isa-kain closed 11 months ago

isa-kain commented 11 months ago

Hello! I'm trying to use dark_twilight_day(), but it seems like it's reporting values that are 3 hours off (plus some change -- maybe this indicates a time zones problem?). Here is the script I am running, mostly taken from the docs:

import datetime as dt
from pytz import timezone
from skyfield import almanac
from skyfield.api import N, W, wgs84, load

zone = timezone('HST') # i've tried multiple time zones, still get incorrect results
now = zone.localize(dt.datetime.now())
tomorrow = now + dt.timedelta(days=1)

ts = load.timescale()
t0 = ts.from_datetime(now)
t1 = ts.from_datetime(tomorrow)

eph = load('de421.bsp')
keckcoords = wgs84.latlon(19.8260 * N, -155.4747 * W)

f = almanac.dark_twilight_day(eph, keckcoords)
times, events = almanac.find_discrete(t0, t1, f) 

for t, e in zip(times, events):

    timestring = str( t.astimezone(zone) )[:16] # Terrestrial Time as Julian date -> convert by time zone
    event = almanac.TWILIGHTS[e]                # day; civil / nautical / astronomical twilight; night

    print( timestring, event )

Here are the outputs:

2023-07-17 22:20 Civil twilight
2023-07-17 22:44 Nautical twilight
2023-07-17 23:12 Astronomical twilight
2023-07-17 23:41 Night
2023-07-18 07:47 Astronomical twilight
2023-07-18 08:16 Nautical twilight
2023-07-18 08:44 Civil twilight
2023-07-18 09:08 Day

…and here are the correct times:

2023-07-17 19:24 Civil twilight
2023-07-17 19:49 Nautical twilight
2023-07-17 20:18 Astronomical twilight
2023-07-17 20:48 Night
2023-07-18 04:41 Astronomical twilight
2023-07-18 05:11 Nautical twilight
2023-07-18 05:41 Civil twilight
2023-07-18 06:05 Day
brandon-rhodes commented 11 months ago

Good morning! The script gave me different output, but by the time I ran it now() was of course a half-day later, so I adjusted to:

now = zone.localize(dt.datetime(2023, 7, 17, 12))

— and was able to get the same output you did.

One problem, I think, is that you are both using the W constant, whose value is -1, and also putting a negative sign in front of 155.4747 to make it -155.4747. Which means that your literal negative sign and the negative one inside of W are canceling each other out, providing the positive value +155.4747 to .latlon(), which means ‘east’ rather than ‘west’. Which, alas, places the position somewhere out in the mid-Pacific between Wake Island and Guam.

So I tried changing the line to:

keckcoords = wgs84.latlon(19.8260 * N, 155.4747 * W)

And got something more vaguely similar to the output you want:

2023-07-17 19:03 Civil twilight
2023-07-17 19:27 Nautical twilight
2023-07-17 19:56 Astronomical twilight
2023-07-17 20:25 Night
2023-07-18 04:30 Astronomical twilight
2023-07-18 05:00 Nautical twilight
2023-07-18 05:28 Civil twilight
2023-07-18 05:52 Day

But, alas, that's still off by quite a bit from the times you want! Could you share how you calculated the list of desired times? Maybe there's a difference in how they define twilight, and how Skyfield does.

brandon-rhodes commented 11 months ago

As I've not heard back after three weeks, I'm going to (perhaps too optimistically?) assume that the fix for your longitude has mostly resolved your problem. If you'd like us to pursue any remaining questions about the times you are getting back, please feel free to make further comments on this issue, and we'll keep investigating! But for the moment I'm going to click "Close" to keep Skyfield's dashboard clear.