For those that are interested.... this is what I am working on. It is for moonphases, moonrise, moonset, etc. This is only part of it. It was originally written by Michael Friedrich for excel, and I am adapting it for .net.
Option Explicit On
Imports System.Math
Module MoonPhases
Public Enum MoonPhases
NewMoon = 0 'New Moon
WaxingMoon = 1 'Waxing Moon
FullMoon = 2 'Full Moon
WaningMoon = 3 'Waning Moon
End Enum
Sub Example_MoonPhase()
Const TimeZone = 1 'means GMT+1
'Calculate the moon phases around today
Debug.Print(MoonPhase(Now, MoonPhases.NewMoon, TimeZone))
Debug.Print(MoonPhase(Now, MoonPhases.WaxingMoon, TimeZone))
Debug.Print(MoonPhase(Now, MoonPhases.FullMoon, TimeZone))
Debug.Print(MoonPhase(Now, MoonPhases.WaningMoon, TimeZone))
End Sub
Function MoonPhase(ByVal ThisDate As Date, ByVal Phase As MoonPhases, Optional ByVal TimeZone As Integer = 0) As Date
'From a code by Michael Friedrich
Dim T, JDE, E, M, MS, F, Omega
Dim A1, A2, A3, A4, A5, A6, A7, A8, A9, A10, A11, A12, A13, A14
Dim PT, W, PK
Dim K, Y
'Original line of code I am trying to rework so it will function
'Y = Year(ThisDate) + CDbl(ThisDate - DateSerial(Year(ThisDate), 1, 1)) / 365.25
Dim JanuaryFirst As DateTime = New DateTime(ThisDate.Year, 1, 1)
Dim tdiff As TimeSpan = ThisDate - JanuaryFirst
Dim ddiff As Double = tdiff.TotalDays / 365.25 '188 / 365.25
Dim resultDate As DateTime = JanuaryFirst.AddDays(ddiff)
Dim strDDiff As String = ddiff.ToString
Dim intIdx As Int16
intIdx = strDDiff.IndexOf(".")
strDDiff = strDDiff.Substring(intIdx, strDDiff.Length - 1)
Y = ThisDate.Year.ToString & strDDiff
K = Int((Y - 2000) * 12.3685)
K = K + Phase * 0.25
T = K / 1236.85
JDE = 2451550.09765 + 29.530588853 * K + 0.0001337 * T ^ 2 - 0.00000015 * T ^ 3 + 0.00000000073 * T ^ 4
E = 1 - 0.002516 * T - 0.0000074 * T ^ 2
M = 2.5534 + 29.10535669 * K - 0.0000218 * T ^ 2 - 0.00000011 * T ^ 3
MS = 201.5643 + 385.81693528 * K + 0.0107438 * T ^ 2 + 0.00001239 * T ^ 3 - 0.000000058 * T ^ 4
F = 160.7108 + 390.67050274 * K - 0.0016341 * T ^ 2 - 0.00000227 * T ^ 3 + 0.000000011 * T ^ 4
Omega = 124.7746 - 1.5637558 * K + 0.0020691 * T ^ 2 + 0.00000215 * T ^ 3
M = ARCMASS(ANGULARREDUCTION(M))
MS = ARCMASS(ANGULARREDUCTION(MS))
F = ARCMASS(ANGULARREDUCTION(F))
Omega = ARCMASS(ANGULARREDUCTION(Omega))
'Planet Arguments
A1 = ARCMASS(ANGULARREDUCTION(299.77 + 0.107408 * K - 0.009173 * T ^ 2))
A2 = ARCMASS(ANGULARREDUCTION(251.88 + 0.016321 * K))
A3 = ARCMASS(ANGULARREDUCTION(251.83 + 26.651886 * K))
A4 = ARCMASS(ANGULARREDUCTION(349.42 + 36.412478 * K))
A5 = ARCMASS(ANGULARREDUCTION(84.66 + 18.206239 * K))
A6 = ARCMASS(ANGULARREDUCTION(141.74 + 53.303771 * K))
A7 = ARCMASS(ANGULARREDUCTION(207.14 + 2.453732 * K))
A8 = ARCMASS(ANGULARREDUCTION(154.84 + 7.30686 * K))
A9 = ARCMASS(ANGULARREDUCTION(34.52 + 27.261239 * K))
A10 = ARCMASS(ANGULARREDUCTION(207.19 + 0.121824 * K))
A11 = ARCMASS(ANGULARREDUCTION(291.34 + 1.844379 * K))
A12 = ARCMASS(ANGULARREDUCTION(161.72 + 24.198154 * K))
A13 = ARCMASS(ANGULARREDUCTION(239.56 + 25.513099 * K))
A14 = ARCMASS(ANGULARREDUCTION(331.55 + 3.592518 * K))
'Periodic Terms New Moon
If Phase = 0 Then
PT = -0.4072 * Sin(MS) + 0.17241 * E * Sin(M) + 0.01608 * Sin(2 * MS) _
+ 0.01039 * Sin(2 * F) + 0.00739 * E * Sin(MS - M) - 0.00514 * E * Sin(MS + M) _
+ 0.00208 * E ^ 2 * Sin(2 * M) - 0.00111 * Sin(MS - 2 * F) - 0.00057 * Sin(MS + 2 * F) _
+ 0.00056 * E * Sin(2 * MS + M) - 0.00042 * Sin(3 * MS) + 0.00042 * E * Sin(M + 2 * F) _
+ 0.00038 * E * Sin(M - 2 * F) - 0.00024 * E * Sin(2 * MS - M) _
- 0.00017 * Sin(Omega) - 0.00007 * Sin(MS + 2 * M) + 0.00004 * Sin(2 * MS - 2 * F) _
+ 0.00004 * Sin(3 * M) + 0.00003 * Sin(MS + M - 2 * F) + 0.00003 * Sin(2 * MS + 2 * F) _
- 0.00003 * Sin(MS + M + 2 * F) + 0.00003 * Sin(MS - M + 2 * F) - 0.00002 * Sin(MS - M - 2 * F) _
- 0.00002 * Sin(3 * MS + M) + 0.00002 * Sin(4 * MS)
End If
'Periodic Terms New Moon
If Phase = 2 Then
PT = -0.40614 * Sin(MS) + 0.17302 * E * Sin(M) + 0.01614 * Sin(2 * MS) _
+ 0.01043 * Sin(2 * F) + 0.00734 * E * Sin(MS - M) - 0.00515 * E * Sin(MS + M) _
+ 0.00209 * E ^ 2 * Sin(2 * M) - 0.00111 * Sin(MS - 2 * F) - 0.00057 * Sin(MS + 2 * F) _
+ 0.00056 * E * Sin(2 * MS + M) - 0.00042 * Sin(3 * MS) + 0.00042 * E * Sin(M + 2 * F) _
+ 0.00038 * E * Sin(M - 2 * F) - 0.00024 * E * Sin(2 * MS - M) _
- 0.00017 * Sin(Omega) - 0.00007 * Sin(MS + 2 * M) + 0.00004 * Sin(2 * MS - 2 * F) _
+ 0.00004 * Sin(3 * M) + 0.00003 * Sin(MS + M - 2 * F) + 0.00003 * Sin(2 * MS + 2 * F) _
- 0.00003 * Sin(MS + M + 2 * F) + 0.00003 * Sin(MS - M + 2 * F) - 0.00002 * Sin(MS - M - 2 * F) _
- 0.00002 * Sin(3 * MS + M) + 0.00002 * Sin(4 * MS)
End If
'Periodic terms first and last phase
If Phase = 1 Or Phase = 3 Then
PT = -0.62801 * Sin(MS) + 0.17172 * E * Sin(M) - 0.01183 * E * Sin(MS + M) _
+ 0.00862 * Sin(2 * MS) + 0.00804 * Sin(2 * F) + 0.00454 * E * Sin(MS - M) _
+ 0.00204 * E ^ 2 * Sin(2 * M) - 0.0018 * Sin(MS - 2 * F) - 0.0007 * Sin(MS + 2 * F) _
- 0.0004 * Sin(3 * MS) - 0.00034 * E * Sin(2 * MS - M) + 0.00032 * E * Sin(M + 2 * F) _
+ 0.00032 * E * Sin(M - 2 * F) - 0.00028 * E ^ 2 * Sin(MS + 2 * M) + 0.00027 * E * Sin(2 * MS + M) _
- 0.00017 * Sin(Omega) - 0.00005 * Sin(MS - M - 2 * F) + 0.00004 * Sin(2 * MS + 2 * F) _
- 0.00004 * Sin(MS + M + 2 * F) + 0.00004 * Sin(MS - 2 * M) + 0.00003 * Sin(MS + M - 2 * F) _
+ 0.00003 * Sin(3 * M) + 0.00002 * Sin(2 * MS - 2 * F) + 0.00002 * Sin(MS - M + 2 * F) _
- 0.00002 * Sin(3 * MS + M)
W = 0.00306 - 0.00038 * E * Cos(M) + 0.00026 * Cos(MS) _
- 0.00002 * Cos(MS - M) + 0.00002 * Cos(MS + M) + 0.00002 * Cos(2 * F)
If Phase = 3 Then W = -W
End If
PK = 0.000325 * Sin(A1) + 0.000165 * Sin(A2) + 0.000164 * Sin(A3) _
+ 0.000126 * Sin(A4) + 0.00011 * Sin(A5) + 0.000062 * Sin(A6) _
+ 0.00006 * Sin(A7) + 0.000056 * Sin(A8) + 0.000047 * Sin(A9) _
+ 0.000042 * Sin(A10) + 0.00004 * Sin(A11) + 0.000037 * Sin(A12) _
+ 0.000035 * Sin(A13) + 0.000023 * Sin(A14)
MoonPhase = CALENDARDATE(JDE + PK + PT + W + (TimeZone / 24) - DELTAT(Year(ThisDate)) / 86400)
End Function
Private Function CALENDARDATE(JD) As Date
Dim Z, F, Alpha, A, B, C, D, E, M, Y, H, Mi, S
JD = JD + 0.5
Z = Int(JD)
F = JD - Z
If Z >= 2299161 Then
Alpha = Int((Z - 1867216.25) / 36524.25)
A = Z + 1 + Alpha - Int(Alpha / 4)
Else
A = Z
End If
B = A + 1524
C = Int((B - 122.1) / 365.25)
D = Int(365.25 * C)
E = Int((B - D) / 30.6001)
D = B - D - Int(30.6001 * E) + F
If E < 14 Then
M = E - 1
Else
M = E - 13
End If
If M > 2 Then
Y = C - 4716
Else
Y = C - 4715
End If
H = D - Int(D)
H = H * 24
Mi = H - Int(H)
Mi = Mi * 60
S = Mi - Int(Mi)
S = S * 60
H = Int(H)
Mi = Int(Mi)
CALENDARDATE = DateSerial(Y, M, Int(D)) + " " + TimeSerial(H, Mi, S)
End Function
Private Function DELTAT(Y)
Dim T, Tage
T = (Y - 2000) / 100
If Y < 948 Then DELTAT = 2715.6 + 573.36 * T + 46.5 * T ^ 2
If Y >= 948 And Y <= 1600 Then DELTAT = 50.6 + 67.5 * T + 22.5 * T ^ 2
If Y >= 1800 And Y < 1900 Then
T = (Y - 1900) / 100
Tage = -0.000009 + 0.003844 * T + 0.083563 * T ^ 2 + 0.865736 * T ^ 3 _
+ 4.867575 * T ^ 4 + 15.845535 * T ^ 5 + 31.332267 * T ^ 6 _
+ 38.291999 * T ^ 7 + 28.316289 * T ^ 8 + 11.636204 * T ^ 9 _
+ 2.043794 * T ^ 10
DELTAT = Tage * 86400
End If
If Y >= 1900 And Y < 1988 Then
T = (Y - 1900) / 100
Tage = -0.00002 + 0.000297 * T + 0.025184 * T ^ 2 - 0.181133 * T ^ 3 _
+ 0.55304 * T ^ 4 - 0.861938 * T ^ 5 + 0.677066 * T ^ 6 _
- 0.212591 * T ^ 7
DELTAT = Tage * 86400
End If
If Y >= 1988 And Y < 2051 Then
DELTAT = (Y - 1990) * 6.6 / 9 + 56.86
End If
End Function
Private Function ARCMASS(Angles)
ARCMASS = Angles * Atan(1) * 4 / 180
End Function
Private Function ANGULARREDUCTION(Angles)
ANGULARREDUCTION = Angles - Int(Angles / 360) * 360
End Function
Private Function JulianDate(ByVal d As Integer, ByVal m As Integer, ByVal y As Integer) As Integer
Dim mm, yy As Integer
Dim k1, k2, k3 As Integer
Dim julian As Integer
yy = y - Math.Floor((12 - m) / 10)
mm = m + 9
If (mm >= 12) Then
mm = mm - 12
End If
k1 = Math.Floor(365.25 * (yy + 4712))
k2 = Math.Floor(30.6001 * mm + 0.5)
k3 = Math.Floor(Math.Floor((yy / 100) + 49) * 0.75) - 38
' "j" for dates in Julian calendar:
julian = k1 + k2 + d + 59
If (julian > 2299160) Then
' For Gregorian calendar:
julian = julian - k3 ' "j" is the Julian date at 12h UT (Universal Time)
End If
Return julian
End Function
End Module