By: Mosè Giordano

Re-posted from: https://giordano.github.io/blog/2017-11-12-analemma/

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.

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 (°)")

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