Problem with dark_twilight_day

[OH6BG]

I am really excited about this new feature, dark_twilight_day. Thank you so much!

However, I think it's slightly misbehaving, only giving one "category 4" date-time, namely sunrise. There should also be another "category 4 event" for sunset, as both sunrise and sunset happen at -0.8333 degrees.

I tested this with this script:

from datetime import datetime, timedelta
from skyfield.api import load, Topos
from skyfield.almanac import find_discrete, sunrise_sunset, dark_twilight_day
import pytz

ts = load.timescale()
e = load('de421.bsp')

lat, lon = 63.1460, 21.5420
start_date = datetime.strptime("2020-12-26", "%Y-%m-%d")
end_date = datetime.strptime("2020-12-27", "%Y-%m-%d")
date_list = [start_date + timedelta(days=x) for x in range((end_date-start_date).days)]

tzn, elv = 'UTC', 0
tz = pytz.timezone(tzn)
qth = Topos(lat, lon, elevation_m=elv)

for date in date_list:

    ldt = tz.localize(date)
    t0 = ts.utc(ldt)
    t1 = ts.utc(tz.normalize(ldt + timedelta(1)))

    print("SUNRISE & SUNSET")
    t, y = find_discrete(t0, t1, sunrise_sunset(e, qth))
    for tt, yt in zip(t, y):
        sun = "Sunrise" if yt else "Sunset"
        print(sun, tt.utc_iso())

    print()
    """
    TWILIGHTS
    0: Dark of night.
    1: Astronomical twilight.
    2: Nautical twilight.
    3: Civil twilight.
    4: Daytime.
    """
    print("TWILIGHTS")
    s, z = find_discrete(t0, t1, dark_twilight_day(e, qth))
    for ti, yi in zip(s, z):
        print(yi, ti.utc_iso())

    print()

The result is this:

SUNRISE & SUNSET
Sunrise 2020-12-26T08:13:50Z
Sunset 2020-12-26T12:55:30Z

TWILIGHTS
1 2020-12-26T04:59:11Z
2 2020-12-26T05:56:46Z
3 2020-12-26T07:02:08Z
4 2020-12-26T08:13:50Z <-- this is sunrise
3 2020-12-26T12:55:30Z <-- this is sunset = 4 ?
2 2020-12-26T14:07:12Z <-- this should be 3 ?
1 2020-12-26T15:12:35Z <-- this should be 2 ?
0 2020-12-26T16:10:11Z <-- this should be 1 ?

Thanks for your consideration!

[brandon-rhodes]

Good question! The event integer like 4 or 3 indicates the new state of affairs that is being entered, so sunset is 3 because it is the moment when the Sun enters the region of -6° to -0.8333° altitude.

How can the documentation be rephrased to make that behavior more clear? I'd like to update it to prevent confusion.

(And, thanks for providing a sample script!)

[OH6BG]

Thanks for the prompt reply! I think that, to prevent confusion in this matter, we may even need better granularity for the zones. I think this can be achieved if the direction of the Sun's movement is known, i.e. are we going towards Day or towards Night within the zone?

Let me illustrate.

Currently, for a "regular" day, the sequence or cycle of the zones can be as follows:

0 1 2 3 4 3 2 1

However, what if we first indicate the direction with + (= "towards Day") and - (= "towards Night") signs as follows:

+0 +1 +2 +3 +4 -3 -2 -1

Now, we know that the change of direction typically happens in the zones 0 and 4, so consequently we could represent the zones like this:

+0 +1 +2 +3 +4 -4 -3 -2 -1 -0

where +0 would mean "the Night Zone after midnight" and -0 "the Night Zone before midnight". Similarly, +4 would mean "the Day Zone before noon" and -4 "the Day Zone after noon".

As you know, there are, however, irregular cycles close to polar areas so the change of direction is not limited to happen in the zones 0 and 4 only. For example, during the polar night period, I think the following cycles may occur:

+0 +1 +2 +3 -3 -2 -1 -0

or perhaps even this:

+0 +1 +2 -2 -1 -0

During a polar day period, the daily cycle (i.e. the changes in the direction of the Sun's movement towards Day or Night) may even happen in one zone only:

+4 -4 +4

So, to recap, for practical purposes I would appreciate it very much if the direction of movement could be indicated, which would then effectively lead to a more granular representation of the day, twilight and night zones.

I hope this makes any sense at all :)

[brandon-rhodes]

That's an interesting proposal!

But I see two issues.

First, the value being measured — altitude — cannot by itself be used to make the distinction “before noon / after noon” so your proposal (I think?) would require two separate almanac searches that then would need to be combined. That's beyond the scope of a single almanac.py routine, which are each designed to search for one simple thing.

Second, it would no longer be possible for users to easily look up the names of phases; an abs() call would be necessary now in all cases, breaking previous code that was working with the old routines. So the proposal looks incompatible with software that is already up and running on Skyfield?

Have you tried computing the sign yourself?

NumPy has a routine np.diff() which, I think, will tell you for each new solar position whether it is on its way up or down. Could you try it out and see if it gives you the result you are after? If it works, I would be happy for you to suggest a new paragraph for the documentation that would help show other users how np.diff() could let them solve this same problem you have run into. Try it out and let me know what happens, thanks!

[ghost]

Not sure if this is helpful, but the Sun's azimuth value, combined with which hemisphere you're in, can tell if the sun is getting higher (azimuth between 0 and 180 in the northern hemisphere) or lower (azimuth between 180 and 360 in the northern hemisphere).

Because the sun's declination changes, there are picayune cases where this won't work, but I think you can safely ignore them-- it's more of an issue for the moon: https://astronomy.stackexchange.com/questions/937/can-moon-set-after-being-up-when-due-north-pseudo-circumpolarity

[foehnwall]

@brandon: "How can the documentation be rephrased to make that behavior more clear? I'd like to update it to prevent confusion."

I think following naming would make the situation clearer: LEGEND = {0: 'begin of dark night ', 1: 'begin of astr. twilight ', 2: 'begin of naut. twilight ', 3: 'begin of civ. twilight ', 4: 'begin of daytime '}

Did I get anything wrong?

Kind regards, Foehn

[brandon-rhodes]

So you do not currently have any problems with the documentation itself, but you would also like a dictionary of names to be provided right in the code?

The names should probably just be (for example) Astronomical twilight rather than "begin of", because the function that takes a time t does not answer 1 only for the beginning of astronomical twilight, but answers 1 over and over again for the whole range of times from the very start of astronomical twilight to its very end.

[foehnwall]

Thank you for your reply.

".. but you would also like a dictionary of names to be provided right in the code?" No, of course not. I just did a lazy and sloppy copy/paste from my own source code.

If I understand you right, 0,1, 3,... are values of a state variable that indicate discrete states of sun altitudes. I did not realize that before. I thought the values marks a "phase change". So my proposal for an improvement of the documentation is obsolete.

Sorry for the inconvenience!

[brandon-rhodes]

No problem! I'm happy to have added the extra docs and the dictionary, and I think your questions helped us make Skyfield a little more useful for its other users.