Astronomical Rendering

TerraLab achieves sub-arcminute accuracy in celestial object positioning by integrating:

Skyfield: Python astronomy library with JPL ephemerides
Gaia DR3: European Space Agency star catalog with parallax and proper motion
DE421: JPL Development Ephemeris (planetary positions 1900–2050)
Atmospheric refraction: Corrects for light bending near the horizon

This page explains the astronomical calculations that power TerraLab’s precision sky rendering.

Skyfield Integration

What is Skyfield?

Skyfield is a pure-Python astronomy library that computes positions of celestial objects using the same algorithms as professional observatories. It provides:

High-precision ephemerides: Planet/moon positions from JPL DE421
Coordinate transformations: ICRS (RA/Dec) ↔ Altitude/Azimuth
Time scales: UTC, UT1, TT (Terrestrial Time), TAI
Observational corrections: Precession, nutation, aberration, refraction

Skyfield is developed by Brandon Rhodes and used by NASA, ESA, and research institutions worldwide.Website: rhodesmill.org/skyfield

Installation

TerraLab loads Skyfield asynchronously to avoid UI freezing:

class SkyfieldLoaderWorker(QObject):
    skyfield_ready = pyqtSignal(object, object)  # (ts, eph)
    
    def load(self):
        try:
            ts = load.timescale()       # Load time scale data
            eph = load('de421.bsp')     # Load JPL ephemeris
            self.skyfield_ready.emit(ts, eph)
        except Exception as e:
            print(f"Skyfield Error: {e}")
            self.skyfield_ready.emit(None, None)

Source: TerraLab/widgets/sky_widget.py:775-790

DE421 Ephemeris

DE421 is a binary SPK (Spacecraft Planet Kernel) file containing:

Sun, Moon, 8 planets: Mercury through Neptune
Pluto, major asteroids: Ceres, Vesta, Pallas
Earth-Moon barycenter: For precise lunar parallax
Time range: 1900–2050 (±75 years from J2000.0)

File location: TerraLab/data/stars/de421.bsp Accuracy: Sub-arcsecond for planets, ~1” for the Moon (after parallax correction) Source: README.md line 67, TerraLab/widgets/sky_widget.py:785

DE421 becomes inaccurate before 1900 or after 2050. For historical/future dates, use DE422 (3000 BCE – 3000 CE) or DE440 (current JPL production ephemeris).

Gaia Star Catalog

Gaia Mission Overview

Gaia is an ESA space observatory launched in 2013 that has measured positions, distances, and motions of 1.8 billion stars with unprecedented precision. Data Release 3 (DR3) includes:

Positions: RA/Dec accurate to 0.1–1 milliarcseconds (mas)
Parallax: Distances via trigonometric parallax (±0.1 mas)
Proper motion: μ_RA, μ_Dec in mas/year
Photometry: G magnitude (broad band), BP (blue), RP (red)
Radial velocity: For bright stars (V < 12 mag)

Source: README.md line 16, 84

TerraLab’s Gaia Subset

TerraLab uses a filtered catalog containing:

~100,000 stars brighter than G = 10.0 mag
All named stars: Proper names + Bayer/Flamsteed designations
Key deep-sky objects: M31, M45, Double Cluster, etc.

File format: JSON array

{
  "data": [
    [
      "5853498713190525696",  // Gaia source_id
      "Sirius",               // Proper name
      101.287,                // RA (degrees)
      -16.716,                // Dec (degrees)
      -1.46,                  // G magnitude
      null,                   // Parallax (not used)
      null,                   // μ_RA
      null,                   // μ_Dec
      0.009                   // BP-RP color index
    ],
    ...
  ]
}

Location: TerraLab/data/stars/gaia_stars.json Source: TerraLab/widgets/sky_widget.py:670-707

Catalog Loading (Optimized)

Parsing JSON is slow (~2 seconds for 100k stars). TerraLab uses binary caching:

class CatalogLoaderWorker(QObject):
    catalog_ready = pyqtSignal(list, object, ...)
    
    def load(self, json_path):
        # Check for .npy cache
        cache_dir = os.path.dirname(json_path)
        npy_ra = os.path.join(cache_dir, 'gaia_cache_ra.npy')
        
        if os.path.exists(npy_ra):
            # Fast path: load NumPy arrays directly
            ra = np.load(npy_ra)
            dec = np.load(npy_dec_path)
            mag = np.load(npy_mag_path)
            # ... (~10ms total)
        else:
            # Slow path: parse JSON, save .npy cache
            with open(json_path) as f:
                data = json.load(f)
            # ... extract arrays, save for next startup

Source: TerraLab/widgets/sky_widget.py:604-772 Performance:

First load (JSON parsing): ~2000ms
Cached load (NumPy): ~10ms (200× faster)

Color Calculation from BP-RP

The BP-RP color index measures a star’s temperature:

BP-RP	Color	Spectral Type	RGB (approx)
Below 0.0	Blue	O, B	(160, 190, 255)
0.0–0.5	Blue-white	A	(200, 220, 255)
0.5–1.0	White	F, G	(255, 255, 220)
1.0–2.0	Orange	K	(255, 200, 140)
Above 2.0	Red	M	(255, 175, 100)

Implementation:

def get_rgb(bp_rp: float) -> tuple:
    if bp_rp < 0.0:
        return (160, 190, 255)  # Hot blue stars
    elif bp_rp < 0.5:
        t = bp_rp / 0.5
        return (160 + int(95*t), 190 + int(65*t), 255)
    elif bp_rp < 1.0:
        t = (bp_rp - 0.5) / 0.5
        return (255, 255, 255 - int(55*t))
    elif bp_rp < 2.0:
        t = (bp_rp - 1.0) / 1.0
        return (255, 255 - int(80*t), 200 - int(100*t))
    else:
        return (255, 175, 100)  # Cool red stars

Source: TerraLab/widgets/sky_widget.py:730-742

Coordinate Transformations

Equatorial to Horizontal

Stars are stored in equatorial coordinates (RA/Dec), but rendering requires horizontal coordinates (Altitude/Azimuth) for the observer’s location and time.

Step 1: Calculate Hour Angle

HA = LST - RA

Where:

$LST$ = Local Sidereal Time (degrees)
$RA$ = Right Ascension (degrees)

Step 2: Compute Altitude

\sin(\text{Alt}) = \sin(\text{Dec}) \cdot \sin(\text{Lat}) + \cos(\text{Dec}) \cdot \cos(\text{Lat}) \cdot \cos(HA)

Step 3: Compute Azimuth

\cos(\text{Az}) = \frac{\sin(\text{Dec}) - \sin(\text{Alt}) \cdot \sin(\text{Lat})}{\cos(\text{Alt}) \cdot \cos(\text{Lat})}

With quadrant correction:

\text{Az} = \begin{cases} \arccos(\cos(\text{Az})) & \text{if } \sin(HA) \leq 0 \\ 360° - \arccos(\cos(\text{Az})) & \text{if } \sin(HA) > 0 \end{cases}

Source: TerraLab/widgets/sky_widget.py:963-997

NumPy Vectorization

TerraLab computes all stars simultaneously using vectorized operations:

# Input arrays (100k stars)
ra = np.array([...])
dec = np.array([...])

# Physics (broadcasted)
ha = lst - ra
ha_rad = np.radians(ha)
dec_rad = np.radians(dec)

sin_alt = (np.sin(dec_rad) * sin_lat + 
           np.cos(dec_rad) * cos_lat * np.cos(ha_rad))
alt_rad = np.arcsin(np.clip(sin_alt, -1.0, 1.0))

# ... azimuth calculation ...

Source: TerraLab/widgets/sky_widget.py:963-997 Performance: 100k stars transformed in ~5ms (vs ~500ms with Python loops)

Atmospheric Effects

Refraction Correction

Earth’s atmosphere bends light rays, making objects appear higher than their geometric position. The effect is strongest near the horizon. Bennett’s formula (accurate to 0.1’):

R = \frac{1.02'}{\tan\left(h + \frac{10.3'}{h + 5.11'}\right)}

Where:

$R$ = refraction correction (arcminutes)
$h$ = apparent altitude (degrees)

Examples:

h = 90° (zenith): R = 0.0’ (no refraction)
h = 45°: R = 1.0’ (1 arcminute)
h = 10°: R = 5.3’ (noticeable shift)
h = 0° (horizon): R = 34’ (Sun’s diameter!)

Implementation:

h_capped = np.maximum(0.1, alt_deg)  # Avoid singularity
refraction_deg = (1.02 / 60.0) / np.tan(np.radians(
    h_capped + 10.3 / (h_capped + 5.11)
))
alt_refined = alt_deg + refraction_deg

Source: TerraLab/widgets/sky_widget.py:985-987

The Sun/Moon at sunset appear ~0.5° higher than their true position due to refraction. This is why you can see the Sun even when it’s geometrically below the horizon.

Atmospheric Extinction

Stars appear fainter near the horizon due to longer atmospheric path length. Airmass formula (Hardie 1962):

X = \frac{1}{\sin(h) + 0.15 \cdot (h + 3.885)^{-1.253}}

Where:

$X$ = airmass (1.0 at zenith, ~2.0 at 30°, ~38 at horizon)
$h$ = altitude (degrees)

Magnitude penalty:

\Delta m = k_{\text{ext}} \cdot (X - 1)

Where

k_{\text{ext}} = 0.15

–

0.45

mag/airmass (default 0.25). Examples:

h = 90°: X = 1.0, Δm = 0.0 (no extinction)
h = 30°: X = 2.0, Δm = 0.25 mag (star ~22% dimmer)
h = 10°: X = 5.6, Δm = 1.15 mag (star ~66% dimmer)

Source: TerraLab/widgets/sky_widget.py:978-983

Planetary Calculations

Moon Position (ELP 2000-82)

TerraLab includes a fallback lunar ephemeris when Skyfield is unavailable:

class AstroEngine:
    @staticmethod
    def get_moon_position_elp(T):
        """
        ELP 2000-82 truncated series for Moon position.
        
        Args:
            T: Julian centuries from J2000.0
        
        Returns:
            (lon, lat, dist): Ecliptic longitude/latitude (deg), distance (km)
        """
        # Fundamental arguments (mean longitudes/anomalies)
        L_prime = 218.3164477 + 481267.8812542 * T  # Mean longitude
        D = 297.8501921 + 445267.1114034 * T        # Mean elongation
        M = 357.5291092 + 35999.0502909 * T         # Sun mean anomaly
        M_prime = 134.9633964 + 477198.8675055 * T  # Moon mean anomaly
        F = 93.2720950 + 483202.0175381 * T         # Argument of latitude
        
        # Periodic terms (truncated to major terms)
        Σl = 6288774 * sin(2*D - M_prime)   # Evection
        Σl -= 1274027 * sin(2*D - 2*M_prime) # Variation
        # ... (7 major terms kept)
        
        lon = L_prime + Σl / 1e6  # Convert from microarcsec
        # ... latitude and distance ...
        
        return lon, lat, dist

Source: TerraLab/widgets/sky_widget.py:107-189 Accuracy: ±2 arcminutes (sufficient for visual rendering, but not for eclipse predictions)

Topocentric Parallax

The Moon is close enough that observer position on Earth matters. Topocentric correction shifts the Moon’s position by up to 1° from its geocentric position.

\Delta RA = \arctan\left(\frac{-\rho \cos(\phi') \cdot \sin(\pi) \cdot \sin(HA)}{\cos(\delta) - \rho \cos(\phi') \cdot \sin(\pi) \cdot \cos(HA)}\right)

Where:

$\rho$ = observer’s distance from Earth center (≈1.0 for sea level)
$\phi'$ = geocentric latitude
$\pi = \arcsin(R_{\oplus} / d)$ = horizontal parallax
$HA$ = hour angle

Source: TerraLab/widgets/sky_widget.py:235-276

Parallax is critical for lunar eclipse predictions. Without correction, predicted eclipse times can be off by several minutes.

Stereographic Projection

TerraLab renders the sky using stereographic projection centered on the horizon:

k = \frac{2}{1 + \cos(\text{Alt}) \cdot \cos(\text{Az}_{\text{rel}})}

x = k \cdot \cos(\text{Alt}) \cdot \sin(\text{Az}_{\text{rel}})

y = k \cdot \sin(\text{Alt})

Where

\text{Az}_{\text{rel}} = \text{Az}_{\text{star}} - \text{Az}_{\text{camera}}

Properties:

Preserves angles (conformal)
Horizon → circle at infinity
Zenith → origin
North/South poles preserved

Source: TerraLab/widgets/sky_widget.py:1036-1074

Visual Magnitude Engine

TerraLab computes effective magnitude accounting for:

Intrinsic magnitude: Catalog G mag
Atmospheric extinction: $\Delta m_{\text{ext}} = k \cdot (X - 1)$
Light pollution: $\Delta m_{\text{LP}}$ from local Bortle
Light domes: $\Delta m_{\text{dome}}(\text{Az}, \text{Alt})$ directional

m_{\text{eff}} = m_{\text{cat}} + \Delta m_{\text{ext}} + \Delta m_{\text{LP}} + \Delta m_{\text{dome}}

Stars with

m_{\text{eff}} > m_{\text{lim}}

are culled. Source: TerraLab/widgets/visual_magnitude_engine.py

Rendering Optimizations

Asynchronous Star Rendering

Stars are rendered in a background thread to avoid UI freezing:

class StarRenderWorker(QObject):
    result_ready = pyqtSignal(QImage, list)
    
    @pyqtSlot(dict)
    def render(self, params):
        # Create QImage
        img = QImage(width, height, QImage.Format_ARGB32_Premultiplied)
        painter = QPainter(img)
        
        # ... coordinate transformations ...
        # ... magnitude filtering ...
        # ... draw stars with bloom ...
        
        painter.end()
        self.result_ready.emit(img, visible_stars)

Source: TerraLab/widgets/sky_widget.py:893-1246

Realistic Star Rendering

Bright stars use radial gradient bloom + diffraction spikes:

if mag < 4.0:  # Bright star
    # Halo (bloom)
    halo_size = core_radius * (5.5 - mag) * 1.8 * zoom
    grad = QRadialGradient(x, y, halo_size)
    grad.setColorAt(0.0, QColor(r, g, b, 90))
    grad.setColorAt(0.3, QColor(r, g, b, 5))
    grad.setColorAt(1.0, Qt.transparent)
    painter.drawEllipse(x, y, halo_size)
    
    # Core
    core_grad = QRadialGradient(x, y, core_radius)
    core_grad.setColorAt(0.0, Qt.white)  # White center
    core_grad.setColorAt(0.6, QColor(r, g, b, 255))
    painter.drawEllipse(x, y, core_radius)

Source: TerraLab/widgets/sky_widget.py:1166-1188

Diffraction Spikes

For stars brighter than mag 2.0:

spike_length = (2.0 - mag) * 20 * zoom
spike_width = (2.0 - mag) * 0.4 * zoom

for angle in [0°, 90°, 180°, 270°]:
    # Draw diamond-shaped ray with gradient
    grad = QLinearGradient(0, 0, spike_length, 0)
    grad.setColorAt(0.0, QColor(r, g, b, 150))
    grad.setColorAt(1.0, Qt.transparent)
    painter.drawPath(diamond_path)

Source: TerraLab/widgets/sky_widget.py:1193-1233

Precision Validation

Test Cases

TerraLab’s astronomical accuracy has been validated against:

Stellarium: Open-source planetarium (VSOP87/ELP2000)
SkySafari: Professional mobile app
JPL Horizons: NASA’s online ephemeris system

Typical Errors

Object	TerraLab Error	Acceptable Limit
Bright stars	Under 0.1°	0.5° (naked eye resolution)
Planets (Skyfield)	Under 1”	10” (telescopic)
Moon (Skyfield)	Under 30”	1’ (visual)
Moon (ELP fallback)	Under 2’	5’ (acceptable)

Next Steps

Architecture - How astronomical rendering fits into the pipeline
Horizon Engine - Terrain occlusion for celestial objects
Light Pollution - Computing limiting magnitude

Getting Started

Core Concepts

User Guide

CLI Tools

Astronomical Rendering

Skyfield Integration

What is Skyfield?

Installation

DE421 Ephemeris

Gaia Star Catalog

Gaia Mission Overview

TerraLab’s Gaia Subset

Catalog Loading (Optimized)

Color Calculation from BP-RP

Coordinate Transformations

Equatorial to Horizontal

Step 1: Calculate Hour Angle

Step 2: Compute Altitude

Step 3: Compute Azimuth

NumPy Vectorization

Atmospheric Effects

Refraction Correction

Atmospheric Extinction

Planetary Calculations

Moon Position (ELP 2000-82)

Topocentric Parallax

Stereographic Projection

Visual Magnitude Engine

Rendering Optimizations

Asynchronous Star Rendering

Realistic Star Rendering

Diffraction Spikes

Precision Validation

Test Cases

Typical Errors

Next Steps

Build docs developers (and LLMs) love

Getting Started

Core Concepts

User Guide

CLI Tools

​Skyfield Integration

​What is Skyfield?

​Installation

​DE421 Ephemeris

​Gaia Star Catalog

​Gaia Mission Overview

​TerraLab’s Gaia Subset

​Catalog Loading (Optimized)

​Color Calculation from BP-RP

​Coordinate Transformations

​Equatorial to Horizontal

​Step 1: Calculate Hour Angle

​Step 2: Compute Altitude

​Step 3: Compute Azimuth

​NumPy Vectorization

​Atmospheric Effects

​Refraction Correction

​Atmospheric Extinction

​Planetary Calculations

​Moon Position (ELP 2000-82)

​Topocentric Parallax

​Stereographic Projection

​Visual Magnitude Engine

​Rendering Optimizations

​Asynchronous Star Rendering

​Realistic Star Rendering

​Diffraction Spikes

​Precision Validation

​Test Cases

​Typical Errors

​Next Steps

Build docs developers (and LLMs) love

Skyfield Integration

What is Skyfield?

Installation

DE421 Ephemeris

Gaia Star Catalog

Gaia Mission Overview

TerraLab’s Gaia Subset

Catalog Loading (Optimized)

Color Calculation from BP-RP

Coordinate Transformations

Equatorial to Horizontal

Step 1: Calculate Hour Angle

Step 2: Compute Altitude

Step 3: Compute Azimuth

NumPy Vectorization

Atmospheric Effects

Refraction Correction

Atmospheric Extinction

Planetary Calculations

Moon Position (ELP 2000-82)

Topocentric Parallax

Stereographic Projection

Visual Magnitude Engine

Rendering Optimizations

Asynchronous Star Rendering

Realistic Star Rendering

Diffraction Spikes

Precision Validation

Test Cases

Typical Errors

Next Steps