// Calculate Durin's Day, from The Hobbit. // // "The first day of the dwarves' New Year, said Thorin, is as all should // know the first day of the last moon of Autumn on the threshold of Winter. // We still call it Durin's Day when the last moon of Autumn and the sun are // in the sky together. But this will not help us much, I fear, for it passes // our skill in these days to guess when such a time will come again." // // Many of these definitions are vague, and we use some conventions to put // the dates earlier than a strict astronomical reading might yield. The // astronomical figures are rigorous, but the clarity of Thorin's definition // is in some question. use sun.frink // How many years should be checked? years = 25 // How chatty should we be? verbose[] := { // true: report everything // false: report New Year's Days and Durin's Days return true } y = parseInt[input["Enter year: ",2000]] limit = y + years // GPS-surveyed position for east side of MIT's Infinite corridor, Cambridge, // Massachusetts. lat = 42.36002 degrees North long = 71.09332 degrees West // NWS for Denver, CO //lat = 39.78 degrees North //long = 104.88 degrees West // Greenwich, England //lat = DMS[51,28,38] North //long = 0 degrees West // Loop until we find a New Year's Day with: // a sunrise or sunset when the moon is entirely above the horizon // AND the separation between sun and moon is at least 7.5 degrees, // which has been proven to be minimal for seeing a crescent. // See http://www.skyandtelescope.com/observing/objects/projects/3308686.html // and http://www.earthsky.org/faq/young-moon-visibility do { found = false dawn = false nyd = NewYearsDay[y] print[" New Year's Day occurs on "] if (verbose[]) println["\$nyd."] else println[(nyd -> ###yyyy-MM-dd###) + "."] if (verbose[]) printSunMoonPos[nyd, lat, long] sunrise = sunrise[nyd, lat, long] sunset = sunset[nyd, lat, long] if (verbose[]) println[" Sunrise is at \$sunrise"] [moonalt, sep] = printSunMoonPos[sunrise, lat, long] if (moonalt > moonRadiusAngle[sunrise]) and (sep > 7.5 degrees) { found = true dawn = true } if (verbose[]) println[" Sunset is at \$sunset"] [moonalt, sep] = printSunMoonPos[sunset, lat, long] if (moonalt > moonRadiusAngle[sunset]) and (sep > 7.5 degrees) { found = true dawn = false } if (! found) and (verbose[]) { println[" This is not Durin's Day."] println[] } if (found) { print["Durin's day occurs at "] if (dawn) print["dawn"] else print["dusk"] println[" on " + (nyd -> ### MMMM d, yyyy ###) + "."] if (verbose[]) println[] } y = y + 1 } while (y < limit) // Print the position of sun and moon at a specified time. printSunMoonPos[nyd, lat, long] := { [sunaz, sunalt] = refractedSunAzimuthAltitude[nyd, lat, long] [moonaz, moonalt] = refractedMoonAzimuthAltitude[nyd, lat, long] sep = sunMoonAngularSeparation[nyd, lat, long] if (verbose[]) { illum = moonIlluminatedFraction[nyd] illum1 = moonIlluminatedFraction[nyd + 1 minute] if (illum1 > illum) indicator = "waxing" else indicator = "waning" println[" At this time, the altitude of the moon is: " + format[moonalt, "degrees" , 5]] println[" Angular separation: " + format[sep,"degrees",3]] println[" Moon illuminated fraction: " + format[illum, percent, 3] + "% (\$indicator)"] println[] } return [moonalt, sep] } // Find New Year's Day, the first day of the last month before the // cross-quarter day between the autumn equinox and winter solstice. NewYearsDay[y] := { autumn = autumnEquinox[y] winter = winterSolstice[y] endautumn = autumn + (winter-autumn)/2 // Find the last new moon of "autumn." nextmoon = endautumn do { newmoon = newMoon[nextmoon] nextmoon = newmoon - lunarmonth } while newmoon > endautumn // Return the new moon return newmoon }