Calculating the date of a full moon

Given the year and month number, I'm trying to create an Excel function (in VBA) that will return the date (and if possible, also the time) when a full moon will occur.

As my 2nd hobby is astrophysics, I have spent some time to search the web for some calculations and code.

After a while I've found an Excel file, which calculates the moon phases relative precisely, just a few seconds difference to the exact moon phases.

I've adapted the code for your needs a little.

Here is the link to the file I've found: http://members.aon.at/excelapps/excelapps.htm#kalender

Andreas.

Option Explicit

Public Enum MondPhases
  NewMoon = 0    'Neumond
  WaxingMoon = 1 'Zunehmender Mond
  FullMoon = 2   'Vollmond
  WaningMoon = 3 'Abnehmender Mond
End Enum

Sub Example_MONDPHASE()
  Const TimeZone = 1 'means GMT+1
  'Calculate the moon phases around today
  Debug.Print MONDPHASE(Now, NewMoon, TimeZone)
  Debug.Print MONDPHASE(Now, WaxingMoon, TimeZone)
  Debug.Print MONDPHASE(Now, FullMoon, TimeZone)
  Debug.Print MONDPHASE(Now, WaningMoon, TimeZone)
End Sub

Function MONDPHASE(ByVal ThisDate As Date, ByVal Phase As MondPhases, _
    Optional ByVal TimeZone As Integer) As Date
  'From a code by Michael Friedrich

  'found at http://members.aon.at/excelapps/excelapps.htm#kalender
  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

  Y = Year(ThisDate) + CDbl(ThisDate - DateSerial(Year(ThisDate), 1, 1)) / 365.25

  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 = BOGENMASS(WINKELREDUZIEREN(M))
  MS = BOGENMASS(WINKELREDUZIEREN(MS))
  F = BOGENMASS(WINKELREDUZIEREN(F))
  Omega = BOGENMASS(WINKELREDUZIEREN(Omega))

  'Argumente der Planeten

  A1 = BOGENMASS(WINKELREDUZIEREN(299.77 + 0.107408 * K - 0.009173 * T ^ 2))
  A2 = BOGENMASS(WINKELREDUZIEREN(251.88 + 0.016321 * K))
  A3 = BOGENMASS(WINKELREDUZIEREN(251.83 + 26.651886 * K))
  A4 = BOGENMASS(WINKELREDUZIEREN(349.42 + 36.412478 * K))
  A5 = BOGENMASS(WINKELREDUZIEREN(84.66 + 18.206239 * K))
  A6 = BOGENMASS(WINKELREDUZIEREN(141.74 + 53.303771 * K))
  A7 = BOGENMASS(WINKELREDUZIEREN(207.14 + 2.453732 * K))
  A8 = BOGENMASS(WINKELREDUZIEREN(154.84 + 7.30686 * K))
  A9 = BOGENMASS(WINKELREDUZIEREN(34.52 + 27.261239 * K))
  A10 = BOGENMASS(WINKELREDUZIEREN(207.19 + 0.121824 * K))
  A11 = BOGENMASS(WINKELREDUZIEREN(291.34 + 1.844379 * K))
  A12 = BOGENMASS(WINKELREDUZIEREN(161.72 + 24.198154 * K))
  A13 = BOGENMASS(WINKELREDUZIEREN(239.56 + 25.513099 * K))
  A14 = BOGENMASS(WINKELREDUZIEREN(331.55 + 3.592518 * K))

  'Periodische Terme Neumond
  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

  'Periodische Terme Vollmond
  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

  'Periodische Terme erstes und letztes 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)

  MONDPHASE = KALENDERDATUM(JDE + PK + PT + W + (TimeZone / 24) - DELTAT(Year(ThisDate)) / 86400)
End Function

Private Function KALENDERDATUM(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)
  KALENDERDATUM = 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 BOGENMASS(Winkel)
  BOGENMASS = Winkel * Atn(1) * 4 / 180
End Function
Private Function WINKELREDUZIEREN(Winkel)
  WINKELREDUZIEREN = Winkel - Int(Winkel / 360) * 360
End Function