Moon Phases

//@version=6
indicator("Moon Phases", "MOON", overlay = true)

waxingColor = input.color(color.blue, "Waxing Moon")
waningColor = input.color(color.white, "Waning Moon")

if timeframe.ismonthly
    runtime.error("Unsupported timeframe.")

dayOfYear(barTime) =>
    yearValue = year(barTime)
    leapYear = yearValue % 400 == 0 or (yearValue % 4 == 0 and yearValue % 100 != 0)
    dayCount = array.from(31, leapYear ? 29 : 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31)
    monthValue = month(barTime)
    days = 0
    if monthValue > 1
        for i = 0 to monthValue - 2
            days += dayCount.get(i)
    days + dayofmonth(barTime)

unixTimeFromJulianDay(julianDay) =>
    z = math.floor(julianDay + 0.5)
    f = (julianDay + 0.5) % 1
    alpha = math.floor((z - 1867216.25) / 36524.25)
    a = z < 2299161 ? z : z + 1 + alpha - math.floor(alpha / 4.0)
    b = a + 1524
    c = math.floor((b - 122.1) / 365.25)
    d = math.floor(365.25 * c)
    e = math.floor((b - d) / 30.6001)
    dayFloat = b - d - math.floor(30.6001 * e) + f
    dayValue = int(dayFloat)
    monthValue = e < 13.5 ? int(e - 1) : int(e - 13)
    yearValue = monthValue > 2.5 ? int(c - 4716) : int(c - 4715)
    secondTotal = int(dayFloat % 1 * 60 * 60 * 24)
    hourValue = math.floor(secondTotal / 3600)
    minuteValue = math.floor((secondTotal - hourValue * 3600) / 60)
    secondValue = secondTotal % 60
    timestamp(yearValue, monthValue, dayValue, hourValue, minuteValue, secondValue)

lastMoonPhase(barTime) =>
    k = (dayOfYear(barTime) / 364.25 + year(barTime) - 1900) * 12.3685
    int moonPhase = na
    if k % 1.0 < 0.5
        k := math.floor(k)
        moonPhase := 1
    else
        k := math.floor(k) + 0.5
        moonPhase := -1
    t = k / 1236.85
    anomalySun = math.toradians(359.2242 + (29.10535608 * k) - (0.0000333 * math.pow(t, 2)) - (0.00000347 * math.pow(t, 3)))
    anomalyMoon = math.toradians(306.0253 + (385.81691806 * k) + (0.0107306 * math.pow(t, 2)) + (0.00001236 * math.pow(t, 3)))
    f = math.toradians(21.2964 + ((390.67050646 * k) - (0.0016528 * math.pow(t, 2)) - (0.00000239 * math.pow(t, 3))))
    dev = (0.1734 - (0.000393 * t)) * math.sin(anomalySun)
          + 0.0021 * math.sin(2 * anomalySun)
          - 0.4068 * math.sin(anomalyMoon)
          + 0.0161 * math.sin(2 * anomalyMoon)
          - 0.0004 * math.sin(3 * anomalyMoon)
          + 0.0104 * math.sin(2 * f)
          - 0.0051 * math.sin(anomalySun + anomalyMoon)
          - 0.0074 * math.sin(anomalySun - anomalyMoon)
          + 0.0004 * math.sin(2 * f + anomalySun)
          - 0.0004 * math.sin(2 * f - anomalySun)
          - 0.0006 * math.sin(2 * f + anomalyMoon)
          + 0.0010 * math.sin(2 * f - anomalyMoon)
          + 0.0005 * math.sin(anomalySun + 2 * anomalyMoon)
    julianDate = 2415020.75933 + (29.53058868 * k) + (0.0001178 * math.pow(t, 2)) - (0.000000155 * math.pow(t, 3)) +
          (0.00033 * math.sin(math.toradians(166.56) + (math.toradians(132.87) * t) - (math.toradians(0.009173) * math.pow(t, 2)))) + dev
    timeLastMoonPhase = unixTimeFromJulianDay(julianDate)
    timeLastMoonPhase := timeLastMoonPhase >= barTime ? timeLastMoonPhase[1] : timeLastMoonPhase
    phaseTimeChanged = ta.change(timeLastMoonPhase) != 0
    var int moonPhaseAdjusted = moonPhase
    moonPhaseAdjusted := barstate.isfirst or phaseTimeChanged ? moonPhase : moonPhaseAdjusted[1]
    moonType = ta.change(moonPhaseAdjusted) != 0 ? moonPhase : na
    [moonPhaseAdjusted, moonType]

[moonPhase, moonType] = lastMoonPhase(time_close)

plotshape(moonType == 1, "New Moon", shape.circle, location.abovebar, color.new(waxingColor, 50), size = size.normal, display = display.pane)
plotshape(moonType == -1, "Full Moon", shape.circle, location.belowbar, color.new(waningColor, 50), size = size.normal, display = display.pane)
bgcolor(color.new(moonPhase == 1 ? waxingColor : waningColor, 95))