Update docs for v0.10.0: benchmarks, Skyfield parity, LLM references

- Add timing numbers for equatorial, aberration, angular distance,
  refraction, and star proper motion+parallax to benchmarks page
- Update From Skyfield page: v0.10.0 now has light-time + aberration
  parity; remaining gap narrowed to nutation (~9 arcsec) and polar motion
- Update llms.txt and llms-full.txt for 114 functions, new DE apparent
  variants, equatorial spatial operators, and aberration/parallax notes

docs/public/llms-full.txt

@@ -1,6 +1,6 @@
 # pg_orrery — Complete LLM Reference
 
-> Celestial mechanics types and functions for PostgreSQL. Native C extension (v0.9.0) with 106 SQL functions, 9 custom types + 1 composite, GiST/SP-GiST indexing. All functions PARALLEL SAFE.
+> Celestial mechanics types and functions for PostgreSQL. Native C extension (v0.10.0) with 114 SQL functions, 9 custom types + 1 composite, GiST/SP-GiST indexing. All functions PARALLEL SAFE.
 
 - Source: https://git.supported.systems/warehack.ing/pg_orrery
 - Docs: https://pg-orrery.warehack.ing
@@ -214,6 +214,8 @@ planet_observe_apparent(body_id int4, observer, timestamptz) → topocentric  IM
 sun_observe_apparent(observer, timestamptz) → topocentric              IMMUTABLE
 planet_equatorial_apparent(body_id int4, timestamptz) → equatorial     IMMUTABLE
 moon_equatorial_apparent(timestamptz) → equatorial                     IMMUTABLE
+
+-- All _apparent() functions include annual aberration correction (~20 arcsec) + light-time
 ```
 
 ### Planetary Moons (4 functions)
@@ -235,6 +237,8 @@ star_equatorial(ra_hours, dec_degrees, timestamptz) → equatorial       IMMUTAB
 -- Proper motion (Hipparcos/Gaia convention: pm_ra = mu_alpha * cos(delta) in mas/yr)
 star_observe_pm(ra_h, dec_deg, pm_ra_masyr, pm_dec_masyr, parallax_mas, rv_kms, observer, timestamptz) → topocentric  IMMUTABLE
 star_equatorial_pm(ra_h, dec_deg, pm_ra_masyr, pm_dec_masyr, parallax_mas, rv_kms, timestamptz) → equatorial  IMMUTABLE
+
+-- When parallax_mas > 0, annual stellar parallax is applied using Earth's heliocentric position (Green 1985)
 ```
 
 RA in hours [0,24), Dec in degrees [-90,90]. Range returned as 0 (infinite distance) unless parallax > 0 in _pm variants.
@@ -283,7 +287,16 @@ predict_passes_refracted(tle, observer, start, end, min_el DEFAULT 0.0) → SETO
 
 Bennett formula: `R = 1/tan(h + 7.31/(h + 4.4))` arcminutes. Domain guard: clamps at -1°, returns 0.0 below. At horizon (0°) refraction is ~0.57°, meaning satellites become visible ~35 seconds earlier.
 
-### DE Ephemeris — Optional High-Precision (13 functions)
+### Equatorial Spatial — Angular Separation (2 functions)
+
+```
+eq_angular_distance(equatorial, equatorial) → float8                   IMMUTABLE  -- degrees, Vincenty formula (stable at 0° and 180°)
+eq_within_cone(equatorial, equatorial, float8) → bool                  IMMUTABLE  -- true if within radius_deg, cosine shortcut
+```
+
+Operator: `equatorial <-> equatorial → float8` (angular separation in degrees, commutative).
+
+### DE Ephemeris — Optional High-Precision (19 functions)
 
 All _de() functions fall back to VSOP87/ELP2000-82B when DE is unavailable. All STABLE (external file dependency).
 
@@ -301,6 +314,14 @@ mars_moon_observe_de(moon_id, observer, timestamptz) → topocentric     STABLE
 planet_equatorial_de(body_id int4, timestamptz) → equatorial           STABLE  -- geocentric RA/Dec via DE
 moon_equatorial_de(timestamptz) → equatorial                           STABLE  -- geocentric RA/Dec via DE
 pg_orrery_ephemeris_info() → (provider, file_path, start_jd, end_jd, version, au_km)  STABLE
+
+-- Apparent DE variants (light-time + aberration, falls back to VSOP87)
+planet_observe_apparent_de(body_id int4, observer, timestamptz) → topocentric  STABLE
+sun_observe_apparent_de(observer, timestamptz) → topocentric                   STABLE
+moon_observe_apparent_de(observer, timestamptz) → topocentric                  STABLE
+planet_equatorial_apparent_de(body_id int4, timestamptz) → equatorial          STABLE
+moon_equatorial_apparent_de(timestamptz) → equatorial                          STABLE
+small_body_observe_apparent_de(orbital_elements, observer, timestamptz) → topocentric  STABLE
 ```
 
 Configure: `ALTER SYSTEM SET pg_orrery.ephemeris_path = '/path/to/de441.bin'; SELECT pg_reload_conf();`
@@ -358,6 +379,12 @@ CREATE INDEX ON satellites USING spgist (elements tle_spgist_ops);
 
 SP-GiST is a 2-level orbital trie (SMA → inclination) with query-time RAAN filter. Returns a conservative superset — survivors need `predict_passes()` for ground truth.
 
+### Equatorial distance
+
+| Operator | Meaning | Usage |
+|----------|---------|-------|
+| `<->` (equatorial) | Angular separation in degrees (Vincenty formula) | `ORDER BY pos1 <-> pos2` or `WHERE pos1 <-> pos2 < 5.0` |
+
 ## Common Query Patterns
 
 ### Observe a satellite
@@ -513,6 +540,21 @@ FROM satellites,
 WHERE name = 'ISS';
 ```
 
+### Angular separation and cone search
+
+```sql
+-- Find angular separation between two objects
+SELECT planet_equatorial(5, NOW()) <-> moon_equatorial(NOW()) AS jupiter_moon_sep_deg;
+
+-- Objects within 10 degrees of Jupiter
+SELECT eq_within_cone(
+    star_equatorial(ra_h, dec_deg, NOW()),
+    planet_equatorial(5, NOW()),
+    10.0
+) AS near_jupiter
+FROM star_catalog;
+```
+
 ## Error Handling
 
 ### _safe() variants

docs/public/llms.txt

@@ -1,6 +1,6 @@
 # pg_orrery
 
-> Celestial mechanics types and functions for PostgreSQL. Native C extension with 106 SQL functions, 9 custom types, GiST/SP-GiST indexing. Covers satellites (SGP4/SDP4), planets (VSOP87), Moon (ELP2000-82B), 19 planetary moons, stars (with proper motion), comets, asteroids (MPC catalog), Jupiter radio bursts, orbit determination, interplanetary Lambert transfers, equatorial RA/Dec coordinates, atmospheric refraction, and light-time correction. Optional JPL DE440/441 ephemeris for sub-arcsecond accuracy.
+> Celestial mechanics types and functions for PostgreSQL. Native C extension with 114 SQL functions, 9 custom types, GiST/SP-GiST indexing. Covers satellites (SGP4/SDP4), planets (VSOP87), Moon (ELP2000-82B), 19 planetary moons, stars (with proper motion), comets, asteroids (MPC catalog), Jupiter radio bursts, orbit determination, interplanetary Lambert transfers, equatorial RA/Dec coordinates, atmospheric refraction, light-time correction, annual stellar aberration, and equatorial angular separation. Optional JPL DE440/441 ephemeris for sub-arcsecond accuracy.
 
 - [Source code](https://git.supported.systems/warehack.ing/pg_orrery)
 - [Full LLM reference](https://pg-orrery.warehack.ing/llms-full.txt): All function signatures, types, body IDs, operators, and query patterns inline
@@ -47,7 +47,8 @@
 - [Functions: Radio](https://pg-orrery.warehack.ing/reference/functions-radio/): Jupiter decametric radio burst prediction — Io phase, CML, burst probability
 - [Functions: Transfers](https://pg-orrery.warehack.ing/reference/functions-transfers/): Lambert transfer solver for interplanetary trajectory design
 - [Functions: Refraction](https://pg-orrery.warehack.ing/reference/functions-refraction/): Bennett (1982) atmospheric refraction, P/T correction, apparent elevation, refracted pass prediction
-- [Functions: DE Ephemeris](https://pg-orrery.warehack.ing/reference/functions-de/): Optional JPL DE440/441 variants of observation and equatorial functions
+- [Functions: Equatorial Spatial](https://pg-orrery.warehack.ing/reference/functions-equatorial/): Angular separation (Vincenty formula), cone search, `<->` operator on equatorial type
+- [Functions: DE Ephemeris](https://pg-orrery.warehack.ing/reference/functions-de/): Optional JPL DE440/441 variants of observation, equatorial, and apparent functions
 - [Functions: Orbit Determination](https://pg-orrery.warehack.ing/reference/functions-od/): TLE fitting from ECI, topocentric, and angles-only observations
 - [Operators & Indexes](https://pg-orrery.warehack.ing/reference/operators-gist/): GiST (&&, <->) and SP-GiST (&?) operator classes for orbital indexing
 - [Body ID Reference](https://pg-orrery.warehack.ing/reference/body-ids/): Planet IDs 0–10, Galilean 0–3, Saturn 0–7, Uranus 0–4, Mars 0–1

docs/src/content/docs/performance/benchmarks.mdx

@@ -24,6 +24,10 @@ All benchmarks use PostgreSQL's `EXPLAIN (ANALYZE, BUFFERS)` for timing. The num
 | Galilean moon observation | 1,000 | 63 ms | 15.9K/sec | L1.2 + VSOP87 pipeline |
 | Saturn moon observation | 800 | 53 ms | 15.1K/sec | TASS17 + VSOP87 |
 | Star observation | 500 | 0.7 ms | 714K/sec | Precession + az/el only |
+| Star with proper motion + parallax | 8,761 | 343 ms | 25.5K/sec | Proper motion + annual parallax (VSOP87 Earth) |
+| Planet equatorial + apparent | 1,414 | 107 ms | 13.2K/sec | RA/Dec geometric + light-time + aberration |
+| Angular distance (`<->` equatorial) | 10,000 | 7 ms | 1.43M/sec | Vincenty formula on equatorial pairs |
+| Atmospheric refraction | 18,200 | 3.5 ms | 5.2M/sec | Bennett (1982) + P/T correction |
 | Lambert transfer solve | 100 | 0.1 ms | 800K/sec | Single-rev prograde |
 | Pork chop plot (150 x 150) | 22,500 | 8.3 s | 2.7K/sec | Full VSOP87 + Lambert pipeline |
 
@@ -145,6 +149,87 @@ LIMIT 500;
 
 This is nearly as fast as SGP4 propagation because the only computation is matrix multiplication (precession) and a trigonometric transform (az/el). No series evaluation, no iteration.
 
+## Star observation with proper motion and parallax
+
+The full stellar pipeline: proper motion correction (RA/Dec rates), annual parallax displacement from VSOP87 Earth position, IAU 1976 precession, sidereal time, and equatorial output.
+
+```sql
+-- Benchmark: proper motion + parallax for stars across a year of epochs
+EXPLAIN (ANALYZE)
+SELECT star_equatorial_pm(
+    mod(a * 1.5, 24.0), mod(a * 0.7, 180.0) - 90.0,
+    a * 0.1, a * 0.05, a * 10.0, a * 0.01,
+    '2024-01-01'::timestamptz + (a || ' hours')::interval
+)
+FROM generate_series(1, 8761) AS a;
+```
+
+**8,761 evaluations in 343 ms --- 25,500 per second.**
+
+When `parallax_mas > 0`, each call adds a `GetVsop87Coor()` evaluation for Earth's heliocentric position to compute the annual parallax displacement. This makes it ~28x slower than basic star observation (which skips VSOP87 entirely). The cost is dominated by the single VSOP87 Earth call per star --- the parallax displacement math itself is negligible.
+
+For catalogs where parallax is zero or unknown (most survey catalogs), `star_equatorial()` without proper motion runs at the full 714K/sec rate.
+
+## Planet equatorial and apparent positions
+
+Equatorial output (RA/Dec) for planets, including the `_apparent()` pipeline that adds light-time correction and annual stellar aberration.
+
+```sql
+-- Benchmark: equatorial + apparent for all non-Earth planets
+EXPLAIN (ANALYZE)
+SELECT planet_equatorial(body_id, t),
+       planet_equatorial_apparent(body_id, t)
+FROM generate_series(1, 8) AS body_id,
+     generate_series(
+       '2024-01-01'::timestamptz,
+       '2024-01-01'::timestamptz + interval '100 hours',
+       interval '1 hour'
+     ) AS t
+WHERE body_id != 3;
+```
+
+**1,414 evaluations (707 geometric + 707 apparent) in 107 ms --- 13,200 per second.**
+
+The `_apparent()` variant adds two operations beyond the geometric pipeline: iterating on light-travel time (1--2 iterations to converge) and applying annual stellar aberration from Earth's VSOP87 velocity vector. The aberration calculation itself is cheap (~20 ns) but the extra VSOP87 evaluation for the retarded position roughly doubles the cost per call.
+
+## Angular distance (`<->` on equatorial)
+
+The `<->` operator on the `equatorial` type computes angular separation in degrees using the Vincenty formula, which is numerically stable at all separations including 0° and 180°.
+
+```sql
+-- Benchmark: angular distance between 10,000 star pairs
+EXPLAIN (ANALYZE)
+SELECT eq_angular_distance(
+    star_equatorial(mod(a * 1.3, 24.0), mod(a * 0.7, 180.0) - 90.0, now()),
+    star_equatorial(mod(a * 2.1, 24.0), mod(a * 0.9, 180.0) - 90.0, now())
+)
+FROM generate_series(1, 10000) AS a;
+```
+
+**10,000 angular distances in 7 ms --- 1.43 million per second.**
+
+The Vincenty formula involves four trigonometric evaluations and an `atan2` --- comparable cost to a single coordinate transform. The `eq_within_cone()` predicate is even faster because it uses a cosine shortcut (`cos(dist) >= cos(radius)`) that avoids the `atan2`.
+
+<Aside type="tip" title="Future: GiST index for cone search">
+Currently `eq_angular_distance()` and `eq_within_cone()` require sequential evaluation. A future GiST or SP-GiST index on the `equatorial` type would enable indexed cone search, turning "what's within 10° of Jupiter?" into an index scan rather than a full table scan over a star catalog.
+</Aside>
+
+## Atmospheric refraction
+
+Bennett's (1982) empirical formula for atmospheric refraction, with optional pressure/temperature correction.
+
+```sql
+-- Benchmark: refraction + P/T-corrected refraction for elevation range
+EXPLAIN (ANALYZE)
+SELECT atmospheric_refraction(a * 0.01),
+       atmospheric_refraction_ext(a * 0.01, 1013.25, 15.0)
+FROM generate_series(1, 9100) AS a;
+```
+
+**18,200 evaluations (9,100 base + 9,100 extended) in 3.5 ms --- 5.2 million per second.**
+
+Refraction is pure arithmetic --- a single `cot()` approximation with a polynomial correction term. No series evaluation, no iteration. The `_ext()` variant adds a pressure/temperature scaling factor (one division). This makes refraction essentially free when layered onto the observation pipeline.
+
 ## Lambert transfer
 
 A single Lambert solve: given two planet positions and a time of flight, find the transfer orbit.
@@ -412,6 +497,8 @@ The CTE pattern works correctly but forces PostgreSQL to compute all distances a
 
 The benchmarks demonstrate that pg_orrery's computation cost is low enough to treat orbital mechanics as a SQL primitive. Propagating an entire satellite catalog takes less time than a typical index scan on a moderately-sized table. Planet observation is fast enough to generate ephemeris tables with `generate_series`. Pork chop plots are feasible as interactive queries rather than batch jobs.
 
+The v0.10.0 additions --- aberration, angular distance, refraction --- range from negligible overhead (refraction at 5.2M/sec) to moderate (apparent positions at 13.2K/sec, roughly half the geometric rate due to the extra VSOP87 call for light-time iteration). Angular distance at 1.43M/sec means cone-search predicates over star catalogs are fast even without index support.
+
 The visibility cone filter (`&?`) is the fastest operation per evaluation --- three floating-point comparisons vs. the full SGP4 pipeline --- and its 84--90% pruning rate means the most expensive operation in a pass prediction pipeline (SGP4 propagation) only runs on the small fraction of the catalog that could actually produce a visible pass.
 
 The GiST index provides the clearest speedup for conjunction screening: 5.8x faster than sequential scan for ISS `&&` queries, with 0 false positives or negatives verified against exhaustive sequential evaluation. KNN queries find the nearest orbits in 2 ms via index-ordered traversal using 2-D orbital distance (altitude + inclination), which would otherwise require computing and sorting all 66,440 distances.

docs/src/content/docs/workflow/from-skyfield.mdx

@@ -273,11 +273,11 @@ Predict when a satellite will be visible from a location. This is where Skyfield
 pg_orrery does not replace Skyfield for all use cases. Be clear about where the trade-offs fall.
 </Aside>
 
-**Apparent-position corrections.** Skyfield uses the full IAU 2000A nutation model, polar motion corrections, delta-T from IERS data, and iterates for light-time and stellar aberration. Since v0.3.0, pg_orrery can [optionally use DE441](/guides/de-ephemeris/) for the same underlying geometric accuracy (~0.1 milliarcsecond), but Skyfield still applies corrections that pg_orrery does not — corrections that matter for precision apparent-coordinate work like occultation timing or sub-arcsecond astrometry.
+**Apparent-position corrections.** Since v0.9.0, pg_orrery applies light-time correction, and since v0.10.0, annual stellar aberration (~20 arcsec) in all `_apparent()` functions — closing the two largest gaps with Skyfield. However, Skyfield still uses the full IAU 2000A nutation model (~9 arcsec more accurate than pg_orrery's IAU 1976 precession), polar motion corrections, and delta-T from IERS data. These matter for precision apparent-coordinate work like occultation timing or sub-arcsecond astrometry.
 
 **Rise/set finding.** `find_events()` uses numerical root-finding to pinpoint the exact moment a body crosses an elevation threshold. pg_orrery's `predict_passes` uses a step-and-refine approach that's fast for batches but less precise for individual events.
 
-**Aberration and light-time.** Skyfield iterates to correct for light travel time and applies stellar aberration. pg_orrery uses geometric positions without light-time iteration — the difference is under 20 arcseconds for planets and irrelevant for satellite tracking, but it matters for some applications.
+**Nutation and polar motion.** pg_orrery uses IAU 1976 precession, which introduces a ~9 arcsecond gap compared to Skyfield's IAU 2000A nutation model. Skyfield also applies polar motion corrections from IERS data. For work requiring sub-arcsecond accuracy (occultation timing, precise astrometric reductions), this gap still matters.
 
 **Visualization integration.** Skyfield works directly with matplotlib, numpy, and the rest of the Python scientific stack. pg_orrery produces rows and columns — you export to CSV or JSON and then plot separately.
 
@@ -336,7 +336,7 @@ FROM fit,
 You don't have to choose one or the other. A practical migration path:
 
 <Steps>
-1. **Keep Skyfield for apparent-position work.** Anything requiring aberration corrections, polar motion, nutation at IAU 2000A level, or custom BSP kernels stays in Python. For raw geometric position accuracy, pg_orrery with [DE enabled](/guides/de-ephemeris/) matches Skyfield.
+1. **Keep Skyfield for high-precision apparent positions.** Anything requiring polar motion, nutation at IAU 2000A level (~9 arcsec beyond pg_orrery's IAU 1976 precession), or custom BSP kernels stays in Python. For geometric accuracy and apparent-position work including light-time and aberration, pg_orrery v0.10.0 with [DE enabled](/guides/de-ephemeris/) is now competitive.
 
 2. **Move batch observation to SQL.** If you're computing positions for hundreds of objects to filter or correlate with database records, pg_orrery eliminates the Python-to-PostgreSQL round trip.