// This program calculates and plots the date of the winter solstice // over the course of a few centuries. I made this because I'm old // and I remember the equinoxes usually being on the 21st or the // 22nd (probably in UTC.) This shows how the date of the equinox drifts over // time (and is corrected back by leap year rules.) use sun.frink use planets.frink use secant.frink longfunc = {|date| highAccuracySunApparentLongitude[date]} // Let's calculate the longitude angle that corresponds to 1 millisecond of // time angleres = ms circle/solaryear tz = "US/Mountain" startyear = 1800 endyear = 2150 g = new graphics g.font["Monospaced", .1 day] p = new polyline for y = startyear to endyear { winterSolstice = secantInvert[longfunc, 270 degrees, winterSolstice[y] - 1 hour, winterSolstice[y] + 1 hour, angleres] d = winterSolstice - parseDate["$y $tz"] println["$y\t" + (winterSolstice -> tz) + "\t" + format[d, "days", 2]] p.addPoint[y,-d] g.fillEllipseCenter[y, -d, 1, .01 day] } g.color[0,0,0,.3] g.add[p] for y = startyear to endyear step 50 { g.color[0,0,0,.3] g.line[y, -353 days, y, -357 days] g.color[0,0,0] g.text["$y", y, -353 days, "center", "top"] } df = ### MM-dd ### for d = 353 days to 357 days step day { g.color[0,0,0,.3] g.line[startyear, -d, endyear, -d] usualDate = (parseDate["2015 $tz"] + d) -> [df, tz] leapDate = (parseDate["2016 $tz"] + d) -> [df, tz] g.color[0,0,0] g.text[format[d,"days",0] + " \n(usually " + usualDate + " \n on leap years, $leapDate) ", startyear, -d, "right", "center"] } g.color[0,0,0,0] g.line[startyear-100, -353 days, startyear-100, -357 days] g.color[0,0,0] g.font["SansSerif", .15 day] g.text["Winter Solstice (days after beginning of year, $tz)", 1975, -353.1 days, "center", "bottom"] g.show[] stz = tz =~ %s/\//./g g.write["wintersolstice$stz.png",2000,1000] g.write["wintersolstice$stz.html",2000,1000] g.write["wintersolstice$stz.svg",2000,1000] browse["wintersolstice$stz.html"] browse["wintersolstice$stz.svg"]