(ANSWERED) apparent asymmetrical effect of season on morning vs evening monthly overall civil twilight periods

[jrtuenge]

I've adapted your When will it get dark tonight? code, and was surprised to find that for any given location, if you calculate a month's overall civil twilight periods, its morning MOCTP can differ substantially from its evening MOCTP:

• (latest sunrise in month) - (earliest civil dawn in month)
• (latest civil dusk in month) - (earliest sunset in month)

For example, I get a 20-minute difference in duration for Bristol Twp, PA, in Jan 2014:

• morning MOCTP = (7:21 on Jan 1) - (6:40 on Jan 31) = 41 minutes
• evening MOCTP = (17:45 on Jan 31) - (16:44 on Jan 1) = 61 minutes

I confirmed that I get essentially identical results on timeanddate.com for January 2014 — Sun in Township of Bristol, but have no easy means of determining whether any bug might be shared by their code and Skyfield (e.g., they could be using Skyfield or identical algorithm under the hood).

I found an equal difference (with same sign) for Sep, and absolute differences of 1 to 13 minutes for eight other months (excluding the DST-change months of Nov and Mar here for simplicity). My expectation was that morning and evening MOCTPs would be essentially identical, given the symmetry in sun path charts. Which of these four explanations applies here?

1. This is a known/understood phenomenon, ideally documented somewhere (e.g., by Meeus) or easily explained (hopefully I'm not just missing something obvious).
2. There's a rounding error somewhere, which could be resolved locally via improved resolution.
3. There's a bug in Skyfield (and similar algorithms).
4. Something else.

P.S. The identical sunrise times from timeanddate.com for Jan 1 through Jan 11 might offer a possible clue?

[brandon-rhodes]

If we could back up for a moment first — what does MOCTP stand for? If you explained the acronym, I haven't found it yet in your question. Second: why are you subtracting the minutes between an event on January 1 and on January 31? That doesn't make sense to me yet.

[jrtuenge]

Sorry, I thought it'd be clear from the first sentence that "month's overall civil twilight periods" = MOCTP.

I'm looking at MOCTP because some datasets (e.g., PennDOT Public Crash Datasets) do not give date, only month and time of day.

But if you think about it, it shouldn't really matter whether/why I'm looking at MOCTP -- I seem to have identified an underlying discrepancy in the process.

[brandon-rhodes]

Sorry, I thought it'd be clear from the first sentence that "month's overall civil twilight periods" = MOCTP.

Thanks!

I'm looking at MOCTP because some datasets (e.g., [PennDOT Public Crash Datasets] https://crashinfo.penndot.gov/PCIT/welcome.html)) do not give date, only month and time of day. But if you think about it, it shouldn't really matter whether/why I'm looking at MOCTP -- I seem to have identified an underlying discrepancy in the process.

I agree that the reason is unimportant — though it’s interesting to hear what you’re working on, so thanks for sharing.

Let’s move on to the second part of my question: why are you subtracting the end of civil twilight of one date from the beginning of civil twilight on another date? I can’t quite make sense of what the meaning of the resulting number would be.

To know the number of minutes of, say, morning civil twilight over a month, wouldn’t you need to go through every single day, subtract that particular day’s beginning-of-civil-twilight from its end-of-civil-twilight, and then add up all the resulting minutes?

[jrtuenge]

Take a crash for which we have month and time of day (but not date). Did it happen at night, during civil twilight, or in daylight? The strictly correct answer is that we cannot say. But we can define civil twilight more widely, such that its morning and evening periods span the times I defined in my initial post here. Although this looser definition should clearly be expected to result in some daytime and nighttime crashes being misclassified as occurring during civil twilight, it should result in no crashes being misclassified as occurring during day or night.

So I'm not counting minutes spent in civil twilight over the month -- I'm instead tracing the morning and evening periods of civil twilight across the month, and using the boundaries of these two traces to set unusually stringent (non-standard) bounds for daytime and nighttime (with everything else loosely denoted "civil twilight").

[brandon-rhodes]

Although this looser definition should clearly be expected to result in some daytime and nighttime crashes being misclassified as occurring during civil twilight, it should result in no crashes being misclassified as occurring during day or night.

I think I have to disagree (at least slightly?), because January 1 is not the January date of the latest sunrise — and thus not of the latest end-of-civil-twilight. That occurs, typically, on January 4.

I think what you want for your algorithm is min(all_civil_twilight_ends_for_january) and max(all_civil_twilight_beginnings_for_january) — a calculation that checks every day of the month, not just the month’s endpoints?

[jrtuenge]

I think I have to disagree (at least slightly?), because January 1 is not the January date of the latest sunrise — and thus not of the latest end-of-civil-twilight. That occurs, typically, on January 4.

I had actually noticed the "Latest sunrise is on January 4" text (which still makes no sense to me) on the linked timeanddate.com page. Except for the following four months, I'd expect the relevant dates to be at start and end of month, given that those dates are closest to the two solstices (which by definition have the longest and shortest daytime periods of the year).

• Mar (DST change)
• Jun (solstice)
• Nov (DST change)
• Dec (solstice)

But as you can see, it makes no meaningful difference for Jan, so let's please put that detail aside for the time being. Moving on to your response here:

I think what you want for your algorithm is min(all_civil_twilight_ends_for_january) and max(all_civil_twilight_beginnings_for_january)

No, there are two issues with your proposed algorithm:

1. Taken literally, you'd get just two periods (which you might loosely call day and night), rather than the four periods (roughly corresponding to day, night, morning CT, evening CT) that I'm seeking. In addition, it'd draw one boundary at earliest sunrise of 7:09, and the other at latest sunset of 17:17 (which is not at all my intent).
2. And if I instead interpret your algorithm as being applied to morning and evening CT periods separately (yielding something more akin to the four periods I'm seeking), you get narrowed CT periods (whereas I had been seeking widened CT periods as I had explained here).

Using interpretation 2 of your approach above, I actually get a negative duration for the evening period:

• morning = (7:09 on Jan 31) - (6:51 on Jan 4) = 18 minutes
• evening = (17:15 on Jan 4) - (17:17 on Jan 31) = -2 minutes

The definitions I gave in my initial initial post reflect my intent, which is for wide (rather than narrow) windows for MOCTP:

• morning MOCTP = (latest sunrise in month) - (earliest civil dawn in month)
• evening MOCTP = (latest civil dusk in month) - (earliest sunset in month)

So to return to my original point: Is the following statement not true?

• Durations of morning MOCTP and evening MOCTP (as I've defined them above) should be equal for Jan 2014 in Bristol Twp, PA.

The diagram on this page, which I had linked to previously, might facilitate visualization. For example, note that night by definition includes nautical and astronomical (but not civil) twilight.

P.S. I've abbreviated civil twilight as CT, if that wasn't clear.

[brandon-rhodes]

Your interpretation 2 is the one I intended. I do appear to have swapped min() and max() — my intention (and this time I'll spell out the result so there's a double-check on my meaning!) was to grab the earliest CT start and the latest CT end, separately for both morning and evening. Thus, for morning, 6:40am through 7:21am (taking those numbers from the timeanddate.com link); and for evening, 4:44pm through 5:45pm.

But, as you point out, the error in this case is only a couple of minutes, since Jan 4 is so very close to Jan 1, and therefore stands in as a good approximation for it.

But my answer to your overall question is still that: No, I would not expect morning MOCTP to match evening MOCTP, because the symmetry you see in CT morning and CT evening durations only holds for particular days, but not between different days.

I completely agree that on any particular day the two periods of twilight should be very close in duration. As indeed is the case:

from skyfield.api import load, wgs84
from skyfield import almanac

ts = load.timescale()
t0 = ts.utc(2014, 1)
t1 = ts.utc(2014, 2)
bristol = wgs84.latlon(40.1259, -74.8683)
eph = load('de421.bsp')
f = almanac.dark_twilight_day(eph, bristol)
times, events = almanac.find_discrete(t0, t1, f)

# Looking at each day, morning minutes of CT always almost exactly
# equals the evening minutes of CT: exactly the symmetry you correctly
# observed in the website diagrams.

prev_e = None
for ti, ei in zip(times, events):
    if prev_e == 3 and ei == 4:
        morning_ct_minutes = (ti - prev_t) * 24 * 60
    if prev_e == 3 and ei == 2:
        evening_ct_minutes = (ti - prev_t) * 24 * 60
        print(f'{ti.utc_strftime("%b %d")}  '
              f'Morning: {morning_ct_minutes:.2f} minutes  '
              f'Evening: {evening_ct_minutes:.2f} minutes')
    prev_t = ti
    prev_e = ei

Which outputs:

Jan 01  Morning: 30.48 minutes  Evening: 30.48 minutes
Jan 02  Morning: 30.45 minutes  Evening: 30.45 minutes
Jan 03  Morning: 30.41 minutes  Evening: 30.41 minutes
Jan 04  Morning: 30.37 minutes  Evening: 30.38 minutes
Jan 05  Morning: 30.33 minutes  Evening: 30.33 minutes
Jan 06  Morning: 30.29 minutes  Evening: 30.29 minutes
Jan 07  Morning: 30.24 minutes  Evening: 30.24 minutes
Jan 08  Morning: 30.19 minutes  Evening: 30.20 minutes
Jan 09  Morning: 30.14 minutes  Evening: 30.15 minutes
Jan 10  Morning: 30.09 minutes  Evening: 30.09 minutes
Jan 11  Morning: 30.04 minutes  Evening: 30.04 minutes
Jan 12  Morning: 29.98 minutes  Evening: 29.99 minutes
Jan 13  Morning: 29.92 minutes  Evening: 29.93 minutes
Jan 14  Morning: 29.87 minutes  Evening: 29.87 minutes
Jan 15  Morning: 29.81 minutes  Evening: 29.81 minutes
Jan 16  Morning: 29.74 minutes  Evening: 29.75 minutes
Jan 17  Morning: 29.68 minutes  Evening: 29.69 minutes
Jan 18  Morning: 29.62 minutes  Evening: 29.62 minutes
Jan 19  Morning: 29.55 minutes  Evening: 29.56 minutes
Jan 20  Morning: 29.49 minutes  Evening: 29.49 minutes
Jan 21  Morning: 29.42 minutes  Evening: 29.42 minutes
Jan 22  Morning: 29.35 minutes  Evening: 29.36 minutes
Jan 23  Morning: 29.28 minutes  Evening: 29.29 minutes
Jan 24  Morning: 29.21 minutes  Evening: 29.22 minutes
Jan 25  Morning: 29.15 minutes  Evening: 29.15 minutes
Jan 26  Morning: 29.08 minutes  Evening: 29.09 minutes
Jan 27  Morning: 29.01 minutes  Evening: 29.02 minutes
Jan 28  Morning: 28.94 minutes  Evening: 28.95 minutes
Jan 29  Morning: 28.87 minutes  Evening: 28.88 minutes
Jan 30  Morning: 28.80 minutes  Evening: 28.81 minutes
Jan 31  Morning: 28.73 minutes  Evening: 28.74 minutes

The length of morning CT and evening CT are very close, within seconds of each other, every day, which confirms the correct intuition you derived from the clear symmetry in the day/twilight/dark diagram.

But I would expect that symmetry to disappear once you start comparing values between dates, because morning and evening move around so asymmetrically over a month. The Earth's elliptical orbit will, I expect, throw off any attempt to find symmetries between how much morning CT moves during the month compared to evening CT.

Returning to your original question, then, I think we can answer “This is a known/understood phenomenon, ideally documented somewhere (e.g., by Meeus) or easily explained (hopefully I'm not just missing something obvious)” and that you can read about it here:

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

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

You may have been assuming that the Earth’s orbit is circular and that therefore each solar day would be of the same length, with 24 UTC hours between each solar noon? That might be the assumption that led you to expect a meaningful comparison between how far morning CT shifts over a month and how much evening CT shifts. If that’s the case, hopefully the above articles will clear things up!

(Drat, I didn’t mean to split this discussion into two threads; I must have clicked the wrong Reply box. Thanks for patiently following me over to this new thread; I'm still learning to be wary of which text box I click in Discussions; in normal Issues there's only one Reply box, down at the bottom, and I'm used to jumping down to it.)

[jrtuenge]

I do appear to have swapped min() and max() — my intention (and this time I'll spell out the result so there's a double-check on my meaning!) was to grab the earliest CT start and the latest CT end, separately for both morning and evening. Thus, for morning, 6:40am through 7:21am (taking those numbers from the timeanddate.com link); and for evening, 4:44pm through 5:45pm.

Max and min in your proposal mean the same thing as latest and earliest (respectively) in my original definition, so you've reproduced my algorithm -- hence the four times in your results excerpted above matching what I had shared initially:

• morning MOCTP = (7:21 on Jan 1) - (6:40 on Jan 31) = 41 minutes
• evening MOCTP = (17:45 on Jan 31) - (16:44 on Jan 1) = 61 minutes

In any case, we seem to be on the same page now regarding my intent. So to return to my original question: I had seen that Skyfield calculated morning and evening CT periods of roughly equal size for any given day, as expected. What still didn't make sense to me is that day length would increase by shifting morning CT only slightly, while shifting evening CT by a greater amount -- I had expected both CT periods to shift equally (symmetrically about solar noon) over the month.

Although I knew that our orbit wasn't perfectly circular (and that axis of rotation is tilted relative to ecliptic), I had overlooked the fact that the solstices don't coincide with the apsides. The perihelion seemed to potentially clear up the relatively minor Jan 4 issue (regarding latest sunrise), but then I realized I'd need to pinpoint a similar/related extreme (for earliest sunset) in early December -- for example, from Weather Underground:

In the Northern Hemisphere midlatitudes, the earliest sunrises occur a few days before the summer solstice, and the latest sunsets occur a few days after the solstice. The displacement is even more pronounced in winter: the earliest sunsets (as measured in standard time) can occur several weeks before the solstice, and the latest sunrises can happen well into the new year.

The author cites a webpage maintained by Larry Denenberg, which explains the roles of tilt and rotation and orbit, and is excerpted below:

near the solstices, solar noon occurs at a slightly later clock time each day because the solar day is more than 24 hours and the clock day is always exactly 24 hours. Next: Since solar noon is a little later each day, sunset is also a little later each day. (This is true because the time from solar noon to sunset doesn't change very much near the winter solstice.) But if sunset is a little later each day, that means the earliest sunset has already happened! By the same reasoning, later solar noons at the solstice imply later sunrises, hence sunrise is getting later and the latest sunrise is yet to occur.

Using the code below, adapted from yours above, I found that the earliest sunset occurred on Dec 7. I think this confirms that explanation 1 from my initial post is indeed correct (as you suspected), so it seems clear that I need to abandon any expectation of symmetry between morning and evening MOCTPs for any given month. It also clarifies that for January (for example), I cannot assume that its CT boundary extremes will occur on the first or last day of the month. Thanks for pointing me in the right direction -- case closed!

from skyfield.api import load, wgs84
from skyfield import almanac

tz = timezone('US/Eastern')
ts = load.timescale()
t0 = ts.from_datetime(tz.localize(dt.datetime(2013, 12, 1)))
t1 = ts.from_datetime(tz.localize(dt.datetime(2014, 2, 1)))
eph = load('de421.bsp')
bristol = wgs84.latlon(40.1259, -74.8683)
f = almanac.dark_twilight_day(eph, bristol)
times, events = almanac.find_discrete(t0, t1, f)

prev_e = None
for ti, ei in zip(times, events):
    if prev_e == 3 and ei == 4:
        am_ct_start = prev_t
        am_ct_end = ti
    elif prev_e == 3 and ei == 2:
        pm_ct_start = prev_t
        pm_ct_end = ti
        print(f'{ti.utc_strftime("%b %d")}  '
              f'am_start = {am_ct_start.astimezone(tz).strftime("%H:%M:%S")}  '
              f'am_end = {am_ct_end.astimezone(tz).strftime("%H:%M:%S")}  '
              f'pm_start = {pm_ct_start.astimezone(tz).strftime("%H:%M:%S")}  '
              f'pm_end = {pm_ct_end.astimezone(tz).strftime("%H:%M:%S")}')
    prev_t = ti
    prev_e = ei