Формула Харви для возраста Луны

Если вы когда-нибудь решите самостоятельно вычислить фазу Луны для какой-нибудь даты в прошлом или будущем, скажем, жуткое второе февраля 1959 года, то обязательно найдёте ссылку на формулу Харви. С комментариями тех, кто её публикует. И ошибками, как правило.

А также глубокомысленными рассуждениями. Вот, мол, формула неточна, может давать ошибку в один день. А всё из-за того, что деление в ней на 30, а надо бы на 29,53 - такова продолжительность лунного месяца (синодического, нас как землян это интересует, марсиане должны взять месяц сидерический, разумеется).

Ещё могут сказать, что формула Харви описывается словами и не существует в привычном виде, зато прекрасно ложится в компьютерный алгоритм. (Только результаты не радуют: Луна на небе чихать хотела на эти астрологические умствования).

А ещё скажут, что лунный год начинается 1 марта, поэтому для января и февраля год уменьшим на единицу. И чуть позже: для января и февраля добавим 12 месяцев, мы же год уменьшили!

Очаровательно. Просто поток положительных эмоций! Человек, то есть, сидит за компьютером и стенает, что не может умножить на 29,53 и вынужден умножать на 30. Это просто праздник какой-то! не говоря уж про отнятый год и добавленные 12 месяцев.

Но помилуйте! А чем вас не устраивает примитивный прямой счёт? Период обращения Луны известен с точностью до 10-го знака включительно, значит, та же ошибка в один день, которая постоянно путается под ногами по формуле Харви, при прямом счёте никак не может появиться раньше, чем через 27 миллионов лет.

Занесённый в глуши метровыми снегами, я набрался терпения причесать формулу Харви. А также сравнить её результаты с результатами прямого счёта. И с фактической Луной на небе. Строчки можно прямо копировать в Visual Studio, всё работает должным образом. Формулу Харви цитирую (в том числе комментарии), мои правки и примечания оговорены отдельно.

Программа ниже выводит сравнение значений, полученных при помощи формулы Харви, с прямым счётом. Кстати, в причёсанном виде формула Харви не так уж и плоха. Я добирался до 4092 года - результаты были всё так же ноздря в ноздрю с прямым счётом: плюс-минус день. Гений этот Харви, я вам скажу: его формула до сих пор не имеет научного обоснования. Особенно загадочно выглядит число 209. Но если вспомнить, что 209 - это 11 раз 19-летний Метонов цикл, то не особенно загадочно.

Отсутствие обоснования, видимо, вызвано тем, что заняться некому, а вовсе не антинаучностью формулы Харви: её первые аккорды до ужаса напоминают алгоритм перевода григорианской даты в юлианскую. Конечно, старорежимный юлианский календарь удобнее для ручных лунных вычислений, но функция DateDiff (интервал между двумя датами, можно даже в секундах), календарные тонкости учитывает автоматически. Одной строкой кода она сводит на нет все преимущества громоздкого юлианства. Тем более, что количество дней, как количество пней: безразлично, в каком порядке считать. У меня оба алгоритма совпали полностью, минута в минуту.

Прямой счёт точен, прост и очевиден: из общего времени, прошедшего от известного новолуния, надо убрать все полные лунные месяцы - остаток и есть возраст Луны от последнего новолуния. Только _не путать "третий день" и "три дня"_. А если кто заметит расхождение в 4093 году - введите в программу свежую дату фактического новолуния, и опять не будете знать забот ещё 4092 года как минимум.

Принципиальную ошибку округления (+/- полдня), если это важно, можно уменьшить до секунды, задав в нижеприведённой программе для переменной Age тип Single или Double. Эротизма добавляет то обстоятельство, что, во-первых, лунный день легко может смениться следующим прямо средь бела дня земного, во-вторых, 2000-й год - ещё не 2000 лет, число полных столетий в 2000-м году только 19.

В качестве домашнего задания: попробуйте модифицировать программу для формулы Харви с учётом последнего обстоятельства, это ровно одна строка кода.

Итак, в бой!

Private Sub Lunar_Day()
Dim MyDate As Date, kStol as integer, Resl As Double

' С первого числа текущего месяца вычислим на 39 дней вперёд:

    MyDate = 1 & "." & Month(Now) & "." & Year(Now)
    For i = 1 To 39   

' Формула Харви для определения возраста Луны
' Исходные данные: год (четыре цифры), месяц (1 - 12), день. Если месяц - январь или февраль, то из года вычесть единицу (январь и февраль считаются относящимися к предыдущему году).
' Вычислим коэффициент столетия:
' Число полных столетий разделить на 3 и оставить целую часть.
' Число столетий разделить на 4 и тоже оставить целую часть.
' Полученные два числа сложить. Прибавить 6. Вычесть число полных столетий.
    kStol = Int(Int(Year(MyDate) / 100) / 3) + Int(Int(Year(MyDate) / 100) / 4) + 6 - Int(Year(MyDate) / 100)
' Для XX и XXI веков kStol = -3.

' Четыре цифры года разделить на 19, учесть, что лунный год год начинается 1 марта
    Resl = Year(MyDate) / 19
    If Month(MyDate) < 3 Then Resl = (Year(MyDate) - 1) / 19
' дробную часть результата умножить на 209, результат округлить до ближайшего целого.
    Resl = (Resl - Int(Resl)) * 209
    If Resl - Int(Resl) > 0.5 Then Resl = Resl + 1
    Resl = Int(Resl)
' Прибавить месяц. Если это январь или февраль, то прибавить еще 12

' Примечание 1: и получим расшатай на границе февраль-март.

    Resl = Resl + Month(MyDate)
    If Month(MyDate) < 3 Then Resl = Resl + 12
' Прибавить коэффициент столетия и день, разделить на 30 и оставить дробную часть.
    Resl = (Resl + kStol + Day(MyDate)) / 30
    Resl = Resl - Int(Resl)
' Умножить на 30 и округлить до ближайшего целого.

' Примечание 2: лучше умножать не на 30, а на синодический период обращения.

    Resl = Resl * 29.530588
    If Resl - Int(Resl) >= 0.5 Then Resl = Resl + 1
    Resl = Int(Resl)
' Результат является возрастом Луны в заданную дату. День новолуния - нулевой.

' Примечание 3:
' Классическая глупость! в новолуние возраст Луны нулевой, а день-то уже первый! Неужели это так трудно - _не путать "третий день" и "три дня"_?

' Примечание 4:
' Легко видеть, что 01 марта 1900 года было новолуние. (Фактически в 14:25 Msk).

' Примечание 5:
' Так вот зачем лунный год начинается 1 марта!
' Это набежавшие за год ошибки спрятаны в несуществующие числа в конце февраля.

' Примечание 6:
' Эпизодически наблюдается рассинхронизация с Луной
' (+- день из-за разного числа дней в месяцах), затем снова синхронизация.

' Между тем продолжительность лунного месяца известна: 29.5305888531 земных дней.
' Прямой счёт от фактического новолуния в принципе не может дать ошибку в один день
' в ближайшие 27 000 000 лет.

' Для этого из числа дней от новолуния 01.03.1900 до нужной даты вычтем все полные лунные месяцы.
' Оставшееся целое число дней, меньше лунного месяца, и есть возраст Луны на эту дату.

Dim Age As Long  ' тип переменной Long округляет результат автоматически

    Age = DateDiff("d", "01.03.1900", MyDate)
    Age = Age - (Int(Age / 29.5305888531)) * 29.5305888531

' Напечатаем результаты
    List1.AddItem Str(MyDate) + " Harvey = " + Str(Resl) + " DC: " + Str(Age)

' Перейдём к следующему дню, посчитаем и его
    MyDate = MyDate + 1
    Next i

End Sub

' Программа замечательно подходит для экспериментов над ней.
' Небольшое усложнение позволяет вычислить, что в первый день нашей эры было новолуние,
' тогда как формула Харви для ранних веков не годится: календарь с тех пор
' неоднократно корректировали, вырезая то 11 дней, то 14.
' Функция DateDiff учтёт и это.

Что же касается формулы Харви, после преобразований к цивилизованному рабочему виду с учётом январско-февральской коррекции она выглядит так:

Function LunarDay(MyDate As Date)
Dim LunarDay As Integer

    If Month(MyDate) > 2 Then
    LunarDay = (((Year(MyDate) - 1900) Mod 19) * 11 + Month(MyDate) + Day(MyDate) + Int(Year(MyDate) / 300) + Int(Year(MyDate) / 400) + 6 - Int(Year(MyDate) / 100)) Mod 30
    Else
    LunarDay = (((Year(MyDate) - 1901) Mod 19) * 11 + Month(MyDate) + 12 + Day(MyDate) + Int(Year(MyDate) / 300) + Int(Year(MyDate) / 400) + 6 - Int(Year(MyDate) / 100)) Mod 30
    End If

End Function

' Однако всё изложенное справедливо лишь в первом приближении. По большому счёту, жуткая профанация.
' Из-за притяжения Солнца и планет, эллиптичности орбиты, даже светового давления - ошибка более полусуток (плюс-минус 14 часов с минутами!) всё равно остаётся.
' Реальная Луна движется настолько причудливо, что каждое следующее уточнение даётся ценой всё более сложной математики.
' Исключительно интересное занятие.

Меня раз упрекнули, мол, размещаешь программы на литературном сайте. Ну, во-первых, это таки "проза на других языках", в смысле прога, любой язык обладает избыточностью средств выражения, а владеют им не все; во-вторых, это просто красиво!

Если Вы дочитали досюда, Вам это и(ли) интересно, и(ли) нужно. Тогда продолжим.

Итак. Алгоритм Жана Меёуса работает с точностью до минуты, но выглядит устрашающе и занимает тысячу строк кода. Это точно не для здесь, спросите Алису. Внизу моей основной страницы есть ссылка: "Программа Zeitgeist", алгоритм Меёуса реализован в ней.

По пути на основе формул Меёуса мне удалось написать компактный, легко читаемый и неожиданно хорошо работающий фрагмент кода. Среднее отклонение на 16-летнем горизонте - 12 минут, это не сутки по формуле Харви и не 14 с лишним часов, как при прямом счёте.

На сайте НАСА есть мемориальный раздел, посвящённый Фреду Эспенаку, увлечённому астроному, умершему в 2025 году. Там есть таблицы новолуний на 6.000 лет, хорошая проверка для моего алгоритма.

Собственно, вот он; в нарушение астрономических традиций я считаю время Т не веками, а годами, что даёт две дополнительные значащие цифры к точности вычислений, а также использую необычные константы, имеющие физический смысл:

Option Explicit
Const SynodicMonth As Double = 29.5305888531

Function LastNewMoonUTC(ByVal YourDate As Date) As Double
Dim T As Double, N As Long
Dim M As Double, M_sun As Double, D As Double
Dim Res As Double
Dim Y As Double, DeltaT As Double, totalCorrDay As Double
Dim daysSinceRef As Double, CandiDate As Double
Const refDate As Date = "6.1.2000 14:15:00" '
Const PeriHelium As Date = "3.1.2000 12:00:00"
Const AnomalisticMonth As Double = 27.554551
Const TropicYear As Double = 365.259636
Const DegToRad As Double = 3.14159265358979 / 180#
   
    ' число дней между реперной и выбранной датами
    daysSinceRef = DateDiff("n", refDate, YourDate) / 1440#
    ' Число полных синодических месяцев от реперной даты до выбранной
    N = Int(daysSinceRef / SynodicMonth)
    ' кандидат на новолуние по "прямой сетке"
    CandiDate = N * SynodicMonth + refDate
    ' шкала времени, годы текущего столетия
    T = DateDiff("s", PeriHelium, CandiDate) / (TropicYear * 86400#)
 
   ' Аргументы и нормализация в [0,360°)
    M = Norm360(291.29 + 360# * T * TropicYear / AnomalisticMonth)
    M_sun = Norm360(359.53 + 360# * T)
    D = Norm360(272.01 + 360# * T * TropicYear / SynodicMonth)

    ' поправка на эллиптичность орбиты, градусов, основная гармоника
    Res = 6.28875 * Sin(DegToRad * (M + 64#))

    Res = Res + 1.499 * Sin(DegToRad * (2# * D - M))      'эвекция
    Res = Res - 0.21 * Sin(DegToRad * 2# * D)        'вариация
    Res = Res + 2.102 * Sin(DegToRad * M_sun)        'годовое уравнение
   
   ' высшие гармоники
    Res = Res - 0.197 * Sin(DegToRad * (2# * D - 2# * M))
    Res = Res + 0.139 * Sin(DegToRad * (2# * D + M))
    Res = Res + 0.116 * Sin(DegToRad * (D - M))

    Res = Res + 0.32 * Sin(DegToRad * (3# * D - M))
    Res = Res + 0.049 * Sin(DegToRad * (3# * D))

    ' Прецессия узлов орбиты
    Res = Res + 0.0045 * Sin(2# * 3.14159265358979 * Y / 18.6 + 0.7)

    ' Перевод отклонения из градусов в дни
    totalCorrDay = SynodicMonth * Res / 360

   ' Метод Espenak-Meeus:
    Y = (Year(CandiDate) + (Month(CandiDate) - 0.5) / 12#) - 2000#

    ' уход всемирного времени от астрономического
    DeltaT = 63.85 - 12.25 * Y + 1.05 * Y * Y
   
    LastNewMoonUTC = CandiDate + totalCorrDay - DeltaT / 86400#

Debug.Print "LastNewMoonUTC: "; Format(CDate(LastNewMoonUTC), "d:mm:yyyy h:nn:ss")
End Function

Private Function Norm360(ByVal x As Double) As Double
    x = x - 360# * Int(x / 360#)
    If x < 0 Then x = x + 360#
    Norm360 = x
End Function

' Теперь, зная время последнего новолуния, легко определить возраст Луны:

Private Sub MoonAge()
Dim Age As Long

'' Метод Espenak-Kayotkin:
'Как уже говорилось, лунные месяцы имеют переменную длительность. Поэтому:
'1 - Может статься, что фактическое новолуние уже настало, а Вы всё ещё считаете возраст Луны от предыдущего;
'2 - длительность "резиновых" лунных фаз не соответствует среднестатистической.
'
'Чтобы разом убить двух зайцев, я применяю "способ трёх лун": вычисляю не только прошедшее, но и два будущих новолуния. Если второе уже наступило, прошедшим считается оно, а следующим становится третье:
'
Dim LastNL, NextNL
    LastNL = LastNewMoonUTC(YourDate)
    NextNL = LastNewMoonUTC(YourDate + SynodicMonth)
    If Now > NextNL Then
    LastNL = NextNL
    NextNL = LastNewMoonUTC(YourDate + 2 * SynodicMonth)
    End If
   
    Age = DateDiff("h", LastNewMoonUTC(YourDate), YourDate) ' возраст Луны в часах на выбранный момент

Debug.Print Age, CDate(LastNL), CDate(NextNL)

End Sub

Остаётся временной интервал между прошлым LastNL и будущим NextNL разделить на число фаз и в нужные интервалы времени использовать нужные фазы, это тривиальная задача.

Например, 12 сентября 2026 года вышеприведённый код ошибся на 5 минут для прошлого и на 4 минуты для будущего новолуний. Такая точность вполне достаточна для отображения правильной фазы и возраста Луны. Более того: ещё шажок - и точность станет абсолютной!

Можно, конечно, просто создать базу данных для новолуний, но для выбора актуального всё равно понадобится похожий код, плюс сама база займёт мегабайты.

А можно пойти и дальше: сначала с помощью Debug.Print и упомянутых 6.000-летних таблиц один раз создать сводку минутных поправок (от 0 до 31), затем в вышеприведённой программе непосредственно перед "способом трёх лун" добавлять к результату в нужный месяц нужную поправку. Как ни смешно, итог - точное время новолуния и точный возраст Луны ценой крошечной, несколько килобайт, добавки к программному коду!

Успехов!