Fun with Python: Phase of the Moon Emoji

Hey everyone,

I went grocery shopping yesterday, which for me is usually a soul-crushing experience. This time, however, was a pleasant surprise. For the first time in months, my wife and I left the store without feeling defeated by the cost. We got back, I put things away, brewed a fresh pot of coffee, and had a feeling it was going to be a good day.

Then, slowly at first, and then all at once, I felt ill. After just two sips of coffee, I was back in bed, shivering and nauseous. I ended up sleeping for most of the day, which meant I was wide awake in the middle of the night. While catching up on emails, I saw a link in PyCoder's Weekly to a neat little script that generated ASCII art of the moon's phases. I loved the idea but felt the ASCII art was a bit much. So, I decided it would be a perfect "Fun with Python" project to make my own version using emojis. There is an emoji for each of the 8 phases of the moon!

On the surface, it seems simple, but it led me down a fascinating rabbit hole I wanted to explore: calculating the Julian Day.

The Julian Day

For astronomical calculations like tracking the moon's phase, using a standard calendar is a nightmare. To solve this, astronomers use the Julian Day, which is simply a continuous count of days that have passed since noon on January 1, 4713 BC. This turns complex date arithmetic into simple subtraction.

I could have easily used a library like ephem to do this, but I wanted to understand the math myself. So, I decided to write my own function to perform the calculation.

def _julian_day(dt_utc: datetime) -> float:
    """Convert a UTC datetime to Julian Day."""
    y, m = dt_utc.year, dt_utc.month
    d = (
        dt_utc.day
        + (
            dt_utc.hour
            + (dt_utc.minute + (dt_utc.second + dt_utc.microsecond / 1_000_000) / 60)
            / 60
        )
        / 24
    )
    if m <= 2:
        y -= 1
        m += 12
    A = y // 100
    if dt_utc >= datetime(1582, 10, 15, tzinfo=timezone.utc):
        B = 2 - A + A // 4
    else:
        B = 0
    return int(365.25 * (y + 4716)) + int(30.6001 * (m + 1)) + d + B - 1524.5

What the Ever-Living Hell Is Going On Here?

That function looks like a mess of magic numbers, but each part serves a specific purpose in a well-established algorithm. Let's break it down.

d = dt_utc.day + ...: This part is straightforward. It takes the day of the month and adds the time (hours, minutes, seconds) as a fractional part of a day.
if m <= 2: y -= 1; m += 12: This is a classic programmer's trick. It treats January and February as months 13 and 14 of the previous year. This simplifies the math for leap years, as the leap day (Feb 29th) now always falls at the end of this adjusted calendar year.
int(365.25 * (y + 4716)): This is the core of the day-counting. The formula's internal starting point is March 1, 4716 BC. The y + 4716 calculates the number of years from that point. Multiplying by 365.25 accounts for the average length of a year, including leap years. We start with this day because of the Jan/Feb shift and it's the first leap year.
int(30.6001 * (m + 1)): This calculates the number of days that have passed since March 1st of your target year. The constant 30.6001 is a clever mathematical trick that approximates the pattern of 30- and 31-day months in our shifted calendar.
B = 2 - A + A // 4: This is the Gregorian calendar correction. The Julian calendar, which was used before October 15, 1582, had a slightly different leap year rule. This part of the formula applies a correction for all modern dates to account for that historical change. (years divisible by 100 are not leap years unless they are also divisible by 400).
- 1524.5: This is the final and most important correction. It's the precise number of days between the formula's internal starting point (March 1, 4716 BC) and the official Julian Day starting point (noon on January 1, 4713 BC). By subtracting this value, the formula shifts its result to align perfectly with the standard astronomical convention. The .5 is there because Julian Days start at noon, not midnight.

By combining these steps, this function can accurately convert any date into its corresponding Julian Day number, which we can then use to calculate the moon's phase with simple math.

Trust me, it would have been easier to use ephem.

import ephem

date = ephem.Date("2025/10/15")
julian_day = ephem.julian_date(date)
print("Julian Day:", julian_day)

Full Code

The complete script, including the part that converts the final phase into a cool emoji, is available on Gist.

Gist: https://gist.github.com/TheCrazyGM/365b3d208d099feefbc3de8d79367736

As always,
Michael Garcia a.k.a. TheCrazyGM