exactSeasons.frink

View or download exactSeasons.frink in plain text format

/*
   These algorithms reference Jean Meeus book _Astronomical Algorithms_

   This program uses the full VSOP87 theory to calculate the exact dates
   of the equinoxes and solstices to the nearest millisecond.  While the
   sun.frink library has routines to calculate the solstices and equinoxes
   rather simply, these can be inaccurate up to several minutes or so.

   This program uses very high precision routines to calculate the solar
   longitude of earth for each season, and uses a secant method solver to
   obtain the exact time of the solstices and equinoxes to the nearest
   millisecond.

   The definition of the equinoxes and solstices are when earth's heliocentric
   *longitude* is equal to 180 degrees (spring), 270 degrees (summer),
   0 degrees (autumn), and 90 degrees (winter.)

   The solstices are defined in terms of longitude, not latitude.  Thus, you
   should be warned that the sun should be quite close to being over the
   equator at equinoxes, and close to being at max/min latitude at solstices,
   but these will not be exact.  Beware anyone who tells you differently.
*/


use planets.frink
use secant.frink

// Let's calculate the longitude angle that corresponds to 1 millisecond of
// time
angleres = ms circle/solaryear

timezone = "US/Mountain"
// Function to use planets.frink to calculate the heliocentric longitude
// at the specified time.
longfunc = {|date| highAccuracySunApparentLongitude[date]}

for yearnum = 1996 to 2020
{
   println[]

   // These functions use the low-precision equinox/solstice functions from
   // sun.frink to estimate a good starting point.  The low-accuracy guess is
   // within 51 seconds for the years 1951-2050 (see Meeus table 27.D).
   // We then use the secantInvert function to use the secant method to invert
   // the higher-accuracy sun longitude function fron planets.frink using the
   // full VSOP-87 theory.
   // The sun longitude function gives the apparent longitude, given a
   // date/time.  We invert this to find the date/time that the sun's apparent
   // longitude matches a given value.
   // We bound the guesses of the secant method to one hour before or after
   // the initial guess. 
   springEquinox = secantInvert[longfunc,
                                180 degrees,
                                springEquinox[yearnum] - 1 hour,
                                springEquinox[yearnum] + 1 hour,
                                angleres]
   println["Spring equinox is:   " + (springEquinox -> timezone)]


   summerSolstice = secantInvert[longfunc,
                                 270 degrees,
                                 summerSolstice[yearnum] - 1 hour,
                                 summerSolstice[yearnum] + 1 hour,
                                 angleres]
   println["Summer solstice is:  " + (summerSolstice -> timezone)]

   autumnEquinox = secantInvert[longfunc,
                                0 degrees,
                                autumnEquinox[yearnum] - 1 hour,
                                autumnEquinox[yearnum] + 1 hour,
                                angleres]
   println["Autumn equinox is:   " + (autumnEquinox -> timezone)]

   winterSolstice = secantInvert[longfunc,
                                 90 degrees,
                                 winterSolstice[yearnum] - 1 hour,
                                 winterSolstice[yearnum] + 1 hour,
                                 angleres]

   println["Winter solstice is:  " + (winterSolstice -> timezone)]
}


View or download exactSeasons.frink in plain text format


This is a program written in the programming language Frink.
For more information, view the Frink Documentation or see More Sample Frink Programs.

Alan Eliasen was born 18666 days, 1 hours, 40 minutes ago.