Accuracy claimed is plus or minus one minute, more than sufficient for Pagan ceremonial purposes.
You will need to specify the calendar year Greg_year, the required event xx, and your time zone TZ.
Code: Select all
'=======================================================================
' Equinoxes and Solstices, to nearest minute.
' Chapter 26. Astronomical Algorithms. Jean Meeus. 1991.
'=======================================================================
' user parameters
Dim As Integer Greg_Year = 2013
Dim As Integer xx = 0 ' select event, see Function event() at line 20
'-------------------------------------------------------------------
Dim Shared As Double dt = +35 ' TAI - UTC
Dim Shared As Integer TZ = + 10 ' time zone in hours, AEST Australia = +10
If (xx = 0) Or (xx = 3) Then TZ += 1 ' adjust for summertime
' The tracking of dynamical or ephemeris time with UTC is effected by the
' variation of Atmospheric Angular Momentum (AAM) requires leap seconds
' to be inserted into UTC. On Feb 15, 2013, TAI-UTC = 35.0 seconds.
' Get the new correction from here; ftp://maia.usno.navy.mil/ser7/ser7.dat
'=======================================================================
' x = 0 march equinox
' x = 1 june solstice
' x = 2 september equinox
' x = 3 december solstice
'-----------------------------------------------------------------------
Function event(Byval x As Integer, Byval greg_year As Integer) As Double
'-------------------------------------------------------------------
Dim As Integer e = 0 ' default to early era, 1000 BC to 1000 AD
Dim As Double y = (Greg_Year / 1000.) ' Greg_Year must be an integer
If Greg_Year > 1000 Then
e = 1 ' the period is from 1000 AD to 3000 AD
y = y - 2 ' this adjusts y to; y = ( Greg_year - 2000 ) / 1000
End If
Dim As Double y2 = y * y, y3 = y2 * y, y4 = y2 * y2 ' powers of y
'-------------------------------------------------------------------
' polynomial coefficients d( era , event , term )
Static As Double d(0 To 1, 0 To 3, 0 To 4) => _
{{ _ ' earlier era 1000 BC to 1000 AD Table 26.A page 166
{+1721139.29189, +365242.13740, +0.06134, +0.00111, -0.00071}, _ ' march
{+1721233.25401, +365241.72562, -0.05323, +0.00907, +0.00025}, _ ' june
{+1721325.70455, +365242.49558, -0.11677, -0.00297, +0.00074}, _ ' sept
{+1721414.39987, +365242.88257, -0.00769, -0.00933, -0.00006} _ ' dec
},{ _ ' later era 1000 AD to 3000 AD Table 26.B page 166
{+2451623.80984, +365242.37404, +0.05169, -0.00411, -0.00057}, _ ' march
{+2451716.56767, +365241.62603, +0.00325, +0.00888, -0.00030}, _ ' june
{+2451810.21715, +365242.01767, -0.11575, +0.00337, +0.00078}, _ ' sept
{+2451900.05952, +365242.74049, -0.06223, -0.00823, +0.00032}}} ' dec
'-------------------------------------------------------------------
' this gives the instant of the "mean" equinox or solstice
Dim As Double JDEo = d(e,x,0) + y*d(e,x,1) + y2*d(e,x,2) + y3*d(e,x,3) + y4*d(e,x,4)
' which can now be adjusted with the periodic terms
'-------------------------------------------------------------------
Type periodic ' table of periodic terms
As Double a, b, c
End Type
#define abc Type ' hide the Type keyword from the formating
Static As periodic term(0 To 23) => { _ ' Table 26.C page 167
_ ' a b c a b c
abc(485, 324.96, 1934.136), abc( 45, 247.54, 29929.562), _
abc(203, 337.23, 32964.467), abc( 44, 325.15, 31555.956), _
abc(199, 342.08, 20.186), abc( 29, 60.93, 4443.417), _
abc(182, 27.85, 445267.112), abc( 18, 155.12, 67555.328), _
abc(156, 73.14, 45036.886), abc( 17, 288.79, 4562.452), _
abc(136, 171.52, 22518.443), abc( 16, 198.04, 62894.029), _
abc( 77, 222.54, 65928.934), abc( 14, 199.76, 31436.921), _
abc( 74, 296.72, 3034.906), abc( 12, 95.39, 14577.848), _
abc( 70, 243.58, 9037.513), abc( 12, 287.11, 31931.756), _
abc( 58, 119.81, 33718.147), abc( 12, 320.81, 34777.259), _
abc( 52, 297.17, 150.678), abc( 9, 227.73, 1222.114), _
abc( 50, 21.02, 2281.226), abc( 8, 15.45, 16859.074)}
#undef abc
'-------------------------------------------------------------------
#define deg *(Atn(1)/45) ' convert degrees to radians
Dim As Double t = ( JDEo - 2451545.0 ) / 36525
Dim As Double w = t * 35999.373 - 2.47
Dim As Double delta_lambda = 1 + 0.0334 * Cos( w deg ) + 0.0007 * Cos( 2 * w deg )
Dim As Double sum = 0
For i As Integer = 0 To 23 ' evaluate periodic terms, order not important
With term(i)
sum += ( .a * Cos(( .b + .c * t) deg ) )
End With
Next i
#undef deg
'-------------------------------------------------------------------
Return JDEo + (0.00001 * sum / delta_lambda ) ' JDE
'-------------------------------------------------------------------
End Function
'=======================================================================
' convert JDE back to Gregorian calendar, year, month, day and time
Sub JD2GC(Byval JD As Double,_
Byref gy As Integer,_
Byref gm As Integer,_
Byref gd As Double) ' includes time as fraction
Dim As Double t = JD + 0.5
'-------------------------------------------------------------------
t += (TZ / 24) ' rotate JDE to the local time zone
t -= (dt / (24 * 3600)) ' apply the current TAI-UTC correction
'-------------------------------------------------------------------
Dim As Double z = Int(t) ' integer part
Dim As Double f = t - z ' fractional part
Dim As Integer a, b, c, d, e
If z < 2299161 Then
a = z
Else
Dim As Integer temp = Int((z - 1867216.25) / 36524.25 )
a = z + 1 + temp - Int(temp / 4)
End If
b = a + 1524
c = Int( ( b - 122.1 ) / 365.25 )
d = Int( 365.25 * c )
e = Int( ( b - d ) / 30.6001 )
gd = b - d - Int( 30.6001 * e ) + f
If e < 14 Then gm = e - 1 Else gm = e - 13
If gm > 2 Then gy = c - 4716 Else gy = c - 4715
End Sub
'=======================================================================
' test the function
Dim As Double jde = event(xx, Greg_Year)
Dim As Integer gy, gm, i
Dim As Double gd, f
JD2GC( jde, gy, gm, gd)
Print " year "; gy
Print "month "; gm
i = Int(gd)
f = gd - i
Print " day "; i
f = f * 24
i = Int(f)
f = f - i
Print i; " hour,";
f = f * 60
i = Int(f)
f = f - i
Print i; " min,";
f = f * 60
Print Int(f + 0.5); " sec. Local time."
'=======================================================================
Print
Print "Approximate test data for the year 2013."
Print "Event "; xx, "Timezone "; TZ
Print "0 Autumn equinox Mar 20, 22:01 daylight saving"
Print "1 Winter solstice Jun 21, 15:04 standard time"
Print "2 Spring equinox Sep 23, 06:44 standard time"
Print "3 Summer solstice Dec 22, 04:11 daylight saving"
'=======================================================================
Sleep
'=======================================================================