VB.NET - Equinoxes and Solstices

This formula calculates dates for Spring Equinoxes, Summer Solstices, Autumn Equinoxes and Winter Solstices. However, as with any Calendrical Calculation, programmers need to be aware of certain limitations.

The main function is based on the Meeus books on Calendrecal Calculations. Take care to avoid infringing the books' Copyright etc.
It is not practical for a short function to predict future Equinox and Solstice times with absolute precision.
- There are too many "periodic" variables due to planetary motions, planetory "wobbles" etc. for a short program
- At the time of writing, it is difficult to find precision algorhithms for selected factors (lunar tidal drift, leap-seconds etc.)
- The Meeus formula appears to have an accuracy that could be measured in minutes, this is near-enough for most applications.
The modern (Gregorian) calendar was not used until 15 Oct 1582, and many countries retained the old (Julian) calander for centuries afterwards.
This verion of the formula calculates Gregorian dates - The UK switched calendars in 1752 and the formula can only be used in the UK for dates after 1752.

The following code requires two variables - year (an integer indicating the Gregorian year for required event) and k. k is an integer (always 0, 1, 2 or 3) used to tell the function which event is required. k=0 - Spring Equinox ... k=1 - Summer Solstice ... k=2 - Autumn(Fall) Equinox ... k=3 - Winter Solstice. The returned result is a Date. You can extract years, months, days, hours and minutes from the returned Date.

This function returns the date (and time) in UT (universal time) sometimes called GMT (Greenwich Mean Time) or Time-Zone 0. VB.NET has a built-in function to convert GMT to your local time (adjusting for yor computer's time-zone and seasonal clock adjustments) There is an explanation of how to use this feature at the end of this page.

Note ... Private ReadOnly TBL() ...
This is one long line, the " _" at end of line is used to tell the VB.Net compiler that the line of code continues on next line of editor
Depending on how your .NET version reacts to "pasting" this line into your code, you may have to manually add the " _" afterwards. " _" does not have any quotes, I am just trying to indicate that it is space followed by _ at end of line and the space is important!
The first block of code does not go inside any function or subroutine - we are defining variables which will be needed by the main function


' this goes outside of any functions or sub routines and before the main functions
Private Shared ReadOnly TBL() As Double = {121, 112, 103, 95, 88, 82, _
 77, 72, 68, 63, 60, 56, 53, 51, 48, 46, 44, 42, 40, 38, _
 35, 33, 31, 29, 26, 24, 22, 20, 18, 16, 14, 12, 11, 10, 9, 8, 7, 7, 7, 7, _
 7, 7, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 10, 10, 10, 11, 11, 11, _
 11, 11, 12, 12, 12, 12, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 15, 15, 16, 16, _
 16, 16, 16, 16, 16, 16, 15, 15, 14, 13, _
 13.1, 12.5, 12.2, 12.0, 12.0, 12.0, 12.0, 12.0, 12.0, 11.9, 11.6, 11.0, 10.2, 9.2, 8.2, _
 7.1, 6.2, 5.6, 5.4, 5.3, 5.4, 5.6, 5.9, 6.2, 6.5, 6.8, 7.1, 7.3, 7.5, 7.6, _
 7.7, 7.3, 6.2, 5.2, 2.7, 1.4, -1.2, -2.8, -3.8, -4.8, -5.5, -5.3, -5.6, -5.7, -5.9, _
 -6.0, -6.3, -6.5, -6.2, -4.7, -2.8, -0.1, 2.6, 5.3, 7.7, 10.4, 13.3, 16.0, 18.2, 20.2, _
 21.1, 22.4, 23.5, 23.8, 24.3, 24.0, 23.9, 23.9, 23.7, 24.0, 24.3, 25.3, 26.2, 27.3, 28.2, _
 29.1, 30.0, 30.7, 31.4, 32.2, 33.1, 34.0, 35.0, 36.5, 38.3, 40.2, 42.2, 44.5, 46.5, 48.5, _
 50.5, 52.5, 53.8, 54.9, 55.8, 56.9, 58.3, 60.0, 61.6, 63.0, 63.8, 64.3, 64.6, 64.8, 65.5, _
 66.1, 66.6, 67.3, 68.1, 69.0, 69.4, 69.4, 69.2}

Private Shared ReadOnly A() As Integer = {485, 203, 199, 182, 156, 136, 77, 74, 70, 58, 52, _
 50, 45, 44, 29, 18, 17, 16, 14, 12, 12, 12, 9, 8}
Private Shared ReadOnly B() As Double = {324.96, 337.23, 342.08, 27.85, 73.14, 171.52, 222.54, _
 296.72, 243.58, 119.81, 297.17, 21.02, 247.54, 325.15, 60.93, 155.12, 288.79, 198.04, 199.76, _
 95.39, 287.11, 320.81, 227.73, 15.45}
Private Shared ReadOnly C() As Double = {1934.136, 32964.467, 20.186, 445267.112, 45036.886, _
 22518.443, 65928.934, 3034.906, 9037.513, 33718.147, 150.678, 2281.226, 29929.562, 31555.956, _
 4443.417, 67555.328, 4562.452, 62894.029, 31436.921, 14577.848, 31931.756, 34777.259, _
 1222.114, 16859.074}

' this is the main function that should be called when calculating Equinoxes/Solstices
Private Shared Function Equinox(year As Integer, k As Integer) As Date
 ' This function is based on Calendrical Calculations (1991) by Meeus - Chapter 26 Page 165
 ' Many of teh comments relate to sections within this book.
 ' The formula has been simplified because VB.Net has "built-in" features which offer
 ' simplified alternatives to some aspects of the calculation.
 ' Find a JD "near" to expected date of Equinox/Solstice
 Dim JDE0 As Double : Dim y As Double = (year - 2000) / 1000 ';
 ' k=0...March Equinox k=1...June_Solstice
 Select Case (k)' {
  Case 0 : JDE0 = 2451623.80984 + 365242.37404 * y + 0.05169 * y * y - 0.00411 * y * y * y - 0.00057 * y * y * y * y
  Case 1 : JDE0 = 2451716.56767 + 365241.62603 * y + 0.00325 * y * y + 0.00888 * y * y * y - 0.0003 * y * y * y * y
  Case 2 : JDE0 = 2451810.21715 + 365242.01767 * y - 0.11575 * y * y + 0.00337 * y * y * y + 0.00078 * y * y * y * y
  Case 3 : JDE0 = 2451900.05952 + 365242.74049 * y - 0.06223 * y * y - 0.00823 * y * y * y + 0.00032 * y * y * y * y
 End Select
 ' Calculate sum of periodic terms in table 26.C
 Dim T As Double = (JDE0 - 2451545.0) / 36525 '
 Dim W As Double = (35999.373 * T - 2.47) Mod 360
 Dim dl As Double = 1 + 0.0334 * CosDeg(W) + 0.0007 * CosDeg(2 * W)
 Dim S As Double = 0
 For i As Integer = 0 To 23 ' {
  S += A(i) * CosDeg(B(i) + (C(i) * T)) '; }
 Next
 ' We are now on Page 167
 Dim JDE As Double = JDE0 + ((0.00001 * S) / dl)
 ' Convert JDE to a Gregorian Date (VB.Net AddDays simplifies this task)
 Dim JDE_Greg As New Date(1858, 11, 17, 0, 0, 0, 0, 0)
 JDE_Greg = JDE_Greg.AddDays(JDE - 2400000.5)
 Dim GMT As Date = FromTDTtoUTC(JDE_Greg) ' Tiny adjustment for leap seconds
 Return GMT
End Function ' // End calcEquiSol

Friend Shared Function FromTDTtoUTC(dat As Date) As Date
 ' We only need to consider years after the start of Gregorian Calender
 ' this is trying to calculate number of leap seconds and use that as UTC offset
 ' This function based on formula from NASA Eclipse website
 ' https://eclipse.gsfc.nasa.gov/SEhelp/deltatpoly2004.html
 ' Values for Delta T for 2000 thru 2002 from NASA
 ' Values for Delta T for 2004 thru 2024 from US Naval Observatory
 Dim TBLfirst As Integer = 1620 : Dim TBLlast As Integer = 2024 ' Entries in Lookup Table
 Dim deltaT, t As Double
 Dim daysInYear As Integer = 365
 If Date.IsLeapYear(dat.Year) Then daysInYear = 366
 Dim y As Double = dat.Year + ((dat.DayOfYear - 1) / daysInYear)
 Dim C As Double = (y - 2000) / 100
 deltaT = 0
 ' Check Lookup Table (if entry present it will provide best approximation (usually within seconds!)
 If (dat.Year >= TBLfirst And dat.Year <= TBLlast) Then '{ // Find correction In table
  If (dat.Year Mod 2) Then ' { // Odd year - interpolate
   deltaT = (TBL((dat.Year - TBLfirst - 1) / 2) + TBL((dat.Year - TBLfirst + 1) / 2)) / 2
  Else ' { // Even year - direct table lookup
   deltaT = TBL((dat.Year - TBLfirst) / 2)
  End If
  Return dat.AddSeconds(0 - deltaT)
 End If

 ' If year not in lookup table we use a slightly less accurate estimate with averaged/expected deltaT trends
 Select Case dat.Year
  Case < -500
   t = (y - 1820) / 100
   deltaT = -20 + 32 * t * t
  Case -500 To 499
   t = y / 100
   deltaT = 10583.6 - 1014.41 * t + 33.78311 * t * t - 5.952053 * t * t * t
   deltaT -= 0.1798452 * Math.Pow(t, 4) + 0.022174192 * Math.Pow(t, 5)
   deltaT += 0.0090316521 * Math.Pow(t, 6)
  Case 500 To 1599
   t = (y - 1000) / 100
   deltaT = 1574.2 - 556.01 * t + 71.23472 * t * t + 0.319781 * t * t * t
   deltaT -= 0.8503463 * Math.Pow(t, 4) - 0.005050998 * Math.Pow(t, 5)
   deltaT += 0.0083572073 * Math.Pow(t, 6)
  Case 1600 To 1699
   t = y - 1600
   deltaT = 120 - 0.9808 * t - 0.01532 * t * t + t * t * t / 7129
  Case 1700 To 1799
   t = y - 1700
   deltaT = 8.83 + 0.1603 * t - 0.0059285 * t * t + 0.00013336 * t * t * t
   deltaT -= Math.Pow(t, 4) / 1174000
  Case 1800 To 1859
   t = y - 1800
   deltaT = 13.72 - 0.332447 * t + 0.0068612 * t * t + 0.0041116 * t * t * t
   deltaT -= 0.00037436 * Math.Pow(t, 4) + 0.0000121272 * Math.Pow(t, 5)
   deltaT -= 0.0000001699 * Math.Pow(t, 6) + 0.000000000875 * Math.Pow(t, 7)
  Case 1860 To 1899
   t = y - 1860
   deltaT = 7.62 + 0.5737 * t - 0.251754 * t * t + 0.01680668 * t * t * t
   deltaT -= 0.0004473624 * Math.Pow(t, 4) + Math.Pow(t, 5) / 233174
  Case 1900 To 1919
   t = y - 1900
   deltaT = -2.79 + 1.494119 * t - 0.0598939 * t * t + 0.0061966 * t * t * t
   deltaT -= 0.000197 * Math.Pow(t, 4)
  Case 1920 To 1940
   t = y - 1920
   deltaT = 21.2 + 0.84493 * t - 0.0761 * t * t + 0.0020936 * t * t * t
  Case 1941 To 1960
   t = y - 1950
   deltaT = 29.07 + 0.407 * t - t * t / 233 + t * t * t / 2547
  Case 1961 To 1985
   t = y - 1975
   deltaT = 45.45 + 1.067 * t - t * t / 260 - t * t * t / 718
  Case 1986 To 2004
   t = y - 2000
   deltaT = 63.86 + 0.3345 * t - 0.060374 * t * t + 0.0017275 * t * t * t
   deltaT += 0.000651814 * Math.Pow(t, 4) + 0.00002373599 * Math.Pow(t, 5)
  Case 2005 To 2049
   t = y - 2000
   deltaT = 62.92 + 0.32217 * t + 0.005589 * t * t
  Case 2050 To 2149
   t = y
   deltaT = -20 + 32 * ((t - 1820) / 100) * ((t - 1820) / 100) - 0.5628 * (2150 - t)
  Case >= 2150
   t = (y - 1820) / 100
   deltaT = -20 + 32 * t * t
 End Select

 ' correction linked to how data was gathered
 Dim cor As Double = 0
 If dat.Year < 1955 OrElse dat.Year > 2005 Then
  cor = -0.000012932 * (y - 1955) * (y - 1955)
 End If
 deltaT += cor
 Return dat.AddSeconds(-deltaT)
End Function

Private Shared Function CosDeg(a As Double) As Double
 Return Math.Cos(a * Math.PI / 180)
End Function

GMT or Local Time?

Calling the function will return the date of en equinox or solstice as "Universal Time," this time setting is also known as UT, GMT, Greenwich Mean Time or Time-Zone Zero. The UK operates on this time during the winter months.


  Dim EquinoxDate as date = Equinox(2026,3)

VB.Net has a built in feature which converts Univarsal Time to local time Simply add


   .ToLocalTime()

to a date requiring conversion. This call will return the Equinox time/date using the the local time using timezone/summertime settings from the host computer.


  Dim EquinoxDate as date = Equinox(2026,3).ToLocalTime()

DigitalDan.co.uk