## Calculate Astronomical Seasons

ALX 2018-09-04

# with a deviation less than 15 minutes

proc AstronomicalSeasons {year} {
set Y [expr {(\$year - 2000) / 1000.0}]

# Julian to Gregorian
proc J2G { jda } {
set jd [expr {floor(\$jda + 0.5)}]
set z [expr {\$jda - \$jd + 0.500001}]
set tz [expr {\$jd - 1721117}]
if {\$jd > 2299160} {
set tz [expr {\$tz + (floor((\$tz + 2) / 36524.25) - floor((\$tz + 2) / 146097) - 2)}]
}
set y [expr {floor((\$tz - 0.2) / 365.25)}]
set r [expr {\$tz - floor(\$y * 365.25)}]
set m [expr {floor((\$r - 0.5) / 30.6)}]
set d [expr {\$r - floor(\$m * 30.6 + 0.5)}]
set m [expr {\$m + 3}]
if {\$m > 12} {
set m [expr {\$m - 12}]
set y [expr {\$y + 1}]
}

set h1 [expr {\$z * 24}]
set ho [expr {floor(\$h1)}]
set mi [expr {floor((\$h1 - \$ho) * 60)}]

return [format {%04d-%02d-%02d %02d:%02d} [expr {int(\$y)}] [expr {int(\$m)}] [expr {int(\$d)}] [expr {int(\$ho)}] [expr {int(\$mi)}]]
}

# March-equinox (Beginning of the astronomical spring):
set jda [expr {2451623.80984 + 365242.37404 * \$Y + 0.05169 * pow(\$Y,2) - 0.00411 * pow(\$Y,3) - 0.00057 * pow(\$Y,4)}]
puts "Spring [J2G \$jda]"

# June-solstice (Beginning of the astronomical summer):
set jda [expr {2451716.56767 + 365241.62603 * \$Y + 0.00325 * pow(\$Y,2) + 0.00888 * pow(\$Y,3) - 0.00030 * pow(\$Y,4)}]
puts "Summer [J2G \$jda]"

# September-equinox (Beginning of the astronomical autumn)
set jda [expr {2451810.21715 + 365242.01767 * \$Y - 0.11575 * pow(\$Y,2) + 0.00337 * pow(\$Y,3) + 0.00078 * pow(\$Y,4)}]
puts "Autumn [J2G \$jda]"

# December-solstice (Beginning of the astronomical winter):
set jda [expr {2451900.05952 + 365242.74049 * \$Y - 0.06223 * pow(\$Y,2) - 0.00823 * pow(\$Y,3) + 0.00032 * pow(\$Y,4)}]
puts "Winter [J2G \$jda]"
}

AstronomicalSeasons 2018

Result:

Spring 2018-03-20 16:08
Summer 2018-06-21 10:00
Autumn 2018-09-23 01:45
Winter 2018-12-21 22:17