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

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

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

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

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

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

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

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

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

Прямой же счёт точен, прост и очевиден: из общего числа дней надо убрать все полные лунные месяцы - вот она и лунная дата! Только не путать "третий день" и "три дня". А если кто заметит расхождение в 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:
' и получим расшатай на границе февраль-март. Прибавлять надо синхронизирующий коэффициент 13, а не 12.

    Resl = Resl + Month(MyDate)
    If Month(MyDate) < 3 Then Resl = Resl + 13
' Прибавить коэффициент столетия и день, разделить на 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:
' Легко видеть, что 01 марта 1900 года было новолуние.

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

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

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

' Даты новолуний в 2026 году (мск):
'18 января: 22:52 <= хорошая точка отсчёта, полночь посреди Евразии
'17 февраля: 15:01
'19 марта: 04:23
'17 апреля: 14:52
'16 мая: 23:01  <= ещё полночь в средней России
'15 июня: 05:54
'14 июля: 12:43
'12 августа: 20:37
'11 сентября: 06:27
'10 октября: 18:50
'9 ноября: 10:02    <= полдень над Уралом, тоже интересная точка
'9 декабря: 03:52

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

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

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

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

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

End Sub

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