TinhAmLich.bas

Attribute VB_Name = "LunarSolar"
Option Explicit
Const PI As Double = 3.14159265358979 ' Atn(1) * 4

Nguồn: https://www.informatik.uni-leipzig.de/~duc/amlich/LunarSolar.bas

' Compute the (integral) Julian day number of day dd/mm/yyyy, i.e., the number
' of days between 1/1/4713 BC (Julian calendar) and dd/mm/yyyy.
' Formula from http://www.tondering.dk/claus/calendar.html
Function jdFromDate(ByVal dd As Long, ByVal mm As Long, ByVal yy As Long) As Long
    Dim a As Double, y As Long, m As Long, jd As Long
    a = Fix((14 - mm) / 12)
    y = yy + 4800 - a
    m = mm + 12 * a - 3
    jd = dd + Fix((153 * m + 2) / 5) + 365 * y _
        + Fix(y / 4) - Fix(y / 100) + Fix(y / 400) - 32045
    If jd < 2299161 Then
        jd = dd + Fix((153 * m + 2) / 5) + 365 * y + Fix(y / 4) - 32083
    End If
    jdFromDate = jd
End Function

' Convert a Julian day number to day/month/year. Parameter jd is an integer
Function jdToDate(ByVal jd As Long)
    Dim a As Long, b As Long, c As Long, d As Long, e As Long, m As Long
    Dim Day As Long, Month As Long, Year As Long
    If (jd > 2299160) Then ' After 5/10/1582, Gregorian calendar
        a = jd + 32044
        b = Fix((4 * a + 3) / 146097)
        c = a - Fix((b * 146097) / 4)
    Else
        b = 0
        c = jd + 32082
    End If
    d = Fix((4 * c + 3) / 1461)
    e = c - Fix((1461 * d) / 4)
    m = Fix((5 * e + 2) / 153)
    Day = e - Fix((153 * m + 2) / 5) + 1
    Month = m + 3 - 12 * Fix(m / 10)
    Year = b * 100 + d - 4800 + Fix(m / 10)
    jdToDate = Array(Day, Month, Year)
End Function

' Compute the time of the k-th new moon after the new moon of 1/1/1900 13:52 UCT
' (measured as the number of days since 1/1/4713 BC noon UCT,
' e.g., 2451545.125 is 1/1/2000 15:00 UTC).
' Returns a floating number, e.g.,
' 2415079.9758617813 for k=2 or 2414961.935157746 for k=-2

Function NewMoon(ByVal k As Long) As Double
    Dim T As Double, T2 As Double, T3 As Double, dr As Double
    Dim Jd1 As Double, m As Double, Mpr As Double
    Dim F As Double, C1 As Double, deltat As Double, JdNew As Double
    T = k / 1236.85 ' Time in Julian centuries from 1900 January 0.5
    T2 = T * T
    T3 = T2 * T
    dr = PI / 180
    Jd1 = 2415020.75933 + 29.53058868 * k + 0.0001178 * T2 - 0.000000155 * T3
    Jd1 = Jd1 + 0.00033 * Sin((166.56 + 132.87 * T - 0.009173 * T2) * dr)
        ' Mean new moon
    m = 359.2242 + 29.10535608 * k - 0.0000333 * T2 - 0.00000347 * T3
        ' Sun's mean anomaly
    Mpr = 306.0253 + 385.81691806 * k + 0.0107306 * T2 + 0.00001236 * T3
        ' Moon's mean anomaly
    F = 21.2964 + 390.67050646 * k - 0.0016528 * T2 - 0.00000239 * T3
        ' Moon's argument of latitude
    C1 = (0.1734 - 0.000393 * T) * Sin(m * dr) + 0.0021 * Sin(2 * dr * m)
    C1 = C1 - 0.4068 * Sin(Mpr * dr) + 0.0161 * Sin(dr * 2 * Mpr)
    C1 = C1 - 0.0004 * Sin(dr * 3 * Mpr)
    C1 = C1 + 0.0104 * Sin(dr * 2 * F) - 0.0051 * Sin(dr * (m + Mpr))
    C1 = C1 - 0.0074 * Sin(dr * (m - Mpr)) + 0.0004 * Sin(dr * (2 * F + m))
    C1 = C1 - 0.0004 * Sin(dr * (2 * F - m)) - 0.0006 * Sin(dr * (2 * F + Mpr))
    C1 = C1 + 0.001 * Sin(dr * (2 * F - Mpr)) + 0.0005 * Sin(dr * (2 * Mpr + m))
    If (T < -11) Then
        deltat = 0.001 + 0.000839 * T + 0.0002261 * T2 _
                - 0.00000845 * T3 - 0.000000081 * T * T3
    Else
        deltat = -0.000278 + 0.000265 * T + 0.000262 * T2
    End If
    JdNew = Jd1 + C1 - deltat
    NewMoon = JdNew
End Function

' Compute the longitude of the sun at any time.
' Parameter: floating number jdn, the number of days since 1/1/4713 BC noon

Function SunLongitude(ByVal jdn As Double) As Double
    Dim T As Double, T2 As Double, dr As Double, m As Double
    Dim L0 As Double, DL As Double, L As Double
    T = (jdn - 2451545) / 36525
        ' Time in Julian centuries from 2000-01-01 12:00:00 GMT
    T2 = T * T
    dr = PI / 180 ' degree to radian
    m = 357.5291 + 35999.0503 * T - 0.0001559 * T2 - 0.00000048 * T * T2
        ' mean anomaly, degree
    L0 = 280.46645 + 36000.76983 * T + 0.0003032 * T2
        ' mean longitude, degree
    DL = (1.9146 - 0.004817 * T - 0.000014 * T2) * Sin(dr * m)
    DL = DL + (0.019993 - 0.000101 * T) * Sin(dr * 2 * m) _
        + 0.00029 * Sin(dr * 3 * m)
    L = L0 + DL ' true longitude, degree
    L = L * dr
    L = L - PI * 2 * (Fix(L / (PI * 2))) ' Normalize to (0, 2*PI)
    SunLongitude = L
End Function

' Compute sun position at midnight of the day with the given Julian day number.
' The time zone if the time difference between local time and UTC: 7.0 for UTC+7:00.
' The function returns a number between 0 and 11.
' From the day after March equinox and the 1st major term after March equinox,
' 0 is returned. After that, return 1, 2, 3 ...
Function getSunLongitude(ByVal dayNumber As Double, ByVal timeZone As Byte) As Long
    getSunLongitude = Fix(SunLongitude(dayNumber - 0.5 - timeZone / 24) / PI * 6)
End Function

' Compute the day of the k-th new moon in the given time zone.
' The time zone if the time difference between local time and UTC: 7.0 for UTC+7:00
Function getNewMoonDay(ByVal k As Long, ByVal timeZone As Long) As Long
    getNewMoonDay = Fix(NewMoon(k) + 0.5 + timeZone / 24)
End Function

' Find the day that starts the luner month 11 of the given year
' for the given time zone
Function getLunarMonth11(ByVal yy As Long, ByVal timeZone As Long) As Long
    Dim k As Long, off As Double, nm As Long, sunLong As Double
    '' off = jdFromDate(31, 12, yy) - 2415021.076998695
    off = jdFromDate(31, 12, yy) - 2415021
    k = Fix(off / 29.530588853)
    nm = getNewMoonDay(k, timeZone)
    sunLong = getSunLongitude(nm, timeZone) ' sun longitude at local midnight
    If (sunLong >= 9) Then
        nm = getNewMoonDay(k - 1, timeZone)
    End If
    getLunarMonth11 = nm
End Function

' Find the index of the leap month after the month starting on the day a11.
Function getLeapMonthOffset(ByVal a11 As Double, ByVal timeZone As Long) As Long
    Dim k As Long, last As Long, Arc As Long, I As Long
    k = Fix((a11 - 2415021.07699869) / 29.530588853 + 0.5)
    last = 0
    I = 1 ' We start with the month following lunar month 11
    Arc = getSunLongitude(getNewMoonDay(k + I, timeZone), timeZone)
    Do
        last = Arc
        I = I + 1
        Arc = getSunLongitude(getNewMoonDay(k + I, timeZone), timeZone)
    Loop While (Arc <> last And I < 14)
    getLeapMonthOffset = I - 1
End Function

' Comvert solar date dd/mm/yyyy to the corresponding lunar date
Function Solar2Lunar( _
        ByVal dd As Long, _
        ByVal mm As Long, _
        Optional ByVal yy As Long = 0, _
        Optional ByVal timeZone As Long = 7) As String
    Dim k As Long, diff As Long, leapMonthDiff As Long, dayNumber As Long
    Dim monthStart As Double, a11 As Long, b11 As Long
    Dim lunarDay As Double, lunarMonth As Long, lunarYear As Long, lunarLeap As Long
    '~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    If yy = 0 Then yy = Year(Date)
    dayNumber = jdFromDate(dd, mm, yy)
    k = Fix((dayNumber - 2415021.07699869) / 29.530588853)
    monthStart = getNewMoonDay(k + 1, timeZone)
    If (monthStart > dayNumber) Then
        monthStart = getNewMoonDay(k, timeZone)
    End If
    ' alert(dayNumber + " -> " + monthStart)
    a11 = getLunarMonth11(yy, timeZone)
    b11 = a11
    If (a11 >= monthStart) Then
        lunarYear = yy
        a11 = getLunarMonth11(yy - 1, timeZone)
    Else
        lunarYear = yy + 1
        b11 = getLunarMonth11(yy + 1, timeZone)
    End If
    lunarDay = dayNumber - monthStart + 1
    diff = Fix((monthStart - a11) / 29)
    lunarLeap = 0
    lunarMonth = diff + 11
    If (b11 - a11 > 365) Then
        leapMonthDiff = getLeapMonthOffset(a11, timeZone)
        If (diff >= leapMonthDiff) Then
            lunarMonth = diff + 10
            If (diff = leapMonthDiff) Then lunarLeap = 1
        End If
    End If
    If (lunarMonth > 12) Then lunarMonth = lunarMonth - 12
    If (lunarMonth >= 11 And diff < 4) Then lunarYear = lunarYear - 1
    Solar2Lunar = Format(lunarDay, "00") & _
                "/" & Format(lunarMonth, "00") & _
                "/" & Format(lunarYear, "0000 \A\L") & IIf(lunarLeap, " (" & lunarMonth & " N)", "")
End Function

' Convert a lunar date to the corresponding solar date
Function Lunar2Solar( _
        ByVal lunarDay As Long, _
        ByVal lunarMonth As Long, _
        Optional ByVal lunarYear As Long = 0, _
        Optional ByVal lunarLeap As Long = 0, _
        Optional ByVal timeZone As Long = 7) As Date
    Dim k As Long, a11 As Long, b11 As Long, off As Long, leapOff As Long
    Dim LeapMonth As Long, monthStart As Long
    '~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    If lunarYear = 0 Then lunarYear = Year(Date)
    If (lunarMonth < 11) Then
        a11 = getLunarMonth11(lunarYear - 1, timeZone)
        b11 = getLunarMonth11(lunarYear, timeZone)
    Else
        a11 = getLunarMonth11(lunarYear, timeZone)
        b11 = getLunarMonth11(lunarYear + 1, timeZone)
    End If
    k = Fix(0.5 + (a11 - 2415021.07699869) / 29.530588853)
    off = lunarMonth - 11
    If (off < 0) Then off = off + 12
    If (b11 - a11 > 365) Then
        leapOff = getLeapMonthOffset(a11, timeZone)
        LeapMonth = leapOff - 2
        If (LeapMonth < 0) Then LeapMonth = LeapMonth + 12
        If (lunarLeap <> 0 And lunarMonth <> LeapMonth) Then
            Lunar2Solar = Array(0, 0, 0)
            Exit Function
        ElseIf (lunarLeap <> 0 Or off >= leapOff) Then
            off = off + 1
        End If
    End If
    monthStart = getNewMoonDay(k + off, timeZone)
    Dim R
    R = jdToDate(monthStart + lunarDay - 1)
    Lunar2Solar = DateSerial(R(2), R(1), R(0))
End Function
                                    
Function AmLich(d as Date)
    AmLich = Solar2Lunar(day(d),month(d),year(d))
end Function