The roofers are here, with very large and impressive equipment.
They are starting at my end of the building, and I am on the top floor.
The old roof is being torn off and it is dramatic. All this to say I will share the
main reference points for the math I was planning to blog about today but
leave the discussion for tomorrow. No puttering around the place for me today!
Yesterday I looked at Longitude, the East/West position of Montreal and Edmonton,
Alberta.
360 degrees, divided over 24 hours, gives 15 degrees. So knowing the Longitude
of both cities allows us to see they are 2 hours and 40 minutes apart.
But to discuss patterns of enlightenment, we need Latitude, the North/South component.
And that is more complicated because the great sphere that represents Earth is itself
bent by 23 degrees and its position with respect to the Sun changes with time of the
year and time of day.
Below, the formula proposed by Copilot:
Gemini's explanation:
https://share.gemini.google/Hr9g5bHDbZXD
* * *
Gemini also gave me a Python script which allowed me to generate daily enlightenment
for a city from its position numbers. The daylight hours are blue curve minus the orange
for any given time.
The code:
import matplotlib.pyplot as plt
import numpy as np
# --- 1. CONFIGURATION & LOCATION PARAMETERS ---
latitude_deg = 45.31 # Saint-Jean-sur-Richelieu latitude
longitude_deg = -73.25 # Saint-Jean-sur-Richelieu longitude (West is negative)
lat_rad = np.radians(latitude_deg)
tz_offset = -4.0 # Forced EDT (UTC-4) year-round
# --- 2. ANNUAL SOLAR CALCULATIONS (365 Days) ---
n = np.arange(1, 366) # Day of year (1 to 365)
gamma = (2 * np.pi / 365) * (n - 1) # Fractional year angle
# Solar Declination (radians) - Spencer (1971) approximation
delta = (
0.006918
- 0.399912 * np.cos(gamma)
+ 0.070257 * np.sin(gamma)
- 0.006758 * np.cos(2 * gamma)
+ 0.000907 * np.sin(2 * gamma)
- 0.002697 * np.cos(3 * gamma)
+ 0.00148 * np.sin(3 * gamma)
)
# Equation of Time (minutes)
eot = 229.18 * (
0.000075
+ 0.001868 * np.cos(gamma)
- 0.032077 * np.sin(gamma)
- 0.014615 * np.cos(2 * gamma)
- 0.040849 * np.sin(2 * gamma)
)
# Solar Noon (UTC)
# Corrected: West longitude adds time to UTC because Montreal is behind Greenwich
solar_noon_utc = 12.0 + (-longitude_deg / 15.0) - (eot / 60.0)
# Hour Angle ($\omega_s$) accounting for refraction (-0.833 degrees)
alt_rad = np.radians(-0.833)
cos_omega = (np.sin(alt_rad) - np.sin(lat_rad) * np.sin(delta)) / (
np.cos(lat_rad) * np.cos(delta)
)
cos_omega = np.clip(cos_omega, -1.0, 1.0)
omega_s = np.degrees(np.arccos(cos_omega))
# --- 3. TIME CONVERSION ---
sunrise_utc = solar_noon_utc - (omega_s / 15.0)
sunset_utc = solar_noon_utc + (omega_s / 15.0)
# Convert UTC to uniform EDT (UTC-4)
sunrise_local = sunrise_utc + tz_offset
sunset_local = sunset_utc + tz_offset
# --- 4. PLOTTING ---
fig, ax = plt.subplots(figsize=(10, 6), dpi=100)
ax.plot(
n, sunrise_local, color="#d97706", linewidth=2, label="Sunrise (EDT Year-Round)"
)
ax.plot(
n, sunset_local, color="#4f46e5", linewidth=2, label="Sunset (EDT Year-Round)"
)
ax.fill_between(
n,
sunrise_local,
sunset_local,
color="#fef3c7",
alpha=0.5,
label="Daylight Span",
)
ax.set_title(
f"Annual Sunrise & Sunset Curves — EDT Year-Round ({latitude_deg}°N)",
fontsize=12,
pad=15,
)
ax.set_xlabel("Day of Year", fontsize=10)
ax.set_ylabel("Local Time (Hours)", fontsize=10)
month_days = [1, 32, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335]
month_labels = [
"Jan",
"Feb",
"Mar",
"Apr",
"May",
"Jun",
"Jul",
"Aug",
"Sep",
"Oct",
"Nov",
"Dec",
]
ax.set_xticks(month_days)
ax.set_xticklabels(month_labels)
ax.set_xlim(1, 365)
ax.set_ylim(4, 23)
ax.set_yticks(np.arange(4, 24, 2))
ax.grid(True, linestyle="--", alpha=0.5)
ax.legend(loc="upper center", frameon=True)
plt.tight_layout()
plt.show()
For St-Jean-sur-Richelieu:
* * *
And because a bit of fun is always welcome as one trudges through Math, the Wikipedia
page on 'Terminator', a concept in this scientific field:
Terminator (solar) - Wikipedia
🎬🍿
