Catching the moon at the right time with Julia

The moon on a night sky is surprisingly bright, especially for a camera with limited dynamic range. And there’s only so many compositions one can do of the moon on a naked sky. For those reasons taking moon pictures at sunrise or sunset, when the moon is close to the horizon and its relative brightness is comparable to ambient light often brings more interesting (or at least varied) results.

Maybe I haven’t searched well enough but I haven’t found a tool giving you the days of the year at which moon rise/set coincide with the sun rise/set for a given location. In any case that’s an excuse for a fun Julia project.

Disclaimer: I’m not a physicist or professional astronomer, if you’re planning something critical (like a rocket launch), maybe double-check with actual experts!

The Plan

The idea is simple: find dates when the sun and moon rise or set within a short time window of each other. This gives us those perfect conditions where you get a nicely lit sky and a visible moon. The main functionalities, namely the moon and sun position at a given date, are available via the AstroLib package.

Setup and converting times

First, let’s define where we’re observing from. I’m using Lyon, France as the example location:

using AstroLib, TimeZones
using Dates, Plots

# Location setup: latitude, longitude in degrees
latitude = 45.7640     # Lyon
longitude = 4.8357
obs_altitude = 212
local_tz = TimeZone("Europe/Paris")

start_date = Date(2025, 8, 4)

We also need our timezone to convert between universal time and our local observations.

AstroLib routines take time in Julian days, which surprisingly aren’t Julia specific units, but units used by astronomers. We can convert from a Julia date to Julian days with jdcnv :

date_to_jd(date) = jdcnv(DateTime(date))

Computing Altitudes

The get_altitude function is doing the real work here. It takes a Julian date and a position function (sunpos or moonpos), gets the celestial coordinates, and converts them to altitude above our horizon. When altitude crosses zero, we have a rise or set event.

function get_altitude(jd, posfun=sunpos)
    sun_ra_dec = posfun(jd)
    sun_alt, sun_az, _ = eq2hor(sun_ra_dec[1], sun_ra_dec[2], jd, latitude, longitude, obs_altitude)
    sun_alt
end

The sunpos function returns the right ascension and declination of the Sun at a given date, right ascension and declination being angles defining a direction on the celestial sphere.

The eq2hor function converts these two angles into an altitude (angle to the horizon) for our local position.

We can make a plot of the altitude for a day. As cou can see on that day the mon rises as the sun is setting, that’s the type of events we are loooking for!

start_date = Date(2025, 8, 8)

jd_start = date_to_jd(start_date)
jd_end   = date_to_jd(start_date + Day(1))

t = LinRange(jd_start, jd_end, 100)

p = plot(t, get_altitude.(t, sunpos), label = "Sun altitude", xlabel = "Julian time")
plot!(t, get_altitude.(t, moonpos), label = "Moon altitude", ylabel = "Altitude")
hline!([0], label = "Horizon",c="gray", dpi=150)

Finding the Crossings

To find the horizon crossings I’m using a basic zero-crossing algorithm, with a 1 minute time resolution. I’ve tried root finding with Roots.jl but it surprisingly was getting stuck sometimes. I don’t need more precision anyway.

When altitude goes from negative to positive, that’s a rise. Positive to negative? That’s a set. Simple but effective.

function find_crossing(start_date, num_days, posfun)
    jd_start = date_to_jd(start_date)
    jd_end = date_to_jd(start_date + Day(num_days))

    times = LinRange(jd_start, jd_end, 24*60*num_days) # 1 minute resolution
    altitudes = get_altitude.(times, posfun)

    crossings =  @NamedTuple{time::DateTime, event::Symbol}[]

    for i in 2:length(altitudes)
        if altitudes[i-1] < 0 && altitudes[i] ≥ 0
            push!(crossings, (
                time = daycnv(times[i]),
                event = :rise,
            ))
        elseif altitudes[i-1] ≥ 0 && altitudes[i] < 0
            push!(crossings, (
                time = daycnv(times[i]),
                event = :set,
            ))
        end
    end
    crossings
end

The same plot as before with the crossings added :

Moon & Sun Coincidences

No need for sophistication either here, we’ll just test every possible pairs of crossing for a year; if they occur within 30 minutes of each other, we keep them.

# Compute crossings for a full year
d = 360
sun_crossings  = find_crossing(start_date, d, sunpos)
moon_crossings = find_crossing(start_date, d, moonpos)

# Find coincidences within a 30-minute window
coincidences = Tuple{T,T}[]

for s in sun_crossings
    for m in moon_crossings
        if abs(s.time - m.time) < Minute(30)
            push!(coincidences, (s,m))
        end
    end
end

Pretty Printing the Results

Finally, let’s format the results nicely, including the moon phase (illuminated fraction) and local times:

We can double check that a full moon will indeed around 17h the 3 Jan, 2026 in Lyon :

https://www.timeanddate.com/moon/france/lyon?month=1&year=2026

Code

You can generate similar tables for your location and timezone. It would be cool to provide an interactive page doing that.

https://github.com/jonathanBieler/sun_and_moon