9
$\begingroup$

I have in mind to do a yearlong photo project of the solar analemma.

Can someone suggest a software program that will provide (& chart) the altitude and azimuth of the sun, for a given location, at a given time of day throughout the year?

I’d like a graph, and perhaps a table of values.

For example, for Dennis di Cicco’s famous analemma photograph, taken from Boston, at 8:30 AM (throughout the course of a year), the parameters would be

Boston (42.3601° N, 71.0589° W), and 8:30 AM (13:30 UT).

I would like to be able to use a different location and time of day than the di Cicco photograph.

(I have Stellarium on my iphone, if there’s a way to use it).

Analemma, Dennis di Cicco

$\endgroup$
1
  • 2
    $\begingroup$ Another project that's similar that might interest you is a solarigraph. It's a continuous 6 month exposure. It has a higher success rate as one cloudy day won't ruin a year's worth of work. $\endgroup$ Jun 27, 2023 at 15:40

3 Answers 3

12
$\begingroup$

If you're open to using a bit of Python, you can use the following script to generate an analemma for arbitrary times and locations.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.projections
import astropy.units as u
import astropy.time
import astropy.coordinates
import astropy.visualization
import cartopy

time = astropy.time.Time("2023-06-26T13:30") + np.arange(365) * u.day
location = astropy.coordinates.EarthLocation.of_address("Boston")
frame = astropy.coordinates.AltAz(obstime=time, location=location)

sun = astropy.coordinates.get_body("sun", time=time).transform_to(frame)

alt = sun.alt.deg
az = sun.az.deg

projection = cartopy.crs.Orthographic(az.mean(), alt.mean())
fig, ax = plt.subplots(subplot_kw=dict(projection=projection))
ax.plot(az, alt, transform=cartopy.crs.PlateCarree());
ax.gridlines(draw_labels=True)
plt.show()

which outputs:

analemma for Boston

$\endgroup$
9
  • $\begingroup$ Nice. Thanks. Can you provide a tweak to your code for a projection where the meridian lines converge to the pole, and the shape of the analemma doesn’t “droop” as much relative to a more symmetrical figure-8? This might make the graph look more like the di Cicco photograph. Perhaps polar coordinates, or maybe there’s a projection that models the field of view of a camera? (It’s been quite a while since I coded graphs in python.) $\endgroup$ Jun 27, 2023 at 13:54
  • 1
    $\begingroup$ @BruceSimonson, I've updated the answer with an improved projection. Let me know if that's closer to what you're looking for. $\endgroup$
    – Roy Smart
    Jun 27, 2023 at 15:58
  • 1
    $\begingroup$ What are those horizontal coordinates (110deg E, 120deg E) on the output diagram? They seem wrong for the time. $\endgroup$ Jan 14 at 4:06
  • 1
    $\begingroup$ @JacobMiller, the horizontal coordinate is the local azimuthal coordinate in Boston, MA. Why do they seem wrong? Stellarium produces roughly the same result for 13:30 UTC using MikeG's answer. $\endgroup$
    – Roy Smart
    Jan 14 at 18:08
  • 1
    $\begingroup$ Ah yes, it is a continual annoyance to me that Astropy does not support timezones, but oh well. I meant to add the plt.show(), but I must've forgotten. In Jupyter notebooks the plt.show() is implicit which is why I forgot. $\endgroup$
    – Roy Smart
    Jan 14 at 23:13
9
$\begingroup$

In the desktop version of Stellarium, you can do this on the Ephemeris tab of the Astronomical calculations window. Set the beginning of the date range with the desired time of day, and make the time step a whole number of solar days. Check the box H.C. (horizontal coordinates), otherwise you get a ring around the ecliptic. You can write a CSV file using the Export ephemeris... dialog.

Analemma generated by Stellarium

$\endgroup$
3
$\begingroup$

Here's a Sage / Python program which creates analemma diagrams using data from JPL Horizons. The azimuth + elevation data is displayed in 3D on a globe, which can be rotated & zoomed. (Sage uses three.js for interactive 3D graphics).

Please note that this globe displays the azimuth and elevation, so its poles are the zenith and nadir, the horizontal circles are lines of equal elevation (corresponding to parallels of latitude on a normal globe), and the vertical semicircles are lines of equal azimuth (corresponding to meridians of longitude on a normal globe).

You can supply the location using coordinates in longitude, latitude, altitude form, or you can supply the site ID using the MPC code of an observatory. The altitude, given in kilometres, is relative to the WGS-84 reference ellipsoid (as used by GPS), so it's very close to altitude relative to mean sea level.

If you supply a site ID the data in the lon, lat, alt field is ignored. To specify lon, lat, alt, put coord in the site field, or leave it blank. If you put * in the site field, Horizons will list all the MPC sites. You can also fetch the list with this URL.

To create a correct analemma you need to specify the time zone as part of the start time field. Eg, to make a noon analemma for longitude 30°E, the time should be 12:00 UT+2 or 12:00 UT+2:00.

Horizons can estimate the effects of atmospheric refraction, if desired. (The true refraction depends on the weather). The effects of refraction are fairly small around noon, but are significant near dawn and dusk.

Here's the relevant entry from the Horizons manual:

4. Apparent azimuth & elevation (AZ-EL) :

Apparent azimuth and elevation of target. Adjusted for light-time, the gravitational deflection of light, stellar aberration, precession and nutation. There is an optional (approximate) adjustment for atmospheric refraction (Earth only).

Azimuth is measured clockwise from north: North(0) -> East(90) -> South(180) -> West(270) Elevation angle is with respect to plane perpendicular to local zenith direction.
TOPOCENTRIC ONLY. Units: DEGREES

Here's a noon analemma for Greenwich observatory, using a 7 day timestep, and a datestep of 4 (so every 4th dot is larger). This plot is for 2023. It doesn't use refraction.

Analemma globe

The program can display numeric labels near the datestep points, but they tend to clutter the diagram, so I toggled the labels off for that example.

Here's the code.

""" Retrieve solar azimuth & elevation from Horizons for a specified Earth location
    and plot it in 3D on a globe.

    Written by PM 2Ring 2024.01.13
"""

import re, requests
from functools import lru_cache

url = "https://ssd.jpl.nasa.gov/api/horizons_file.api"
api_version = "1.0"

# Greenwich Lon, Lat, Alt (km)
greenwich = '0.0, 51.4773207, 0.06707'

base_cmd = """
MAKE_EPHEM=YES
EPHEM_TYPE=OBSERVER
QUANTITIES=4
CAL_FORMAT=CAL
CAL_TYPE=GREGORIAN
CSV_FORMAT=YES
OBJ_DATA=NO
COMMAND=10
"""

@lru_cache(maxsize=int(3))
def fetch_data(site, coord, start, stop, step, refracted, verbose=True):       
    if site == "":
        site = "coord"
    apparent = "REFRACTED" if refracted else "AIRLESS"
    
    cmd = f"""
!$$SOF
{base_cmd}
CENTER='{site}@399'
SITE_COORD='{coord}'
START_TIME='{start}'
STOP_TIME='{stop}'
STEP_SIZE='{step}'
APPARENT='{apparent}'
!$$EOF
"""
    #print(cmd)
    req = requests.post(url, data={'format': 'text'}, files={'input': ('cmd', cmd)})
    version = re.search(r"API VERSION:\s*(\S*)", req.text).group(1)
    if version != api_version:
        print(f"Warning: API version is {version}, but this script is for {api_version}")

    m = re.search(r"(?s)\\\$\\\$SOE(.*)\\\$\\\$EOE", req.text)
    if m is None:
        print("NO EPHEMERIS")
        print(req.text)
        return None

    if verbose:
        print(req.text)
    else:
        lines = req.text.splitlines()
        print("\n".join(lines[5:19]), "\n")

    return m.group(1)[1:]

def extract_data(data, datestep):
    data = data.splitlines()
    data = [s.split(',')[:-1] for s in data]
    #for row in data: print(row)
    #return

    pos, dates = [], []
    for i, row in enumerate(data):
        # Extract azimuth & elevation
        az, el = float(row[3]), float(row[4])
        pos.append((az, el))

        # Make timestamps
        if datestep and i % datestep == 0:
            # Extract calendar date
            ymd, hms = row[0].split()
            dates.append(ymd)

    return pos, dates

# Globe with lat-lon grid
var('u,v')
pp = parametric_plot3d
deg = pi / 180
dth = 15 * deg

def xyz(lat, lon):
    r = cos(lat)
    return vector([r*cos(lon), r*sin(lon), sin(lat)])

udom, vdom = (u, -pi/2, pi/2), (v, -pi, pi)
G = pp(xyz(u, v), udom, vdom, color='#007', plot_points=75)
G += sum(pp(xyz(u, lon), udom, color='#00c' if lon != pi else '#0ac') for lon in srange(0, pi*2, dth))
G += sum(pp(xyz(lat, v), vdom, color='#00c' if lat else '#0ac') for lat in srange(dth-pi/2, pi/2, dth))

@interact
def main(txt=HtmlBox("<h3>Analemma globe</h3>"),
  site='', coord=('lon, lat, alt (km)', greenwich),
  start="2023-Jan-01 12:00", stop="2024-Jan-01", step="7d",
  datestep=4, show_labels=True, refracted=False, color="red",
  verbose=False, auto_update=False):

    result = fetch_data(site.strip(), coord.strip(),
      start.strip(), stop.strip(), step.strip(), refracted, verbose=verbose)
    if result is None:
        print("No ephemeris data found!")
        return

    pos, dates = extract_data(result, datestep)
    #for row in pos: print(row)

    for i, s in enumerate(dates):
        print(i, s)

    # Make a blank globe
    P = G

    for i, (az, el) in enumerate(pos):
        do_label = datestep and i % datestep == 0
        v = xyz(el * deg, -az * deg)
        P += point3d(v, size=5 if do_label else 3, color=color)
        if show_labels and do_label:
            label = str(i // datestep)
            vv = v * 1.1
            P += text3d(label, vv, fontsize=8)

    P.show(frame=False, theme='dark', projection='orthographic', online=True)

The data-fetching part of the code is plain Python. The Sage features are only needed to do the 3D display.

Here's a live version of the program, running on the SageMathCell server.

The program maintains a cache of its last 3 queries of Horizons data, so if you only change "cosmetic" options the plot is generated without re-fetching the data from Horizons.

Here are the 3D interface controls.

  • Orbit - right mouse, or left mouse + ctrl/meta/shiftKey
    • touch: one-finger rotate
  • Zoom - middle mouse, or mousewheel
    • touch: two-finger spread or squish
  • Pan - left mouse, or arrow keys
    • touch: two-finger move
$\endgroup$
7
  • $\begingroup$ Awesome. It can be called, the AIO Earth Analemma generator :). But can you Label the coordinate system and directions in the diagram? Or only numbers allowed? Because it would be easy for layman to understand the diagram. $\endgroup$ Jan 13 at 16:41
  • $\begingroup$ @Jacob I suppose I could add coordinates, but they would get in the way, especially for noon analemmas. The grid steps are 15°. The "equator" is 0° elevation. The prime "meridian" is the azimuth=180° (south) line. They're in a lighter, greener shade than the other azimuth & elevation circles, but I guess that may be not so easy to see. $\endgroup$
    – PM 2Ring
    Jan 13 at 16:53
  • $\begingroup$ So this view is outside the celestial sphere. Then if you made the globe transparent and looked from the opposite side (in the direction of the north at noon), view correspond to how observers on Earth see? $\endgroup$ Jan 13 at 17:11
  • 1
    $\begingroup$ @Jacob Right, but you can't actually zoom the camera into the sphere. I guess it's possible by playing with the viewpoint parameter... But it's easy to make the globe transparent: just add something like opacity=0.5 to the function call that creates the globe: G = pp(xyz(u, v), ... $\endgroup$
    – PM 2Ring
    Jan 13 at 17:40
  • 1
    $\begingroup$ @Jacob Right, it's not the celestial equator, which is why I wrote "equator" in quotes. My globe displays azimuth & elevation, so the poles are the zenith and nadir, the horizontal circles are lines of equal elevation (corresponding to parallels of latitude on a normal globe), and the vertical semicircles are lines of equal azimuth (corresponding to meridians of longitude on a normal globe). I suppose I should clarify that in my answer... $\endgroup$
    – PM 2Ring
    Jan 16 at 11:15

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .