Solar Module — Day/Night Terminator & Tissot Indicatrix

The cleopatra.basemap.solar module is the flat-2D counterpart to the 3-D globe's directional lighting (TexturedGlobeGlyph, which shades a sphere and produces no geometry). It has two layers:

1. CRS-free solar geometry — pure lon/lat maths that depends on nothing but the clock, with no ephemeris dependency:

   • subsolar_point(when) — the point where the sun is overhead (declination + the equation of time), via the low-precision NOAA/Meeus formulae (sub-degree over the relevant centuries).
   • terminator(when, ...) — the day/night terminator as a small circle at angular distance 90 - refraction from the subsolar point (a great circle only when refraction == 0); pass refraction=-6/-12/-18 for the civil/nautical/astronomical twilight lines.
   • night_polygon(when, ...) — the filled night region as lon/lat rings, closed along the dark pole's edge when a pole is in darkness and otherwise split at the antimeridian so a flat map does not get a band smeared across the whole world.
   • tissot_circles(lons, lats, radius_m, ...) — geodesic circles of a fixed ground radius, in lon/lat, for building a Tissot indicatrix.

2. Artists that draw in data coordinates — add_nightshade builds night_polygon and draws it as a PolyCollection, preserving the current axis limits and returning the artist; add_tissot draws a set of rings the caller supplies (already in the axes' coordinates).

Scope boundary. cleopatra computes lon/lat and draws — it never resolves a CRS or decides the axes' projection. On a plain lon/lat axes the rings are drawn directly. On a projected axes the consumer either pre-projects the vertices and passes them, or passes a transform callable (lon, lat) -> (x, y). The one convenience is an optional crs= shortcut (mutually exclusive with transform=) that reprojects from EPSG:4326 through pyproj — the same optional cleopatra[tiles] extra add_features uses.

Usage

import matplotlib
matplotlib.use("Agg")  # any backend; Agg shown for headless rendering
from datetime import UTC, datetime

import matplotlib.pyplot as plt

from cleopatra.basemap.solar import add_nightshade, add_tissot, tissot_circles

when = datetime(2026, 6, 21, 12, tzinfo=UTC)

fig, ax = plt.subplots()
ax.set_xlim(-180, 180)
ax.set_ylim(-90, 90)

# Shade the night region on a plain lon/lat axes.
add_nightshade(ax, when, alpha=0.35, color="black", zorder=5)

# Tissot indicatrix: geodesic circles in lon/lat. On a lon/lat axes they can be
# drawn directly; on a projected axes, project each ring first (that distortion
# is exactly what the indicatrix shows).
circles = tissot_circles(lons=[-120, 0, 120], lats=[-45, 0, 45], radius_m=5e5)
add_tissot(ax, circles, facecolor="none", edgecolor="crimson")

fig.savefig("nightshade.png")

Note

add_nightshade and add_tissot draw in data coordinates and reproject nothing unless you ask: pass a transform callable, or the optional crs= shortcut (which needs the cleopatra[tiles] extra and raises an actionable ImportError without it). transform= and crs= are mutually exclusive. The night region always spans about half the globe, so a non-global projection (orthographic, azimuthal, a regional CRS) maps the far side outside its domain, to non-finite coordinates; add_nightshade drops those vertices, so the drawn region is then only the part that lies inside the projection — choose a projection whose domain covers the area you are shading. If an out-of-domain mapping drops every vertex, the returned artist is empty (nothing is drawn). A conformal projection such as Web Mercator (crs=3857) instead maps the poles to very large but finite coordinates: nothing is dropped, and a pole-covering fill is simply drawn far outside the usual map extent (the axis limits, which add_nightshade preserves, crop it).

Module Documentation

cleopatra.basemap.solar

Day/night terminator and Tissot-indicatrix artists for flat 2-D axes.

cleopatra can already shade a day/night terminator, but only as lighting on a 3-D sphere (cleopatra.glyphs.globe.textured_globe_glyph.TexturedGlobeGlyph): the input is a world-space light direction, the output is per-face shaded facecolors on an Axes3D, and no lon/lat geometry is ever produced. This module is the flat-map counterpart: it computes the terminator as a small circle about the subsolar point in lon/lat and draws it (and the filled night region, and Tissot distortion circles) on an ordinary matplotlib.axes.Axes.

Scope boundary -- the same split the rest of basemap keeps:

• Solar geometry is CRS-free maths, so it lives here. The subsolar point for a datetime (solar declination plus the Greenwich hour angle via the equation of time) and the terminator small circle depend on nothing but the clock; they are not CRS-dependent. Likewise the geodesic circles behind a Tissot indicatrix are generic spherical geometry.
• Any CRS/projection transform is the consumer's. These helpers compute lon/lat and draw in data coordinates. A consumer whose axes is not lon/lat either pre-projects the vertices and passes them, or passes a transform callable (lon, lat) -> (x, y). This module never resolves a CRS, never decides the axes' projection, and never grows a projection engine -- the same contract cleopatra.basemap.projection.apply_projection_frame states.
• The one convenience is an optional crs=, mirroring cleopatra.basemap.reference.add_features: it reprojects from EPSG:4326 through pyproj (the [tiles] extra) for the common EPSG case, and is mutually exclusive with transform=. It is a shortcut, not a licence to own projections.

Accuracy is deliberately low-precision: the NOAA/Meeus solar-position formulae are sub-degree over the relevant centuries, which is far finer than a shaded overlay needs. No ephemeris dependency.

Example (lon/lat axes)::

from datetime import UTC, datetime
from cleopatra.basemap.solar import add_nightshade

when = datetime(2026, 6, 21, 12, tzinfo=UTC)
add_nightshade(ax, when, alpha=0.35, color="black", zorder=5)   # shade the night

Consumers whose axes is not lon/lat pass a transform callable, or the optional crs= pyproj shortcut, so cleopatra never resolves a projection.

See also cleopatra.glyphs.globe.textured_globe_glyph for the 3-D globe's directional lighting, which answers a different question (shading on a sphere, not geometry on a flat map).

add_nightshade(ax, when, *, refraction=DEFAULT_REFRACTION, n=DEFAULT_TERMINATOR_SAMPLES, transform=None, crs=None, **style)

Shade the night region on ax and return the artist.

Builds the night_polygon rings, maps them into the axes' coordinates, draws them, preserves the current axis limits (the way cleopatra.basemap.reference.add_features does), and returns the artist so the caller can restyle or remove it. Several calls on one axes are fine; each returns its own artist.

Coordinate mapping is the consumer's responsibility:

• Neither transform nor crs: the rings are drawn as lon/lat, i.e. the axes is assumed to be in EPSG:4326 data coordinates.
• transform: a callable applied to each (m, 2) lon/lat ring, returning (m, 2) axes coordinates.
• crs: the optional pyproj shortcut -- reproject from EPSG:4326 to crs (reuses cleopatra.basemap.reference._make_transformer, so it needs the [tiles] extra).

Vertices that a transform or crs maps to non-finite values -- the night region always spans ~half the globe, so a non-global projection sends its far side outside the projection's domain -- are dropped so the fill stays valid (as cleopatra.basemap.reference.add_features does). If a mapping drops every vertex, the returned artist has no polygons (nothing is drawn).

Parameters:

Name	Type	Description	Default
ax	Any	A matplotlib ~matplotlib.axes.Axes with data already plotted.	required
when	datetime	An aware datetime.datetime (naive is treated as UTC).	required
refraction	float	Solar altitude in degrees; see DEFAULT_REFRACTION.	DEFAULT_REFRACTION
n	int	Number of samples along the terminator.	DEFAULT_TERMINATOR_SAMPLES
transform	Callable[[ndarray], ndarray] | None	Optional (lon, lat) -> (x, y) callable operating on (m, 2) arrays. Mutually exclusive with crs.	None
crs	int | str | None	Optional target CRS (EPSG int or CRS string) for the pyproj shortcut. Mutually exclusive with transform.	None
**style	Any	Style overrides forwarded to the PolyCollection (facecolor/color, alpha, zorder, ...).	{}

Returns:

Type	Description
PolyCollection	matplotlib.collections.PolyCollection: The night-shade artist.

Raises:

Type	Description
TypeError	If ax is not a matplotlib Axes.
ValueError	If both transform and crs are given; if refraction/n are out of range (propagated from night_polygon); or if crs cannot be parsed (propagated from cleopatra.basemap.reference._make_transformer).
ImportError	If crs requires reprojection but pyproj (the [tiles] extra) is not installed.

Examples:

• Shade the night region on a lon/lat axes and keep the artist; the axis limits are preserved:

  >>> import matplotlib
  >>> matplotlib.use("Agg")
  >>> import matplotlib.pyplot as plt
  >>> from datetime import UTC, datetime
  >>> fig, ax = plt.subplots()
  >>> _ = ax.set_xlim(-180, 180)
  >>> _ = ax.set_ylim(-90, 90)
  >>> art = add_nightshade(ax, datetime(2026, 6, 21, 12, tzinfo=UTC), alpha=0.3)
  >>> art in ax.collections
  True
  >>> bool(ax.get_xlim() == (-180.0, 180.0))
  True
  

• A consumer whose axes is not lon/lat passes a transform callable (here a trivial scale), so cleopatra never resolves the projection:

  >>> import matplotlib
  >>> matplotlib.use("Agg")
  >>> import matplotlib.pyplot as plt
  >>> from datetime import UTC, datetime
  >>> fig, ax = plt.subplots()
  >>> art = add_nightshade(
  ...     ax, datetime(2026, 3, 20, 12, tzinfo=UTC), transform=lambda a: a * 2.0
  ... )
  >>> len(art.get_paths()) >= 1
  True
  

Source code in src/cleopatra/basemap/solar.py

def add_nightshade(
    ax: Any,
    when: datetime,
    *,
    refraction: float = DEFAULT_REFRACTION,
    n: int = DEFAULT_TERMINATOR_SAMPLES,
    transform: Callable[[np.ndarray], np.ndarray] | None = None,
    crs: int | str | None = None,
    **style: Any,
) -> PolyCollection:
    """Shade the night region on ``ax`` and return the artist.

    Builds the `night_polygon` rings, maps them into the axes' coordinates,
    draws them, preserves the current axis limits (the way
    `cleopatra.basemap.reference.add_features` does), and returns the artist so
    the caller can restyle or remove it. Several calls on one axes are fine;
    each returns its own artist.

    Coordinate mapping is the consumer's responsibility:

    - Neither ``transform`` nor ``crs``: the rings are drawn as lon/lat, i.e.
      the axes is assumed to be in EPSG:4326 data coordinates.
    - ``transform``: a callable applied to each ``(m, 2)`` lon/lat ring,
      returning ``(m, 2)`` axes coordinates.
    - ``crs``: the optional pyproj shortcut -- reproject from EPSG:4326 to
      ``crs`` (reuses `cleopatra.basemap.reference._make_transformer`, so it
      needs the ``[tiles]`` extra).

    Vertices that a ``transform`` or ``crs`` maps to non-finite values -- the
    night region always spans ~half the globe, so a non-global projection sends
    its far side outside the projection's domain -- are dropped so the fill stays
    valid (as `cleopatra.basemap.reference.add_features` does). If a mapping drops
    every vertex, the returned artist has no polygons (nothing is drawn).

    Args:
        ax: A matplotlib `~matplotlib.axes.Axes` with data already plotted.
        when: An aware `datetime.datetime` (naive is treated as UTC).
        refraction: Solar altitude in degrees; see `DEFAULT_REFRACTION`.
        n: Number of samples along the terminator.
        transform: Optional ``(lon, lat) -> (x, y)`` callable operating on
            ``(m, 2)`` arrays. Mutually exclusive with ``crs``.
        crs: Optional target CRS (EPSG int or CRS string) for the pyproj
            shortcut. Mutually exclusive with ``transform``.
        **style: Style overrides forwarded to the `PolyCollection`
            (``facecolor``/``color``, ``alpha``, ``zorder``, ...).

    Returns:
        matplotlib.collections.PolyCollection: The night-shade artist.

    Raises:
        TypeError: If ``ax`` is not a matplotlib Axes.
        ValueError: If both ``transform`` and ``crs`` are given; if
            ``refraction``/``n`` are out of range (propagated from
            `night_polygon`); or if ``crs`` cannot be parsed (propagated from
            `cleopatra.basemap.reference._make_transformer`).
        ImportError: If ``crs`` requires reprojection but ``pyproj`` (the
            ``[tiles]`` extra) is not installed.

    Examples:
        - Shade the night region on a lon/lat axes and keep the artist; the axis
          limits are preserved:
            ```python
            >>> import matplotlib
            >>> matplotlib.use("Agg")
            >>> import matplotlib.pyplot as plt
            >>> from datetime import UTC, datetime
            >>> fig, ax = plt.subplots()
            >>> _ = ax.set_xlim(-180, 180)
            >>> _ = ax.set_ylim(-90, 90)
            >>> art = add_nightshade(ax, datetime(2026, 6, 21, 12, tzinfo=UTC), alpha=0.3)
            >>> art in ax.collections
            True
            >>> bool(ax.get_xlim() == (-180.0, 180.0))
            True

            ```
        - A consumer whose axes is not lon/lat passes a ``transform`` callable
          (here a trivial scale), so cleopatra never resolves the projection:
            ```python
            >>> import matplotlib
            >>> matplotlib.use("Agg")
            >>> import matplotlib.pyplot as plt
            >>> from datetime import UTC, datetime
            >>> fig, ax = plt.subplots()
            >>> art = add_nightshade(
            ...     ax, datetime(2026, 3, 20, 12, tzinfo=UTC), transform=lambda a: a * 2.0
            ... )
            >>> len(art.get_paths()) >= 1
            True

            ```
    """
    if transform is not None and crs is not None:
        raise ValueError("Pass at most one of transform= or crs=, not both.")
    _validate_axes(ax)

    rings = night_polygon(when, refraction=refraction, n=n)
    if transform is not None:
        rings = [np.asarray(transform(ring), dtype=float) for ring in rings]
    elif crs is not None and not _is_4326(crs):
        transformer = _make_transformer(crs)
        rings = [_reproject_arr(ring, transformer) for ring in rings]

    # A transform or reprojection can map points outside the projection's domain
    # (e.g. a non-global projection's far side) to non-finite values; the night
    # region always spans ~half the globe, so drop those rows per ring -- and
    # rings left with < 3 vertices -- to keep the fill valid, mirroring
    # cleopatra.basemap.reference.add_features. If every vertex is dropped the
    # artist is simply empty.
    rings = [ring[np.isfinite(ring).all(axis=1)] for ring in rings]
    rings = [ring for ring in rings if len(ring) >= 3]

    opts: dict[str, Any] = {"facecolor": "black", "edgecolor": "none", "alpha": 0.35}
    if "color" in style:
        # `color` sets both face and edge; drop the split defaults to avoid a clash.
        opts.pop("facecolor")
        opts.pop("edgecolor")
    opts.update(style)

    xlim, ylim = ax.get_xlim(), ax.get_ylim()
    artist = PolyCollection(rings, **opts)
    artist.set_transform(ax.transData)
    ax.add_collection(artist)
    ax.set_xlim(xlim)
    ax.set_ylim(ylim)
    return artist

add_tissot(ax, ellipses, **style)

Draw pre-projected Tissot rings on ax and return the artist.

cleopatra draws what it is given and computes no distortion: the shape a circle takes is a property of the consumer's projection, so the consumer projects tissot_circles output through its own transform and passes the result here. It draws the rings as an unfilled PolyCollection (outline only by default) -- the same "sequence of polygons with differing vertex counts" shape cleopatra.glyphs.primitives.polygon_glyph.PolygonGlyph handles -- and preserves the current axis limits.

Parameters:

Name	Type	Description	Default
ax	Any	A matplotlib ~matplotlib.axes.Axes.	required
ellipses	Sequence[ndarray]	A sequence of (m, 2) vertex arrays in the axes' own coordinates (already projected by the consumer).	required
**style	Any	Style overrides forwarded to the underlying collection (facecolor, edgecolor, linewidth, ...).	{}

Returns:

Type	Description
PolyCollection	matplotlib.collections.PolyCollection: The Tissot artist.

Raises:

Type	Description
TypeError	If ax is not a matplotlib Axes.
ValueError	If an ellipses element cannot be coerced to a float (m, 2) array (e.g. a ragged ring), propagated from numpy.asarray.

Examples:

• Draw two supplied rings (already in axes coordinates) and keep them:

  >>> import matplotlib
  >>> matplotlib.use("Agg")
  >>> import matplotlib.pyplot as plt
  >>> import numpy as np
  >>> theta = np.linspace(0, 2 * np.pi, 16)
  >>> circle = np.column_stack([np.cos(theta), np.sin(theta)])
  >>> fig, ax = plt.subplots()
  >>> art = add_tissot(ax, [circle, circle + 3.0], edgecolor="crimson")
  >>> len(art.get_paths())
  2
  >>> art in ax.collections
  True
  

• Generate geodesic circles with tissot_circles and draw them on a lon/lat axes:

  >>> import matplotlib
  >>> matplotlib.use("Agg")
  >>> import matplotlib.pyplot as plt
  >>> from cleopatra.basemap.solar import tissot_circles
  >>> circles = tissot_circles([-90.0, 0.0, 90.0], [0.0, 0.0, 0.0], 5e5)
  >>> fig, ax = plt.subplots()
  >>> art = add_tissot(ax, circles, edgecolor="navy")
  >>> len(art.get_paths())
  3
  

Source code in src/cleopatra/basemap/solar.py

def add_tissot(ax: Any, ellipses: Sequence[np.ndarray], **style: Any) -> PolyCollection:
    """Draw pre-projected Tissot rings on ``ax`` and return the artist.

    **cleopatra draws what it is given and computes no distortion**: the shape a
    circle takes is a property of the consumer's projection, so the consumer
    projects `tissot_circles` output through its own transform and passes the
    result here. It draws the rings as an unfilled `PolyCollection` (outline only
    by default) -- the same "sequence of polygons with differing vertex counts"
    shape `cleopatra.glyphs.primitives.polygon_glyph.PolygonGlyph` handles -- and
    preserves the current axis limits.

    Args:
        ax: A matplotlib `~matplotlib.axes.Axes`.
        ellipses: A sequence of ``(m, 2)`` vertex arrays in the axes' own
            coordinates (already projected by the consumer).
        **style: Style overrides forwarded to the underlying collection
            (``facecolor``, ``edgecolor``, ``linewidth``, ...).

    Returns:
        matplotlib.collections.PolyCollection: The Tissot artist.

    Raises:
        TypeError: If ``ax`` is not a matplotlib Axes.
        ValueError: If an ``ellipses`` element cannot be coerced to a float
            ``(m, 2)`` array (e.g. a ragged ring), propagated from `numpy.asarray`.

    Examples:
        - Draw two supplied rings (already in axes coordinates) and keep them:
            ```python
            >>> import matplotlib
            >>> matplotlib.use("Agg")
            >>> import matplotlib.pyplot as plt
            >>> import numpy as np
            >>> theta = np.linspace(0, 2 * np.pi, 16)
            >>> circle = np.column_stack([np.cos(theta), np.sin(theta)])
            >>> fig, ax = plt.subplots()
            >>> art = add_tissot(ax, [circle, circle + 3.0], edgecolor="crimson")
            >>> len(art.get_paths())
            2
            >>> art in ax.collections
            True

            ```
        - Generate geodesic circles with `tissot_circles` and draw them on a
          lon/lat axes:
            ```python
            >>> import matplotlib
            >>> matplotlib.use("Agg")
            >>> import matplotlib.pyplot as plt
            >>> from cleopatra.basemap.solar import tissot_circles
            >>> circles = tissot_circles([-90.0, 0.0, 90.0], [0.0, 0.0, 0.0], 5e5)
            >>> fig, ax = plt.subplots()
            >>> art = add_tissot(ax, circles, edgecolor="navy")
            >>> len(art.get_paths())
            3

            ```
    """
    _validate_axes(ax)
    verts = [np.asarray(ring, dtype=float) for ring in ellipses]

    opts: dict[str, Any] = {"facecolor": "none", "edgecolor": "black", "linewidth": 0.8}
    if "color" in style:
        # `color` sets both face and edge; drop the split defaults to avoid a clash.
        opts.pop("facecolor")
        opts.pop("edgecolor")
    opts.update(style)

    xlim, ylim = ax.get_xlim(), ax.get_ylim()
    artist = PolyCollection(verts, **opts)
    artist.set_transform(ax.transData)
    ax.add_collection(artist)
    ax.set_xlim(xlim)
    ax.set_ylim(ylim)
    return artist

night_polygon(when, *, refraction=DEFAULT_REFRACTION, n=DEFAULT_TERMINATOR_SAMPLES)

Return the filled night region as lon/lat rings.

The rings are pre-split at the antimeridian so that a flat (lon/lat) map does not get a band smeared across the whole world when the night region wraps.

Parameters:

Name	Type	Description	Default
when	datetime	An aware datetime.datetime (naive is treated as UTC).	required
refraction	float	Solar altitude in degrees; see DEFAULT_REFRACTION.	DEFAULT_REFRACTION
n	int	Number of samples along the terminator used to build the region.	DEFAULT_TERMINATOR_SAMPLES

Returns:

Type	Description
list[ndarray]	list[numpy.ndarray]: One or more (m, 2) lon/lat rings covering the
list[ndarray]	night side; more than one where the region crosses the antimeridian.
list[ndarray]	The rings are open (first vertex != last); matplotlib closes polygons
list[ndarray]	on fill.

Raises:

Type	Description
ValueError	If refraction is outside (-90, 0] or n < 4 (via terminator).

Examples:

• At a solstice one pole is in darkness, so the night region is a single ring that runs to that pole (the south pole, lat -90, in June):

  >>> from datetime import UTC, datetime
  >>> rings = night_polygon(datetime(2026, 6, 21, 12, 0, tzinfo=UTC))
  >>> len(rings)
  1
  >>> float(rings[0][:, 1].min())
  -90.0
  

• When the night region straddles the antimeridian it is split into two rings so a flat map does not smear a band across the world:

  >>> from datetime import UTC, datetime
  >>> rings = night_polygon(datetime(2026, 3, 20, 12, 0, tzinfo=UTC))
  >>> len(rings)
  2
  

Source code in src/cleopatra/basemap/solar.py

def night_polygon(
    when: datetime,
    *,
    refraction: float = DEFAULT_REFRACTION,
    n: int = DEFAULT_TERMINATOR_SAMPLES,
) -> list[np.ndarray]:
    """Return the filled night region as lon/lat rings.

    The rings are pre-split at the antimeridian so that a flat (lon/lat) map
    does not get a band smeared across the whole world when the night region
    wraps.

    Args:
        when: An aware `datetime.datetime` (naive is treated as UTC).
        refraction: Solar altitude in degrees; see `DEFAULT_REFRACTION`.
        n: Number of samples along the terminator used to build the region.

    Returns:
        list[numpy.ndarray]: One or more ``(m, 2)`` lon/lat rings covering the
        night side; more than one where the region crosses the antimeridian.
        The rings are open (first vertex != last); matplotlib closes polygons
        on fill.

    Raises:
        ValueError: If ``refraction`` is outside ``(-90, 0]`` or ``n < 4`` (via
            `terminator`).

    Examples:
        - At a solstice one pole is in darkness, so the night region is a single
          ring that runs to that pole (the south pole, lat -90, in June):
            ```python
            >>> from datetime import UTC, datetime
            >>> rings = night_polygon(datetime(2026, 6, 21, 12, 0, tzinfo=UTC))
            >>> len(rings)
            1
            >>> float(rings[0][:, 1].min())
            -90.0

            ```
        - When the night region straddles the antimeridian it is split into two
          rings so a flat map does not smear a band across the world:
            ```python
            >>> from datetime import UTC, datetime
            >>> rings = night_polygon(datetime(2026, 3, 20, 12, 0, tzinfo=UTC))
            >>> len(rings)
            2

            ```
    """
    _, lat_s = subsolar_point(when)
    ring = terminator(when, refraction=refraction, n=n)

    # The sun's altitude at the north/south pole is +lat_s / -lat_s. A pole is
    # in night when its altitude drops below ``refraction``.
    north_dark = lat_s < refraction
    south_dark = -lat_s < refraction
    if north_dark != south_dark:
        # A geographic pole is enclosed: the night cap wraps every longitude, so
        # close the (single-valued) terminator along that pole's map edge.
        dark_lat = 90.0 if north_dark else -90.0
        return [_pole_cap_ring(ring, dark_lat)]

    # No pole enclosed: the night cap is a simple sub-hemispheric region bounded
    # by the terminator loop. Split it wherever it crosses the antimeridian.
    return _split_antimeridian(ring)

subsolar_point(when)

Return the subsolar point (lon, lat) in degrees for when.

The subsolar point is where the sun is directly overhead: its latitude is the solar declination, its longitude is derived from the Greenwich hour angle (apparent solar time via the equation of time). Pure maths -- no matplotlib import needed.

Parameters:

Name	Type	Description	Default
when	datetime	An aware datetime.datetime. A naive datetime is treated as UTC.	required

Returns:

Type	Description
float	tuple[float, float]: (lon, lat) in degrees, lon in
float	(-180, 180] and lat in [-90, 90].

Examples:

• At the June solstice the sun is overhead near the Tropic of Cancer (~23.4 deg N), close to Greenwich at 12:00 UTC:

  >>> from datetime import UTC, datetime
  >>> lon, lat = subsolar_point(datetime(2026, 6, 21, 12, 0, tzinfo=UTC))
  >>> 23.0 < lat < 24.0
  True
  >>> abs(lon) < 5.0
  True
  

• The subsolar meridian tracks the sun ~15 deg westward each hour:

  >>> from datetime import UTC, datetime, timedelta
  >>> noon = datetime(2026, 3, 20, 12, 0, tzinfo=UTC)
  >>> before, _ = subsolar_point(noon)
  >>> after, _ = subsolar_point(noon + timedelta(hours=1))
  >>> -15.5 < after - before < -14.5
  True
  

Source code in src/cleopatra/basemap/solar.py

def subsolar_point(when: datetime) -> tuple[float, float]:
    """Return the subsolar point ``(lon, lat)`` in degrees for ``when``.

    The subsolar point is where the sun is directly overhead: its latitude is
    the solar declination, its longitude is derived from the Greenwich hour
    angle (apparent solar time via the equation of time). Pure maths -- no
    matplotlib import needed.

    Args:
        when: An aware `datetime.datetime`. A naive datetime is treated as UTC.

    Returns:
        tuple[float, float]: ``(lon, lat)`` in degrees, ``lon`` in
        ``(-180, 180]`` and ``lat`` in ``[-90, 90]``.

    Examples:
        - At the June solstice the sun is overhead near the Tropic of Cancer
          (~23.4 deg N), close to Greenwich at 12:00 UTC:
            ```python
            >>> from datetime import UTC, datetime
            >>> lon, lat = subsolar_point(datetime(2026, 6, 21, 12, 0, tzinfo=UTC))
            >>> 23.0 < lat < 24.0
            True
            >>> abs(lon) < 5.0
            True

            ```
        - The subsolar meridian tracks the sun ~15 deg westward each hour:
            ```python
            >>> from datetime import UTC, datetime, timedelta
            >>> noon = datetime(2026, 3, 20, 12, 0, tzinfo=UTC)
            >>> before, _ = subsolar_point(noon)
            >>> after, _ = subsolar_point(noon + timedelta(hours=1))
            >>> -15.5 < after - before < -14.5
            True

            ```
    """
    dt = when.replace(tzinfo=UTC) if when.tzinfo is None else when.astimezone(UTC)
    hours = dt.hour + dt.minute / 60.0 + dt.second / 3600.0 + dt.microsecond / 3.6e9

    # Julian Day (Gregorian) then Julian centuries since J2000.0.
    a = (14 - dt.month) // 12
    y = dt.year + 4800 - a
    m = dt.month + 12 * a - 3
    jdn = dt.day + (153 * m + 2) // 5 + 365 * y + y // 4 - y // 100 + y // 400 - 32045
    jd = jdn + (hours - 12.0) / 24.0
    t = (jd - 2451545.0) / 36525.0

    # NOAA low-precision solar position (Meeus, Astronomical Algorithms).
    mean_long = np.radians((280.46646 + t * (36000.76983 + t * 0.0003032)) % 360.0)
    mean_anom = np.radians(357.52911 + t * (35999.05029 - 0.0001537 * t))
    eccentricity = 0.016708634 - t * (0.000042037 + 0.0000001267 * t)
    center = (
        np.sin(mean_anom) * (1.914602 - t * (0.004817 + 0.000014 * t))
        + np.sin(2 * mean_anom) * (0.019993 - 0.000101 * t)
        + np.sin(3 * mean_anom) * 0.000289
    )
    true_long = np.degrees(mean_long) + center
    omega = np.radians(125.04 - 1934.136 * t)
    app_long = np.radians(true_long - 0.00569 - 0.00478 * np.sin(omega))
    obliquity = np.radians(
        23.0
        + (26.0 + (21.448 - t * (46.8150 + t * (0.00059 - t * 0.001813))) / 60.0) / 60.0
        + 0.00256 * np.cos(omega)
    )

    declination = np.degrees(np.arcsin(np.sin(obliquity) * np.sin(app_long)))

    # Equation of time (minutes) -> subsolar longitude via apparent solar time.
    var_y = np.tan(obliquity / 2.0) ** 2
    eqn_time = 4.0 * np.degrees(
        var_y * np.sin(2 * mean_long)
        - 2 * eccentricity * np.sin(mean_anom)
        + 4 * eccentricity * var_y * np.sin(mean_anom) * np.cos(2 * mean_long)
        - 0.5 * var_y**2 * np.sin(4 * mean_long)
        - 1.25 * eccentricity**2 * np.sin(2 * mean_anom)
    )
    lon = -15.0 * (hours - 12.0 + eqn_time / 60.0)
    return float(_wrap_longitude(lon)), float(declination)

terminator(when, *, refraction=DEFAULT_REFRACTION, n=DEFAULT_TERMINATOR_SAMPLES)

Return the day/night terminator as lon/lat vertices.

The terminator is the small circle at angular distance 90 - refraction from the subsolar point -- a great circle only when refraction == 0.

Parameters:

Name	Type	Description	Default
when	datetime	An aware datetime.datetime (naive is treated as UTC).	required
refraction	float	Solar altitude in degrees defining the terminator; see DEFAULT_REFRACTION. Use -6 / -12 / -18 for civil / nautical / astronomical twilight.	DEFAULT_REFRACTION
n	int	Number of samples along the circle. The ring is closed, so n includes the duplicated endpoint and the effective resolution is n - 1.	DEFAULT_TERMINATOR_SAMPLES

Returns:

Type	Description
ndarray	numpy.ndarray: An (n, 2) array of (lon, lat) degrees, densified
ndarray	for a smooth curve. The ring is closed (last vertex equals the first).

Raises:

Type	Description
ValueError	If refraction is outside (-90, 0] (0 is the geometric terminator; a positive value is not a terminator and a value of -90 or below is degenerate), or if n < 4 (the ring's duplicated endpoint means 4 samples are the fewest that give 3 distinct vertices).

Examples:

• The default terminator is a closed ring of 720 lon/lat vertices:

  >>> import numpy as np
  >>> from datetime import UTC, datetime
  >>> ring = terminator(datetime(2026, 6, 21, 12, 0, tzinfo=UTC))
  >>> ring.shape
  (720, 2)
  >>> bool(np.allclose(ring[0], ring[-1]))
  True
  

• A coarser ring stays within the lon/lat bounds; refraction shifts the ring's angular radius (-6 gives the civil-twilight line):

  >>> import numpy as np
  >>> from datetime import UTC, datetime
  >>> ring = terminator(datetime(2026, 6, 21, 12, 0, tzinfo=UTC), refraction=-6.0, n=180)
  >>> ring.shape
  (180, 2)
  >>> bool(np.all(np.abs(ring[:, 0]) <= 180.0) and np.all(np.abs(ring[:, 1]) <= 90.0))
  True
  

Source code in src/cleopatra/basemap/solar.py

def terminator(
    when: datetime,
    *,
    refraction: float = DEFAULT_REFRACTION,
    n: int = DEFAULT_TERMINATOR_SAMPLES,
) -> np.ndarray:
    """Return the day/night terminator as lon/lat vertices.

    The terminator is the small circle at angular distance ``90 - refraction``
    from the subsolar point -- a great circle only when ``refraction == 0``.

    Args:
        when: An aware `datetime.datetime` (naive is treated as UTC).
        refraction: Solar altitude in degrees defining the terminator; see
            `DEFAULT_REFRACTION`. Use ``-6`` / ``-12`` / ``-18`` for civil /
            nautical / astronomical twilight.
        n: Number of samples along the circle. The ring is closed, so ``n``
            includes the duplicated endpoint and the effective resolution is
            ``n - 1``.

    Returns:
        numpy.ndarray: An ``(n, 2)`` array of ``(lon, lat)`` degrees, densified
        for a smooth curve. The ring is closed (last vertex equals the first).

    Raises:
        ValueError: If ``refraction`` is outside ``(-90, 0]`` (0 is the
            geometric terminator; a positive value is not a terminator and a
            value of -90 or below is degenerate), or if ``n < 4`` (the ring's
            duplicated endpoint means 4 samples are the fewest that give 3
            distinct vertices).

    Examples:
        - The default terminator is a closed ring of 720 lon/lat vertices:
            ```python
            >>> import numpy as np
            >>> from datetime import UTC, datetime
            >>> ring = terminator(datetime(2026, 6, 21, 12, 0, tzinfo=UTC))
            >>> ring.shape
            (720, 2)
            >>> bool(np.allclose(ring[0], ring[-1]))
            True

            ```
        - A coarser ring stays within the lon/lat bounds; ``refraction`` shifts
          the ring's angular radius (``-6`` gives the civil-twilight line):
            ```python
            >>> import numpy as np
            >>> from datetime import UTC, datetime
            >>> ring = terminator(datetime(2026, 6, 21, 12, 0, tzinfo=UTC), refraction=-6.0, n=180)
            >>> ring.shape
            (180, 2)
            >>> bool(np.all(np.abs(ring[:, 0]) <= 180.0) and np.all(np.abs(ring[:, 1]) <= 90.0))
            True

            ```
    """
    if not -90.0 < refraction <= 0.0:
        raise ValueError(
            "refraction must be in (-90, 0] degrees (0 is the geometric "
            f"terminator; -6/-12/-18 are the twilight lines); got {refraction}."
        )
    if n < 4:
        # linspace(0, 2*pi, n) duplicates the endpoint, so n counts one closing
        # vertex; a non-degenerate ring needs 3 distinct vertices, i.e. n >= 4.
        raise ValueError(f"n must be at least 4 to form a ring; got {n}.")

    lon_s, lat_s = subsolar_point(when)
    # The terminator is the small circle at angular distance 90 - refraction from
    # the subsolar point (a great circle only when refraction == 0).
    return _small_circle(lon_s, lat_s, np.radians(90.0 - refraction), n)

tissot_circles(lons, lats, radius_m, *, n=DEFAULT_TISSOT_SAMPLES)

Return geodesic circles of a fixed ground radius as lon/lat rings.

Generic spherical geometry: each circle is the locus of points a fixed great-circle distance (radius_m on a sphere of MEAN_EARTH_RADIUS_M) from its centre. No CRS is involved; the consumer projects these rings to show how a projection distorts them (that is the Tissot indicatrix).

Parameters:

Name	Type	Description	Default
lons	Sequence[float] | ndarray	Circle-centre longitudes in degrees.	required
lats	Sequence[float] | ndarray	Circle-centre latitudes in degrees. Same length as lons.	required
radius_m	float	Ground radius of each circle in metres.	required
n	int	Number of vertices per circle.	DEFAULT_TISSOT_SAMPLES

Returns:

Type	Description
list[ndarray]	list[numpy.ndarray]: One (n, 2) lon/lat ring per centre, in input
list[ndarray]	order. Each ring is closed (the last vertex duplicates the first).

Raises:

Type	Description
ValueError	If lons or lats is not 1-D or they differ in shape, if radius_m is not a positive sub-antipodal radius (0 < radius_m < pi * R), or if n < 4.

Examples:

• One ~500 km circle around the origin, sampled coarsely:

  >>> import numpy as np
  >>> rings = tissot_circles([0.0], [0.0], 5e5, n=8)
  >>> len(rings)
  1
  >>> rings[0].shape
  (8, 2)
  >>> bool(np.isfinite(rings[0]).all())
  True
  

• The ground radius maps to a fixed angular radius (~4.5 deg at 500 km):

  >>> rings = tissot_circles([0.0], [0.0], 5e5, n=4)
  >>> round(float(rings[0][:, 1].max()), 1)
  4.5
  

Source code in src/cleopatra/basemap/solar.py

def tissot_circles(
    lons: Sequence[float] | np.ndarray,
    lats: Sequence[float] | np.ndarray,
    radius_m: float,
    *,
    n: int = DEFAULT_TISSOT_SAMPLES,
) -> list[np.ndarray]:
    """Return geodesic circles of a fixed ground radius as lon/lat rings.

    Generic spherical geometry: each circle is the locus of points a fixed
    great-circle distance (``radius_m`` on a sphere of `MEAN_EARTH_RADIUS_M`)
    from its centre. No CRS is involved; the consumer projects these rings to
    show how a projection distorts them (that is the Tissot indicatrix).

    Args:
        lons: Circle-centre longitudes in degrees.
        lats: Circle-centre latitudes in degrees. Same length as ``lons``.
        radius_m: Ground radius of each circle in metres.
        n: Number of vertices per circle.

    Returns:
        list[numpy.ndarray]: One ``(n, 2)`` lon/lat ring per centre, in input
        order. Each ring is closed (the last vertex duplicates the first).

    Raises:
        ValueError: If ``lons`` or ``lats`` is not 1-D or they differ in shape,
            if ``radius_m`` is not a positive sub-antipodal radius
            (``0 < radius_m < pi * R``), or if ``n < 4``.

    Examples:
        - One ~500 km circle around the origin, sampled coarsely:
            ```python
            >>> import numpy as np
            >>> rings = tissot_circles([0.0], [0.0], 5e5, n=8)
            >>> len(rings)
            1
            >>> rings[0].shape
            (8, 2)
            >>> bool(np.isfinite(rings[0]).all())
            True

            ```
        - The ground radius maps to a fixed angular radius (~4.5 deg at 500 km):
            ```python
            >>> rings = tissot_circles([0.0], [0.0], 5e5, n=4)
            >>> round(float(rings[0][:, 1].max()), 1)
            4.5

            ```
    """
    lon_arr = np.atleast_1d(np.asarray(lons, dtype=float))
    lat_arr = np.atleast_1d(np.asarray(lats, dtype=float))
    if lon_arr.ndim != 1 or lat_arr.ndim != 1:
        raise ValueError(
            f"lons and lats must be 1-D sequences of centres; got {lon_arr.ndim}-D "
            f"and {lat_arr.ndim}-D."
        )
    if lon_arr.shape != lat_arr.shape:
        raise ValueError(
            f"lons and lats must have the same shape; got {lon_arr.shape} and {lat_arr.shape}."
        )
    if not 0.0 < radius_m < np.pi * MEAN_EARTH_RADIUS_M:
        raise ValueError(
            "radius_m must be a positive, sub-antipodal ground radius in metres "
            f"(0, {np.pi * MEAN_EARTH_RADIUS_M:.0f}); got {radius_m}."
        )
    if n < 4:
        raise ValueError(f"n must be at least 4 to form a ring; got {n}.")

    radius_rad = radius_m / MEAN_EARTH_RADIUS_M
    return [
        _small_circle(float(lon), float(lat), radius_rad, n)
        for lon, lat in zip(lon_arr, lat_arr, strict=True)
    ]