You may know that if you check the position of the Sun every day in the same place at the same time (accounting for daylight saving time if necessary), you’ll find that it slightly moves. This is a combination of the tilt of the Earth’s axis and the Earth’s orbital eccentricity. The path traced out by the position in the sky of the Sun during its wandering is called analemma.

Analemma Murray Hill

Afternoon analemma taken in 1998–99 by Jack Fishburn in Murray Hill, New Jersey, USA. Image credit: Jfishburn, Wikimedia Commons, GFDL 1.2+ and CC-BY-SA 3.0.

We can use Julia to plot the analemma. In particular, we’ll employ AstroLib.jl to do the needed calculations. Throughout this post I’ll assume you have installed the latest stable version of Julia and the necessary packages with the built-in package manager.

What we want to do is to determine the position of the Sun for a specific time every day in a year, say at noon for the whole 2018. This is the recipe:

  1. compute the Julian dates of all the wanted times
  2. calculate the equatorial coordinates for the given Julian dates
  3. convert the equatorial coordinates to horizontal coordinates in the desired place. For example, we choose Heidelberg, in Germany, which has coordinates 49°25′N 08°43′E and elevation of 114 m.

The trickiest part is to get the right Julian dates. The jdcnv function in AstroLib.jl assumes that times are given in UTC standard, but Heidelberg is one hour ahead of Greenwich. In order to work around this issue we can use the TimeZones.zdt2julian provided by the TimeZones.jl package which takes care of the time zones. In addition, Germany adopts daylight saving time from March to October, thus noon on May 15th is not actually the same time of day as noon on November 7th. However, noon on January 1st is the same time of day as noon on December 31st, so we can create a range between these two times with step one (Julian) day.

using AstroLib, TimeZones

function analemma(start_date, end_date,
                  latitude, longitude, elevation)
    julian_dates = TimeZones.zdt2julian(start_date):TimeZones.zdt2julian(end_date)
    right_ascension, declination = sunpos(julian_dates)
    altaz = eq2hor.(right_ascension, declination,
	                julian_dates, latitude, longitude, elevation)
    altitude = getindex.(altaz, 1)
    azimuth  = getindex.(altaz, 2)
    return azimuth, altitude
end

We have used sunpos to get the position of the Sun in equatorial coordinates and converted them with eq2hor to horizontal coordinates, specifying the coordinates of Heidelberg. The broadcast version of this function returns an array of 2-tuples, being the first element the altitude of the Sun and the second element its azimuth. We’ve used getindex.(altaz, i) to obtain the arrays with the i-th elements of the tuples. Now we can draw the analemma. I recommend using the Plots.jl package, which provides a single interface to several different back-ends (GR, PyPlot, PGFPlots, etc…).

using Plots, Base.Dates

azimuth, altitude =
    analemma(ZonedDateTime(2018,  1,  1, 12, tz"Europe/Berlin"),
             ZonedDateTime(2018, 12, 31, 12, tz"Europe/Berlin"),
             ten(49, 25), ten(8, 43), 114)

scatter(azimuth, altitude, aspect_ratio = :equal,
        xlabel = "Azimuth (°)", ylabel = "Altitude (°)")

Analemma Heidelberg

You can check with the JPL HORIZONS System that this is accurate within a few arcminutes.