Observation Pipeline

When a user calls planet_observe(5, '40.0N 105.3W 1655m'::observer, now()) to find Jupiter’s position in the sky, seven coordinate transformations execute in sequence. Each stage consumes the output of the previous stage and produces input for the next. Understanding this pipeline is a prerequisite for modifying any part of it.

The full pipeline

flowchart TD
    A["VSOP87: Target heliocentric<br/>ecliptic J2000 (AU)"] --> C["Geocentric ecliptic<br/>target - Earth"]
    B["VSOP87: Earth heliocentric<br/>ecliptic J2000 (AU)"] --> C
    C --> D["Ecliptic → Equatorial J2000<br/>obliquity rotation"]
    D --> E["RA/Dec J2000<br/>Cartesian → spherical"]
    E --> F["Precession J2000 → date<br/>IAU 1976"]
    F --> G["Sidereal time → hour angle<br/>GMST (Vallado Eq. 3-47)"]
    G --> H["Equatorial → Horizontal<br/>az/el for observer"]
    H --> I["pg_topocentric result<br/>(az, el, range, range_rate)"]

Stage-by-stage breakdown

1. Heliocentric ecliptic position of the target

   VSOP87 (Bretagnon & Francou, 1988) computes the target planet’s position in the heliocentric ecliptic J2000 frame. The output is three Cartesian coordinates in AU.

   VSOP87 is a semi-analytical theory: it expands each coordinate as a sum of trigonometric series with polynomial time arguments. The truncated series used in pg_orbit provides ~1 arcsecond accuracy for the inner planets and ~1-2 arcseconds for the outer planets over the period 2000 BCE to 6000 CE.

   For the Sun, this stage returns $(0, 0, 0)$ --- the Sun is at the origin of heliocentric coordinates. The Sun’s apparent position is computed by inverting Earth’s heliocentric position.

   For the Moon, this stage is replaced by ELP2000-82B (Chapront-Touze & Chapront, 1988), which computes geocentric ecliptic coordinates directly, skipping stage 2.

   Code: src/vsop87.c:GetVsop87Coor(), src/elp82b.c:GetElp82bCoor()

2. Heliocentric ecliptic position of Earth

   The same VSOP87 call, but for body ID 3 (Earth). This gives Earth’s position in the same frame as the target. Without this, there is no way to compute where the target appears from Earth’s perspective.

   Code: src/vsop87.c:GetVsop87Coor() with body = 2 (VSOP87 uses 0-indexed: Venus=1, Earth=2, Mars=3, etc.)

3. Geocentric ecliptic vector

   Subtract Earth’s heliocentric position from the target’s:

   $$\vec{r}_\text{geo} = \vec{r}_\text{target} - \vec{r}_\text{Earth}$$

   The result is the target’s position as seen from Earth’s center, still in the ecliptic J2000 frame. For the Moon, ELP2000-82B provides this directly.

   Code: src/planet_funcs.c:planet_observe(), lines computing geo_ecl_au[3]

4. Ecliptic to equatorial rotation

   The ecliptic and equatorial planes are tilted relative to each other by the obliquity of the ecliptic. At J2000, this angle is:

   $$\varepsilon_0 = 23.4392911\degree = 0.40909280422232897\ \text{rad}$$

   The rotation is around the X-axis (the vernal equinox direction, which is shared between both frames):

   $$\begin{bmatrix} x_\text{equ} \\ y_\text{equ} \\ z_\text{equ} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos\varepsilon_0 & -\sin\varepsilon_0 \\ 0 & \sin\varepsilon_0 & \cos\varepsilon_0 \end{bmatrix} \begin{bmatrix} x_\text{ecl} \\ y_\text{ecl} \\ z_\text{ecl} \end{bmatrix}$$

   This uses the J2000 obliquity, not the obliquity of date, because both VSOP87 and the subsequent precession model are referenced to J2000.

   Code: src/planet_funcs.c:ecliptic_to_equatorial() (via astro_math.h)

5. Precession from J2000 to date

   The equatorial coordinate system itself rotates slowly due to lunisolar and planetary precession. Right ascension and declination at J2000 must be precessed to the epoch of observation.

   pg_orbit uses the IAU 1976 precession model (Lieske et al., 1977), which expresses the three Euler angles $\zeta_A$, $z_A$, and $\theta_A$ as cubic polynomials in centuries from J2000:

   $$\zeta_A = 0\overset{\prime\prime}{.}6406161 \cdot T + 0\overset{\prime\prime}{.}0000839 \cdot T^2 + 0\overset{\prime\prime}{.}0000050 \cdot T^3$$

   The full rotation matrix $R_3(-z_A) \cdot R_2(\theta_A) \cdot R_3(-\zeta_A)$ transforms J2000 equatorial to equatorial of date.

   Code: src/precession.c:precess_j2000_to_date()

6. Sidereal time and hour angle

   To convert from the celestial sphere to the observer’s local sky, we need the observer’s relationship to the vernal equinox. This requires two quantities:

   • GMST (Greenwich Mean Sidereal Time): the hour angle of the vernal equinox at Greenwich, computed from the Julian date using Vallado Eq. 3-47 (IAU 1982 model).
   • LST (Local Sidereal Time): $\text{LST} = \text{GMST} + \lambda_\text{observer}$, where $\lambda$ is the observer’s east longitude in radians.
   • Hour angle: $h = \text{LST} - \alpha$, where $\alpha$ is the right ascension of date from stage 5.

   The GMST formula:

   $$\text{GMST}_\text{sec} = 67310.54841 + (876600 \times 3600 + 8640184.812866) \cdot T + 0.093104 \cdot T^2 - 6.2 \times 10^{-6} \cdot T^3$$

   where $T = (JD - 2451545.0) / 36525.0$. The result in seconds is converted to radians: $\text{GMST}_\text{rad} = \text{GMST}_\text{sec} \times \pi / 43200$.

   Code: src/sidereal_time.c:gmst_from_jd()

7. Equatorial to horizontal

   Given hour angle $h$, declination $\delta$, and observer latitude $\phi$, compute azimuth $A$ and elevation $a$:

   $$\sin a = \sin\phi \sin\delta + \cos\phi \cos\delta \cos h$$ $$\tan A = \frac{-\sin h}{\cos\phi \tan\delta - \sin\phi \cos h}$$

   Azimuth is measured from north through east (0 = north, 90 = east, 180 = south, 270 = west).

   The result is packed into a pg_topocentric struct with range computed from the geocentric distance ($\text{range}_\text{km} = d_\text{AU} \times 149597870.7$). Range rate is set to 0.0 for planetary observations (velocity computation is not yet implemented for the VSOP87 pipeline).

   Code: src/planet_funcs.c:observe_from_geocentric()

Pipeline variants

The seven-stage pipeline applies to planets observed via VSOP87. Other observation targets use modified versions.

Moon (ELP2000-82B)

Stages 1-3 are replaced by a single call to GetElp82bCoor(), which returns geocentric ecliptic coordinates directly. The remaining stages (4-7) are identical.

ELP2000-82B is a semi-analytical lunar theory with ~10 arcsecond accuracy. It accounts for the principal perturbations from the Sun, but not the full set of planetary perturbations included in the DE series.

Planetary moons

Each moon theory (L1.2 for Galilean moons, TASS17 for Saturn, GUST86 for Uranus, MarsSat for Mars) computes the moon’s position relative to its parent planet. The pipeline:

1. Compute parent planet heliocentric position via VSOP87
2. Compute moon position relative to parent in parent-equatorial frame
3. Transform to heliocentric ecliptic J2000 (parent offset + frame rotation)
4. Proceed from stage 3 of the standard pipeline (geocentric ecliptic)

Satellites (SGP4)

Satellites use a fundamentally different pipeline. SGP4 outputs TEME (True Equator, Mean Equinox) positions, not heliocentric ecliptic. The satellite pipeline:

1. SGP4/SDP4 propagation to TEME position/velocity (WGS-72)
2. GMST rotation to ECEF
3. ECEF to geodetic (WGS-84) or topocentric via SEZ transform

Stars

Stars are effectively at infinite distance. The pipeline:

1. J2000 catalog coordinates (RA, Dec)
2. Precession to date (IAU 1976)
3. Sidereal time and hour angle
4. Equatorial to horizontal

No VSOP87 call, no geocentric vector, no distance computation. This makes star observation the fastest pipeline --- 714K observations per second.

Why this pipeline and not another

Several simplifications are deliberate.

No nutation

The pipeline uses precession but not nutation. For 1-arcsecond VSOP87 positions, the ~9-arcsecond nutation correction is below the noise floor of the ephemeris. Adding nutation would increase computation cost without improving practical accuracy.

No aberration correction. Annual aberration shifts apparent positions by up to 20 arcseconds, but for observation planning (which quadrant of the sky is Jupiter in tonight?) this is irrelevant. Sub-arcsecond work should use SPICE or Skyfield with DE441.

No light-time iteration. The positions returned are geometric, not apparent. Light-time corrections of a few minutes for the outer planets shift the apparent position by a fraction of an arcsecond at most --- again, below the VSOP87 accuracy floor.

No atmospheric refraction. Refraction near the horizon can shift apparent elevation by half a degree. pg_orbit reports geometric elevation; the user must apply refraction corrections for their local conditions if needed. This is a deliberate choice --- refraction depends on temperature, pressure, and humidity that pg_orbit does not model.

Extending the pipeline

To add a new observation target, identify which stages change:

New target	Stages replaced	Stages reused
New moon theory	1-3 (new theory + parent VSOP87)	4-7
New planet ephemeris	1 (new theory replaces VSOP87)	2-7 (Earth still VSOP87)
Near-Earth asteroid	1 (Keplerian propagation)	2-7
Distant star	1-3 (catalog lookup, no distance)	5-7 (skip obliquity)

The modularity of the pipeline means new targets require implementing only the first few stages. The coordinate transformation machinery from stage 4 onward is shared across all targets.