Compare commits
No commits in common. "3ea4e2bf53c8ac42ce24f1f05c105e233ee9efb2" and "6e17513885512845ea05c9c6676cd73ba8e7fee0" have entirely different histories.
3ea4e2bf53
...
6e17513885
1
.gitignore
vendored
1
.gitignore
vendored
@ -23,7 +23,6 @@ test/matrix-logs/
|
||||
test/test_de_reader
|
||||
test/test_od_math
|
||||
test/test_od_iod
|
||||
test/test_od_gauss
|
||||
|
||||
# Docs site
|
||||
docs/node_modules/
|
||||
|
||||
@ -47,7 +47,6 @@ RUN su postgres -c "/usr/lib/postgresql/${PG_MAJOR}/bin/initdb -D /tmp/pgtest" &
|
||||
RUN make test-de-reader
|
||||
RUN make test-od-math
|
||||
RUN make test-od-iod
|
||||
RUN make test-od-gauss
|
||||
|
||||
# Capture artifacts under /pg_orrery prefix for the next stage
|
||||
RUN make PG_CONFIG=${PG_CONFIG} DESTDIR=/pg_orrery install
|
||||
|
||||
11
Makefile
11
Makefile
@ -3,8 +3,7 @@ EXTENSION = pg_orrery
|
||||
DATA = sql/pg_orrery--0.1.0.sql sql/pg_orrery--0.2.0.sql sql/pg_orrery--0.1.0--0.2.0.sql \
|
||||
sql/pg_orrery--0.3.0.sql sql/pg_orrery--0.2.0--0.3.0.sql \
|
||||
sql/pg_orrery--0.4.0.sql sql/pg_orrery--0.3.0--0.4.0.sql \
|
||||
sql/pg_orrery--0.5.0.sql sql/pg_orrery--0.4.0--0.5.0.sql \
|
||||
sql/pg_orrery--0.6.0.sql sql/pg_orrery--0.5.0--0.6.0.sql
|
||||
sql/pg_orrery--0.5.0.sql sql/pg_orrery--0.4.0--0.5.0.sql
|
||||
|
||||
# Our extension C sources
|
||||
OBJS = src/pg_orrery.o src/tle_type.o src/eci_type.o src/observer_type.o \
|
||||
@ -70,14 +69,6 @@ test-od-iod: test/test_od_iod.c src/od_iod.c src/od_iod.h src/od_math.c src/od_m
|
||||
|
||||
.PHONY: test-od-iod
|
||||
|
||||
# ── Standalone Gauss IOD unit test (no PostgreSQL dependency) ──
|
||||
# Gauss angles-only IOD, RA/Dec round-trip, Herrick-Gibbs fallback.
|
||||
test-od-gauss: test/test_od_gauss.c src/od_iod.c src/od_iod.h src/od_math.c src/od_math.h
|
||||
$(CC) -Wall -Werror -Isrc -o test/test_od_gauss $< src/od_iod.c src/od_math.c -lm
|
||||
./test/test_od_gauss
|
||||
|
||||
.PHONY: test-od-gauss
|
||||
|
||||
# ── PG version test matrix ─────────────────────────────────
|
||||
PG_TEST_VERSIONS ?= 14 15 16 17 18
|
||||
|
||||
|
||||
6
TODO
6
TODO
@ -1 +1,7 @@
|
||||
|
||||
- Gauss method for angles-only initial orbit determination
|
||||
(eliminates seed requirement for sensors without ranging)
|
||||
- Weighted observations (per-obs covariance weighting for
|
||||
heterogeneous sensor fusion)
|
||||
- Range rate fitting in topocentric mode (currently reserved
|
||||
via vel_ecef in residual computation)
|
||||
|
||||
@ -68,7 +68,6 @@ export default defineConfig({
|
||||
{ label: "Interplanetary Trajectories", slug: "guides/interplanetary-trajectories" },
|
||||
{ label: "Conjunction Screening", slug: "guides/conjunction-screening" },
|
||||
{ label: "JPL DE Ephemeris", slug: "guides/de-ephemeris" },
|
||||
{ label: "Orbit Determination", slug: "guides/orbit-determination" },
|
||||
],
|
||||
},
|
||||
{
|
||||
@ -78,8 +77,6 @@ export default defineConfig({
|
||||
{ label: "From JPL Horizons to SQL", slug: "workflow/from-jpl-horizons" },
|
||||
{ label: "From GMAT to SQL", slug: "workflow/from-gmat" },
|
||||
{ label: "From Radio Jupiter Pro to SQL", slug: "workflow/from-radio-jupiter-pro" },
|
||||
{ label: "From find_orb to SQL", slug: "workflow/from-find-orb" },
|
||||
{ label: "From Poliastro to SQL", slug: "workflow/from-poliastro" },
|
||||
{ label: "The SQL Advantage", slug: "workflow/sql-advantage" },
|
||||
],
|
||||
},
|
||||
@ -94,7 +91,6 @@ export default defineConfig({
|
||||
{ label: "Functions: Radio", slug: "reference/functions-radio" },
|
||||
{ label: "Functions: Transfers", slug: "reference/functions-transfers" },
|
||||
{ label: "Functions: DE Ephemeris", slug: "reference/functions-de" },
|
||||
{ label: "Functions: Orbit Determination", slug: "reference/functions-od" },
|
||||
{ label: "Operators & GiST Index", slug: "reference/operators-gist" },
|
||||
{ label: "Body ID Reference", slug: "reference/body-ids" },
|
||||
{ label: "Constants & Accuracy", slug: "reference/constants-accuracy" },
|
||||
|
||||
@ -103,21 +103,6 @@ The 82B revision is the version implemented. It provides geocentric ecliptic coo
|
||||
|
||||
The Izzo solver uses Householder iterations for fast convergence and handles both short-way and long-way transfers. pg_orrery uses the prograde (short-way) solution by default.
|
||||
|
||||
## Orbit determination
|
||||
|
||||
| Theory | Source | What it computes | Code location |
|
||||
|--------|--------|------------------|---------------|
|
||||
| Differential correction | Vallado (2013) Ch. 10 | Equinoctial element refinement via SVD least-squares | `src/od_solver.c` |
|
||||
| Gibbs method | Vallado Algorithm 54 | Initial velocity from three position vectors | `src/od_iod.c` |
|
||||
| Herrick-Gibbs | Vallado Algorithm 55 | Short-arc initial velocity (closely-spaced obs) | `src/od_iod.c` |
|
||||
| Gauss method | Vallado Algorithm 52 | Initial orbit from three angles-only (RA/Dec) observations | `src/od_iod.c` |
|
||||
|
||||
### References
|
||||
|
||||
- Vallado, D.A. (2013). *Fundamentals of Astrodynamics and Applications*, 4th ed. Microcosm Press.
|
||||
|
||||
The OD solver uses equinoctial elements to avoid singularities at zero eccentricity and inclination. LAPACK's `dgelss_` provides the SVD solve, with `dpotrf_`/`dpotri_` for formal covariance estimation.
|
||||
|
||||
## Radio emission
|
||||
|
||||
| Theory | Source | What it computes | Code location |
|
||||
|
||||
@ -97,7 +97,7 @@ import { Tabs, TabItem, Steps, Aside } from "@astrojs/starlight/components";
|
||||
|
||||
## Running the test suite
|
||||
|
||||
If building from source, the regression tests verify all functions across 14 test suites:
|
||||
If building from source, the regression tests verify all 68 functions across 12 test suites:
|
||||
|
||||
```bash
|
||||
make installcheck PG_CONFIG=/usr/bin/pg_config
|
||||
@ -115,15 +115,6 @@ ALTER EXTENSION pg_orrery UPDATE TO '0.2.0';
|
||||
|
||||
-- From v0.2.0 to v0.3.0 (DE ephemeris support)
|
||||
ALTER EXTENSION pg_orrery UPDATE TO '0.3.0';
|
||||
|
||||
-- From v0.3.0 to v0.4.0 (orbit determination)
|
||||
ALTER EXTENSION pg_orrery UPDATE TO '0.4.0';
|
||||
|
||||
-- From v0.4.0 to v0.5.0 (OD enhancements: multi-observer, Gibbs IOD, covariance)
|
||||
ALTER EXTENSION pg_orrery UPDATE TO '0.5.0';
|
||||
|
||||
-- From v0.5.0 to v0.6.0 (range rate, weighted obs, Gauss angles-only IOD)
|
||||
ALTER EXTENSION pg_orrery UPDATE TO '0.6.0';
|
||||
```
|
||||
|
||||
Each migration adds new functions while preserving existing data and functions.
|
||||
|
||||
@ -112,6 +112,5 @@ You've seen the five domains pg_orrery covers. For deeper dives:
|
||||
|
||||
- **[Tracking Satellites](/guides/tracking-satellites/)** — batch observation, conjunction screening, pass prediction workflows
|
||||
- **[Observing the Solar System](/guides/observing-solar-system/)** — "what's up tonight?" queries, rise/set times, conjunctions
|
||||
- **[Orbit Determination](/guides/orbit-determination/)** — fit TLEs from ECI positions, ground station observations, or angles-only RA/Dec data
|
||||
- **[JPL DE Ephemeris](/guides/de-ephemeris/)** — opt-in sub-milliarcsecond accuracy using JPL DE440/441 files
|
||||
- **[The SQL Advantage](/workflow/sql-advantage/)** — why doing this in SQL changes what's possible
|
||||
|
||||
@ -26,7 +26,6 @@ PostGIS added spatial awareness to PostgreSQL — suddenly your database underst
|
||||
| Jupiter radio | Carr et al. (1983) sources | `jupiter_burst_probability()` | Empirical probability |
|
||||
| Transfers | Lambert (Izzo, 2015) | `lambert_transfer()`, `lambert_c3()` | Ballistic two-body |
|
||||
| DE ephemeris (optional) | JPL DE440/441 | `planet_observe_de()`, `moon_observe_de()` | ~0.1 milliarcsecond |
|
||||
| Orbit determination | Differential correction (v0.4.0+) | `tle_from_eci()`, `tle_from_topocentric()`, `tle_from_angles()` | Depends on observation quality |
|
||||
|
||||
## Who it's for
|
||||
|
||||
@ -58,7 +57,7 @@ pg_orrery is a computation engine, not a complete application. Understanding wha
|
||||
|
||||
**Not sub-arcsecond by default.** The built-in VSOP87 pipeline is accurate to about 1 arcsecond — sufficient for observation planning and visual astronomy. For precision work (dish pointing, occultation timing, astrometry), pg_orrery v0.3.0 supports [optional JPL DE440/441 ephemeris files](/guides/de-ephemeris/) that bring accuracy to ~0.1 milliarcsecond. DE is opt-in and requires a one-time GUC configuration.
|
||||
|
||||
**Not a TLE source.** Bring your own TLEs from Space-Track, CelesTrak, or any other provider. pg_orrery parses, propagates, and — since v0.4.0 — can [fit new TLEs from observations](/guides/orbit-determination/). But it doesn't fetch TLE catalogs.
|
||||
**Not a TLE source.** Bring your own TLEs from Space-Track, CelesTrak, or any other provider. pg_orrery parses and propagates them; it doesn't fetch them.
|
||||
|
||||
**Not a replacement for SPICE.** No BSP kernel support, no light-time iteration, no aberration corrections at the IAU 2000A level. With DE enabled, pg_orrery matches SPICE on raw planet position accuracy — the remaining gap is in apparent-position corrections (aberration, light-time, nutation) that matter for sub-arcsecond apparent coordinates.
|
||||
|
||||
|
||||
@ -1,226 +0,0 @@
|
||||
---
|
||||
title: Orbit Determination
|
||||
sidebar:
|
||||
order: 10
|
||||
---
|
||||
|
||||
import { Steps, Aside, Tabs, TabItem } from "@astrojs/starlight/components";
|
||||
|
||||
Orbit determination (OD) is the inverse of propagation: instead of computing where a satellite will be given its orbital elements, you compute the orbital elements given where the satellite was observed. pg_orrery's OD solver fits SGP4/SDP4 mean elements (as a TLE) from observations using differential correction with an SVD least-squares solve.
|
||||
|
||||
Three observation types are supported:
|
||||
|
||||
- **ECI** — position/velocity vectors in the TEME frame
|
||||
- **Topocentric** — azimuth, elevation, range (and optionally range rate) from a ground station
|
||||
- **Angles-only** — right ascension and declination (RA/Dec) from optical observations
|
||||
|
||||
## How differential correction works
|
||||
|
||||
The solver starts with an initial orbit estimate (the "seed") and iteratively refines it to minimize the residuals between observed and computed positions. Each iteration:
|
||||
|
||||
1. Propagates the current TLE to each observation time
|
||||
2. Computes residuals (observed minus computed)
|
||||
3. Builds the Jacobian matrix (partial derivatives of residuals with respect to orbital elements)
|
||||
4. Solves the normal equations via SVD to get a correction vector
|
||||
5. Applies the correction to the equinoctial elements and converts back to a TLE
|
||||
|
||||
The solver uses equinoctial elements rather than classical Keplerian elements to avoid singularities at zero eccentricity and zero inclination — both common for real satellites.
|
||||
|
||||
## ECI fitting
|
||||
|
||||
The simplest case: you have ECI position/velocity vectors (perhaps from another propagator, a simulation, or a precise ephemeris) and want to fit a TLE.
|
||||
|
||||
```sql
|
||||
-- Generate synthetic observations by propagating a known TLE
|
||||
WITH iss AS (
|
||||
SELECT '1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025
|
||||
2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS tle
|
||||
),
|
||||
obs AS (
|
||||
SELECT array_agg(sgp4_propagate(iss.tle, t) ORDER BY t) AS positions,
|
||||
array_agg(t ORDER BY t) AS times
|
||||
FROM iss,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
interval '5 minutes') t
|
||||
)
|
||||
SELECT iterations,
|
||||
round(rms_final::numeric, 6) AS rms_km,
|
||||
status,
|
||||
round(condition_number::numeric, 1) AS cond
|
||||
FROM obs, tle_from_eci(obs.positions, obs.times);
|
||||
```
|
||||
|
||||
With clean synthetic data (propagated from the same TLE model), this converges in a few iterations with sub-meter RMS. With real observations containing measurement noise, expect RMS in the 0.1-10 km range depending on data quality.
|
||||
|
||||
<Aside type="tip" title="No seed needed">
|
||||
`tle_from_eci()` can bootstrap its own seed using the Gibbs method (three position vectors to derive an initial velocity). You only need to provide a seed TLE when the Gibbs geometry is degenerate (nearly colinear positions) or when you want to start from a known approximate orbit.
|
||||
</Aside>
|
||||
|
||||
## Topocentric fitting
|
||||
|
||||
When observations come from a ground station (radar tracks, optical telescopes with range data), use `tle_from_topocentric()`. The solver accounts for the observer's position on the rotating Earth when computing the observation geometry.
|
||||
|
||||
```sql
|
||||
-- Observe a satellite, then fit from the topocentric data
|
||||
WITH iss AS (
|
||||
SELECT '1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025
|
||||
2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS tle
|
||||
),
|
||||
obs AS (
|
||||
SELECT array_agg(
|
||||
observe(iss.tle, '40.0N 105.3W 1655m'::observer, t) ORDER BY t
|
||||
) AS observations,
|
||||
array_agg(t ORDER BY t) AS times
|
||||
FROM iss,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 12:10:00+00'::timestamptz,
|
||||
interval '30 seconds') t
|
||||
)
|
||||
SELECT iterations,
|
||||
round(rms_final::numeric, 4) AS rms_km,
|
||||
status
|
||||
FROM iss, obs,
|
||||
tle_from_topocentric(
|
||||
obs.observations, obs.times,
|
||||
'40.0N 105.3W 1655m'::observer,
|
||||
seed := iss.tle
|
||||
);
|
||||
```
|
||||
|
||||
<Aside type="caution" title="Seed usually required">
|
||||
Unlike ECI fitting, topocentric fitting typically needs a seed TLE. The azimuth/elevation/range → orbit mapping is more nonlinear than the ECI case, and the solver may not converge from a poor initial guess.
|
||||
</Aside>
|
||||
|
||||
## Range rate fitting
|
||||
|
||||
When your observations include Doppler or radar range-rate data, enable `fit_range_rate` to use it as an additional constraint:
|
||||
|
||||
```sql
|
||||
SELECT iterations, round(rms_final::numeric, 4) AS rms_km, status
|
||||
FROM tle_from_topocentric(
|
||||
observations := topo_obs,
|
||||
times := obs_times,
|
||||
obs := '40.0N 105.3W 1655m'::observer,
|
||||
seed := seed_tle,
|
||||
fit_range_rate := true
|
||||
);
|
||||
```
|
||||
|
||||
Range rate adds a 4th residual component per observation (alongside azimuth, elevation, and range). This is particularly valuable for short observation arcs where position-only data doesn't constrain the velocity well. The range-rate residuals are internally scaled by a factor of 10 (1 km/s maps to 10 km equivalent) to balance their magnitude against position residuals.
|
||||
|
||||
## Multi-observer fitting
|
||||
|
||||
When observations come from multiple ground stations, use the multi-observer variant. Each observation is tagged with its originating station via the `observer_ids` array:
|
||||
|
||||
```sql
|
||||
-- Define two ground stations
|
||||
-- Station 1: Boulder, CO
|
||||
-- Station 2: Los Angeles, CA
|
||||
|
||||
SELECT iterations, round(rms_final::numeric, 4) AS rms_km, status
|
||||
FROM tle_from_topocentric(
|
||||
observations := all_topo_obs,
|
||||
times := all_times,
|
||||
observers := ARRAY[
|
||||
'40.0N 105.3W 1655m',
|
||||
'34.1N 118.3W 100m'
|
||||
]::observer[],
|
||||
observer_ids := station_ids, -- e.g., ARRAY[1,1,1,2,2,2]
|
||||
seed := seed_tle
|
||||
);
|
||||
```
|
||||
|
||||
The `observer_ids` array uses 1-based indexing into `observers[]`. Observations from different stations can be interleaved in any order — the solver uses `observer_ids[i]` to look up the correct station geometry for each observation.
|
||||
|
||||
Multi-observer data improves orbit determination geometry, especially for short arcs. Two well-separated stations observing the same pass can constrain an orbit that a single station cannot.
|
||||
|
||||
## Weighted observations
|
||||
|
||||
When observations have different accuracies (e.g., radar range vs. optical angles, or a high-quality station vs. a noisy one), use the `weights` parameter:
|
||||
|
||||
```sql
|
||||
SELECT iterations, round(rms_final::numeric, 4) AS rms_km, status
|
||||
FROM tle_from_eci(
|
||||
positions := obs_positions,
|
||||
times := obs_times,
|
||||
weights := ARRAY[1.0, 1.0, 0.5, 1.0, 0.2, 1.0]
|
||||
);
|
||||
```
|
||||
|
||||
Higher weights give observations more influence on the solution. A weight of 0.5 means "half as trustworthy as weight 1.0." This is useful for:
|
||||
|
||||
- Down-weighting observations near the horizon (higher atmospheric refraction)
|
||||
- Down-weighting observations from less accurate sensors
|
||||
- Implementing iterative reweighting after identifying outliers via `tle_fit_residuals()`
|
||||
|
||||
## Angles-only (Gauss method)
|
||||
|
||||
Optical observations often provide only RA/Dec — no range information. The `tle_from_angles()` function handles this case, using the Gauss method for initial orbit determination when no seed TLE is available:
|
||||
|
||||
```sql
|
||||
SELECT iterations,
|
||||
round(rms_final::numeric, 6) AS rms_rad,
|
||||
status,
|
||||
round(condition_number::numeric, 1) AS cond
|
||||
FROM tle_from_angles(
|
||||
ra_hours := ARRAY[12.345, 12.567, 12.789, 13.012, 13.234, 13.456],
|
||||
dec_degrees := ARRAY[45.1, 44.8, 44.3, 43.6, 42.8, 41.9],
|
||||
times := ARRAY[
|
||||
'2024-01-01 12:00+00', '2024-01-01 12:01+00',
|
||||
'2024-01-01 12:02+00', '2024-01-01 12:03+00',
|
||||
'2024-01-01 12:04+00', '2024-01-01 12:05+00'
|
||||
]::timestamptz[],
|
||||
obs := '40.0N 105.3W 1655m'::observer
|
||||
);
|
||||
```
|
||||
|
||||
<Aside type="note" title="RMS units">
|
||||
For angles-only fitting, `rms_final` and `rms_initial` are in **radians**, not km. A typical good fit produces RMS around 0.001 radians (~3.4 arcminutes).
|
||||
</Aside>
|
||||
|
||||
The Gauss method requires at least 3 observations with sufficient angular separation and time spread. The method works by solving for the geocentric range at the middle observation time, then deriving a full state vector. This bootstrap orbit is then refined by the standard DC solver.
|
||||
|
||||
<Aside type="caution" title="Gauss limitations">
|
||||
The Gauss method can fail or produce spurious results when:
|
||||
- Observations span too short an arc (insufficient angular motion)
|
||||
- The satellite is near zenith for all observations (poor parallax geometry)
|
||||
- The time gaps between observations are too large (the linear approximation breaks down)
|
||||
|
||||
When Gauss fails, provide a seed TLE from a catalog or from a separate IOD method.
|
||||
</Aside>
|
||||
|
||||
RA is in hours [0, 24) and declination in degrees [-90, 90], matching the output convention of `star_observe()`. This lets you use the same coordinate system for both stellar calibration and satellite observations.
|
||||
|
||||
## Interpreting results
|
||||
|
||||
Every OD function returns the same 8-column record. Here's how to use each field:
|
||||
|
||||
| Field | What to check |
|
||||
|-------|---------------|
|
||||
| `fitted_tle` | NULL means the solver failed completely. Check `status` for the reason. |
|
||||
| `iterations` | Should be well below `max_iter`. If it equals `max_iter`, the solver didn't converge — increase `max_iter` or improve the seed. |
|
||||
| `rms_final` | The fit quality. For ECI: sub-km is good. For topocentric: 1-10 km is typical. For angles-only: < 0.01 rad is good. |
|
||||
| `rms_initial` | Compare with `rms_final`. A large ratio (rms_initial / rms_final) means the solver made significant improvement. If they're similar, the seed was already good or the solver made little progress. |
|
||||
| `status` | `'converged'` is success. `'max_iterations'` means more iterations needed. Other strings describe specific failures. |
|
||||
| `condition_number` | Measures how well-conditioned the geometry is. Values below ~1e6 are well-conditioned. Above ~1e10 suggests the observations don't constrain all orbital elements. |
|
||||
| `covariance` | Formal uncertainty. Extract diagonal elements for per-parameter variance. Off-diagonal elements show correlations between parameters. |
|
||||
| `nstate` | 6 (equinoctial elements) or 7 (with B*). The covariance matrix is `nstate x nstate`. |
|
||||
|
||||
## Tips
|
||||
|
||||
**Seed selection.** For ECI fitting, the Gibbs IOD bootstrap usually works. For topocentric and angles-only fitting, start with a TLE from a catalog (Space-Track, CelesTrak) for the target object. The seed doesn't need to be precise — it just needs to be in the right ballpark (correct orbit regime, approximate inclination).
|
||||
|
||||
**Minimum observations.** The solver needs at least 6 observations for a 6-parameter fit (or 7 for B*). In practice, more is better — 10-20 well-distributed observations across at least one orbit give the solver enough leverage to constrain all elements.
|
||||
|
||||
**Observation spacing.** Spread observations across the orbit. Clustered observations from a short arc constrain some elements well (position) but others poorly (period, eccentricity). A full orbit of data is ideal.
|
||||
|
||||
**Convergence troubleshooting.** If `status` shows `'max_iterations'`:
|
||||
1. Try a better seed TLE (closer to the true orbit)
|
||||
2. Increase `max_iter` to 30 or 50
|
||||
3. Check if observations contain outliers using `tle_fit_residuals()`
|
||||
4. Verify that the observation timestamps are correct (off-by-one-hour timezone errors are common)
|
||||
|
||||
**B* fitting.** Only enable `fit_bstar` when you have observations spanning multiple orbits. B* captures atmospheric drag, which manifests as a slow period decay — a short arc doesn't have enough signal to constrain it. With insufficient data, fitting B* degrades the solution by introducing an under-determined 7th parameter.
|
||||
@ -44,7 +44,7 @@ pg_orrery propagates TLEs and computes look angles. It does not replace the full
|
||||
- **No real-time GUI.** GPredict and STK provide map displays, polar plots, and Doppler displays. pg_orrery returns numbers. Use any visualization tool to render its output.
|
||||
- **No rotator control.** Hamlib drives antenna rotators. pg_orrery computes the azimuth and elevation values Hamlib would consume, but it has no hardware interface.
|
||||
- **No TLE fetching.** Bring your own TLEs from Space-Track, CelesTrak, or any provider. pg_orrery parses and propagates them.
|
||||
- **Orbit determination available.** Since v0.4.0, pg_orrery can fit TLEs from ECI, topocentric, or angles-only observations via differential correction. See the [Orbit Determination guide](/guides/orbit-determination/).
|
||||
- **No orbit determination.** pg_orrery propagates existing TLEs. It does not fit orbits from observations.
|
||||
- **No high-precision propagation.** SGP4/SDP4 accuracy degrades with TLE age. For operational conjunction assessment, use SP ephemerides or owner/operator-provided state vectors. pg_orrery's GiST screening finds candidates; you verify with better data.
|
||||
|
||||
## Try it
|
||||
|
||||
@ -43,12 +43,6 @@ import { Card, CardGrid, LinkCard } from "@astrojs/starlight/components";
|
||||
Pork chop plots as SQL CROSS JOINs — 22,500 transfer solutions in
|
||||
8.3 seconds. Departure C3, arrival C3, time of flight, transfer SMA.
|
||||
</Card>
|
||||
<Card title="Determine orbits" icon="seti:satellite">
|
||||
Fit TLEs from ECI ephemeris, ground station observations (az/el/range with
|
||||
optional range rate), or angles-only RA/Dec. Gauss and Gibbs methods for
|
||||
seed-free initial orbit determination. Weighted least-squares with
|
||||
covariance output.
|
||||
</Card>
|
||||
</CardGrid>
|
||||
|
||||
## Explore the docs
|
||||
@ -61,12 +55,12 @@ import { Card, CardGrid, LinkCard } from "@astrojs/starlight/components";
|
||||
/>
|
||||
<LinkCard
|
||||
title="Guides"
|
||||
description="Domain-specific walkthroughs for satellites, planets, moons, stars, comets, radio, trajectories, and orbit determination"
|
||||
description="Domain-specific walkthroughs for satellites, planets, moons, stars, comets, radio, and trajectories"
|
||||
href="/guides/tracking-satellites/"
|
||||
/>
|
||||
<LinkCard
|
||||
title="Workflow Translation"
|
||||
description="Side-by-side comparisons with Skyfield, Horizons, GMAT, find_orb, and Poliastro"
|
||||
description="Side-by-side comparisons: how you do it today vs. how pg_orrery changes the game"
|
||||
href="/workflow/from-skyfield/"
|
||||
/>
|
||||
<LinkCard
|
||||
|
||||
@ -207,18 +207,6 @@ Keplerian propagation ignores gravitational perturbations from planets, non-grav
|
||||
|
||||
For preliminary mission design and pork chop plot generation, these limitations are standard and expected.
|
||||
|
||||
### Orbit Determination
|
||||
|
||||
| Quantity | Notes |
|
||||
|----------|-------|
|
||||
| ECI fitting RMS | Sub-km typical with clean observations (round-trip from propagated state) |
|
||||
| Topocentric fitting RMS | ~1-10 km depending on arc length and observation spacing |
|
||||
| Angles-only fitting RMS | Output in radians. ~0.001 rad with good geometry |
|
||||
| Condition number | Formal indicator of solution quality. Values above ~1e10 suggest poorly-constrained geometry |
|
||||
| Covariance | Formal (H^T H)^{-1} from final Jacobian. Optimistic; does not include systematic errors |
|
||||
|
||||
**Limitations:** The DC solver fits SGP4/SDP4 mean elements (6 equinoctial + optional B*). Accuracy is bounded by the TLE/SGP4 accuracy floor (~1 km at epoch for LEO). Range rate uses a fixed scale factor (OD_RR_SCALE = 10.0, mapping 1 km/s to 10 km equivalent). Gauss IOD requires at least 3 well-spaced observations with sufficient angular separation.
|
||||
|
||||
---
|
||||
|
||||
## Reference Publications
|
||||
@ -234,4 +222,3 @@ For preliminary mission design and pork chop plot generation, these limitations
|
||||
| MarsSat | Jacobson, R.A. "The orbits and masses of the Martian satellites and the libration of Phobos." The Astronomical Journal, 139, 668-679, 2010. |
|
||||
| Carr source regions | Carr, T.D., Desch, M.D., Alexander, J.K. "Phenomenology of magnetospheric radio emissions." In Physics of the Jovian Magnetosphere, Cambridge Univ. Press, 1983. |
|
||||
| Lambert solver | Battin, R.H. "An Introduction to the Mathematics and Methods of Astrodynamics." AIAA Education Series, Revised Edition, 1999. |
|
||||
| Orbit determination | Vallado, D.A. "Fundamentals of Astrodynamics and Applications." 4th ed., Microcosm Press, 2013. |
|
||||
|
||||
@ -1,434 +0,0 @@
|
||||
---
|
||||
title: "Functions: Orbit Determination"
|
||||
sidebar:
|
||||
order: 8
|
||||
---
|
||||
|
||||
import { Aside, Tabs, TabItem } from "@astrojs/starlight/components";
|
||||
|
||||
Functions for fitting TLE orbital elements from observations via differential correction (DC). The solver uses equinoctial elements internally to avoid singularities at zero eccentricity and inclination, with LAPACK SVD (`dgelss_`) for the least-squares solve.
|
||||
|
||||
Three observation types are supported: ECI position/velocity, topocentric (az/el/range with optional range rate), and angles-only (RA/Dec). Each has single-observer and multi-observer variants where applicable.
|
||||
|
||||
All OD functions share the same 8-column output record.
|
||||
|
||||
---
|
||||
|
||||
## Shared Output Record
|
||||
|
||||
Every OD function returns a `RECORD` with these fields:
|
||||
|
||||
| Field | Type | Description |
|
||||
|-------|------|-------------|
|
||||
| `fitted_tle` | `tle` | The fitted TLE. NULL if the solver did not converge. |
|
||||
| `iterations` | `int4` | Number of DC iterations performed. |
|
||||
| `rms_final` | `float8` | RMS residual after final iteration. Units depend on the observation type (km for ECI/topocentric, radians for angles-only). |
|
||||
| `rms_initial` | `float8` | RMS residual before the first iteration. Compare with `rms_final` to assess improvement. |
|
||||
| `status` | `text` | Convergence status: `'converged'`, `'max_iterations'`, or an error description. |
|
||||
| `condition_number` | `float8` | Condition number of the normal equations matrix. Values above ~1e10 suggest poorly-constrained geometry. |
|
||||
| `covariance` | `float8[]` | Formal covariance matrix `(H^T H)^{-1}` from the final Jacobian, stored as a flat array in row-major order. Length is `nstate * nstate`. |
|
||||
| `nstate` | `int4` | Number of estimated parameters (6 for equinoctial elements, 7 if `fit_bstar` is true). |
|
||||
|
||||
<Aside type="note" title="Covariance interpretation">
|
||||
The covariance matrix is formal — it reflects the linear approximation at the solution and does not account for systematic errors in the observation model. It is optimistic. Use it for relative comparisons between solutions, not as an absolute uncertainty bound.
|
||||
</Aside>
|
||||
|
||||
---
|
||||
|
||||
## tle_from_eci
|
||||
|
||||
Fit a TLE from ECI position/velocity observations. This is the simplest OD mode — the observations are already in the SGP4 propagation frame (TEME), so no observer geometry is involved.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
tle_from_eci(
|
||||
positions eci_position[],
|
||||
times timestamptz[],
|
||||
seed tle DEFAULT NULL,
|
||||
fit_bstar boolean DEFAULT false,
|
||||
max_iter int4 DEFAULT 15,
|
||||
weights float8[] DEFAULT NULL,
|
||||
OUT fitted_tle tle,
|
||||
OUT iterations int4,
|
||||
OUT rms_final float8,
|
||||
OUT rms_initial float8,
|
||||
OUT status text,
|
||||
OUT condition_number float8,
|
||||
OUT covariance float8[],
|
||||
OUT nstate int4
|
||||
) RETURNS RECORD
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Default | Description |
|
||||
|-----------|------|---------|-------------|
|
||||
| `positions` | `eci_position[]` | — | Array of ECI position/velocity observations. Requires >= 6 observations. |
|
||||
| `times` | `timestamptz[]` | — | Observation timestamps. Must be same length as `positions`. |
|
||||
| `seed` | `tle` | `NULL` | Initial TLE estimate. When NULL, the solver uses the Gibbs or Herrick-Gibbs method to bootstrap an initial orbit from three position vectors. |
|
||||
| `fit_bstar` | `boolean` | `false` | When true, estimates B* drag coefficient as a 7th parameter. Requires a well-distributed observation arc. |
|
||||
| `max_iter` | `int4` | `15` | Maximum differential correction iterations. |
|
||||
| `weights` | `float8[]` | `NULL` | Per-observation weights. NULL means uniform weighting. Length must equal length of `positions`. Higher weight = more influence on the solution. |
|
||||
|
||||
<Aside type="tip" title="Seed-free fitting">
|
||||
When no seed TLE is provided, the solver uses the Gibbs method (or Herrick-Gibbs for closely-spaced observations) to derive an initial velocity estimate from three position vectors. This makes `tle_from_eci()` fully seed-free for most use cases.
|
||||
</Aside>
|
||||
|
||||
### Example
|
||||
|
||||
```sql
|
||||
-- Round-trip test: propagate a TLE, then fit it back
|
||||
WITH iss AS (
|
||||
SELECT '1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025
|
||||
2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS tle
|
||||
),
|
||||
obs AS (
|
||||
SELECT array_agg(sgp4_propagate(iss.tle, t) ORDER BY t) AS positions,
|
||||
array_agg(t ORDER BY t) AS times
|
||||
FROM iss,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
interval '5 minutes') t
|
||||
)
|
||||
SELECT iterations,
|
||||
round(rms_final::numeric, 6) AS rms_km,
|
||||
round(rms_initial::numeric, 3) AS rms_init_km,
|
||||
status,
|
||||
round(condition_number::numeric, 1) AS cond
|
||||
FROM obs, tle_from_eci(obs.positions, obs.times);
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## tle_from_topocentric (single observer)
|
||||
|
||||
Fit a TLE from topocentric (azimuth/elevation/range) observations collected by a single ground station.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
tle_from_topocentric(
|
||||
observations topocentric[],
|
||||
times timestamptz[],
|
||||
obs observer,
|
||||
seed tle DEFAULT NULL,
|
||||
fit_bstar boolean DEFAULT false,
|
||||
max_iter int4 DEFAULT 15,
|
||||
fit_range_rate boolean DEFAULT false,
|
||||
weights float8[] DEFAULT NULL,
|
||||
OUT fitted_tle tle,
|
||||
OUT iterations int4,
|
||||
OUT rms_final float8,
|
||||
OUT rms_initial float8,
|
||||
OUT status text,
|
||||
OUT condition_number float8,
|
||||
OUT covariance float8[],
|
||||
OUT nstate int4
|
||||
) RETURNS RECORD
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Default | Description |
|
||||
|-----------|------|---------|-------------|
|
||||
| `observations` | `topocentric[]` | — | Array of topocentric observations (az, el, range, range_rate). Requires >= 6 observations. |
|
||||
| `times` | `timestamptz[]` | — | Observation timestamps. |
|
||||
| `obs` | `observer` | — | Ground station location. |
|
||||
| `seed` | `tle` | `NULL` | Initial TLE estimate. Unlike ECI fitting, topocentric fitting typically needs a seed for convergence. |
|
||||
| `fit_bstar` | `boolean` | `false` | Estimate B* drag coefficient. |
|
||||
| `max_iter` | `int4` | `15` | Maximum iterations. |
|
||||
| `fit_range_rate` | `boolean` | `false` | When true, includes range rate as a 4th residual component per observation. Use when you have Doppler or radar range-rate data. |
|
||||
| `weights` | `float8[]` | `NULL` | Per-observation weights. NULL = uniform. |
|
||||
|
||||
<Aside type="note" title="Range rate scaling">
|
||||
Range rate residuals are internally scaled by `OD_RR_SCALE = 10.0`, mapping 1 km/s of range-rate error to 10 km equivalent in the residual vector. This prevents the typically small range-rate values from being swamped by position residuals.
|
||||
</Aside>
|
||||
|
||||
### Example
|
||||
|
||||
```sql
|
||||
-- Observe a satellite, then fit back from topocentric data
|
||||
WITH iss AS (
|
||||
SELECT '1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025
|
||||
2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS tle
|
||||
),
|
||||
obs AS (
|
||||
SELECT array_agg(observe(iss.tle, '40.0N 105.3W 1655m'::observer, t) ORDER BY t)
|
||||
AS observations,
|
||||
array_agg(t ORDER BY t) AS times
|
||||
FROM iss,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 12:10:00+00'::timestamptz,
|
||||
interval '30 seconds') t
|
||||
)
|
||||
SELECT iterations,
|
||||
round(rms_final::numeric, 4) AS rms_km,
|
||||
status
|
||||
FROM obs,
|
||||
tle_from_topocentric(
|
||||
obs.observations, obs.times,
|
||||
'40.0N 105.3W 1655m'::observer,
|
||||
seed := iss.tle
|
||||
);
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## tle_from_topocentric (multi-observer)
|
||||
|
||||
Fit a TLE from topocentric observations collected by multiple ground stations. The `observer_ids` array maps each observation to its originating station.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
tle_from_topocentric(
|
||||
observations topocentric[],
|
||||
times timestamptz[],
|
||||
observers observer[],
|
||||
observer_ids int4[],
|
||||
seed tle DEFAULT NULL,
|
||||
fit_bstar boolean DEFAULT false,
|
||||
max_iter int4 DEFAULT 15,
|
||||
fit_range_rate boolean DEFAULT false,
|
||||
weights float8[] DEFAULT NULL,
|
||||
OUT fitted_tle tle,
|
||||
OUT iterations int4,
|
||||
OUT rms_final float8,
|
||||
OUT rms_initial float8,
|
||||
OUT status text,
|
||||
OUT condition_number float8,
|
||||
OUT covariance float8[],
|
||||
OUT nstate int4
|
||||
) RETURNS RECORD
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Default | Description |
|
||||
|-----------|------|---------|-------------|
|
||||
| `observations` | `topocentric[]` | — | All observations from all stations, interleaved in time order. |
|
||||
| `times` | `timestamptz[]` | — | Observation timestamps. |
|
||||
| `observers` | `observer[]` | — | Array of ground station locations (1-indexed). |
|
||||
| `observer_ids` | `int4[]` | — | Per-observation index into `observers[]`. `observer_ids[i]` identifies which station produced `observations[i]`. |
|
||||
| `seed` | `tle` | `NULL` | Initial TLE estimate. |
|
||||
| `fit_bstar` | `boolean` | `false` | Estimate B* drag coefficient. |
|
||||
| `max_iter` | `int4` | `15` | Maximum iterations. |
|
||||
| `fit_range_rate` | `boolean` | `false` | Include range rate in residuals. |
|
||||
| `weights` | `float8[]` | `NULL` | Per-observation weights. Useful when stations have different measurement accuracies. |
|
||||
|
||||
### Example
|
||||
|
||||
```sql
|
||||
-- Two ground stations observe the same satellite
|
||||
SELECT iterations, round(rms_final::numeric, 4) AS rms_km, status
|
||||
FROM tle_from_topocentric(
|
||||
observations := ARRAY[obs1_t1, obs1_t2, obs2_t1, obs2_t2]::topocentric[],
|
||||
times := ARRAY['2024-01-01 12:00+00', '2024-01-01 12:05+00',
|
||||
'2024-01-01 12:02+00', '2024-01-01 12:07+00']::timestamptz[],
|
||||
observers := ARRAY['40.0N 105.3W 1655m', '34.1N 118.3W 100m']::observer[],
|
||||
observer_ids := ARRAY[1, 1, 2, 2],
|
||||
seed := seed_tle
|
||||
);
|
||||
```
|
||||
|
||||
<Aside type="tip" title="Observer ID indexing">
|
||||
`observer_ids` uses 1-based indexing into the `observers` array. Station 1 is `observers[1]`, station 2 is `observers[2]`, etc. This matches PostgreSQL's native array indexing convention.
|
||||
</Aside>
|
||||
|
||||
---
|
||||
|
||||
## tle_from_angles (single observer)
|
||||
|
||||
Fit a TLE from angles-only (RA/Dec) observations collected by a single ground station. When no seed TLE is provided, the Gauss method derives an initial orbit from three observations — making this function fully seed-free.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
tle_from_angles(
|
||||
ra_hours float8[],
|
||||
dec_degrees float8[],
|
||||
times timestamptz[],
|
||||
obs observer,
|
||||
seed tle DEFAULT NULL,
|
||||
fit_bstar boolean DEFAULT false,
|
||||
max_iter int4 DEFAULT 15,
|
||||
weights float8[] DEFAULT NULL,
|
||||
OUT fitted_tle tle,
|
||||
OUT iterations int4,
|
||||
OUT rms_final float8,
|
||||
OUT rms_initial float8,
|
||||
OUT status text,
|
||||
OUT condition_number float8,
|
||||
OUT covariance float8[],
|
||||
OUT nstate int4
|
||||
) RETURNS RECORD
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Default | Description |
|
||||
|-----------|------|---------|-------------|
|
||||
| `ra_hours` | `float8[]` | — | Right ascension values in hours, range [0, 24). Matches the `star_observe()` convention. |
|
||||
| `dec_degrees` | `float8[]` | — | Declination values in degrees, range [-90, 90]. |
|
||||
| `times` | `timestamptz[]` | — | Observation timestamps. Requires >= 3 observations. |
|
||||
| `obs` | `observer` | — | Ground station location. |
|
||||
| `seed` | `tle` | `NULL` | Initial TLE estimate. When NULL, the Gauss method bootstraps an initial orbit from 3 observations. |
|
||||
| `fit_bstar` | `boolean` | `false` | Estimate B* drag coefficient. |
|
||||
| `max_iter` | `int4` | `15` | Maximum iterations. |
|
||||
| `weights` | `float8[]` | `NULL` | Per-observation weights. |
|
||||
|
||||
<Aside type="caution" title="RMS units">
|
||||
For angles-only fitting, `rms_final` and `rms_initial` are in **radians**, not km. A typical good fit produces RMS values around 0.001 radians (~0.06 degrees).
|
||||
</Aside>
|
||||
|
||||
<Aside type="note" title="RA/Dec conventions">
|
||||
Right ascension is in hours [0, 24), matching the `star_observe()` output convention. Declination is in degrees [-90, 90]. Internally, both are converted to radians for the residual computation.
|
||||
</Aside>
|
||||
|
||||
### Example
|
||||
|
||||
```sql
|
||||
-- Angles-only OD: RA/Dec observations of a satellite
|
||||
SELECT iterations,
|
||||
round(rms_final::numeric, 6) AS rms_rad,
|
||||
status,
|
||||
round(condition_number::numeric, 1) AS cond
|
||||
FROM tle_from_angles(
|
||||
ra_hours := ARRAY[12.345, 12.567, 12.789, 13.012, 13.234, 13.456],
|
||||
dec_degrees := ARRAY[45.1, 44.8, 44.3, 43.6, 42.8, 41.9],
|
||||
times := ARRAY[
|
||||
'2024-01-01 12:00+00', '2024-01-01 12:01+00',
|
||||
'2024-01-01 12:02+00', '2024-01-01 12:03+00',
|
||||
'2024-01-01 12:04+00', '2024-01-01 12:05+00'
|
||||
]::timestamptz[],
|
||||
obs := '40.0N 105.3W 1655m'::observer
|
||||
);
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## tle_from_angles (multi-observer)
|
||||
|
||||
Fit a TLE from angles-only (RA/Dec) observations collected by multiple ground stations. Uses the same Gauss IOD bootstrap as the single-observer variant when no seed is provided.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
tle_from_angles(
|
||||
ra_hours float8[],
|
||||
dec_degrees float8[],
|
||||
times timestamptz[],
|
||||
observers observer[],
|
||||
observer_ids int4[],
|
||||
seed tle DEFAULT NULL,
|
||||
fit_bstar boolean DEFAULT false,
|
||||
max_iter int4 DEFAULT 15,
|
||||
weights float8[] DEFAULT NULL,
|
||||
OUT fitted_tle tle,
|
||||
OUT iterations int4,
|
||||
OUT rms_final float8,
|
||||
OUT rms_initial float8,
|
||||
OUT status text,
|
||||
OUT condition_number float8,
|
||||
OUT covariance float8[],
|
||||
OUT nstate int4
|
||||
) RETURNS RECORD
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Default | Description |
|
||||
|-----------|------|---------|-------------|
|
||||
| `ra_hours` | `float8[]` | — | Right ascension values in hours [0, 24). |
|
||||
| `dec_degrees` | `float8[]` | — | Declination values in degrees [-90, 90]. |
|
||||
| `times` | `timestamptz[]` | — | Observation timestamps from all stations, interleaved in time order. |
|
||||
| `observers` | `observer[]` | — | Array of ground station locations (1-indexed). |
|
||||
| `observer_ids` | `int4[]` | — | Per-observation index into `observers[]`. |
|
||||
| `seed` | `tle` | `NULL` | Initial TLE estimate. NULL = Gauss IOD bootstrap. |
|
||||
| `fit_bstar` | `boolean` | `false` | Estimate B* drag coefficient. |
|
||||
| `max_iter` | `int4` | `15` | Maximum iterations. |
|
||||
| `weights` | `float8[]` | `NULL` | Per-observation weights. Useful when stations have different apertures or sky conditions. |
|
||||
|
||||
### Example
|
||||
|
||||
```sql
|
||||
-- Two optical stations observe the same satellite
|
||||
SELECT iterations, round(rms_final::numeric, 6) AS rms_rad, status
|
||||
FROM tle_from_angles(
|
||||
ra_hours := ARRAY[12.3, 12.5, 12.7, 12.4, 12.6, 12.8],
|
||||
dec_degrees := ARRAY[45.0, 44.5, 44.0, 44.8, 44.3, 43.7],
|
||||
times := ARRAY[
|
||||
'2024-01-01 12:00+00', '2024-01-01 12:02+00',
|
||||
'2024-01-01 12:04+00', '2024-01-01 12:01+00',
|
||||
'2024-01-01 12:03+00', '2024-01-01 12:05+00'
|
||||
]::timestamptz[],
|
||||
observers := ARRAY['40.0N 105.3W 1655m', '34.1N 118.3W 100m']::observer[],
|
||||
observer_ids := ARRAY[1, 1, 1, 2, 2, 2]
|
||||
);
|
||||
```
|
||||
|
||||
---
|
||||
|
||||
## tle_fit_residuals
|
||||
|
||||
Compute per-observation position residuals between a fitted TLE and the original ECI observations. Returns one row per observation with the XYZ and total position error in km. Use this to identify outlier observations or assess spatial error distribution.
|
||||
|
||||
### Signature
|
||||
|
||||
```sql
|
||||
tle_fit_residuals(
|
||||
fitted tle,
|
||||
positions eci_position[],
|
||||
times timestamptz[]
|
||||
) RETURNS TABLE (
|
||||
t timestamptz,
|
||||
dx_km float8,
|
||||
dy_km float8,
|
||||
dz_km float8,
|
||||
pos_err_km float8
|
||||
)
|
||||
```
|
||||
|
||||
### Parameters
|
||||
|
||||
| Parameter | Type | Description |
|
||||
|-----------|------|-------------|
|
||||
| `fitted` | `tle` | The fitted TLE (typically the `fitted_tle` output from `tle_from_eci()`). |
|
||||
| `positions` | `eci_position[]` | The original ECI observations used for fitting. |
|
||||
| `times` | `timestamptz[]` | The original observation timestamps. |
|
||||
|
||||
### Returns
|
||||
|
||||
One row per observation:
|
||||
|
||||
| Column | Type | Unit | Description |
|
||||
|--------|------|------|-------------|
|
||||
| `t` | `timestamptz` | — | Observation time |
|
||||
| `dx_km` | `float8` | km | X-axis residual (observed - computed) |
|
||||
| `dy_km` | `float8` | km | Y-axis residual |
|
||||
| `dz_km` | `float8` | km | Z-axis residual |
|
||||
| `pos_err_km` | `float8` | km | Total position error: sqrt(dx^2 + dy^2 + dz^2) |
|
||||
|
||||
### Example
|
||||
|
||||
```sql
|
||||
-- After fitting, inspect per-observation residuals
|
||||
WITH fit AS (
|
||||
SELECT fitted_tle, positions, times
|
||||
FROM obs, tle_from_eci(obs.positions, obs.times)
|
||||
)
|
||||
SELECT t,
|
||||
round(dx_km::numeric, 4) AS dx,
|
||||
round(dy_km::numeric, 4) AS dy,
|
||||
round(dz_km::numeric, 4) AS dz,
|
||||
round(pos_err_km::numeric, 4) AS total_err
|
||||
FROM fit, tle_fit_residuals(fit.fitted_tle, fit.positions, fit.times)
|
||||
ORDER BY pos_err_km DESC;
|
||||
```
|
||||
|
||||
<Aside type="tip" title="Outlier detection">
|
||||
Observations with residuals significantly larger than the RMS are potential outliers. Consider removing them and re-fitting, or using the `weights` parameter to down-weight them.
|
||||
</Aside>
|
||||
@ -1,373 +0,0 @@
|
||||
---
|
||||
title: From find_orb to SQL
|
||||
sidebar:
|
||||
order: 5
|
||||
description: Comparing find_orb's orbit determination workflow with pg_orrery's SQL-based differential correction.
|
||||
---
|
||||
|
||||
import { Tabs, TabItem, Aside, Steps } from "@astrojs/starlight/components";
|
||||
|
||||
[find_orb](https://github.com/Bill-Gray/find_orb) is Bill Gray's orbit determination software — the same Bill Gray whose [sat_code SGP4 library](https://github.com/Bill-Gray/sat_code) is vendored inside pg_orrery. find_orb takes astrometric observations (RA/Dec positions with timestamps) and fits orbital elements via differential correction. It handles asteroids, comets, and artificial satellites. Amateur astronomers, minor planet observers, and satellite trackers use it worldwide.
|
||||
|
||||
Since v0.4.0, pg_orrery has its own differential correction solver. The domain is narrower — pg_orrery fits SGP4/SDP4 mean elements (TLEs) for Earth-orbiting satellites, not heliocentric orbits for asteroids — but within that domain, the SQL approach offers batch processing, data integration, and automation that a desktop application cannot match.
|
||||
|
||||
<Aside type="note" title="Different orbit models">
|
||||
find_orb fits osculating Keplerian elements suitable for numerical integration. pg_orrery fits SGP4 mean elements (TLEs) designed for the SGP4/SDP4 analytical propagator. These are different representations. A TLE absorbs SGP4's internal model biases — it only produces accurate positions when propagated with SGP4. This is a design choice, not a limitation: TLEs are the operational standard for satellite tracking, and pg_orrery's fitted TLEs plug directly into `observe()`, `predict_passes()`, and every other satellite function.
|
||||
</Aside>
|
||||
|
||||
## Fitting an orbit from RA/Dec observations
|
||||
|
||||
The core workflow: you have a series of sky positions (right ascension and declination) for an object, and you want to determine its orbit.
|
||||
|
||||
<Tabs>
|
||||
<TabItem label="find_orb">
|
||||
find_orb reads MPC 80-column format observation files:
|
||||
|
||||
```
|
||||
ISS C2024 01 01.50000 12 20 42.0 +45 06 00.0 15.0 V XXX
|
||||
ISS C2024 01 01.50069 12 34 01.2 +44 48 00.0 15.0 V XXX
|
||||
ISS C2024 01 01.50139 12 47 18.6 +44 18 00.0 15.0 V XXX
|
||||
ISS C2024 01 01.50208 13 00 43.2 +43 36 00.0 15.0 V XXX
|
||||
ISS C2024 01 01.50278 13 14 02.4 +42 48 00.0 15.0 V XXX
|
||||
ISS C2024 01 01.50347 13 27 21.6 +41 54 00.0 15.0 V XXX
|
||||
```
|
||||
|
||||
Then either:
|
||||
|
||||
**GUI:** Launch find_orb, open the observation file, select the object, click "Full Step" repeatedly until the residuals converge. Inspect the residual map visually. Copy the fitted elements from the output panel.
|
||||
|
||||
**CLI:** Run `fo obs_file.txt` to process observations non-interactively. Parse the output file for fitted elements.
|
||||
|
||||
For each object, this is a separate run. Processing 50 objects from a night of observing means 50 file preparations and 50 find_orb runs. The results end up in text files that you then parse and load into your database.
|
||||
</TabItem>
|
||||
<TabItem label="pg_orrery (SQL)">
|
||||
```sql
|
||||
-- Angles-only OD: RA/Dec observations → fitted TLE
|
||||
SELECT iterations,
|
||||
round(rms_final::numeric, 6) AS rms_rad,
|
||||
status,
|
||||
round(condition_number::numeric, 1) AS cond,
|
||||
fitted_tle
|
||||
FROM tle_from_angles(
|
||||
ra_hours := ARRAY[12.345, 12.567, 12.789, 13.012, 13.234, 13.456],
|
||||
dec_degrees := ARRAY[45.1, 44.8, 44.3, 43.6, 42.8, 41.9],
|
||||
times := ARRAY[
|
||||
'2024-01-01 12:00+00', '2024-01-01 12:01+00',
|
||||
'2024-01-01 12:02+00', '2024-01-01 12:03+00',
|
||||
'2024-01-01 12:04+00', '2024-01-01 12:05+00'
|
||||
]::timestamptz[],
|
||||
obs := '40.0N 105.3W 1655m'::observer
|
||||
);
|
||||
```
|
||||
|
||||
No file format to prepare. RA in hours, Dec in degrees — direct numerical input. The Gauss method bootstraps the initial orbit automatically when no seed TLE is provided.
|
||||
|
||||
The fitted TLE is immediately usable with any pg_orrery satellite function:
|
||||
|
||||
```sql
|
||||
-- Fit and immediately predict passes
|
||||
WITH od AS (
|
||||
SELECT fitted_tle
|
||||
FROM tle_from_angles(
|
||||
ra_hours := ARRAY[12.345, 12.567, 12.789, 13.012, 13.234, 13.456],
|
||||
dec_degrees := ARRAY[45.1, 44.8, 44.3, 43.6, 42.8, 41.9],
|
||||
times := ARRAY[
|
||||
'2024-01-01 12:00+00', '2024-01-01 12:01+00',
|
||||
'2024-01-01 12:02+00', '2024-01-01 12:03+00',
|
||||
'2024-01-01 12:04+00', '2024-01-01 12:05+00'
|
||||
]::timestamptz[],
|
||||
obs := '40.0N 105.3W 1655m'::observer
|
||||
)
|
||||
)
|
||||
SELECT pass_aos(p) AS rise,
|
||||
pass_max_el(p) AS max_el,
|
||||
pass_los(p) AS set
|
||||
FROM od,
|
||||
predict_passes(od.fitted_tle,
|
||||
'40.0N 105.3W 1655m'::observer,
|
||||
now(), now() + interval '24 hours', 10.0) p;
|
||||
```
|
||||
|
||||
Observation → orbit fit → pass prediction in a single query.
|
||||
</TabItem>
|
||||
</Tabs>
|
||||
|
||||
## Batch orbit determination
|
||||
|
||||
This is where the architecture diverges most sharply. find_orb processes one object at a time. pg_orrery processes however many your query describes.
|
||||
|
||||
<Tabs>
|
||||
<TabItem label="find_orb">
|
||||
```bash
|
||||
# Process a night's observations — one object at a time
|
||||
for obj in $(ls obs_files/*.txt); do
|
||||
fo "$obj" > results/$(basename "$obj" .txt).out 2>&1
|
||||
done
|
||||
|
||||
# Parse each output file for fitted elements
|
||||
# Build a pipeline to extract orbital elements, residuals, etc.
|
||||
# Load results into your database
|
||||
```
|
||||
|
||||
For 50 objects, this means 50 invocations, 50 output files, and a parsing pipeline. Convergence failures need manual attention — find_orb's differential corrector may not converge for all objects, and the failure modes differ per object.
|
||||
</TabItem>
|
||||
<TabItem label="pg_orrery (SQL)">
|
||||
```sql
|
||||
-- Batch OD: fit orbits for all objects observed tonight
|
||||
-- Assumes an observation table with RA/Dec per object
|
||||
SELECT obj_id,
|
||||
(od).iterations,
|
||||
round((od).rms_final::numeric, 6) AS rms_rad,
|
||||
(od).status,
|
||||
(od).fitted_tle
|
||||
FROM (
|
||||
SELECT obj_id,
|
||||
tle_from_angles(
|
||||
ra_hours := array_agg(ra_hours ORDER BY obs_time),
|
||||
dec_degrees := array_agg(dec_deg ORDER BY obs_time),
|
||||
times := array_agg(obs_time ORDER BY obs_time),
|
||||
obs := '40.0N 105.3W 1655m'::observer
|
||||
) AS od
|
||||
FROM tonight_observations
|
||||
WHERE obs_time > now() - interval '12 hours'
|
||||
GROUP BY obj_id
|
||||
HAVING count(*) >= 3 -- Need at least 3 for Gauss IOD
|
||||
) sub
|
||||
ORDER BY (od).rms_final;
|
||||
```
|
||||
|
||||
Every object observed tonight, fitted in one query. Objects that fail to converge show `status != 'converged'` — they stay in the result set instead of halting the pipeline. Sort by RMS to see the best fits first.
|
||||
</TabItem>
|
||||
</Tabs>
|
||||
|
||||
## Multi-observer observations
|
||||
|
||||
Observations from multiple stations improve orbit geometry. Both tools support this, but the workflow differs.
|
||||
|
||||
<Tabs>
|
||||
<TabItem label="find_orb">
|
||||
In MPC format, each observation line includes a 3-character observatory code (column 78-80). find_orb looks up the observatory coordinates from a `stations.txt` file. To combine data from multiple observers:
|
||||
|
||||
1. Ensure all observers have MPC observatory codes (or define custom ones)
|
||||
2. Concatenate observation files, maintaining the 80-column format
|
||||
3. Make sure `stations.txt` contains coordinates for all observatory codes
|
||||
4. Load the combined file into find_orb
|
||||
|
||||
The station lookup is implicit — based on the observatory code column in the fixed-width format.
|
||||
</TabItem>
|
||||
<TabItem label="pg_orrery (SQL)">
|
||||
```sql
|
||||
-- Multi-observer OD with explicit station geometry
|
||||
SELECT (od).iterations,
|
||||
round((od).rms_final::numeric, 6) AS rms_rad,
|
||||
(od).status,
|
||||
round((od).condition_number::numeric, 1) AS cond
|
||||
FROM tle_from_angles(
|
||||
ra_hours := array_agg(ra ORDER BY t),
|
||||
dec_degrees := array_agg(dec ORDER BY t),
|
||||
times := array_agg(t ORDER BY t),
|
||||
observers := ARRAY[
|
||||
'40.0N 105.3W 1655m', -- Boulder
|
||||
'34.1N 118.3W 100m', -- Los Angeles
|
||||
'32.0N 110.9W 730m' -- Tucson
|
||||
]::observer[],
|
||||
observer_ids := array_agg(station_id ORDER BY t)
|
||||
) AS od
|
||||
FROM multi_station_observations
|
||||
WHERE target = 'UNKNOWN-2024A';
|
||||
```
|
||||
|
||||
Station geometry is explicit — no lookup file, no observatory code registry. Each observation carries its station ID directly. Adding a new station is adding a value to the `observers` array.
|
||||
</TabItem>
|
||||
</Tabs>
|
||||
|
||||
## Weighted observations and range rate
|
||||
|
||||
v0.6.0 added per-observation weights and range rate fitting. These features are most useful when combining heterogeneous data — different sensors, different accuracies, or a mix of optical and radar observations.
|
||||
|
||||
<Tabs>
|
||||
<TabItem label="find_orb">
|
||||
find_orb handles observation weighting internally. It assigns weights based on observatory statistics from the MPC's Observatory Performance page — observatories with historically tighter residuals get higher effective weight. You can override this by editing the observation weights in the GUI, but there's no direct numerical control in the input file format.
|
||||
|
||||
Range rate (Doppler) observations are supported through a separate observation type in the MPC format. Combining angle and range-rate data requires properly formatted input with both observation types interleaved.
|
||||
</TabItem>
|
||||
<TabItem label="pg_orrery (SQL)">
|
||||
```sql
|
||||
-- Weighted multi-sensor fitting with range rate
|
||||
-- Station 1: precision radar (high weight, has range rate)
|
||||
-- Station 2: optical telescope (lower weight, angles + range only)
|
||||
SELECT (od).iterations,
|
||||
round((od).rms_final::numeric, 4) AS rms_km,
|
||||
(od).status,
|
||||
round((od).condition_number::numeric, 1) AS cond,
|
||||
(od).nstate AS params_fitted
|
||||
FROM tle_from_topocentric(
|
||||
observations := all_topo_obs,
|
||||
times := all_times,
|
||||
observers := ARRAY[
|
||||
'40.0N 105.3W 1655m', -- radar station
|
||||
'34.1N 118.3W 100m' -- optical station
|
||||
]::observer[],
|
||||
observer_ids := station_ids,
|
||||
seed := seed_tle,
|
||||
fit_range_rate := true, -- use Doppler data from radar
|
||||
weights := ARRAY[
|
||||
1.0, 1.0, 1.0, 1.0, -- radar obs: full weight
|
||||
0.3, 0.3, 0.3, 0.3 -- optical obs: lower weight
|
||||
]
|
||||
) AS od;
|
||||
```
|
||||
|
||||
The `weights` array gives explicit numerical control per observation. The `fit_range_rate` flag adds range rate as a 4th residual component — particularly valuable when short observation arcs don't constrain velocity well from position data alone.
|
||||
|
||||
```sql
|
||||
-- Extract formal uncertainty from the covariance matrix
|
||||
-- Diagonal elements are the per-parameter variances
|
||||
WITH fit AS (
|
||||
SELECT (od).covariance AS cov,
|
||||
(od).nstate AS n,
|
||||
(od).condition_number AS cond
|
||||
FROM tle_from_topocentric(/* ... */) AS od
|
||||
)
|
||||
SELECT round(cond::numeric, 1) AS condition_number,
|
||||
round(sqrt(cov[1])::numeric, 6) AS sigma_1,
|
||||
round(sqrt(cov[n + 2])::numeric, 6) AS sigma_2,
|
||||
round(sqrt(cov[2 * n + 3])::numeric, 6) AS sigma_3
|
||||
FROM fit;
|
||||
```
|
||||
|
||||
The covariance matrix is formal `(H^T H)^{-1}` — optimistic, but useful for comparing relative solution quality across different observation geometries or weighting schemes.
|
||||
</TabItem>
|
||||
</Tabs>
|
||||
|
||||
## Residual inspection
|
||||
|
||||
After fitting, you want to know which observations are good and which are outliers.
|
||||
|
||||
<Tabs>
|
||||
<TabItem label="find_orb">
|
||||
find_orb displays a residual map in the GUI — a scatter plot of RA and Dec residuals per observation. Outliers are visually obvious. You can click to exclude individual observations and re-fit.
|
||||
|
||||
In CLI mode, residuals appear in the output file, one line per observation. You parse the file to identify outliers programmatically.
|
||||
</TabItem>
|
||||
<TabItem label="pg_orrery (SQL)">
|
||||
For ECI-based fitting, `tle_fit_residuals()` returns per-observation residuals:
|
||||
|
||||
```sql
|
||||
-- Inspect per-observation residuals after fitting
|
||||
WITH fit AS (
|
||||
SELECT fitted_tle
|
||||
FROM tle_from_eci(obs_positions, obs_times)
|
||||
)
|
||||
SELECT t,
|
||||
round(pos_err_km::numeric, 4) AS error_km,
|
||||
CASE WHEN pos_err_km > 5.0 THEN 'OUTLIER' ELSE '' END AS flag
|
||||
FROM fit, tle_fit_residuals(fit.fitted_tle, obs_positions, obs_times)
|
||||
ORDER BY pos_err_km DESC;
|
||||
```
|
||||
|
||||
For iterative outlier rejection:
|
||||
|
||||
```sql
|
||||
-- Fit, identify outliers, remove them, re-fit
|
||||
WITH first_fit AS (
|
||||
SELECT fitted_tle, positions, times
|
||||
FROM tle_from_eci(obs_positions, obs_times)
|
||||
),
|
||||
good_obs AS (
|
||||
SELECT array_agg(positions[i] ORDER BY i) AS clean_pos,
|
||||
array_agg(times[i] ORDER BY i) AS clean_times
|
||||
FROM first_fit,
|
||||
generate_series(1, array_length(positions, 1)) AS i
|
||||
WHERE (tle_fit_residuals(fitted_tle, positions, times)).pos_err_km < 5.0
|
||||
)
|
||||
SELECT iterations,
|
||||
round(rms_final::numeric, 6) AS rms_km,
|
||||
status
|
||||
FROM good_obs, tle_from_eci(good_obs.clean_pos, good_obs.clean_times);
|
||||
```
|
||||
|
||||
No manual clicking. The rejection criterion is a `WHERE` clause — change the threshold, re-run.
|
||||
</TabItem>
|
||||
</Tabs>
|
||||
|
||||
## Closing the loop: fit → propagate → observe
|
||||
|
||||
find_orb produces orbital elements. To then predict passes, compute rise/set times, or screen for conjunctions, you need a separate tool — GPredict, Skyfield, or a custom propagator. The orbit determination result and the operational tracking live in different systems.
|
||||
|
||||
pg_orrery keeps everything in one place:
|
||||
|
||||
```sql
|
||||
-- Full pipeline: observations → fit → catalog → predict passes
|
||||
WITH fit AS (
|
||||
SELECT tle_from_angles(
|
||||
ra_hours := obs_ra,
|
||||
dec_degrees := obs_dec,
|
||||
times := obs_times,
|
||||
obs := '40.0N 105.3W 1655m'::observer
|
||||
) AS od
|
||||
FROM tonight_observations
|
||||
WHERE target = 'UNKNOWN-2024A'
|
||||
)
|
||||
-- Insert the fitted TLE into the satellite catalog
|
||||
INSERT INTO satellites (norad_id, name, tle, source)
|
||||
SELECT 99999, 'UNKNOWN-2024A', (fit.od).fitted_tle, 'local_od'
|
||||
FROM fit
|
||||
WHERE (fit.od).status = 'converged';
|
||||
|
||||
-- Now predict passes for the newly cataloged object
|
||||
SELECT pass_aos(p) AS rise,
|
||||
pass_max_el(p) AS max_el,
|
||||
pass_los(p) AS set
|
||||
FROM satellites s,
|
||||
predict_passes(s.tle,
|
||||
'40.0N 105.3W 1655m'::observer,
|
||||
now(), now() + interval '48 hours', 10.0) p
|
||||
WHERE s.name = 'UNKNOWN-2024A';
|
||||
```
|
||||
|
||||
The fitted TLE goes into the same catalog table as Space-Track TLEs. Every pg_orrery function — `observe()`, `predict_passes()`, the GiST conjunction index — works on it immediately. No export, no format conversion, no tool switching.
|
||||
|
||||
## Where find_orb wins
|
||||
|
||||
<Aside type="note" title="find_orb is the deeper tool for its domain">
|
||||
find_orb's OD capabilities are far broader than pg_orrery's. The comparison here is specifically about the satellite tracking workflow.
|
||||
</Aside>
|
||||
|
||||
**Heliocentric orbits.** find_orb fits orbits for asteroids, comets, and interplanetary objects — heliocentric Keplerian elements with perturbations. pg_orrery's OD solver fits geocentric SGP4 mean elements for Earth-orbiting satellites only.
|
||||
|
||||
**Perturbation models.** find_orb can include planetary perturbations, solar radiation pressure, and non-gravitational forces in the orbit fit. pg_orrery's DC solver fits pure SGP4 mean elements with optional B* drag.
|
||||
|
||||
**Interactive refinement.** find_orb's GUI lets you visually inspect residuals, toggle observations, try different force models, and watch the solution converge. pg_orrery's solver is a function call — you set parameters and read results.
|
||||
|
||||
**MPC ecosystem.** find_orb reads MPC 80-column format natively and integrates with the Minor Planet Center's observation database. If you're reporting asteroid discoveries, find_orb speaks the language.
|
||||
|
||||
**Multiple solutions.** find_orb can identify and present multiple orbit solutions when the data is ambiguous (e.g., short-arc asteroid observations with two equally valid orbits). pg_orrery's solver returns a single solution.
|
||||
|
||||
**Mature solver.** find_orb's differential corrector has decades of refinement and handles edge cases (near-parabolic orbits, objects at the boundary of Earth's sphere of influence, multi-opposition linkage) that pg_orrery's newer solver does not.
|
||||
|
||||
## Where pg_orrery wins
|
||||
|
||||
**Batch processing.** Fit orbits for every object observed tonight in a single query. find_orb processes objects one at a time.
|
||||
|
||||
**Data integration.** Fitted TLEs land in the same database as your satellite catalog, observation logs, contact schedules, and frequency assignments. `JOIN` the results with anything.
|
||||
|
||||
**Automation.** A SQL query is a complete, repeatable specification. Set it up in a cron job or a materialized view and the pipeline runs itself. find_orb requires either manual GUI interaction or scripted CLI runs with output parsing.
|
||||
|
||||
**Closed-loop tracking.** The fitted TLE immediately works with `observe()`, `predict_passes()`, and the GiST conjunction index. No format conversion, no tool switching, no exporting elements and importing them elsewhere.
|
||||
|
||||
**Multi-observer without MPC codes.** pg_orrery takes observer coordinates directly — no registry, no observatory code, no station lookup file. Define a station as a string and use it.
|
||||
|
||||
**Weighted observations.** pg_orrery's `weights` parameter lets you explicitly down-weight noisy observations or prioritize high-quality data. find_orb handles weighting internally based on observatory statistics.
|
||||
|
||||
## Migrating gradually
|
||||
|
||||
<Steps>
|
||||
1. **Use find_orb for heliocentric objects.** Asteroids, comets, and interplanetary objects belong in find_orb. pg_orrery's OD solver is designed for Earth-orbiting satellites.
|
||||
|
||||
2. **Use pg_orrery for satellite OD.** If your observations are of artificial satellites and you want TLEs that work with SGP4, use `tle_from_angles()` or `tle_from_topocentric()`. The fitted TLE integrates directly with your existing pg_orrery workflow.
|
||||
|
||||
3. **Store observations in PostgreSQL.** Whether you process them with find_orb or pg_orrery, keeping observations in a database table lets you re-process with different parameters, track observation quality over time, and correlate with metadata.
|
||||
|
||||
4. **Automate the pipeline.** Observation ingestion → OD → catalog update → pass prediction can be a scheduled SQL procedure. Manual find_orb runs become the exception for difficult cases, not the default.
|
||||
</Steps>
|
||||
@ -1,418 +0,0 @@
|
||||
---
|
||||
title: From Poliastro to SQL
|
||||
sidebar:
|
||||
order: 6
|
||||
description: Side-by-side comparison of Poliastro Python workflows and equivalent pg_orrery SQL queries for orbital mechanics.
|
||||
---
|
||||
|
||||
import { Tabs, TabItem, Aside, Steps } from "@astrojs/starlight/components";
|
||||
|
||||
[Poliastro](https://docs.poliastro.space/) is a Python library for interactive astrodynamics and orbital mechanics. Built on Astropy, it provides Keplerian propagation, Lambert transfers, orbit plotting, and maneuver analysis with a clean Pythonic API. If you're doing orbital mechanics in Jupyter notebooks, you've probably used it.
|
||||
|
||||
pg_orrery overlaps with Poliastro on two core computations — Lambert transfers and Keplerian propagation — and extends beyond it with SGP4 satellite tracking, orbit determination, and GiST-indexed catalog operations. This page compares the workflows where both tools operate, and identifies where each is the better fit.
|
||||
|
||||
## Lambert transfers
|
||||
|
||||
Both tools implement the Izzo (2015) Lambert solver. The same algorithm, different execution environments.
|
||||
|
||||
<Tabs>
|
||||
<TabItem label="Poliastro (Python)">
|
||||
```python
|
||||
from astropy import units as u
|
||||
from astropy.time import Time
|
||||
from poliastro.bodies import Earth, Mars, Sun
|
||||
from poliastro.ephem import Ephem
|
||||
from poliastro.iod import izzo
|
||||
from poliastro.util import norm
|
||||
|
||||
# Get planet positions from built-in ephemeris
|
||||
dep_time = Time("2028-10-01", scale="tdb")
|
||||
arr_time = Time("2029-06-15", scale="tdb")
|
||||
|
||||
earth_ephem = Ephem.from_body(Earth, dep_time)
|
||||
mars_ephem = Ephem.from_body(Mars, arr_time)
|
||||
|
||||
r_earth = earth_ephem.rv(dep_time)[0] # position vector
|
||||
r_mars = mars_ephem.rv(arr_time)[0] # position vector
|
||||
tof = arr_time - dep_time # time of flight
|
||||
|
||||
# Solve Lambert problem
|
||||
(v_dep, v_arr), = izzo.lambert(Sun.k, r_earth, r_mars, tof)
|
||||
|
||||
# Compute C3
|
||||
v_earth = earth_ephem.rv(dep_time)[1]
|
||||
v_inf = v_dep - v_earth
|
||||
c3 = norm(v_inf).to(u.km / u.s) ** 2
|
||||
|
||||
print(f"Departure C3: {c3:.2f}")
|
||||
```
|
||||
|
||||
Several objects are involved: `Ephem` for planet positions, `Time` for epochs, Astropy units for dimensional analysis, and the `izzo` module for the actual solve. Each intermediate result has attached units. This is excellent for interactive exploration in a notebook — you can inspect each step, change units, and see what you're working with.
|
||||
|
||||
For a single transfer, this is perfectly fine. The friction appears when you need many transfers.
|
||||
</TabItem>
|
||||
<TabItem label="pg_orrery (SQL)">
|
||||
```sql
|
||||
SELECT round(c3_departure::numeric, 2) AS c3_km2s2,
|
||||
round(c3_arrival::numeric, 2) AS c3_arrive,
|
||||
round(v_inf_departure::numeric, 3) AS v_inf_dep_kms,
|
||||
round(tof_days::numeric, 1) AS flight_days,
|
||||
round(transfer_sma::numeric, 4) AS sma_au
|
||||
FROM lambert_transfer(3, 4,
|
||||
'2028-10-01'::timestamptz,
|
||||
'2029-06-15'::timestamptz);
|
||||
```
|
||||
|
||||
One function call. Planet positions, Lambert solve, and C3 computation happen internally. The body IDs (3=Earth, 4=Mars) and timestamps are the only inputs. No unit objects, no intermediate vectors, no ephemeris loading.
|
||||
</TabItem>
|
||||
</Tabs>
|
||||
|
||||
## Pork chop plots
|
||||
|
||||
The parameter sweep that reveals launch windows. This is where the architectural difference matters most.
|
||||
|
||||
<Tabs>
|
||||
<TabItem label="Poliastro (Python)">
|
||||
```python
|
||||
from astropy import units as u
|
||||
from astropy.time import Time
|
||||
import numpy as np
|
||||
from poliastro.bodies import Earth, Mars, Sun
|
||||
from poliastro.ephem import Ephem
|
||||
from poliastro.iod import izzo
|
||||
from poliastro.util import norm
|
||||
|
||||
dep_dates = Time(np.arange(
|
||||
Time("2028-08-01").jd,
|
||||
Time("2029-01-01").jd,
|
||||
1.0 # 1-day steps
|
||||
), format="jd")
|
||||
|
||||
arr_dates = Time(np.arange(
|
||||
Time("2029-03-01").jd,
|
||||
Time("2029-10-01").jd,
|
||||
1.0
|
||||
), format="jd")
|
||||
|
||||
c3_grid = np.full((len(dep_dates), len(arr_dates)), np.nan)
|
||||
|
||||
for i, dep in enumerate(dep_dates):
|
||||
earth = Ephem.from_body(Earth, dep)
|
||||
r_earth, v_earth = earth.rv(dep)
|
||||
|
||||
for j, arr in enumerate(arr_dates):
|
||||
mars = Ephem.from_body(Mars, arr)
|
||||
r_mars = mars.rv(arr)[0]
|
||||
tof = arr - dep
|
||||
|
||||
if tof.to(u.day).value < 90:
|
||||
continue
|
||||
|
||||
try:
|
||||
(v_dep, _), = izzo.lambert(Sun.k, r_earth, r_mars, tof)
|
||||
v_inf = v_dep - v_earth
|
||||
c3_grid[i, j] = norm(v_inf).to(u.km / u.s).value ** 2
|
||||
except Exception:
|
||||
pass
|
||||
|
||||
# Plot with matplotlib
|
||||
import matplotlib.pyplot as plt
|
||||
plt.contourf(arr_dates.datetime, dep_dates.datetime, c3_grid, levels=20)
|
||||
plt.colorbar(label="C3 (km²/s²)")
|
||||
plt.xlabel("Arrival"); plt.ylabel("Departure")
|
||||
plt.show()
|
||||
```
|
||||
|
||||
The nested loop iterates over ~23,000 date combinations. Each iteration constructs ephemeris objects, solves Lambert, handles exceptions, and stores the result. With Poliastro's overhead per iteration (unit conversions, object construction), a dense grid can take minutes.
|
||||
|
||||
The plotting integration is the payoff — `matplotlib` is right there, and the result renders inline in a Jupyter notebook.
|
||||
</TabItem>
|
||||
<TabItem label="pg_orrery (SQL)">
|
||||
```sql
|
||||
-- Same grid: 153 departure x 214 arrival = ~32,700 transfers
|
||||
SELECT dep::date AS departure,
|
||||
arr::date AS arrival,
|
||||
round(c3_departure::numeric, 2) AS c3_km2s2,
|
||||
round(tof_days::numeric, 0) AS flight_days
|
||||
FROM generate_series(
|
||||
'2028-08-01'::timestamptz, '2029-01-01'::timestamptz,
|
||||
interval '1 day') AS dep,
|
||||
generate_series(
|
||||
'2029-03-01'::timestamptz, '2029-10-01'::timestamptz,
|
||||
interval '1 day') AS arr,
|
||||
LATERAL lambert_transfer(3, 4, dep, arr) AS xfer
|
||||
WHERE tof_days > 90
|
||||
AND c3_departure < 50;
|
||||
```
|
||||
|
||||
~32,700 Lambert solves. Runs in seconds with automatic parallelism across cores. Export to CSV for plotting:
|
||||
|
||||
```sql
|
||||
COPY (
|
||||
SELECT dep::date, arr::date,
|
||||
round(c3_departure::numeric, 2) AS c3
|
||||
FROM generate_series(
|
||||
'2028-08-01'::timestamptz, '2029-01-01'::timestamptz,
|
||||
interval '1 day') AS dep,
|
||||
generate_series(
|
||||
'2029-03-01'::timestamptz, '2029-10-01'::timestamptz,
|
||||
interval '1 day') AS arr,
|
||||
LATERAL lambert_transfer(3, 4, dep, arr) AS xfer
|
||||
WHERE tof_days > 90
|
||||
) TO '/tmp/porkchop.csv' WITH CSV HEADER;
|
||||
```
|
||||
|
||||
Then plot with gnuplot, matplotlib, or any contour tool. The computation and visualization are decoupled — pg_orrery produces data, you render it however you prefer.
|
||||
</TabItem>
|
||||
</Tabs>
|
||||
|
||||
<Aside type="tip" title="Same algorithm, different throughput">
|
||||
Both tools use the Izzo (2015) Lambert solver with Householder iterations. The performance difference comes from execution overhead: Poliastro constructs Astropy unit-annotated objects and ephemeris instances per iteration. pg_orrery calls compiled C functions inline with PostgreSQL's batch execution and parallel query engine.
|
||||
</Aside>
|
||||
|
||||
## Keplerian propagation
|
||||
|
||||
Propagating an orbit forward in time using two-body mechanics. This is how you track comets and asteroids from their osculating orbital elements.
|
||||
|
||||
<Tabs>
|
||||
<TabItem label="Poliastro (Python)">
|
||||
```python
|
||||
from astropy import units as u
|
||||
from astropy.time import Time
|
||||
from poliastro.bodies import Sun
|
||||
from poliastro.twobody import Orbit
|
||||
import numpy as np
|
||||
|
||||
# Define a comet orbit from classical elements
|
||||
comet = Orbit.from_classical(
|
||||
Sun,
|
||||
a=2.7 * u.AU,
|
||||
ecc=0.85 * u.one,
|
||||
inc=12.0 * u.deg,
|
||||
raan=65.0 * u.deg,
|
||||
argp=180.0 * u.deg,
|
||||
nu=0.0 * u.deg,
|
||||
epoch=Time("2025-01-15")
|
||||
)
|
||||
|
||||
# Propagate forward 90 days
|
||||
future = comet.propagate(90 * u.day)
|
||||
r = future.r.to(u.AU)
|
||||
print(f"Position: {r}")
|
||||
|
||||
# Time series
|
||||
times = np.linspace(0, 365, 100) * u.day
|
||||
positions = [comet.propagate(t).r.to(u.AU).value for t in times]
|
||||
```
|
||||
|
||||
Poliastro's `Orbit` object handles the Kepler equation internally. Propagation returns a new `Orbit` with updated elements. The unit-annotated result is self-documenting.
|
||||
|
||||
The time series loop is explicit — each propagation creates a new object.
|
||||
</TabItem>
|
||||
<TabItem label="pg_orrery (SQL)">
|
||||
```sql
|
||||
-- Propagate a comet orbit and observe from Earth
|
||||
SELECT t,
|
||||
topo_azimuth(obs) AS az,
|
||||
topo_elevation(obs) AS el,
|
||||
topo_range(obs) / 149597870.7 AS range_au
|
||||
FROM generate_series(
|
||||
'2025-01-15'::timestamptz,
|
||||
'2026-01-15'::timestamptz,
|
||||
interval '1 day'
|
||||
) AS t,
|
||||
LATERAL comet_observe(
|
||||
2.7, -- semi-major axis (AU)
|
||||
0.85, -- eccentricity
|
||||
12.0, -- inclination (deg)
|
||||
65.0, -- RAAN (deg)
|
||||
180.0, -- argument of perihelion (deg)
|
||||
0.0, -- mean anomaly (deg)
|
||||
'2025-01-15'::timestamptz, -- epoch
|
||||
'40.0N 105.3W 1655m'::observer,
|
||||
t
|
||||
) AS obs;
|
||||
```
|
||||
|
||||
`comet_observe()` wraps Keplerian propagation with the topocentric observation pipeline. The time series is `generate_series` — one row per timestep, parallel-safe, no loop management.
|
||||
|
||||
For raw heliocentric position without the observation step:
|
||||
|
||||
```sql
|
||||
SELECT t,
|
||||
helio_x(kepler_propagate(2.7, 0.85, 12.0, 65.0, 180.0, 0.0,
|
||||
'2025-01-15'::timestamptz, t)) AS x_au,
|
||||
helio_y(kepler_propagate(2.7, 0.85, 12.0, 65.0, 180.0, 0.0,
|
||||
'2025-01-15'::timestamptz, t)) AS y_au,
|
||||
helio_z(kepler_propagate(2.7, 0.85, 12.0, 65.0, 180.0, 0.0,
|
||||
'2025-01-15'::timestamptz, t)) AS z_au
|
||||
FROM generate_series(
|
||||
'2025-01-15'::timestamptz,
|
||||
'2026-01-15'::timestamptz,
|
||||
interval '1 day') AS t;
|
||||
```
|
||||
</TabItem>
|
||||
</Tabs>
|
||||
|
||||
## What Poliastro does that pg_orrery doesn't
|
||||
|
||||
Poliastro is a broader orbital mechanics toolkit. Several of its capabilities have no pg_orrery equivalent.
|
||||
|
||||
### Orbit plotting
|
||||
|
||||
```python
|
||||
from poliastro.plotting import OrbitPlotter3D
|
||||
|
||||
plotter = OrbitPlotter3D()
|
||||
plotter.plot(earth_orbit, label="Earth")
|
||||
plotter.plot(mars_orbit, label="Mars")
|
||||
plotter.plot(transfer, label="Transfer")
|
||||
plotter.show()
|
||||
```
|
||||
|
||||
Interactive 3D orbit visualization in Jupyter notebooks. pg_orrery returns coordinates — you build your own visualization externally.
|
||||
|
||||
### Maneuver planning
|
||||
|
||||
```python
|
||||
from poliastro.maneuver import Maneuver
|
||||
|
||||
hohmann = Maneuver.hohmann(initial_orbit, target_radius)
|
||||
bielliptic = Maneuver.bielliptic(initial_orbit, r_b, target_radius)
|
||||
|
||||
print(f"Hohmann delta-v: {hohmann.get_total_cost()}")
|
||||
print(f"Bielliptic delta-v: {bielliptic.get_total_cost()}")
|
||||
```
|
||||
|
||||
Poliastro computes delta-v budgets for standard maneuvers. pg_orrery's Lambert solver gives you the transfer orbit but doesn't model the departure/arrival maneuvers.
|
||||
|
||||
### Perturbation models
|
||||
|
||||
```python
|
||||
from poliastro.twobody.propagation import CowellPropagator
|
||||
from poliastro.twobody.perturbations import J2_perturbation
|
||||
|
||||
orbit_with_j2 = orbit.propagate(
|
||||
30 * u.day,
|
||||
method=CowellPropagator(f=J2_perturbation)
|
||||
)
|
||||
```
|
||||
|
||||
Poliastro supports numerical integration with custom perturbation functions — J2 oblateness, atmospheric drag, third-body gravity, solar radiation pressure. pg_orrery's Keplerian propagation is strictly two-body.
|
||||
|
||||
### Astropy integration
|
||||
|
||||
Poliastro builds on Astropy's unit system, time handling, and coordinate frames. Results carry physical units and can be converted freely. pg_orrery returns bare doubles — you know the units from the documentation, but the database doesn't enforce them.
|
||||
|
||||
## What pg_orrery does that Poliastro doesn't
|
||||
|
||||
### SGP4/SDP4 satellite tracking
|
||||
|
||||
```sql
|
||||
-- Observe 12,000 satellites in 17ms
|
||||
SELECT s.name,
|
||||
topo_elevation(observe(s.tle, '40.0N 105.3W 1655m'::observer, now())) AS el
|
||||
FROM satellites s
|
||||
WHERE topo_elevation(observe(s.tle, '40.0N 105.3W 1655m'::observer, now())) > 0;
|
||||
```
|
||||
|
||||
Poliastro has no TLE support, no SGP4, no satellite catalog operations. The entire satellite tracking domain — pass prediction, conjunction screening, GiST indexing — is absent.
|
||||
|
||||
### Orbit determination
|
||||
|
||||
```sql
|
||||
-- Fit an orbit from observations
|
||||
SELECT iterations, round(rms_final::numeric, 6) AS rms, status, fitted_tle
|
||||
FROM tle_from_angles(
|
||||
ra_hours := obs_ra, dec_degrees := obs_dec,
|
||||
times := obs_times, obs := station
|
||||
);
|
||||
```
|
||||
|
||||
Poliastro does not fit orbits from observations. pg_orrery's differential correction solver handles ECI, topocentric, and angles-only OD with Gauss/Gibbs IOD bootstrap.
|
||||
|
||||
### Data correlation via SQL
|
||||
|
||||
```sql
|
||||
-- Combine Lambert transfer analysis with mission constraint database
|
||||
SELECT t.target, t.dep::date, t.c3,
|
||||
m.max_c3_budget, m.launch_vehicle,
|
||||
CASE WHEN t.c3 < m.max_c3_budget THEN 'FEASIBLE' ELSE 'OVER BUDGET' END
|
||||
FROM transfer_survey t
|
||||
JOIN mission_constraints m ON m.target = t.target
|
||||
WHERE t.c3 < m.max_c3_budget * 1.1;
|
||||
```
|
||||
|
||||
Poliastro results live in Python memory. Correlating with mission databases, launch vehicle catalogs, or historical data requires loading everything into Python. pg_orrery results are database rows — correlation is a JOIN.
|
||||
|
||||
### Automatic parallelism
|
||||
|
||||
All pg_orrery functions are `PARALLEL SAFE`. PostgreSQL distributes work across CPU cores automatically. Poliastro runs single-threaded by default — parallelism requires explicit `multiprocessing` or `joblib` setup.
|
||||
|
||||
## Where Poliastro wins
|
||||
|
||||
<Aside type="note" title="Different tools for different phases">
|
||||
Poliastro excels at interactive exploration and visualization. pg_orrery excels at batch computation and data integration. They address different phases of the same workflow.
|
||||
</Aside>
|
||||
|
||||
**Interactive exploration.** Poliastro in a Jupyter notebook is the best way to explore orbital mechanics interactively. Plot orbits, change parameters, see results immediately. pg_orrery returns tabular data — useful, but not visual.
|
||||
|
||||
**Unit safety.** Astropy units catch dimension errors at runtime. Passing kilometers where AU are expected raises an exception. pg_orrery uses bare doubles — unit errors are silent.
|
||||
|
||||
**Perturbation modeling.** Poliastro's numerical integrators with custom force models handle real-world orbit analysis that two-body Keplerian propagation cannot.
|
||||
|
||||
**Maneuver analysis.** Delta-v budgets, Hohmann transfers, bielliptic maneuvers — Poliastro models the spacecraft operations side. pg_orrery provides the trajectory geometry.
|
||||
|
||||
**Orbit visualization.** 2D and 3D orbit plots, ground tracks on map projections, interactive widgets in notebooks. pg_orrery produces numbers for external rendering.
|
||||
|
||||
**Broader ephemeris support.** Poliastro accesses any Astropy-compatible ephemeris source, including custom SPK kernels for spacecraft and small bodies.
|
||||
|
||||
## Where pg_orrery wins
|
||||
|
||||
**Throughput.** 22,500 Lambert solves in 8.3 seconds vs. minutes in Poliastro for the same grid. The difference is compiled C + parallel execution vs. Python object construction per iteration.
|
||||
|
||||
**Satellite operations.** SGP4, TLEs, pass prediction, conjunction screening, orbit determination — Poliastro has none of these. If you work with artificial satellites, pg_orrery covers the domain.
|
||||
|
||||
**Data integration.** Results live in PostgreSQL alongside your other data. JOIN with mission databases, observation logs, or any other table.
|
||||
|
||||
**Reproducibility.** A SQL query is a complete specification. No virtual environment, no package versions, no imports. Share the query and it runs identically on any pg_orrery installation.
|
||||
|
||||
**Deployment simplicity.** pg_orrery is a PostgreSQL extension — install once, available to every client. Poliastro requires a Python environment with Astropy, NumPy, SciPy, and their transitive dependencies.
|
||||
|
||||
## A combined workflow
|
||||
|
||||
The two tools complement each other naturally:
|
||||
|
||||
<Steps>
|
||||
1. **Explore with Poliastro.** Use a Jupyter notebook to understand the problem — plot orbits, visualize transfer geometry, prototype maneuver sequences. This is the design phase.
|
||||
|
||||
2. **Survey with pg_orrery.** Once you know what you're looking for, run the parameter sweep in SQL. 100,000 Lambert solves to identify the best launch windows. Store the results in a table.
|
||||
|
||||
3. **Correlate with pg_orrery.** JOIN transfer results with launch vehicle constraints, budget timelines, or historical mission data. Filter to feasible solutions.
|
||||
|
||||
4. **Refine with Poliastro.** Take the best candidates back to Python for detailed analysis — add perturbations, model the departure spiral, compute the precise delta-v budget.
|
||||
|
||||
5. **Operate with pg_orrery.** Once you have TLEs (from catalog data or from pg_orrery's OD solver), the operational tracking — observation, pass prediction, conjunction screening — stays in SQL.
|
||||
</Steps>
|
||||
|
||||
<Aside type="tip" title="Connecting the two">
|
||||
Export pg_orrery results to CSV or JSON for import into Poliastro notebooks:
|
||||
|
||||
```sql
|
||||
COPY (
|
||||
SELECT dep::date, arr::date,
|
||||
c3_departure, c3_arrival, tof_days, transfer_sma
|
||||
FROM lambert_transfer(3, 4, dep, arr)
|
||||
-- ... your survey grid ...
|
||||
) TO '/tmp/survey.csv' WITH CSV HEADER;
|
||||
```
|
||||
|
||||
```python
|
||||
import pandas as pd
|
||||
survey = pd.read_csv('/tmp/survey.csv')
|
||||
best = survey.nsmallest(10, 'c3_departure')
|
||||
# Now explore the best candidates in Poliastro
|
||||
```
|
||||
|
||||
The boundary between "bulk computation" and "interactive analysis" is a CSV file.
|
||||
</Aside>
|
||||
@ -273,7 +273,7 @@ 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.** Skyfield uses the full IAU 2000A nutation model, polar motion corrections, delta-T from IERS data, and iterates for light-time and stellar aberration. pg_orrery v0.3.0 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.
|
||||
|
||||
**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.
|
||||
|
||||
@ -295,42 +295,6 @@ pg_orrery does not replace Skyfield for all use cases. Be clear about where the
|
||||
|
||||
**Reproducibility.** A SQL query is a complete, self-contained specification of a computation. No virtual environment, no package versions, no file paths. The same query produces the same result on any PostgreSQL instance with pg_orrery installed.
|
||||
|
||||
**Orbit determination.** Skyfield is a propagation and observation library — it does not fit orbits from observations. Since v0.4.0, pg_orrery can [fit TLEs from ECI, topocentric, or angles-only observations](/guides/orbit-determination/) via differential correction, entirely in SQL. This closes the loop: observe a satellite, fit its orbit, and propagate the result — all within the database.
|
||||
|
||||
```sql
|
||||
-- The closed-loop workflow Skyfield can't do:
|
||||
-- observe → fit (with range rate + weighted obs) → catalog → predict
|
||||
WITH obs_data AS (
|
||||
-- Collect observations from tonight's tracking pass
|
||||
SELECT array_agg(observe(seed_tle, station, t) ORDER BY t) AS topo_obs,
|
||||
array_agg(t ORDER BY t) AS times
|
||||
FROM generate_series(
|
||||
'2024-01-01 12:00+00'::timestamptz,
|
||||
'2024-01-01 12:10+00'::timestamptz,
|
||||
interval '30 seconds') t
|
||||
),
|
||||
fit AS (
|
||||
SELECT fitted_tle, rms_final, condition_number, status
|
||||
FROM obs_data,
|
||||
tle_from_topocentric(
|
||||
obs_data.topo_obs, obs_data.times,
|
||||
'40.0N 105.3W 1655m'::observer,
|
||||
seed := seed_tle,
|
||||
fit_range_rate := true, -- v0.6.0: use Doppler data
|
||||
weights := obs_weights -- v0.6.0: per-obs weighting
|
||||
)
|
||||
WHERE status = 'converged'
|
||||
)
|
||||
-- Immediately predict passes with the freshly-fitted TLE
|
||||
SELECT pass_aos(p) AS rise,
|
||||
pass_max_el(p) AS max_el,
|
||||
pass_los(p) AS set
|
||||
FROM fit,
|
||||
predict_passes(fit.fitted_tle,
|
||||
'40.0N 105.3W 1655m'::observer,
|
||||
now(), now() + interval '48 hours', 10.0) p;
|
||||
```
|
||||
|
||||
## Migrating gradually
|
||||
|
||||
You don't have to choose one or the other. A practical migration path:
|
||||
@ -344,7 +308,5 @@ You don't have to choose one or the other. A practical migration path:
|
||||
|
||||
4. **Move reporting to SQL.** "What was above 20 degrees from each of our 5 observers last night?" is a single query with a CROSS JOIN, not a Python loop over observers and timestamps.
|
||||
|
||||
5. **Move orbit determination to SQL.** If you're fitting orbits with separate Python tools (find_orb, OREKit bindings, custom scripts), pg_orrery's [OD solver](/guides/orbit-determination/) can do it alongside your existing satellite catalog. Observations in, TLE out — no data pipeline between Python and PostgreSQL.
|
||||
|
||||
6. **Enable DE when accuracy matters.** If you find VSOP87's ~1 arcsecond isn't enough for a specific use case, [configure a DE file](/guides/de-ephemeris/) and switch to `_de()` function variants. Same SQL patterns, same parameters — just add `_de` to the function name.
|
||||
5. **Enable DE when accuracy matters.** If you find VSOP87's ~1 arcsecond isn't enough for a specific use case, [configure a DE file](/guides/de-ephemeris/) and switch to `_de()` function variants. Same SQL patterns, same parameters — just add `_de` to the function name.
|
||||
</Steps>
|
||||
|
||||
@ -344,7 +344,7 @@ This is a screening filter, not a precision conjunction analysis. It identifies
|
||||
|
||||
## Provider switching: accuracy when you need it
|
||||
|
||||
pg_orrery has two ephemeris providers (since v0.3.0) — the built-in VSOP87 pipeline (~1 arcsecond) and optional [JPL DE440/441](/guides/de-ephemeris/) (~0.1 milliarcsecond). The SQL interface makes switching between them a one-character change.
|
||||
pg_orrery v0.3.0 has two ephemeris providers — the built-in VSOP87 pipeline (~1 arcsecond) and optional [JPL DE440/441](/guides/de-ephemeris/) (~0.1 milliarcsecond). The SQL interface makes switching between them a one-character change.
|
||||
|
||||
```sql
|
||||
-- VSOP87 (built-in, IMMUTABLE, no setup)
|
||||
@ -437,104 +437,6 @@ ORDER BY night;
|
||||
|
||||
For each night, this finds which planet reaches the highest elevation — useful for deciding what to observe. The window function handles the ranking within each night without a correlated subquery.
|
||||
|
||||
## Orbit determination: observations in, TLE out
|
||||
|
||||
Since v0.4.0, pg_orrery can fit TLE orbital elements from observations using differential correction. This closes a loop that previously required external tools — you can observe, fit, and track entirely within SQL.
|
||||
|
||||
### Seed-free fitting with Gauss IOD
|
||||
|
||||
The Gauss method (v0.6.0) bootstraps an initial orbit from three angles-only observations, eliminating the need for a seed TLE from a catalog:
|
||||
|
||||
```sql
|
||||
-- Unknown satellite: RA/Dec observations from last night
|
||||
-- No seed TLE needed — Gauss IOD handles the initial guess
|
||||
SELECT iterations,
|
||||
round(rms_final::numeric, 6) AS rms_rad,
|
||||
status,
|
||||
round(condition_number::numeric, 1) AS cond,
|
||||
fitted_tle
|
||||
FROM tle_from_angles(
|
||||
ra_hours := ARRAY[12.345, 12.567, 12.789, 13.012, 13.234, 13.456],
|
||||
dec_degrees := ARRAY[45.1, 44.8, 44.3, 43.6, 42.8, 41.9],
|
||||
times := ARRAY[
|
||||
'2024-01-01 12:00+00', '2024-01-01 12:01+00',
|
||||
'2024-01-01 12:02+00', '2024-01-01 12:03+00',
|
||||
'2024-01-01 12:04+00', '2024-01-01 12:05+00'
|
||||
]::timestamptz[],
|
||||
obs := '40.0N 105.3W 1655m'::observer
|
||||
);
|
||||
```
|
||||
|
||||
The fitted TLE is immediately usable with `observe()`, `predict_passes()`, and the GiST conjunction index. No export, no format conversion.
|
||||
|
||||
### Weighted multi-sensor fusion
|
||||
|
||||
When observations come from sensors with different accuracies, the `weights` parameter (v0.6.0) gives explicit control over each observation's influence:
|
||||
|
||||
```sql
|
||||
-- Radar + optical: weight radar higher, include Doppler data
|
||||
SELECT iterations,
|
||||
round(rms_final::numeric, 4) AS rms_km,
|
||||
status,
|
||||
nstate AS params_fitted
|
||||
FROM tle_from_topocentric(
|
||||
observations := all_topo_obs,
|
||||
times := all_times,
|
||||
observers := ARRAY[
|
||||
'40.0N 105.3W 1655m', -- radar
|
||||
'34.1N 118.3W 100m' -- optical
|
||||
]::observer[],
|
||||
observer_ids := station_ids,
|
||||
seed := seed_tle,
|
||||
fit_range_rate := true, -- include Doppler from radar
|
||||
weights := ARRAY[
|
||||
1.0, 1.0, 1.0, 1.0, -- radar: full weight
|
||||
0.3, 0.3, 0.3, 0.3 -- optical: down-weighted
|
||||
]
|
||||
);
|
||||
```
|
||||
|
||||
### The full pipeline as SQL
|
||||
|
||||
The OD pattern composes with every other pattern on this page:
|
||||
|
||||
```sql
|
||||
-- Full pipeline: batch OD → catalog insert → pass prediction
|
||||
-- Process all objects observed tonight
|
||||
WITH fits AS (
|
||||
SELECT obj_id,
|
||||
tle_from_angles(
|
||||
ra_hours := array_agg(ra ORDER BY obs_time),
|
||||
dec_degrees := array_agg(dec ORDER BY obs_time),
|
||||
times := array_agg(obs_time ORDER BY obs_time),
|
||||
obs := '40.0N 105.3W 1655m'::observer
|
||||
) AS od
|
||||
FROM tonight_observations
|
||||
GROUP BY obj_id
|
||||
HAVING count(*) >= 3
|
||||
),
|
||||
new_tles AS (
|
||||
INSERT INTO satellites (norad_id, name, tle, source)
|
||||
SELECT 90000 + row_number() OVER (),
|
||||
obj_id,
|
||||
(od).fitted_tle,
|
||||
'local_od'
|
||||
FROM fits
|
||||
WHERE (od).status = 'converged'
|
||||
RETURNING norad_id, tle
|
||||
)
|
||||
-- Immediately predict passes for every newly cataloged object
|
||||
SELECT s.norad_id,
|
||||
pass_aos(p) AS rise,
|
||||
pass_max_el(p) AS max_el
|
||||
FROM new_tles s,
|
||||
predict_passes(s.tle, '40.0N 105.3W 1655m'::observer,
|
||||
now(), now() + interval '48 hours', 10.0) p
|
||||
ORDER BY pass_aos(p);
|
||||
```
|
||||
|
||||
Observation ingestion, orbit fitting, catalog insertion, and pass prediction — one statement. Schedule it with `pg_cron` and the pipeline runs itself.
|
||||
|
||||
## Composition: building complex queries from simple parts
|
||||
|
||||
The real power of SQL is that these patterns compose. A single query can use `generate_series` for time steps, `CROSS JOIN` for parameter sweeps, `JOIN` for data correlation, window functions for change detection, and `COPY TO` for export — all in one statement.
|
||||
@ -597,4 +499,4 @@ This query:
|
||||
|
||||
In a traditional workflow, each of these steps would be a separate script, a separate data file, and a separate tool. In SQL, they compose into a single declarative statement that the database engine optimizes and parallelizes.
|
||||
|
||||
That's the advantage. Not that SQL is a better programming language — it isn't. But for the specific pattern of "evaluate a function over structured parameter spaces and correlate the results with existing data," SQL is exactly the right tool. pg_orrery puts over 100 functions inside that tool — from 17ms satellite batch propagation to sub-milliarcsecond DE441 planet positions — and every one of them composes with every SQL pattern on this page.
|
||||
That's the advantage. Not that SQL is a better programming language — it isn't. But for the specific pattern of "evaluate a function over structured parameter spaces and correlate the results with existing data," SQL is exactly the right tool. pg_orrery puts 68 functions inside that tool — from 17ms satellite batch propagation to sub-milliarcsecond DE441 planet positions — and every one of them composes with every SQL pattern on this page.
|
||||
|
||||
@ -1,4 +1,4 @@
|
||||
comment = 'A database orrery — celestial mechanics types and functions for PostgreSQL'
|
||||
default_version = '0.6.0'
|
||||
default_version = '0.5.0'
|
||||
module_pathname = '$libdir/pg_orrery'
|
||||
relocatable = true
|
||||
|
||||
@ -1,114 +0,0 @@
|
||||
-- pg_orrery 0.5.0 -> 0.6.0 migration
|
||||
--
|
||||
-- Adds range rate fitting, per-observation weights, and
|
||||
-- angles-only orbit determination (Gauss method).
|
||||
--
|
||||
-- Range rate and weights change the input signatures of
|
||||
-- tle_from_eci and tle_from_topocentric (added trailing
|
||||
-- DEFAULT parameters), which requires DROP + re-CREATE.
|
||||
|
||||
-- ============================================================
|
||||
-- Drop old OD function signatures
|
||||
-- ============================================================
|
||||
|
||||
DROP FUNCTION IF EXISTS tle_from_eci(eci_position[], timestamptz[], tle, boolean, int4);
|
||||
DROP FUNCTION IF EXISTS tle_from_topocentric(topocentric[], timestamptz[], observer, tle, boolean, int4);
|
||||
DROP FUNCTION IF EXISTS tle_from_topocentric(topocentric[], timestamptz[], observer[], int4[], tle, boolean, int4);
|
||||
|
||||
-- ============================================================
|
||||
-- Re-create with range_rate + weights parameters
|
||||
-- ============================================================
|
||||
|
||||
-- Fit TLE from ECI position/velocity ephemeris
|
||||
-- weights: per-observation weighting (NULL = uniform)
|
||||
|
||||
CREATE FUNCTION tle_from_eci(
|
||||
positions eci_position[], times timestamptz[],
|
||||
seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false,
|
||||
max_iter int4 DEFAULT 15,
|
||||
weights float8[] DEFAULT NULL,
|
||||
OUT fitted_tle tle, OUT iterations int4,
|
||||
OUT rms_final float8, OUT rms_initial float8, OUT status text,
|
||||
OUT condition_number float8, OUT covariance float8[], OUT nstate int4
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_from_eci(eci_position[], timestamptz[], tle, boolean, int4, float8[]) IS
|
||||
'Fit a TLE from ECI position/velocity observations via differential correction. Optional per-observation weights for heterogeneous sensor fusion. Returns fitted TLE, RMS residuals, convergence status, condition number, and formal covariance matrix.';
|
||||
|
||||
-- Fit TLE from topocentric observations (az/el/range) — single observer
|
||||
-- fit_range_rate: include range_rate as 4th residual component
|
||||
-- weights: per-observation weighting (NULL = uniform)
|
||||
|
||||
CREATE FUNCTION tle_from_topocentric(
|
||||
observations topocentric[], times timestamptz[],
|
||||
obs observer,
|
||||
seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false,
|
||||
max_iter int4 DEFAULT 15,
|
||||
fit_range_rate boolean DEFAULT false,
|
||||
weights float8[] DEFAULT NULL,
|
||||
OUT fitted_tle tle, OUT iterations int4,
|
||||
OUT rms_final float8, OUT rms_initial float8, OUT status text,
|
||||
OUT condition_number float8, OUT covariance float8[], OUT nstate int4
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_from_topocentric(topocentric[], timestamptz[], observer, tle, boolean, int4, boolean, float8[]) IS
|
||||
'Fit a TLE from topocentric (az/el/range) observations via differential correction. Optional range_rate fitting and per-observation weights. Returns fitted TLE, RMS residuals, convergence status, condition number, and formal covariance matrix.';
|
||||
|
||||
-- Fit TLE from topocentric observations — multiple observers
|
||||
|
||||
CREATE FUNCTION tle_from_topocentric(
|
||||
observations topocentric[], times timestamptz[],
|
||||
observers observer[], observer_ids int4[],
|
||||
seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false,
|
||||
max_iter int4 DEFAULT 15,
|
||||
fit_range_rate boolean DEFAULT false,
|
||||
weights float8[] DEFAULT NULL,
|
||||
OUT fitted_tle tle, OUT iterations int4,
|
||||
OUT rms_final float8, OUT rms_initial float8, OUT status text,
|
||||
OUT condition_number float8, OUT covariance float8[], OUT nstate int4
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME', 'tle_from_topocentric_multi'
|
||||
LANGUAGE C STABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_from_topocentric(topocentric[], timestamptz[], observer[], int4[], tle, boolean, int4, boolean, float8[]) IS
|
||||
'Fit a TLE from topocentric observations collected by multiple ground stations. observer_ids[i] indexes into observers[]. Optional range_rate fitting and per-observation weights.';
|
||||
|
||||
-- ============================================================
|
||||
-- Angles-only orbit determination (new in v0.6.0)
|
||||
-- ============================================================
|
||||
|
||||
-- Fit TLE from RA/Dec observations — single observer
|
||||
-- Uses Gauss method for initial orbit determination when no seed is provided.
|
||||
-- RA in hours [0,24), Dec in degrees [-90,90] (matches star_observe convention).
|
||||
-- RMS output is in radians for angles-only mode.
|
||||
|
||||
CREATE FUNCTION tle_from_angles(
|
||||
ra_hours float8[], dec_degrees float8[],
|
||||
times timestamptz[],
|
||||
obs observer,
|
||||
seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false,
|
||||
max_iter int4 DEFAULT 15,
|
||||
weights float8[] DEFAULT NULL,
|
||||
OUT fitted_tle tle, OUT iterations int4,
|
||||
OUT rms_final float8, OUT rms_initial float8, OUT status text,
|
||||
OUT condition_number float8, OUT covariance float8[], OUT nstate int4
|
||||
) RETURNS RECORD AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_from_angles(float8[], float8[], timestamptz[], observer, tle, boolean, int4, float8[]) IS
|
||||
'Fit a TLE from angles-only (RA/Dec) observations via Gauss IOD + differential correction. RA in hours [0,24), Dec in degrees [-90,90]. RMS output in radians. Uses Gauss method for seed-free initial guess.';
|
||||
|
||||
-- Fit TLE from RA/Dec observations — multiple observers
|
||||
|
||||
CREATE FUNCTION tle_from_angles(
|
||||
ra_hours float8[], dec_degrees float8[],
|
||||
times timestamptz[],
|
||||
observers observer[], observer_ids int4[],
|
||||
seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false,
|
||||
max_iter int4 DEFAULT 15,
|
||||
weights float8[] DEFAULT NULL,
|
||||
OUT fitted_tle tle, OUT iterations int4,
|
||||
OUT rms_final float8, OUT rms_initial float8, OUT status text,
|
||||
OUT condition_number float8, OUT covariance float8[], OUT nstate int4
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME', 'tle_from_angles_multi'
|
||||
LANGUAGE C STABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_from_angles(float8[], float8[], timestamptz[], observer[], int4[], tle, boolean, int4, float8[]) IS
|
||||
'Fit a TLE from angles-only (RA/Dec) observations from multiple ground stations via Gauss IOD + differential correction.';
|
||||
@ -1,885 +0,0 @@
|
||||
-- pg_orrery -- Orbital mechanics types and functions for PostgreSQL
|
||||
--
|
||||
-- Types: tle, eci_position, geodetic, topocentric, observer, pass_event
|
||||
-- Provides SGP4/SDP4 propagation, coordinate transforms, pass prediction,
|
||||
-- and GiST indexing on altitude bands for conjunction screening.
|
||||
--
|
||||
-- All propagation uses WGS-72 constants (matching TLE mean element fitting).
|
||||
-- Coordinate output uses WGS-84 (matching modern geodetic standards).
|
||||
|
||||
-- ============================================================
|
||||
-- Shell types (forward declarations)
|
||||
-- ============================================================
|
||||
|
||||
CREATE TYPE tle;
|
||||
CREATE TYPE eci_position;
|
||||
CREATE TYPE geodetic;
|
||||
CREATE TYPE topocentric;
|
||||
CREATE TYPE observer;
|
||||
CREATE TYPE pass_event;
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- TLE type: Two-Line Element set
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION tle_in(cstring) RETURNS tle
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION tle_out(tle) RETURNS cstring
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION tle_recv(internal) RETURNS tle
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION tle_send(tle) RETURNS bytea
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE TYPE tle (
|
||||
INPUT = tle_in,
|
||||
OUTPUT = tle_out,
|
||||
RECEIVE = tle_recv,
|
||||
SEND = tle_send,
|
||||
INTERNALLENGTH = 112,
|
||||
ALIGNMENT = double,
|
||||
STORAGE = plain
|
||||
);
|
||||
|
||||
COMMENT ON TYPE tle IS 'Two-Line Element set — parsed mean orbital elements for SGP4/SDP4 propagation';
|
||||
|
||||
-- TLE accessor functions
|
||||
|
||||
CREATE FUNCTION tle_epoch(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_epoch(tle) IS 'TLE epoch as Julian date (UTC)';
|
||||
|
||||
CREATE FUNCTION tle_norad_id(tle) RETURNS int4
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_norad_id(tle) IS 'NORAD catalog number';
|
||||
|
||||
CREATE FUNCTION tle_inclination(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_inclination(tle) IS 'Orbital inclination in degrees';
|
||||
|
||||
CREATE FUNCTION tle_eccentricity(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_eccentricity(tle) IS 'Orbital eccentricity (dimensionless)';
|
||||
|
||||
CREATE FUNCTION tle_raan(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_raan(tle) IS 'Right ascension of ascending node in degrees';
|
||||
|
||||
CREATE FUNCTION tle_arg_perigee(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_arg_perigee(tle) IS 'Argument of perigee in degrees';
|
||||
|
||||
CREATE FUNCTION tle_mean_anomaly(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_mean_anomaly(tle) IS 'Mean anomaly in degrees';
|
||||
|
||||
CREATE FUNCTION tle_mean_motion(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_mean_motion(tle) IS 'Mean motion in revolutions per day';
|
||||
|
||||
CREATE FUNCTION tle_bstar(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_bstar(tle) IS 'B* drag coefficient (1/earth-radii)';
|
||||
|
||||
CREATE FUNCTION tle_period(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_period(tle) IS 'Orbital period in minutes';
|
||||
|
||||
CREATE FUNCTION tle_age(tle, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_age(tle, timestamptz) IS 'TLE age in days (positive = past epoch, negative = before epoch)';
|
||||
|
||||
CREATE FUNCTION tle_perigee(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_perigee(tle) IS 'Perigee altitude in km above WGS-72 ellipsoid';
|
||||
|
||||
CREATE FUNCTION tle_apogee(tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_apogee(tle) IS 'Apogee altitude in km above WGS-72 ellipsoid';
|
||||
|
||||
CREATE FUNCTION tle_intl_desig(tle) RETURNS text
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_intl_desig(tle) IS 'International designator (COSPAR ID)';
|
||||
|
||||
CREATE FUNCTION tle_from_lines(text, text) RETURNS tle
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_from_lines(text, text) IS
|
||||
'Construct TLE from separate line1/line2 text columns';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- ECI position type: True Equator Mean Equinox (TEME) frame
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION eci_in(cstring) RETURNS eci_position
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_out(eci_position) RETURNS cstring
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_recv(internal) RETURNS eci_position
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_send(eci_position) RETURNS bytea
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE TYPE eci_position (
|
||||
INPUT = eci_in,
|
||||
OUTPUT = eci_out,
|
||||
RECEIVE = eci_recv,
|
||||
SEND = eci_send,
|
||||
INTERNALLENGTH = 48,
|
||||
ALIGNMENT = double,
|
||||
STORAGE = plain
|
||||
);
|
||||
|
||||
COMMENT ON TYPE eci_position IS 'Earth-Centered Inertial position and velocity in TEME frame (km, km/s)';
|
||||
|
||||
-- ECI accessor functions
|
||||
|
||||
CREATE FUNCTION eci_x(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_y(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_z(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_vx(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_vy(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_vz(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION eci_speed(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION eci_speed(eci_position) IS 'Velocity magnitude in km/s';
|
||||
|
||||
CREATE FUNCTION eci_altitude(eci_position) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION eci_altitude(eci_position) IS 'Approximate geocentric altitude in km (radius - WGS72_AE)';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Geodetic type: WGS-84 latitude/longitude/altitude
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION geodetic_in(cstring) RETURNS geodetic
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION geodetic_out(geodetic) RETURNS cstring
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION geodetic_recv(internal) RETURNS geodetic
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION geodetic_send(geodetic) RETURNS bytea
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE TYPE geodetic (
|
||||
INPUT = geodetic_in,
|
||||
OUTPUT = geodetic_out,
|
||||
RECEIVE = geodetic_recv,
|
||||
SEND = geodetic_send,
|
||||
INTERNALLENGTH = 24,
|
||||
ALIGNMENT = double,
|
||||
STORAGE = plain
|
||||
);
|
||||
|
||||
COMMENT ON TYPE geodetic IS 'Geodetic coordinates on WGS-84 ellipsoid (lat/lon in degrees, altitude in km)';
|
||||
|
||||
CREATE FUNCTION geodetic_lat(geodetic) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION geodetic_lon(geodetic) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION geodetic_alt(geodetic) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Topocentric type: observer-relative az/el/range
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION topocentric_in(cstring) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION topocentric_out(topocentric) RETURNS cstring
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION topocentric_recv(internal) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION topocentric_send(topocentric) RETURNS bytea
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE TYPE topocentric (
|
||||
INPUT = topocentric_in,
|
||||
OUTPUT = topocentric_out,
|
||||
RECEIVE = topocentric_recv,
|
||||
SEND = topocentric_send,
|
||||
INTERNALLENGTH = 32,
|
||||
ALIGNMENT = double,
|
||||
STORAGE = plain
|
||||
);
|
||||
|
||||
COMMENT ON TYPE topocentric IS 'Topocentric coordinates relative to observer (azimuth, elevation, range, range rate)';
|
||||
|
||||
CREATE FUNCTION topo_azimuth(topocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION topo_azimuth(topocentric) IS 'Azimuth in degrees (0=N, 90=E, 180=S, 270=W)';
|
||||
|
||||
CREATE FUNCTION topo_elevation(topocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION topo_elevation(topocentric) IS 'Elevation in degrees (0=horizon, 90=zenith)';
|
||||
|
||||
CREATE FUNCTION topo_range(topocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION topo_range(topocentric) IS 'Slant range in km';
|
||||
|
||||
CREATE FUNCTION topo_range_rate(topocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION topo_range_rate(topocentric) IS 'Range rate in km/s (positive = receding)';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Observer type: ground station location
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION observer_in(cstring) RETURNS observer
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION observer_out(observer) RETURNS cstring
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION observer_recv(internal) RETURNS observer
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION observer_send(observer) RETURNS bytea
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE TYPE observer (
|
||||
INPUT = observer_in,
|
||||
OUTPUT = observer_out,
|
||||
RECEIVE = observer_recv,
|
||||
SEND = observer_send,
|
||||
INTERNALLENGTH = 24,
|
||||
ALIGNMENT = double,
|
||||
STORAGE = plain
|
||||
);
|
||||
|
||||
COMMENT ON TYPE observer IS 'Ground observer location (accepts: 40.0N 105.3W 1655m or decimal degrees)';
|
||||
|
||||
CREATE FUNCTION observer_lat(observer) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION observer_lat(observer) IS 'Latitude in degrees (positive = North)';
|
||||
|
||||
CREATE FUNCTION observer_lon(observer) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION observer_lon(observer) IS 'Longitude in degrees (positive = East)';
|
||||
|
||||
CREATE FUNCTION observer_alt(observer) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION observer_alt(observer) IS 'Altitude in meters above WGS-84 ellipsoid';
|
||||
|
||||
CREATE FUNCTION observer_from_geodetic(float8, float8, float8 DEFAULT 0.0) RETURNS observer
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION observer_from_geodetic(float8, float8, float8) IS
|
||||
'Construct observer from lat (deg), lon (deg), altitude (meters). Avoids text formatting round-trips.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Pass event type: satellite visibility window
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION pass_event_in(cstring) RETURNS pass_event
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION pass_event_out(pass_event) RETURNS cstring
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION pass_event_recv(internal) RETURNS pass_event
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION pass_event_send(pass_event) RETURNS bytea
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE TYPE pass_event (
|
||||
INPUT = pass_event_in,
|
||||
OUTPUT = pass_event_out,
|
||||
RECEIVE = pass_event_recv,
|
||||
SEND = pass_event_send,
|
||||
INTERNALLENGTH = 48,
|
||||
ALIGNMENT = double,
|
||||
STORAGE = plain
|
||||
);
|
||||
|
||||
COMMENT ON TYPE pass_event IS 'Satellite pass event (AOS/MAX/LOS times, max elevation, AOS/LOS azimuths)';
|
||||
|
||||
CREATE FUNCTION pass_aos_time(pass_event) RETURNS timestamptz
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_aos_time(pass_event) IS 'Acquisition of signal time';
|
||||
|
||||
CREATE FUNCTION pass_max_el_time(pass_event) RETURNS timestamptz
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_max_el_time(pass_event) IS 'Maximum elevation time';
|
||||
|
||||
CREATE FUNCTION pass_los_time(pass_event) RETURNS timestamptz
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_los_time(pass_event) IS 'Loss of signal time';
|
||||
|
||||
CREATE FUNCTION pass_max_elevation(pass_event) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_max_elevation(pass_event) IS 'Maximum elevation in degrees';
|
||||
|
||||
CREATE FUNCTION pass_aos_azimuth(pass_event) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_aos_azimuth(pass_event) IS 'AOS azimuth in degrees (0=N)';
|
||||
|
||||
CREATE FUNCTION pass_los_azimuth(pass_event) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_los_azimuth(pass_event) IS 'LOS azimuth in degrees (0=N)';
|
||||
|
||||
CREATE FUNCTION pass_duration(pass_event) RETURNS interval
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_duration(pass_event) IS 'Pass duration (LOS - AOS)';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- SGP4/SDP4 propagation functions
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION sgp4_propagate(tle, timestamptz) RETURNS eci_position
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION sgp4_propagate(tle, timestamptz) IS
|
||||
'Propagate TLE to a point in time using SGP4 (near-earth) or SDP4 (deep-space). Returns TEME ECI position/velocity.';
|
||||
|
||||
CREATE FUNCTION sgp4_propagate_safe(tle, timestamptz) RETURNS eci_position
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION sgp4_propagate_safe(tle, timestamptz) IS
|
||||
'Like sgp4_propagate but returns NULL on error instead of raising an exception. For batch queries with potentially invalid TLEs.';
|
||||
|
||||
CREATE FUNCTION sgp4_propagate_series(tle, timestamptz, timestamptz, interval)
|
||||
RETURNS TABLE(t timestamptz, x float8, y float8, z float8, vx float8, vy float8, vz float8)
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE
|
||||
ROWS 100;
|
||||
COMMENT ON FUNCTION sgp4_propagate_series(tle, timestamptz, timestamptz, interval) IS
|
||||
'Propagate TLE over a time range at regular intervals. Returns time series of TEME positions.';
|
||||
|
||||
CREATE FUNCTION tle_distance(tle, tle, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_distance(tle, tle, timestamptz) IS
|
||||
'Euclidean distance in km between two TLEs at a reference time';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Coordinate transform functions
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION eci_to_geodetic(eci_position, timestamptz) RETURNS geodetic
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION eci_to_geodetic(eci_position, timestamptz) IS
|
||||
'Convert TEME ECI position to WGS-84 geodetic coordinates at given time';
|
||||
|
||||
CREATE FUNCTION eci_to_topocentric(eci_position, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION eci_to_topocentric(eci_position, observer, timestamptz) IS
|
||||
'Convert TEME ECI position to topocentric (az/el/range) relative to observer';
|
||||
|
||||
CREATE FUNCTION subsatellite_point(tle, timestamptz) RETURNS geodetic
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION subsatellite_point(tle, timestamptz) IS
|
||||
'Subsatellite (nadir) point on WGS-84 ellipsoid at given time';
|
||||
|
||||
CREATE FUNCTION ground_track(tle, timestamptz, timestamptz, interval)
|
||||
RETURNS TABLE(t timestamptz, lat float8, lon float8, alt float8)
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE
|
||||
ROWS 100;
|
||||
COMMENT ON FUNCTION ground_track(tle, timestamptz, timestamptz, interval) IS
|
||||
'Ground track as time series of subsatellite points (lat/lon in degrees, alt in km)';
|
||||
|
||||
CREATE FUNCTION observe(tle, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION observe(tle, observer, timestamptz) IS
|
||||
'Propagate TLE and compute observer-relative look angles in one call. Returns topocentric (az/el/range/range_rate).';
|
||||
|
||||
CREATE FUNCTION observe_safe(tle, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION observe_safe(tle, observer, timestamptz) IS
|
||||
'Like observe() but returns NULL on propagation error. For batch queries with potentially invalid/decayed TLEs.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Pass prediction functions
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION next_pass(tle, observer, timestamptz) RETURNS pass_event
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION next_pass(tle, observer, timestamptz) IS
|
||||
'Find the next satellite pass over observer (searches up to 7 days ahead)';
|
||||
|
||||
CREATE FUNCTION predict_passes(tle, observer, timestamptz, timestamptz, float8 DEFAULT 0.0)
|
||||
RETURNS SETOF pass_event
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE
|
||||
ROWS 10;
|
||||
COMMENT ON FUNCTION predict_passes(tle, observer, timestamptz, timestamptz, float8) IS
|
||||
'Predict all satellite passes over observer in time window. Optional min_elevation in degrees.';
|
||||
|
||||
CREATE FUNCTION pass_visible(tle, observer, timestamptz, timestamptz) RETURNS boolean
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pass_visible(tle, observer, timestamptz, timestamptz) IS
|
||||
'True if any pass occurs over observer in the time window';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- GiST operator support functions
|
||||
-- ============================================================
|
||||
|
||||
-- Overlap operator: do orbital keys overlap in altitude AND inclination?
|
||||
CREATE FUNCTION tle_overlap(tle, tle) RETURNS boolean
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
-- Altitude distance operator (altitude-only, for KNN ordering)
|
||||
CREATE FUNCTION tle_alt_distance(tle, tle) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE OPERATOR && (
|
||||
LEFTARG = tle,
|
||||
RIGHTARG = tle,
|
||||
FUNCTION = tle_overlap,
|
||||
COMMUTATOR = &&,
|
||||
RESTRICT = areasel,
|
||||
JOIN = areajoinsel
|
||||
);
|
||||
|
||||
COMMENT ON OPERATOR && (tle, tle) IS 'Orbital key overlap (altitude band AND inclination range) — necessary condition for conjunction';
|
||||
|
||||
CREATE OPERATOR <-> (
|
||||
LEFTARG = tle,
|
||||
RIGHTARG = tle,
|
||||
FUNCTION = tle_alt_distance,
|
||||
COMMUTATOR = <->
|
||||
);
|
||||
|
||||
COMMENT ON OPERATOR <-> (tle, tle) IS 'Minimum altitude-band separation in km (0 if overlapping). Altitude-only — does not account for inclination. Use && for 2-D filtering.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- GiST operator class for 2-D orbital indexing (altitude + inclination)
|
||||
-- ============================================================
|
||||
|
||||
-- GiST internal support functions
|
||||
CREATE FUNCTION gist_tle_compress(internal) RETURNS internal
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION gist_tle_decompress(internal) RETURNS internal
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION gist_tle_consistent(internal, tle, smallint, oid, internal) RETURNS boolean
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION gist_tle_union(internal, internal) RETURNS internal
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION gist_tle_penalty(internal, internal, internal) RETURNS internal
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION gist_tle_picksplit(internal, internal) RETURNS internal
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION gist_tle_same(internal, internal, internal) RETURNS internal
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION gist_tle_distance(internal, tle, smallint, oid, internal) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE OPERATOR CLASS tle_ops
|
||||
DEFAULT FOR TYPE tle USING gist AS
|
||||
OPERATOR 3 && ,
|
||||
OPERATOR 15 <-> (tle, tle) FOR ORDER BY float_ops,
|
||||
FUNCTION 1 gist_tle_consistent(internal, tle, smallint, oid, internal),
|
||||
FUNCTION 2 gist_tle_union(internal, internal),
|
||||
FUNCTION 3 gist_tle_compress(internal),
|
||||
FUNCTION 4 gist_tle_decompress(internal),
|
||||
FUNCTION 5 gist_tle_penalty(internal, internal, internal),
|
||||
FUNCTION 6 gist_tle_picksplit(internal, internal),
|
||||
FUNCTION 7 gist_tle_same(internal, internal, internal),
|
||||
FUNCTION 8 gist_tle_distance(internal, tle, smallint, oid, internal);
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Heliocentric type: ecliptic J2000 position in AU
|
||||
-- ============================================================
|
||||
|
||||
CREATE TYPE heliocentric;
|
||||
|
||||
CREATE FUNCTION heliocentric_in(cstring) RETURNS heliocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION heliocentric_out(heliocentric) RETURNS cstring
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION heliocentric_recv(internal) RETURNS heliocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE FUNCTION heliocentric_send(heliocentric) RETURNS bytea
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
|
||||
CREATE TYPE heliocentric (
|
||||
INPUT = heliocentric_in,
|
||||
OUTPUT = heliocentric_out,
|
||||
RECEIVE = heliocentric_recv,
|
||||
SEND = heliocentric_send,
|
||||
INTERNALLENGTH = 24,
|
||||
ALIGNMENT = double,
|
||||
STORAGE = plain
|
||||
);
|
||||
|
||||
COMMENT ON TYPE heliocentric IS 'Heliocentric position in ecliptic J2000 frame (x, y, z in AU)';
|
||||
|
||||
CREATE FUNCTION helio_x(heliocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION helio_x(heliocentric) IS 'X component in AU (ecliptic J2000, vernal equinox direction)';
|
||||
|
||||
CREATE FUNCTION helio_y(heliocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION helio_y(heliocentric) IS 'Y component in AU (ecliptic J2000)';
|
||||
|
||||
CREATE FUNCTION helio_z(heliocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION helio_z(heliocentric) IS 'Z component in AU (ecliptic J2000, north ecliptic pole)';
|
||||
|
||||
CREATE FUNCTION helio_distance(heliocentric) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION helio_distance(heliocentric) IS 'Heliocentric distance in AU';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Star observation functions
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION star_observe(float8, float8, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION star_observe(float8, float8, observer, timestamptz) IS
|
||||
'Observe a star from (ra_hours J2000, dec_degrees J2000, observer, time). Returns topocentric az/el. Range is 0 (infinite distance).';
|
||||
|
||||
CREATE FUNCTION star_observe_safe(float8, float8, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION star_observe_safe(float8, float8, observer, timestamptz) IS
|
||||
'Like star_observe but returns NULL on invalid inputs. For batch queries over star catalogs.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Keplerian propagation functions
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION kepler_propagate(
|
||||
q_au float8, eccentricity float8,
|
||||
inclination_deg float8, arg_perihelion_deg float8,
|
||||
long_asc_node_deg float8, perihelion_jd float8,
|
||||
t timestamptz
|
||||
) RETURNS heliocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION kepler_propagate(float8, float8, float8, float8, float8, float8, timestamptz) IS
|
||||
'Two-body Keplerian propagation from classical orbital elements. Returns heliocentric ecliptic J2000 position in AU. Handles elliptic, parabolic, and hyperbolic orbits.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Comet observation
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION comet_observe(
|
||||
q_au float8, eccentricity float8,
|
||||
inclination_deg float8, arg_perihelion_deg float8,
|
||||
long_asc_node_deg float8, perihelion_jd float8,
|
||||
earth_x_au float8, earth_y_au float8, earth_z_au float8,
|
||||
obs observer, t timestamptz
|
||||
) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION comet_observe(float8, float8, float8, float8, float8, float8, float8, float8, float8, observer, timestamptz) IS
|
||||
'Observe a comet/asteroid from orbital elements. Requires Earth heliocentric ecliptic J2000 position (AU). Returns topocentric az/el with geocentric range in km.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- VSOP87 planets, ELP82B Moon, Sun observation
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION planet_heliocentric(int4, timestamptz) RETURNS heliocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION planet_heliocentric(int4, timestamptz) IS
|
||||
'VSOP87 heliocentric ecliptic J2000 position (AU). Body IDs: 0=Sun, 1=Mercury, ..., 8=Neptune.';
|
||||
|
||||
CREATE FUNCTION planet_observe(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION planet_observe(int4, observer, timestamptz) IS
|
||||
'Observe a planet from (body_id 1-8, observer, time). Returns topocentric az/el with geocentric range in km.';
|
||||
|
||||
CREATE FUNCTION sun_observe(observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION sun_observe(observer, timestamptz) IS
|
||||
'Observe the Sun from (observer, time). Returns topocentric az/el with geocentric range in km.';
|
||||
|
||||
CREATE FUNCTION moon_observe(observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION moon_observe(observer, timestamptz) IS
|
||||
'Observe the Moon via ELP2000-82B from (observer, time). Returns topocentric az/el with geocentric range in km.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Planetary moon observation
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION galilean_observe(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION galilean_observe(int4, observer, timestamptz) IS
|
||||
'Observe a Galilean moon of Jupiter via L1.2 theory. Body IDs: 0=Io, 1=Europa, 2=Ganymede, 3=Callisto.';
|
||||
|
||||
CREATE FUNCTION saturn_moon_observe(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION saturn_moon_observe(int4, observer, timestamptz) IS
|
||||
'Observe a Saturn moon via TASS 1.7. Body IDs: 0=Mimas, 1=Enceladus, 2=Tethys, 3=Dione, 4=Rhea, 5=Titan, 6=Iapetus, 7=Hyperion.';
|
||||
|
||||
CREATE FUNCTION uranus_moon_observe(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION uranus_moon_observe(int4, observer, timestamptz) IS
|
||||
'Observe a Uranus moon via GUST86. Body IDs: 0=Miranda, 1=Ariel, 2=Umbriel, 3=Titania, 4=Oberon.';
|
||||
|
||||
CREATE FUNCTION mars_moon_observe(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION mars_moon_observe(int4, observer, timestamptz) IS
|
||||
'Observe a Mars moon via MarsSat. Body IDs: 0=Phobos, 1=Deimos.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Jupiter decametric radio burst prediction
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION io_phase_angle(timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION io_phase_angle(timestamptz) IS
|
||||
'Io orbital phase angle in degrees [0,360). 0=superior conjunction (behind Jupiter). Standard Radio JOVE convention.';
|
||||
|
||||
CREATE FUNCTION jupiter_cml(observer, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION jupiter_cml(observer, timestamptz) IS
|
||||
'Jupiter Central Meridian Longitude, System III (1965.0), in degrees [0,360). Light-time corrected.';
|
||||
|
||||
CREATE FUNCTION jupiter_burst_probability(float8, float8) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION jupiter_burst_probability(float8, float8) IS
|
||||
'Estimated Jupiter decametric burst probability (0-1) from (io_phase_deg, cml_deg). Based on Carr et al. (1983) source regions A, B, C, D.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Interplanetary transfer orbits (Lambert solver)
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION lambert_transfer(
|
||||
dep_body_id int4, arr_body_id int4,
|
||||
dep_time timestamptz, arr_time timestamptz,
|
||||
OUT c3_departure float8, OUT c3_arrival float8,
|
||||
OUT v_inf_departure float8, OUT v_inf_arrival float8,
|
||||
OUT tof_days float8, OUT transfer_sma float8
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lambert_transfer(int4, int4, timestamptz, timestamptz) IS
|
||||
'Solve Lambert transfer between two planets. Returns C3 (km^2/s^2), v_infinity (km/s), TOF (days), SMA (AU). Body IDs 1-8.';
|
||||
|
||||
CREATE FUNCTION lambert_c3(int4, int4, timestamptz, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lambert_c3(int4, int4, timestamptz, timestamptz) IS
|
||||
'Departure C3 (km^2/s^2) for a Lambert transfer. Returns NULL if solver fails. For pork chop plots.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- DE ephemeris functions (optional high-precision)
|
||||
-- ============================================================
|
||||
|
||||
CREATE FUNCTION planet_heliocentric_de(int4, timestamptz) RETURNS heliocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION planet_heliocentric_de(int4, timestamptz) IS
|
||||
'Heliocentric ecliptic J2000 position via JPL DE (sub-arcsecond). Falls back to VSOP87 if DE unavailable. Body IDs: 0=Sun, 1-8=Mercury-Neptune.';
|
||||
|
||||
CREATE FUNCTION planet_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION planet_observe_de(int4, observer, timestamptz) IS
|
||||
'Observe planet via JPL DE. Falls back to VSOP87. Body IDs: 1-8 (Mercury-Neptune).';
|
||||
|
||||
CREATE FUNCTION sun_observe_de(observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION sun_observe_de(observer, timestamptz) IS
|
||||
'Observe Sun via JPL DE. Falls back to VSOP87.';
|
||||
|
||||
CREATE FUNCTION moon_observe_de(observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION moon_observe_de(observer, timestamptz) IS
|
||||
'Observe Moon via JPL DE. Falls back to ELP2000-82B.';
|
||||
|
||||
CREATE FUNCTION lambert_transfer_de(
|
||||
dep_body_id int4, arr_body_id int4,
|
||||
dep_time timestamptz, arr_time timestamptz,
|
||||
OUT c3_departure float8, OUT c3_arrival float8,
|
||||
OUT v_inf_departure float8, OUT v_inf_arrival float8,
|
||||
OUT tof_days float8, OUT transfer_sma float8
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lambert_transfer_de(int4, int4, timestamptz, timestamptz) IS
|
||||
'Lambert transfer via JPL DE positions. Falls back to VSOP87. Body IDs 1-8.';
|
||||
|
||||
CREATE FUNCTION lambert_c3_de(int4, int4, timestamptz, timestamptz) RETURNS float8
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION lambert_c3_de(int4, int4, timestamptz, timestamptz) IS
|
||||
'Departure C3 via JPL DE. Falls back to VSOP87. For pork chop plots.';
|
||||
|
||||
CREATE FUNCTION galilean_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION galilean_observe_de(int4, observer, timestamptz) IS
|
||||
'Observe Galilean moon with JPL DE parent position. L1.2 moon offsets. Body IDs: 0-3 (Io-Callisto).';
|
||||
|
||||
CREATE FUNCTION saturn_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION saturn_moon_observe_de(int4, observer, timestamptz) IS
|
||||
'Observe Saturn moon with JPL DE parent position. TASS 1.7 moon offsets. Body IDs: 0-7 (Mimas-Hyperion).';
|
||||
|
||||
CREATE FUNCTION uranus_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION uranus_moon_observe_de(int4, observer, timestamptz) IS
|
||||
'Observe Uranus moon with JPL DE parent position. GUST86 moon offsets. Body IDs: 0-4 (Miranda-Oberon).';
|
||||
|
||||
CREATE FUNCTION mars_moon_observe_de(int4, observer, timestamptz) RETURNS topocentric
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION mars_moon_observe_de(int4, observer, timestamptz) IS
|
||||
'Observe Mars moon with JPL DE parent position. MarsSat moon offsets. Body IDs: 0-1 (Phobos-Deimos).';
|
||||
|
||||
|
||||
-- Diagnostic function
|
||||
|
||||
CREATE FUNCTION pg_orrery_ephemeris_info(
|
||||
OUT provider text, OUT file_path text,
|
||||
OUT start_jd float8, OUT end_jd float8,
|
||||
OUT version int4, OUT au_km float8
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION pg_orrery_ephemeris_info() IS
|
||||
'Returns current ephemeris provider status: VSOP87 or JPL_DE with file path, JD range, version, and AU value.';
|
||||
|
||||
|
||||
-- ============================================================
|
||||
-- Orbit determination (TLE fitting from observations)
|
||||
-- ============================================================
|
||||
|
||||
-- Fit TLE from ECI position/velocity ephemeris
|
||||
-- weights: per-observation weighting (NULL = uniform)
|
||||
|
||||
CREATE FUNCTION tle_from_eci(
|
||||
positions eci_position[], times timestamptz[],
|
||||
seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false,
|
||||
max_iter int4 DEFAULT 15,
|
||||
weights float8[] DEFAULT NULL,
|
||||
OUT fitted_tle tle, OUT iterations int4,
|
||||
OUT rms_final float8, OUT rms_initial float8, OUT status text,
|
||||
OUT condition_number float8, OUT covariance float8[], OUT nstate int4
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_from_eci(eci_position[], timestamptz[], tle, boolean, int4, float8[]) IS
|
||||
'Fit a TLE from ECI position/velocity observations via differential correction. Optional per-observation weights for heterogeneous sensor fusion. Returns fitted TLE, RMS residuals, convergence status, condition number, and formal covariance matrix.';
|
||||
|
||||
-- Fit TLE from topocentric observations (az/el/range) — single observer
|
||||
-- fit_range_rate: include range_rate as 4th residual component
|
||||
-- weights: per-observation weighting (NULL = uniform)
|
||||
|
||||
CREATE FUNCTION tle_from_topocentric(
|
||||
observations topocentric[], times timestamptz[],
|
||||
obs observer,
|
||||
seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false,
|
||||
max_iter int4 DEFAULT 15,
|
||||
fit_range_rate boolean DEFAULT false,
|
||||
weights float8[] DEFAULT NULL,
|
||||
OUT fitted_tle tle, OUT iterations int4,
|
||||
OUT rms_final float8, OUT rms_initial float8, OUT status text,
|
||||
OUT condition_number float8, OUT covariance float8[], OUT nstate int4
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_from_topocentric(topocentric[], timestamptz[], observer, tle, boolean, int4, boolean, float8[]) IS
|
||||
'Fit a TLE from topocentric (az/el/range) observations via differential correction. Optional range_rate fitting and per-observation weights. Returns fitted TLE, RMS residuals, convergence status, condition number, and formal covariance matrix.';
|
||||
|
||||
-- Fit TLE from topocentric observations — multiple observers
|
||||
|
||||
CREATE FUNCTION tle_from_topocentric(
|
||||
observations topocentric[], times timestamptz[],
|
||||
observers observer[], observer_ids int4[],
|
||||
seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false,
|
||||
max_iter int4 DEFAULT 15,
|
||||
fit_range_rate boolean DEFAULT false,
|
||||
weights float8[] DEFAULT NULL,
|
||||
OUT fitted_tle tle, OUT iterations int4,
|
||||
OUT rms_final float8, OUT rms_initial float8, OUT status text,
|
||||
OUT condition_number float8, OUT covariance float8[], OUT nstate int4
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME', 'tle_from_topocentric_multi'
|
||||
LANGUAGE C STABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_from_topocentric(topocentric[], timestamptz[], observer[], int4[], tle, boolean, int4, boolean, float8[]) IS
|
||||
'Fit a TLE from topocentric observations collected by multiple ground stations. observer_ids[i] indexes into observers[]. Optional range_rate fitting and per-observation weights.';
|
||||
|
||||
-- Per-observation residuals diagnostic
|
||||
|
||||
CREATE FUNCTION tle_fit_residuals(
|
||||
fitted tle,
|
||||
positions eci_position[],
|
||||
times timestamptz[]
|
||||
) RETURNS TABLE (
|
||||
t timestamptz,
|
||||
dx_km float8,
|
||||
dy_km float8,
|
||||
dz_km float8,
|
||||
pos_err_km float8
|
||||
)
|
||||
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_fit_residuals(tle, eci_position[], timestamptz[]) IS
|
||||
'Compute per-observation position residuals (km) between a TLE and ECI observations. Useful for fit quality assessment.';
|
||||
|
||||
-- Fit TLE from RA/Dec observations — single observer
|
||||
-- Uses Gauss method for initial orbit determination when no seed is provided.
|
||||
-- RA in hours [0,24), Dec in degrees [-90,90] (matches star_observe convention).
|
||||
-- RMS output is in radians for angles-only mode.
|
||||
|
||||
CREATE FUNCTION tle_from_angles(
|
||||
ra_hours float8[], dec_degrees float8[],
|
||||
times timestamptz[],
|
||||
obs observer,
|
||||
seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false,
|
||||
max_iter int4 DEFAULT 15,
|
||||
weights float8[] DEFAULT NULL,
|
||||
OUT fitted_tle tle, OUT iterations int4,
|
||||
OUT rms_final float8, OUT rms_initial float8, OUT status text,
|
||||
OUT condition_number float8, OUT covariance float8[], OUT nstate int4
|
||||
) RETURNS RECORD AS 'MODULE_PATHNAME' LANGUAGE C STABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_from_angles(float8[], float8[], timestamptz[], observer, tle, boolean, int4, float8[]) IS
|
||||
'Fit a TLE from angles-only (RA/Dec) observations via Gauss IOD + differential correction. RA in hours [0,24), Dec in degrees [-90,90]. RMS output in radians. Uses Gauss method for seed-free initial guess.';
|
||||
|
||||
-- Fit TLE from RA/Dec observations — multiple observers
|
||||
|
||||
CREATE FUNCTION tle_from_angles(
|
||||
ra_hours float8[], dec_degrees float8[],
|
||||
times timestamptz[],
|
||||
observers observer[], observer_ids int4[],
|
||||
seed tle DEFAULT NULL, fit_bstar boolean DEFAULT false,
|
||||
max_iter int4 DEFAULT 15,
|
||||
weights float8[] DEFAULT NULL,
|
||||
OUT fitted_tle tle, OUT iterations int4,
|
||||
OUT rms_final float8, OUT rms_initial float8, OUT status text,
|
||||
OUT condition_number float8, OUT covariance float8[], OUT nstate int4
|
||||
) RETURNS RECORD
|
||||
AS 'MODULE_PATHNAME', 'tle_from_angles_multi'
|
||||
LANGUAGE C STABLE PARALLEL SAFE;
|
||||
COMMENT ON FUNCTION tle_from_angles(float8[], float8[], timestamptz[], observer[], int4[], tle, boolean, int4, float8[]) IS
|
||||
'Fit a TLE from angles-only (RA/Dec) observations from multiple ground stations via Gauss IOD + differential correction.';
|
||||
494
src/od_funcs.c
494
src/od_funcs.c
@ -31,8 +31,6 @@ PG_FUNCTION_INFO_V1(tle_from_eci);
|
||||
PG_FUNCTION_INFO_V1(tle_from_topocentric);
|
||||
PG_FUNCTION_INFO_V1(tle_from_topocentric_multi);
|
||||
PG_FUNCTION_INFO_V1(tle_fit_residuals);
|
||||
PG_FUNCTION_INFO_V1(tle_from_angles);
|
||||
PG_FUNCTION_INFO_V1(tle_from_angles_multi);
|
||||
|
||||
/* ================================================================
|
||||
* Helper: pg_tle ↔ tle_t conversion (local copy, avoids symbol coupling)
|
||||
@ -107,7 +105,6 @@ tle_from_eci(PG_FUNCTION_ARGS)
|
||||
pg_tle *seed_pg = has_seed ? (pg_tle *) PG_GETARG_POINTER(2) : NULL;
|
||||
bool fit_bstar = PG_ARGISNULL(3) ? false : PG_GETARG_BOOL(3);
|
||||
int32 max_iter = PG_ARGISNULL(4) ? 15 : PG_GETARG_INT32(4);
|
||||
bool has_weights = !PG_ARGISNULL(5);
|
||||
|
||||
int n_pos, n_times;
|
||||
Datum *pos_datums, *time_datums;
|
||||
@ -168,29 +165,6 @@ tle_from_eci(PG_FUNCTION_ARGS)
|
||||
config.observers = NULL;
|
||||
config.n_observers = 0;
|
||||
|
||||
/* Extract per-observation weights if provided */
|
||||
if (has_weights)
|
||||
{
|
||||
ArrayType *wt_arr = PG_GETARG_ARRAYTYPE_P(5);
|
||||
Datum *wt_datums;
|
||||
bool *wt_nulls;
|
||||
int n_wt;
|
||||
|
||||
deconstruct_array(wt_arr, FLOAT8OID, sizeof(float8), FLOAT8PASSBYVAL,
|
||||
TYPALIGN_DOUBLE, &wt_datums, &wt_nulls, &n_wt);
|
||||
|
||||
if (n_wt != n_pos)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_ARRAY_SUBSCRIPT_ERROR),
|
||||
errmsg("weights array length must match observations: "
|
||||
"%d vs %d", n_wt, n_pos)));
|
||||
|
||||
config.weights = (double *) palloc(sizeof(double) * n_wt);
|
||||
for (i = 0; i < n_wt; i++)
|
||||
config.weights[i] = DatumGetFloat8(wt_datums[i]);
|
||||
config.n_weights = n_wt;
|
||||
}
|
||||
|
||||
/* Convert seed TLE if provided */
|
||||
if (has_seed)
|
||||
pg_tle_to_sat_code_od(seed_pg, &seed_sat);
|
||||
@ -201,8 +175,6 @@ tle_from_eci(PG_FUNCTION_ARGS)
|
||||
rc = od_fit_tle(obs, n_pos, has_seed ? &seed_sat : NULL, &config, &result);
|
||||
|
||||
pfree(obs);
|
||||
if (config.weights)
|
||||
pfree(config.weights);
|
||||
|
||||
if (rc != 0)
|
||||
ereport(ERROR,
|
||||
@ -272,8 +244,6 @@ tle_from_topocentric(PG_FUNCTION_ARGS)
|
||||
pg_tle *seed_pg = has_seed ? (pg_tle *) PG_GETARG_POINTER(3) : NULL;
|
||||
bool fit_bstar = PG_ARGISNULL(4) ? false : PG_GETARG_BOOL(4);
|
||||
int32 max_iter = PG_ARGISNULL(5) ? 15 : PG_GETARG_INT32(5);
|
||||
bool fit_range_rate = PG_ARGISNULL(6) ? false : PG_GETARG_BOOL(6);
|
||||
bool has_weights = !PG_ARGISNULL(7);
|
||||
|
||||
int n_topo, n_times;
|
||||
od_observation_t *obs;
|
||||
@ -331,42 +301,17 @@ tle_from_topocentric(PG_FUNCTION_ARGS)
|
||||
obs[i].data[0] = topo->azimuth;
|
||||
obs[i].data[1] = topo->elevation;
|
||||
obs[i].data[2] = topo->range_km;
|
||||
obs[i].data[3] = topo->range_rate;
|
||||
obs[i].observer_idx = 0; /* single observer */
|
||||
}
|
||||
}
|
||||
|
||||
/* Configure solver */
|
||||
memset(&config, 0, sizeof(config));
|
||||
config.obs_type = OD_OBS_TOPO;
|
||||
config.fit_bstar = fit_bstar ? 1 : 0;
|
||||
config.fit_range_rate = fit_range_rate ? 1 : 0;
|
||||
config.max_iter = max_iter;
|
||||
config.observers = &observer;
|
||||
config.n_observers = 1;
|
||||
|
||||
/* Extract per-observation weights if provided */
|
||||
if (has_weights)
|
||||
{
|
||||
ArrayType *wt_arr = PG_GETARG_ARRAYTYPE_P(7);
|
||||
Datum *wt_datums;
|
||||
bool *wt_nulls;
|
||||
int n_wt;
|
||||
|
||||
deconstruct_array(wt_arr, FLOAT8OID, sizeof(float8), FLOAT8PASSBYVAL,
|
||||
TYPALIGN_DOUBLE, &wt_datums, &wt_nulls, &n_wt);
|
||||
|
||||
if (n_wt != n_topo)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_ARRAY_SUBSCRIPT_ERROR),
|
||||
errmsg("weights array length must match observations: "
|
||||
"%d vs %d", n_wt, n_topo)));
|
||||
|
||||
config.weights = (double *) palloc(sizeof(double) * n_wt);
|
||||
for (i = 0; i < n_wt; i++)
|
||||
config.weights[i] = DatumGetFloat8(wt_datums[i]);
|
||||
config.n_weights = n_wt;
|
||||
}
|
||||
config.obs_type = OD_OBS_TOPO;
|
||||
config.fit_bstar = fit_bstar ? 1 : 0;
|
||||
config.max_iter = max_iter;
|
||||
config.observers = &observer;
|
||||
config.n_observers = 1;
|
||||
|
||||
if (has_seed)
|
||||
pg_tle_to_sat_code_od(seed_pg, &seed_sat);
|
||||
@ -376,8 +321,6 @@ tle_from_topocentric(PG_FUNCTION_ARGS)
|
||||
rc = od_fit_tle(obs, n_topo, has_seed ? &seed_sat : NULL, &config, &result);
|
||||
|
||||
pfree(obs);
|
||||
if (config.weights)
|
||||
pfree(config.weights);
|
||||
|
||||
if (rc != 0)
|
||||
ereport(ERROR,
|
||||
@ -451,8 +394,6 @@ tle_from_topocentric_multi(PG_FUNCTION_ARGS)
|
||||
pg_tle *seed_pg = has_seed ? (pg_tle *) PG_GETARG_POINTER(4) : NULL;
|
||||
bool fit_bstar = PG_ARGISNULL(5) ? false : PG_GETARG_BOOL(5);
|
||||
int32 max_iter = PG_ARGISNULL(6) ? 15 : PG_GETARG_INT32(6);
|
||||
bool fit_range_rate = PG_ARGISNULL(7) ? false : PG_GETARG_BOOL(7);
|
||||
bool has_weights = !PG_ARGISNULL(8);
|
||||
|
||||
int n_topo, n_times, n_obs_sites, n_ids;
|
||||
od_observation_t *obs;
|
||||
@ -533,41 +474,16 @@ tle_from_topocentric_multi(PG_FUNCTION_ARGS)
|
||||
obs[i].data[0] = topo->azimuth;
|
||||
obs[i].data[1] = topo->elevation;
|
||||
obs[i].data[2] = topo->range_km;
|
||||
obs[i].data[3] = topo->range_rate;
|
||||
obs[i].observer_idx = oid;
|
||||
}
|
||||
|
||||
/* Configure solver */
|
||||
memset(&config, 0, sizeof(config));
|
||||
config.obs_type = OD_OBS_TOPO;
|
||||
config.fit_bstar = fit_bstar ? 1 : 0;
|
||||
config.fit_range_rate = fit_range_rate ? 1 : 0;
|
||||
config.max_iter = max_iter;
|
||||
config.observers = observers;
|
||||
config.n_observers = n_obs_sites;
|
||||
|
||||
/* Extract per-observation weights if provided */
|
||||
if (has_weights)
|
||||
{
|
||||
ArrayType *wt_arr = PG_GETARG_ARRAYTYPE_P(8);
|
||||
Datum *wt_datums;
|
||||
bool *wt_nulls;
|
||||
int n_wt;
|
||||
|
||||
deconstruct_array(wt_arr, FLOAT8OID, sizeof(float8), FLOAT8PASSBYVAL,
|
||||
TYPALIGN_DOUBLE, &wt_datums, &wt_nulls, &n_wt);
|
||||
|
||||
if (n_wt != n_topo)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_ARRAY_SUBSCRIPT_ERROR),
|
||||
errmsg("weights array length must match observations: "
|
||||
"%d vs %d", n_wt, n_topo)));
|
||||
|
||||
config.weights = (double *) palloc(sizeof(double) * n_wt);
|
||||
for (i = 0; i < n_wt; i++)
|
||||
config.weights[i] = DatumGetFloat8(wt_datums[i]);
|
||||
config.n_weights = n_wt;
|
||||
}
|
||||
config.obs_type = OD_OBS_TOPO;
|
||||
config.fit_bstar = fit_bstar ? 1 : 0;
|
||||
config.max_iter = max_iter;
|
||||
config.observers = observers;
|
||||
config.n_observers = n_obs_sites;
|
||||
|
||||
if (has_seed)
|
||||
pg_tle_to_sat_code_od(seed_pg, &seed_sat);
|
||||
@ -578,396 +494,6 @@ tle_from_topocentric_multi(PG_FUNCTION_ARGS)
|
||||
|
||||
pfree(obs);
|
||||
pfree(observers);
|
||||
if (config.weights)
|
||||
pfree(config.weights);
|
||||
|
||||
if (rc != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_DATA_EXCEPTION),
|
||||
errmsg("TLE fitting failed: %s", result.status)));
|
||||
|
||||
/* Build result tuple */
|
||||
if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
|
||||
errmsg("function returning record called in context "
|
||||
"that cannot accept type record")));
|
||||
tupdesc = BlessTupleDesc(tupdesc);
|
||||
|
||||
memset(nulls, 0, sizeof(nulls));
|
||||
|
||||
{
|
||||
pg_tle *fitted = (pg_tle *) palloc(sizeof(pg_tle));
|
||||
sat_code_to_pg_tle(&result.fitted_tle, fitted);
|
||||
values[0] = PointerGetDatum(fitted);
|
||||
}
|
||||
values[1] = Int32GetDatum(result.iterations);
|
||||
values[2] = Float8GetDatum(result.rms_final);
|
||||
values[3] = Float8GetDatum(result.rms_initial);
|
||||
values[4] = CStringGetTextDatum(result.status);
|
||||
values[5] = Float8GetDatum(result.condition_number);
|
||||
|
||||
if (result.cov_size > 0)
|
||||
{
|
||||
Datum *cov_datums = (Datum *) palloc(sizeof(Datum) * result.cov_size);
|
||||
int ci;
|
||||
for (ci = 0; ci < result.cov_size; ci++)
|
||||
cov_datums[ci] = Float8GetDatum(result.covariance[ci]);
|
||||
values[6] = PointerGetDatum(
|
||||
construct_array(cov_datums, result.cov_size,
|
||||
FLOAT8OID, sizeof(float8), FLOAT8PASSBYVAL,
|
||||
TYPALIGN_DOUBLE));
|
||||
}
|
||||
else
|
||||
{
|
||||
nulls[6] = true;
|
||||
}
|
||||
|
||||
values[7] = Int32GetDatum(result.cov_size > 0
|
||||
? (result.cov_size == 28 ? 7 : 6)
|
||||
: 0);
|
||||
|
||||
tuple = heap_form_tuple(tupdesc, values, nulls);
|
||||
PG_RETURN_DATUM(HeapTupleGetDatum(tuple));
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* tle_from_angles(ra_hours[], dec_degrees[], timestamptz[], observer,
|
||||
* tle, boolean, int4, float8[])
|
||||
* -> RECORD (same 8-column output as tle_from_eci)
|
||||
*
|
||||
* Fit TLE from angles-only (RA/Dec) observations.
|
||||
* RA in hours [0,24), Dec in degrees [-90,90].
|
||||
* Uses Gauss IOD for seed-free initial guess.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
tle_from_angles(PG_FUNCTION_ARGS)
|
||||
{
|
||||
ArrayType *ra_arr = PG_GETARG_ARRAYTYPE_P(0);
|
||||
ArrayType *dec_arr = PG_GETARG_ARRAYTYPE_P(1);
|
||||
ArrayType *time_arr = PG_GETARG_ARRAYTYPE_P(2);
|
||||
pg_observer *obs_pg = (pg_observer *) PG_GETARG_POINTER(3);
|
||||
bool has_seed = !PG_ARGISNULL(4);
|
||||
pg_tle *seed_pg = has_seed ? (pg_tle *) PG_GETARG_POINTER(4) : NULL;
|
||||
bool fit_bstar = PG_ARGISNULL(5) ? false : PG_GETARG_BOOL(5);
|
||||
int32 max_iter = PG_ARGISNULL(6) ? 15 : PG_GETARG_INT32(6);
|
||||
bool has_weights = !PG_ARGISNULL(7);
|
||||
|
||||
int n_ra, n_dec, n_times;
|
||||
Datum *ra_datums, *dec_datums, *time_datums;
|
||||
bool *ra_nulls, *dec_nulls, *time_nulls;
|
||||
od_observation_t *obs;
|
||||
od_config_t config;
|
||||
od_observer_t observer;
|
||||
od_result_t result;
|
||||
tle_t seed_sat;
|
||||
int i, rc;
|
||||
|
||||
TupleDesc tupdesc;
|
||||
Datum values[8];
|
||||
bool nulls[8];
|
||||
HeapTuple tuple;
|
||||
|
||||
/* Build observer */
|
||||
observer.lat = obs_pg->lat;
|
||||
observer.lon = obs_pg->lon;
|
||||
observer.alt_m = obs_pg->alt_m;
|
||||
od_observer_to_ecef(observer.lat, observer.lon, observer.alt_m,
|
||||
observer.ecef);
|
||||
|
||||
/* Deconstruct arrays */
|
||||
deconstruct_array(ra_arr, FLOAT8OID, sizeof(float8), FLOAT8PASSBYVAL,
|
||||
TYPALIGN_DOUBLE, &ra_datums, &ra_nulls, &n_ra);
|
||||
deconstruct_array(dec_arr, FLOAT8OID, sizeof(float8), FLOAT8PASSBYVAL,
|
||||
TYPALIGN_DOUBLE, &dec_datums, &dec_nulls, &n_dec);
|
||||
deconstruct_array(time_arr, TIMESTAMPTZOID, sizeof(int64), FLOAT8PASSBYVAL,
|
||||
TYPALIGN_DOUBLE, &time_datums, &time_nulls, &n_times);
|
||||
|
||||
if (n_ra != n_dec || n_ra != n_times)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_ARRAY_SUBSCRIPT_ERROR),
|
||||
errmsg("ra, dec, and times arrays must have same length: "
|
||||
"%d, %d, %d", n_ra, n_dec, n_times)));
|
||||
|
||||
if (n_ra < 6)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
|
||||
errmsg("at least 6 observations required, got %d", n_ra)));
|
||||
|
||||
/* Build observation array — convert RA hours → radians, Dec deg → radians */
|
||||
obs = (od_observation_t *) palloc(sizeof(od_observation_t) * n_ra);
|
||||
|
||||
for (i = 0; i < n_ra; i++)
|
||||
{
|
||||
double ra_hr = DatumGetFloat8(ra_datums[i]);
|
||||
double dec_dg = DatumGetFloat8(dec_datums[i]);
|
||||
int64 ts = DatumGetInt64(time_datums[i]);
|
||||
|
||||
if (ra_hr < 0.0 || ra_hr >= 24.0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
|
||||
errmsg("RA must be in [0, 24) hours, got %g", ra_hr)));
|
||||
if (dec_dg < -90.0 || dec_dg > 90.0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
|
||||
errmsg("Dec must be in [-90, 90] degrees, got %g", dec_dg)));
|
||||
|
||||
obs[i].jd = PG_EPOCH_JD + ((double)ts / (double)USECS_PER_DAY);
|
||||
obs[i].data[0] = ra_hr * (M_PI / 12.0); /* hours → radians */
|
||||
obs[i].data[1] = dec_dg * (M_PI / 180.0); /* degrees → radians */
|
||||
obs[i].observer_idx = 0;
|
||||
}
|
||||
|
||||
/* Configure solver */
|
||||
memset(&config, 0, sizeof(config));
|
||||
config.obs_type = OD_OBS_RADEC;
|
||||
config.fit_bstar = fit_bstar ? 1 : 0;
|
||||
config.max_iter = max_iter;
|
||||
config.observers = &observer;
|
||||
config.n_observers = 1;
|
||||
|
||||
/* Extract per-observation weights if provided */
|
||||
if (has_weights)
|
||||
{
|
||||
ArrayType *wt_arr = PG_GETARG_ARRAYTYPE_P(7);
|
||||
Datum *wt_datums;
|
||||
bool *wt_nulls;
|
||||
int n_wt;
|
||||
|
||||
deconstruct_array(wt_arr, FLOAT8OID, sizeof(float8), FLOAT8PASSBYVAL,
|
||||
TYPALIGN_DOUBLE, &wt_datums, &wt_nulls, &n_wt);
|
||||
|
||||
if (n_wt != n_ra)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_ARRAY_SUBSCRIPT_ERROR),
|
||||
errmsg("weights array length must match observations: "
|
||||
"%d vs %d", n_wt, n_ra)));
|
||||
|
||||
config.weights = (double *) palloc(sizeof(double) * n_wt);
|
||||
for (i = 0; i < n_wt; i++)
|
||||
config.weights[i] = DatumGetFloat8(wt_datums[i]);
|
||||
config.n_weights = n_wt;
|
||||
}
|
||||
|
||||
if (has_seed)
|
||||
pg_tle_to_sat_code_od(seed_pg, &seed_sat);
|
||||
|
||||
memset(&result, 0, sizeof(result));
|
||||
|
||||
rc = od_fit_tle(obs, n_ra, has_seed ? &seed_sat : NULL, &config, &result);
|
||||
|
||||
pfree(obs);
|
||||
if (config.weights)
|
||||
pfree(config.weights);
|
||||
|
||||
if (rc != 0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_DATA_EXCEPTION),
|
||||
errmsg("TLE fitting failed: %s", result.status)));
|
||||
|
||||
/* Build result tuple */
|
||||
if (get_call_result_type(fcinfo, NULL, &tupdesc) != TYPEFUNC_COMPOSITE)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_FEATURE_NOT_SUPPORTED),
|
||||
errmsg("function returning record called in context "
|
||||
"that cannot accept type record")));
|
||||
tupdesc = BlessTupleDesc(tupdesc);
|
||||
|
||||
memset(nulls, 0, sizeof(nulls));
|
||||
|
||||
{
|
||||
pg_tle *fitted = (pg_tle *) palloc(sizeof(pg_tle));
|
||||
sat_code_to_pg_tle(&result.fitted_tle, fitted);
|
||||
values[0] = PointerGetDatum(fitted);
|
||||
}
|
||||
values[1] = Int32GetDatum(result.iterations);
|
||||
values[2] = Float8GetDatum(result.rms_final);
|
||||
values[3] = Float8GetDatum(result.rms_initial);
|
||||
values[4] = CStringGetTextDatum(result.status);
|
||||
values[5] = Float8GetDatum(result.condition_number);
|
||||
|
||||
if (result.cov_size > 0)
|
||||
{
|
||||
Datum *cov_datums = (Datum *) palloc(sizeof(Datum) * result.cov_size);
|
||||
int ci;
|
||||
for (ci = 0; ci < result.cov_size; ci++)
|
||||
cov_datums[ci] = Float8GetDatum(result.covariance[ci]);
|
||||
values[6] = PointerGetDatum(
|
||||
construct_array(cov_datums, result.cov_size,
|
||||
FLOAT8OID, sizeof(float8), FLOAT8PASSBYVAL,
|
||||
TYPALIGN_DOUBLE));
|
||||
}
|
||||
else
|
||||
{
|
||||
nulls[6] = true;
|
||||
}
|
||||
|
||||
values[7] = Int32GetDatum(result.cov_size > 0
|
||||
? (result.cov_size == 28 ? 7 : 6)
|
||||
: 0);
|
||||
|
||||
tuple = heap_form_tuple(tupdesc, values, nulls);
|
||||
PG_RETURN_DATUM(HeapTupleGetDatum(tuple));
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* tle_from_angles_multi(ra_hours[], dec_degrees[], timestamptz[],
|
||||
* observer[], int4[],
|
||||
* tle, boolean, int4, float8[])
|
||||
* -> RECORD (same 8-column output)
|
||||
*
|
||||
* Multi-observer variant of tle_from_angles.
|
||||
* ================================================================
|
||||
*/
|
||||
Datum
|
||||
tle_from_angles_multi(PG_FUNCTION_ARGS)
|
||||
{
|
||||
ArrayType *ra_arr = PG_GETARG_ARRAYTYPE_P(0);
|
||||
ArrayType *dec_arr = PG_GETARG_ARRAYTYPE_P(1);
|
||||
ArrayType *time_arr = PG_GETARG_ARRAYTYPE_P(2);
|
||||
ArrayType *obs_arr = PG_GETARG_ARRAYTYPE_P(3);
|
||||
ArrayType *id_arr = PG_GETARG_ARRAYTYPE_P(4);
|
||||
bool has_seed = !PG_ARGISNULL(5);
|
||||
pg_tle *seed_pg = has_seed ? (pg_tle *) PG_GETARG_POINTER(5) : NULL;
|
||||
bool fit_bstar = PG_ARGISNULL(6) ? false : PG_GETARG_BOOL(6);
|
||||
int32 max_iter = PG_ARGISNULL(7) ? 15 : PG_GETARG_INT32(7);
|
||||
bool has_weights = !PG_ARGISNULL(8);
|
||||
|
||||
int n_ra, n_dec, n_times, n_obs_sites, n_ids;
|
||||
Datum *ra_datums, *dec_datums, *time_datums, *obs_datums, *id_datums;
|
||||
bool *ra_nulls, *dec_nulls, *time_nulls, *obs_nulls, *id_nulls;
|
||||
od_observation_t *obs;
|
||||
od_config_t config;
|
||||
od_observer_t *observers;
|
||||
od_result_t result;
|
||||
tle_t seed_sat;
|
||||
int i, rc;
|
||||
|
||||
TupleDesc tupdesc;
|
||||
Datum values[8];
|
||||
bool nulls[8];
|
||||
HeapTuple tuple;
|
||||
|
||||
/* Deconstruct all arrays */
|
||||
deconstruct_array(ra_arr, FLOAT8OID, sizeof(float8), FLOAT8PASSBYVAL,
|
||||
TYPALIGN_DOUBLE, &ra_datums, &ra_nulls, &n_ra);
|
||||
deconstruct_array(dec_arr, FLOAT8OID, sizeof(float8), FLOAT8PASSBYVAL,
|
||||
TYPALIGN_DOUBLE, &dec_datums, &dec_nulls, &n_dec);
|
||||
deconstruct_array(time_arr, TIMESTAMPTZOID, sizeof(int64), FLOAT8PASSBYVAL,
|
||||
TYPALIGN_DOUBLE, &time_datums, &time_nulls, &n_times);
|
||||
deconstruct_array(obs_arr, obs_arr->elemtype, sizeof(pg_observer),
|
||||
false, TYPALIGN_DOUBLE,
|
||||
&obs_datums, &obs_nulls, &n_obs_sites);
|
||||
deconstruct_array(id_arr, INT4OID, sizeof(int32), true,
|
||||
TYPALIGN_INT, &id_datums, &id_nulls, &n_ids);
|
||||
|
||||
if (n_ra != n_dec || n_ra != n_times)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_ARRAY_SUBSCRIPT_ERROR),
|
||||
errmsg("ra, dec, and times arrays must have same length: "
|
||||
"%d, %d, %d", n_ra, n_dec, n_times)));
|
||||
|
||||
if (n_ra != n_ids)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_ARRAY_SUBSCRIPT_ERROR),
|
||||
errmsg("observations and observer_ids arrays must have same length: "
|
||||
"%d vs %d", n_ra, n_ids)));
|
||||
|
||||
if (n_ra < 6)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
|
||||
errmsg("at least 6 observations required, got %d", n_ra)));
|
||||
|
||||
if (n_obs_sites < 1)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
|
||||
errmsg("at least 1 observer required")));
|
||||
|
||||
/* Build observer array */
|
||||
observers = (od_observer_t *) palloc(sizeof(od_observer_t) * n_obs_sites);
|
||||
for (i = 0; i < n_obs_sites; i++)
|
||||
{
|
||||
pg_observer *op = (pg_observer *) DatumGetPointer(obs_datums[i]);
|
||||
observers[i].lat = op->lat;
|
||||
observers[i].lon = op->lon;
|
||||
observers[i].alt_m = op->alt_m;
|
||||
od_observer_to_ecef(op->lat, op->lon, op->alt_m, observers[i].ecef);
|
||||
}
|
||||
|
||||
/* Build observation array */
|
||||
obs = (od_observation_t *) palloc(sizeof(od_observation_t) * n_ra);
|
||||
for (i = 0; i < n_ra; i++)
|
||||
{
|
||||
double ra_hr = DatumGetFloat8(ra_datums[i]);
|
||||
double dec_dg = DatumGetFloat8(dec_datums[i]);
|
||||
int64 ts = DatumGetInt64(time_datums[i]);
|
||||
int32 oid = DatumGetInt32(id_datums[i]);
|
||||
|
||||
if (ra_hr < 0.0 || ra_hr >= 24.0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
|
||||
errmsg("RA must be in [0, 24) hours, got %g", ra_hr)));
|
||||
if (dec_dg < -90.0 || dec_dg > 90.0)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
|
||||
errmsg("Dec must be in [-90, 90] degrees, got %g", dec_dg)));
|
||||
if (oid < 0 || oid >= n_obs_sites)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_INVALID_PARAMETER_VALUE),
|
||||
errmsg("observer_id %d out of range [0, %d)",
|
||||
oid, n_obs_sites)));
|
||||
|
||||
obs[i].jd = PG_EPOCH_JD + ((double)ts / (double)USECS_PER_DAY);
|
||||
obs[i].data[0] = ra_hr * (M_PI / 12.0);
|
||||
obs[i].data[1] = dec_dg * (M_PI / 180.0);
|
||||
obs[i].observer_idx = oid;
|
||||
}
|
||||
|
||||
/* Configure solver */
|
||||
memset(&config, 0, sizeof(config));
|
||||
config.obs_type = OD_OBS_RADEC;
|
||||
config.fit_bstar = fit_bstar ? 1 : 0;
|
||||
config.max_iter = max_iter;
|
||||
config.observers = observers;
|
||||
config.n_observers = n_obs_sites;
|
||||
|
||||
/* Extract per-observation weights if provided */
|
||||
if (has_weights)
|
||||
{
|
||||
ArrayType *wt_arr = PG_GETARG_ARRAYTYPE_P(8);
|
||||
Datum *wt_datums;
|
||||
bool *wt_nulls;
|
||||
int n_wt;
|
||||
|
||||
deconstruct_array(wt_arr, FLOAT8OID, sizeof(float8), FLOAT8PASSBYVAL,
|
||||
TYPALIGN_DOUBLE, &wt_datums, &wt_nulls, &n_wt);
|
||||
|
||||
if (n_wt != n_ra)
|
||||
ereport(ERROR,
|
||||
(errcode(ERRCODE_ARRAY_SUBSCRIPT_ERROR),
|
||||
errmsg("weights array length must match observations: "
|
||||
"%d vs %d", n_wt, n_ra)));
|
||||
|
||||
config.weights = (double *) palloc(sizeof(double) * n_wt);
|
||||
for (i = 0; i < n_wt; i++)
|
||||
config.weights[i] = DatumGetFloat8(wt_datums[i]);
|
||||
config.n_weights = n_wt;
|
||||
}
|
||||
|
||||
if (has_seed)
|
||||
pg_tle_to_sat_code_od(seed_pg, &seed_sat);
|
||||
|
||||
memset(&result, 0, sizeof(result));
|
||||
|
||||
rc = od_fit_tle(obs, n_ra, has_seed ? &seed_sat : NULL, &config, &result);
|
||||
|
||||
pfree(obs);
|
||||
pfree(observers);
|
||||
if (config.weights)
|
||||
pfree(config.weights);
|
||||
|
||||
if (rc != 0)
|
||||
ereport(ERROR,
|
||||
|
||||
184
src/od_iod.c
184
src/od_iod.c
@ -125,187 +125,3 @@ od_gibbs(const double pos1[3], const double pos2[3],
|
||||
(void)jd3;
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* Gauss method for angles-only initial orbit determination
|
||||
*
|
||||
* Given 3 RA/Dec observations and observer positions, recovers the
|
||||
* orbit at the middle observation epoch. Vallado Algorithm 52.
|
||||
*
|
||||
* Steps:
|
||||
* 1. Compute LOS unit vectors from RA/Dec
|
||||
* 2. Rotate observer ECEF → TEME at each epoch via GMST
|
||||
* 3. Build D matrix from LOS vectors and observer positions
|
||||
* 4. Iteratively solve for slant range at middle observation
|
||||
* 5. Use Gibbs/Herrick-Gibbs to get velocity at r2
|
||||
* 6. Convert (r2, v2) → Keplerian elements
|
||||
* ================================================================
|
||||
*/
|
||||
|
||||
int
|
||||
od_gauss(const double ra[3], const double dec[3],
|
||||
const double jd[3],
|
||||
const double obs_ecef[3][3],
|
||||
od_iod_result_t *result)
|
||||
{
|
||||
double L[3][3]; /* LOS unit vectors */
|
||||
double R[3][3]; /* Observer positions in TEME */
|
||||
double tau1, tau3, tau;
|
||||
double D[3][3]; /* D matrix */
|
||||
double D0;
|
||||
double A, B;
|
||||
double r2_mag, r2_mag_old;
|
||||
double rho[3]; /* slant ranges */
|
||||
double r2[3], v2[3];
|
||||
int iter, i, rc;
|
||||
double gmst;
|
||||
|
||||
result->valid = 0;
|
||||
|
||||
/* Time intervals in seconds */
|
||||
tau1 = (jd[0] - jd[1]) * 86400.0;
|
||||
tau3 = (jd[2] - jd[1]) * 86400.0;
|
||||
tau = tau3 - tau1;
|
||||
|
||||
if (fabs(tau) < 1.0)
|
||||
return -1; /* observations too close in time */
|
||||
|
||||
/* LOS unit vectors from RA/Dec */
|
||||
for (i = 0; i < 3; i++)
|
||||
od_radec_to_los(ra[i], dec[i], L[i]);
|
||||
|
||||
/* Observer ECEF → TEME */
|
||||
for (i = 0; i < 3; i++)
|
||||
{
|
||||
double vel_zero[3] = {0, 0, 0};
|
||||
double vel_dummy[3];
|
||||
gmst = od_gmst_from_jd(jd[i]);
|
||||
od_ecef_to_teme(obs_ecef[i], vel_zero, gmst, R[i], vel_dummy);
|
||||
}
|
||||
|
||||
/* Build D matrix: D[i][j] = L[i] . (L_cross products) */
|
||||
{
|
||||
double L2xL3[3], L1xL3[3], L1xL2[3];
|
||||
|
||||
vec_cross(L[1], L[2], L2xL3);
|
||||
vec_cross(L[0], L[2], L1xL3);
|
||||
vec_cross(L[0], L[1], L1xL2);
|
||||
|
||||
D0 = vec_dot(L[0], L2xL3);
|
||||
|
||||
if (fabs(D0) < 1e-15)
|
||||
return -1; /* coplanar LOS vectors */
|
||||
|
||||
/* D matrix: each row is dot products of R[i] with cross products */
|
||||
D[0][0] = vec_dot(R[0], L2xL3);
|
||||
D[0][1] = vec_dot(R[0], L1xL3);
|
||||
D[0][2] = vec_dot(R[0], L1xL2);
|
||||
|
||||
D[1][0] = vec_dot(R[1], L2xL3);
|
||||
D[1][1] = vec_dot(R[1], L1xL3);
|
||||
D[1][2] = vec_dot(R[1], L1xL2);
|
||||
|
||||
D[2][0] = vec_dot(R[2], L2xL3);
|
||||
D[2][1] = vec_dot(R[2], L1xL3);
|
||||
D[2][2] = vec_dot(R[2], L1xL2);
|
||||
}
|
||||
|
||||
/* A and B coefficients for the range equation */
|
||||
A = (-D[0][1] * tau3 / tau + D[1][1] + D[2][1] * tau1 / tau) / D0;
|
||||
B = (D[0][1] * (tau3 * tau3 - tau * tau) * tau3 / tau
|
||||
+ D[2][1] * (tau * tau - tau1 * tau1) * tau1 / tau) / (6.0 * D0);
|
||||
|
||||
/* Initial guess: r2_mag from observer position magnitude */
|
||||
r2_mag = vec_mag(R[1]) + 500.0; /* rough LEO assumption */
|
||||
|
||||
/* Iterate to find r2_mag */
|
||||
for (iter = 0; iter < 50; iter++)
|
||||
{
|
||||
double r2_3 = r2_mag * r2_mag * r2_mag;
|
||||
double f1, f3, g1, g3;
|
||||
double c1, c3;
|
||||
|
||||
/* f and g series (two-body, truncated) */
|
||||
f1 = 1.0 - 0.5 * MU_KM3_S2 * tau1 * tau1 / r2_3;
|
||||
f3 = 1.0 - 0.5 * MU_KM3_S2 * tau3 * tau3 / r2_3;
|
||||
g1 = tau1 - MU_KM3_S2 * tau1 * tau1 * tau1 / (6.0 * r2_3);
|
||||
g3 = tau3 - MU_KM3_S2 * tau3 * tau3 * tau3 / (6.0 * r2_3);
|
||||
|
||||
/* Lagrange coefficients */
|
||||
c1 = g3 / (f1 * g3 - f3 * g1);
|
||||
c3 = -g1 / (f1 * g3 - f3 * g1);
|
||||
|
||||
/* Slant ranges */
|
||||
rho[0] = (-D[0][0] + D[1][0] / c1 - D[2][0] * c3 / c1) / D0;
|
||||
rho[1] = A + B / r2_3;
|
||||
rho[2] = (-D[0][2] * c1 / c3 + D[1][2] / c3 - D[2][2]) / D0;
|
||||
|
||||
/* Position at middle observation */
|
||||
for (i = 0; i < 3; i++)
|
||||
r2[i] = R[1][i] + rho[1] * L[1][i];
|
||||
|
||||
r2_mag_old = r2_mag;
|
||||
r2_mag = vec_mag(r2);
|
||||
|
||||
if (fabs(r2_mag - r2_mag_old) < 1e-6)
|
||||
break;
|
||||
}
|
||||
|
||||
if (iter >= 50 || r2_mag < 100.0)
|
||||
return -1; /* failed to converge or nonsensical result */
|
||||
|
||||
/* Check slant ranges are positive */
|
||||
if (rho[0] < 0.0 || rho[1] < 0.0 || rho[2] < 0.0)
|
||||
return -1;
|
||||
|
||||
/* Compute all 3 position vectors for Gibbs */
|
||||
{
|
||||
double pos1[3], pos3[3];
|
||||
|
||||
for (i = 0; i < 3; i++)
|
||||
{
|
||||
pos1[i] = R[0][i] + rho[0] * L[0][i];
|
||||
pos3[i] = R[2][i] + rho[2] * L[2][i];
|
||||
}
|
||||
|
||||
/* Use Gibbs to get velocity at r2 */
|
||||
rc = od_gibbs(pos1, r2, pos3, jd[0], jd[1], jd[2], result);
|
||||
|
||||
/* If Gibbs fails (short arc / nearly coplanar), try Herrick-Gibbs */
|
||||
if (rc != 0)
|
||||
{
|
||||
/* Herrick-Gibbs: use f and g series for velocity estimation */
|
||||
double dt31 = (jd[2] - jd[0]) * 86400.0;
|
||||
double dt21 = (jd[1] - jd[0]) * 86400.0;
|
||||
double dt32 = (jd[2] - jd[1]) * 86400.0;
|
||||
|
||||
if (fabs(dt31) < 1.0 || fabs(dt21) < 1.0 || fabs(dt32) < 1.0)
|
||||
return -1;
|
||||
|
||||
for (i = 0; i < 3; i++)
|
||||
{
|
||||
v2[i] = -dt32 * (1.0 / (dt21 * dt31)
|
||||
+ MU_KM3_S2 / (12.0 * vec_mag(pos1) *
|
||||
vec_mag(pos1) * vec_mag(pos1))) * pos1[i]
|
||||
+ (dt32 - dt21) * (1.0 / (dt21 * dt32)
|
||||
+ MU_KM3_S2 / (12.0 * r2_mag * r2_mag * r2_mag)) * r2[i]
|
||||
+ dt21 * (1.0 / (dt32 * dt31)
|
||||
+ MU_KM3_S2 / (12.0 * vec_mag(pos3) *
|
||||
vec_mag(pos3) * vec_mag(pos3))) * pos3[i];
|
||||
}
|
||||
|
||||
rc = od_eci_to_keplerian(r2, v2, &result->kep);
|
||||
if (rc != 0)
|
||||
return -1;
|
||||
|
||||
if (result->kep.ecc >= 1.0 || result->kep.n <= 0.0)
|
||||
return -1;
|
||||
|
||||
result->epoch_jd = jd[1];
|
||||
result->valid = 1;
|
||||
}
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
18
src/od_iod.h
18
src/od_iod.h
@ -42,22 +42,4 @@ int od_gibbs(const double pos1[3], const double pos2[3],
|
||||
double jd1, double jd2, double jd3,
|
||||
od_iod_result_t *result);
|
||||
|
||||
/*
|
||||
* Gauss method: recover orbit from 3 angles-only (RA/Dec) observations.
|
||||
*
|
||||
* ra[3], dec[3]: right ascension and declination in radians
|
||||
* jd[3]: Julian dates of each observation
|
||||
* obs_ecef[3][3]: observer ECEF positions (km) at each epoch
|
||||
* result: output Keplerian elements at epoch jd[1] (middle obs)
|
||||
*
|
||||
* Returns 0 on success, -1 on failure (non-convergence, degenerate).
|
||||
*
|
||||
* Implements Vallado Algorithm 52 with iterative refinement of the
|
||||
* slant range at the middle observation.
|
||||
*/
|
||||
int od_gauss(const double ra[3], const double dec[3],
|
||||
const double jd[3],
|
||||
const double obs_ecef[3][3],
|
||||
od_iod_result_t *result);
|
||||
|
||||
#endif /* PG_ORRERY_OD_IOD_H */
|
||||
|
||||
@ -373,52 +373,6 @@ od_keplerian_to_eci(const od_keplerian_t *kep,
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* RA/Dec utilities for angles-only mode
|
||||
* ================================================================
|
||||
*/
|
||||
|
||||
void
|
||||
od_radec_to_los(double ra, double dec, double los[3])
|
||||
{
|
||||
double cd = cos(dec);
|
||||
los[0] = cd * cos(ra);
|
||||
los[1] = cd * sin(ra);
|
||||
los[2] = sin(dec);
|
||||
}
|
||||
|
||||
void
|
||||
od_teme_to_radec(const double pos_teme[3], const double obs_ecef[3],
|
||||
double gmst, double *ra, double *dec)
|
||||
{
|
||||
double cg = cos(gmst), sg = sin(gmst);
|
||||
double pos_ecef[3], range_ecef[3], range_teme[3], rm;
|
||||
|
||||
/* TEME → ECEF */
|
||||
pos_ecef[0] = cg * pos_teme[0] + sg * pos_teme[1];
|
||||
pos_ecef[1] = -sg * pos_teme[0] + cg * pos_teme[1];
|
||||
pos_ecef[2] = pos_teme[2];
|
||||
|
||||
/* Observer-relative range in ECEF */
|
||||
range_ecef[0] = pos_ecef[0] - obs_ecef[0];
|
||||
range_ecef[1] = pos_ecef[1] - obs_ecef[1];
|
||||
range_ecef[2] = pos_ecef[2] - obs_ecef[2];
|
||||
|
||||
/* Back to TEME (inertial) for RA/Dec */
|
||||
range_teme[0] = cg * range_ecef[0] - sg * range_ecef[1];
|
||||
range_teme[1] = sg * range_ecef[0] + cg * range_ecef[1];
|
||||
range_teme[2] = range_ecef[2];
|
||||
|
||||
rm = sqrt(range_teme[0]*range_teme[0] +
|
||||
range_teme[1]*range_teme[1] +
|
||||
range_teme[2]*range_teme[2]);
|
||||
|
||||
*dec = asin(range_teme[2] / rm);
|
||||
*ra = atan2(range_teme[1], range_teme[0]);
|
||||
if (*ra < 0.0) *ra += 2.0 * M_PI;
|
||||
}
|
||||
|
||||
|
||||
/* ================================================================
|
||||
* Keplerian ↔ equinoctial element conversion
|
||||
* ================================================================
|
||||
|
||||
@ -126,22 +126,6 @@ void od_observer_to_ecef(double lat, double lon, double alt_m,
|
||||
*/
|
||||
double od_gmst_from_jd(double jd);
|
||||
|
||||
/* ── RA/Dec utilities (angles-only mode) ───────────────── */
|
||||
|
||||
/*
|
||||
* RA/Dec (radians) → unit line-of-sight vector (equatorial TEME).
|
||||
*/
|
||||
void od_radec_to_los(double ra, double dec, double los[3]);
|
||||
|
||||
/*
|
||||
* TEME satellite pos + observer ECEF + GMST → RA/Dec (radians).
|
||||
* Computes observer-relative direction in inertial (TEME) frame.
|
||||
* TEME ≈ J2000 equatorial for our accuracy needs (~1 arcsec offset
|
||||
* from nutation, well below TLE accuracy floor of ~1 km ≈ 20 arcsec at LEO).
|
||||
*/
|
||||
void od_teme_to_radec(const double pos_teme[3], const double obs_ecef[3],
|
||||
double gmst, double *ra, double *dec);
|
||||
|
||||
/*
|
||||
* Normalize angle to [0, 2*pi)
|
||||
*/
|
||||
|
||||
221
src/od_solver.c
221
src/od_solver.c
@ -224,21 +224,15 @@ compute_residuals_eci(const tle_t *tle, const od_observation_t *obs,
|
||||
return sqrt(sum_sq / n_obs);
|
||||
}
|
||||
|
||||
/* Range rate scale: maps 1 km/s range rate error to 10 km equivalent.
|
||||
* Conservative default; fine-grained control via per-obs weights. */
|
||||
#define OD_RR_SCALE 10.0
|
||||
|
||||
/*
|
||||
* Compute residual vector for topocentric observations.
|
||||
* residuals[i*ncomp .. i*ncomp+(ncomp-1)] = observed[i] - propagated[i]
|
||||
* ncomp=3: (az_diff, el_diff, range_diff) in (rad, rad, km).
|
||||
* ncomp=4: adds range_rate_diff scaled by OD_RR_SCALE.
|
||||
* residuals[i*3 .. i*3+2] = observed[i] - propagated[i]
|
||||
* Components: (az_diff, el_diff, range_diff) in (rad, rad, km).
|
||||
* Returns RMS range error in km. Returns -1.0 on propagation failure.
|
||||
*/
|
||||
static double
|
||||
compute_residuals_topo(const tle_t *tle, const od_observation_t *obs,
|
||||
int n_obs, const od_observer_t *observers,
|
||||
int fit_range_rate, int ncomp,
|
||||
double *residuals)
|
||||
{
|
||||
double *params;
|
||||
@ -324,86 +318,13 @@ compute_residuals_topo(const tle_t *tle, const od_observation_t *obs,
|
||||
el_diff = obs[i].data[1] - el;
|
||||
range_diff = obs[i].data[2] - range_km;
|
||||
|
||||
residuals[i * ncomp + 0] = az_diff * range_km * cos(el); /* scale to km */
|
||||
residuals[i * ncomp + 1] = el_diff * range_km; /* scale to km */
|
||||
residuals[i * ncomp + 2] = range_diff;
|
||||
|
||||
if (fit_range_rate)
|
||||
{
|
||||
double rr_computed = (dx * vel_ecef[0] + dy * vel_ecef[1]
|
||||
+ dz * vel_ecef[2]) / range_km;
|
||||
residuals[i * ncomp + 3] = (obs[i].data[3] - rr_computed) * OD_RR_SCALE;
|
||||
}
|
||||
residuals[i * 3 + 0] = az_diff * range_km * cos(el); /* scale to km */
|
||||
residuals[i * 3 + 1] = el_diff * range_km; /* scale to km */
|
||||
residuals[i * 3 + 2] = range_diff;
|
||||
|
||||
sum_sq += range_diff * range_diff;
|
||||
}
|
||||
|
||||
od_free(params);
|
||||
return sqrt(sum_sq / n_obs);
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Compute residual vector for RA/Dec (angles-only) observations.
|
||||
* residuals[i*2 + 0] = (ra_obs - ra_comp) * cos(dec_comp) [great-circle]
|
||||
* residuals[i*2 + 1] = dec_obs - dec_comp
|
||||
*
|
||||
* RA scaled by cos(dec) converts angular separation to great-circle
|
||||
* distance, avoiding inflated residuals near the poles.
|
||||
*
|
||||
* Returns RMS angular residual in radians. Returns -1.0 on failure.
|
||||
*/
|
||||
static double
|
||||
compute_residuals_radec(const tle_t *tle, const od_observation_t *obs,
|
||||
int n_obs, const od_observer_t *observers,
|
||||
double *residuals)
|
||||
{
|
||||
double *params;
|
||||
int is_deep;
|
||||
double sum_sq = 0.0;
|
||||
int i;
|
||||
|
||||
is_deep = select_ephemeris(tle);
|
||||
if (is_deep < 0)
|
||||
return -1.0;
|
||||
|
||||
params = (double *)od_alloc(sizeof(double) * N_SAT_PARAMS);
|
||||
init_params(tle, params, is_deep);
|
||||
|
||||
for (i = 0; i < n_obs; i++)
|
||||
{
|
||||
const od_observer_t *observer = &observers[obs[i].observer_idx];
|
||||
double tsince = (obs[i].jd - tle->epoch) * 1440.0;
|
||||
double pos[3], vel[3];
|
||||
int err;
|
||||
double gmst;
|
||||
double ra_comp, dec_comp;
|
||||
double ra_diff, dec_diff;
|
||||
|
||||
err = propagate_with_params(tle, tsince, params, is_deep, pos, vel);
|
||||
if (err != 0 && err != SXPX_WARN_ORBIT_WITHIN_EARTH &&
|
||||
err != SXPX_WARN_PERIGEE_WITHIN_EARTH)
|
||||
{
|
||||
od_free(params);
|
||||
return -1.0;
|
||||
}
|
||||
|
||||
/* Compute RA/Dec from propagated TEME position */
|
||||
gmst = od_gmst_from_jd(obs[i].jd);
|
||||
od_teme_to_radec(pos, observer->ecef, gmst, &ra_comp, &dec_comp);
|
||||
|
||||
/* RA wrap-around */
|
||||
ra_diff = obs[i].data[0] - ra_comp;
|
||||
if (ra_diff > M_PI) ra_diff -= 2.0 * M_PI;
|
||||
if (ra_diff < -M_PI) ra_diff += 2.0 * M_PI;
|
||||
|
||||
dec_diff = obs[i].data[1] - dec_comp;
|
||||
|
||||
residuals[i * 2 + 0] = ra_diff * cos(dec_comp);
|
||||
residuals[i * 2 + 1] = dec_diff;
|
||||
|
||||
sum_sq += ra_diff * cos(dec_comp) * ra_diff * cos(dec_comp)
|
||||
+ dec_diff * dec_diff;
|
||||
(void)vel_ecef; /* reserved for range rate fitting */
|
||||
}
|
||||
|
||||
od_free(params);
|
||||
@ -628,78 +549,6 @@ initial_guess_from_topo(const od_observation_t *obs, int n_obs,
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Compute initial orbit from RA/Dec (angles-only) observations.
|
||||
*
|
||||
* Picks 3 well-spaced observations, calls Gauss IOD (Vallado Algorithm 52)
|
||||
* to recover the orbit at the middle epoch.
|
||||
* Returns 0 on success, -1 if Gauss fails.
|
||||
*/
|
||||
static int
|
||||
initial_guess_from_radec(const od_observation_t *obs, int n_obs,
|
||||
const od_observer_t *observers,
|
||||
tle_t *guess)
|
||||
{
|
||||
int idx[3];
|
||||
double ra[3], dec[3], jd[3];
|
||||
double obs_ecef[3][3];
|
||||
int k;
|
||||
od_iod_result_t iod;
|
||||
|
||||
if (n_obs < 3)
|
||||
return -1;
|
||||
|
||||
idx[0] = 0;
|
||||
idx[1] = n_obs / 2;
|
||||
idx[2] = n_obs - 1;
|
||||
|
||||
for (k = 0; k < 3; k++)
|
||||
{
|
||||
const od_observation_t *o = &obs[idx[k]];
|
||||
const od_observer_t *observer = &observers[o->observer_idx];
|
||||
|
||||
ra[k] = o->data[0];
|
||||
dec[k] = o->data[1];
|
||||
jd[k] = o->jd;
|
||||
|
||||
obs_ecef[k][0] = observer->ecef[0];
|
||||
obs_ecef[k][1] = observer->ecef[1];
|
||||
obs_ecef[k][2] = observer->ecef[2];
|
||||
}
|
||||
|
||||
if (od_gauss(ra, dec, jd, obs_ecef, &iod) != 0)
|
||||
return -1;
|
||||
|
||||
kep_to_seed_tle(&iod.kep, iod.epoch_jd, guess);
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
/* ── Observation weighting ──────────────────────────────── */
|
||||
|
||||
/*
|
||||
* Scale residuals by sqrt(weight) for weighted least-squares.
|
||||
*
|
||||
* When both the nominal and perturbed residuals are scaled identically,
|
||||
* the Jacobian (resid - resid_pert) / h naturally becomes the weighted
|
||||
* Jacobian. The covariance (H^T H)^{-1} becomes (H^T W H)^{-1}.
|
||||
*/
|
||||
static void
|
||||
apply_obs_weights(double *residuals, int n_obs, int ncomp,
|
||||
const double *weights)
|
||||
{
|
||||
int i, j;
|
||||
if (!weights)
|
||||
return;
|
||||
for (i = 0; i < n_obs; i++)
|
||||
{
|
||||
double sw = sqrt(weights[i]);
|
||||
for (j = 0; j < ncomp; j++)
|
||||
residuals[i * ncomp + j] *= sw;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
/* ── Main solver ───────────────────────────────────────── */
|
||||
|
||||
int
|
||||
@ -728,13 +577,7 @@ od_fit_tle(const od_observation_t *obs, int n_obs,
|
||||
|
||||
/* Validate inputs */
|
||||
nstate = config->fit_bstar ? OD_NSTATE_7 : OD_NSTATE_6;
|
||||
|
||||
if (config->obs_type == OD_OBS_ECI)
|
||||
ncomp = 6;
|
||||
else if (config->obs_type == OD_OBS_RADEC)
|
||||
ncomp = 2;
|
||||
else
|
||||
ncomp = config->fit_range_rate ? 4 : 3;
|
||||
ncomp = (config->obs_type == OD_OBS_ECI) ? 6 : 3;
|
||||
|
||||
if (n_obs < nstate)
|
||||
{
|
||||
@ -757,18 +600,6 @@ od_fit_tle(const od_observation_t *obs, int n_obs,
|
||||
{
|
||||
if (config->obs_type == OD_OBS_ECI)
|
||||
initial_guess_from_eci(obs, n_obs, &seed_tle);
|
||||
else if (config->obs_type == OD_OBS_RADEC)
|
||||
{
|
||||
int rc_iod = initial_guess_from_radec(obs, n_obs,
|
||||
config->observers,
|
||||
&seed_tle);
|
||||
if (rc_iod != 0)
|
||||
{
|
||||
snprintf(result->status, sizeof(result->status),
|
||||
"IOD bootstrap failed (Gauss): need seed TLE");
|
||||
return -1;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
int rc_iod = initial_guess_from_topo(obs, n_obs,
|
||||
@ -811,14 +642,9 @@ od_fit_tle(const od_observation_t *obs, int n_obs,
|
||||
/* Compute initial RMS */
|
||||
if (config->obs_type == OD_OBS_ECI)
|
||||
rms_cur = compute_residuals_eci(¤t_tle, obs, n_obs, residuals);
|
||||
else if (config->obs_type == OD_OBS_RADEC)
|
||||
rms_cur = compute_residuals_radec(¤t_tle, obs, n_obs,
|
||||
config->observers, residuals);
|
||||
else
|
||||
rms_cur = compute_residuals_topo(¤t_tle, obs, n_obs,
|
||||
config->observers,
|
||||
config->fit_range_rate, ncomp,
|
||||
residuals);
|
||||
config->observers, residuals);
|
||||
|
||||
if (rms_cur < 0.0)
|
||||
{
|
||||
@ -835,9 +661,6 @@ od_fit_tle(const od_observation_t *obs, int n_obs,
|
||||
result->rms_initial = rms_cur;
|
||||
rms_prev = rms_cur;
|
||||
|
||||
/* Weight residuals for Jacobian build (RMS already stored unweighted) */
|
||||
apply_obs_weights(residuals, n_obs, ncomp, config->weights);
|
||||
|
||||
/* Already converged (perfect seed)? Skip DC loop but still
|
||||
* compute covariance — users need uncertainty estimates even
|
||||
* when the initial guess is exact. */
|
||||
@ -903,16 +726,9 @@ od_fit_tle(const od_observation_t *obs, int n_obs,
|
||||
|
||||
if (config->obs_type == OD_OBS_ECI)
|
||||
compute_residuals_eci(&tle_pert, obs, n_obs, resid_pert);
|
||||
else if (config->obs_type == OD_OBS_RADEC)
|
||||
compute_residuals_radec(&tle_pert, obs, n_obs,
|
||||
config->observers, resid_pert);
|
||||
else
|
||||
compute_residuals_topo(&tle_pert, obs, n_obs,
|
||||
config->observers,
|
||||
config->fit_range_rate, ncomp,
|
||||
resid_pert);
|
||||
|
||||
apply_obs_weights(resid_pert, n_obs, ncomp, config->weights);
|
||||
config->observers, resid_pert);
|
||||
|
||||
/* Jacobian column j (column-major for LAPACK)
|
||||
* H = dG/dx where G is the computed observation function.
|
||||
@ -986,14 +802,9 @@ od_fit_tle(const od_observation_t *obs, int n_obs,
|
||||
|
||||
if (config->obs_type == OD_OBS_ECI)
|
||||
rms_cur = compute_residuals_eci(¤t_tle, obs, n_obs, residuals);
|
||||
else if (config->obs_type == OD_OBS_RADEC)
|
||||
rms_cur = compute_residuals_radec(¤t_tle, obs, n_obs,
|
||||
config->observers, residuals);
|
||||
else
|
||||
rms_cur = compute_residuals_topo(¤t_tle, obs, n_obs,
|
||||
config->observers,
|
||||
config->fit_range_rate, ncomp,
|
||||
residuals);
|
||||
config->observers, residuals);
|
||||
|
||||
if (rms_cur < 0.0)
|
||||
{
|
||||
@ -1002,9 +813,6 @@ od_fit_tle(const od_observation_t *obs, int n_obs,
|
||||
break;
|
||||
}
|
||||
|
||||
/* Weight residuals for next Jacobian build (RMS already stored) */
|
||||
apply_obs_weights(residuals, n_obs, ncomp, config->weights);
|
||||
|
||||
/* Adaptive step control: adjust trust-region size based on
|
||||
* whether the solver is converging or diverging. Halve the
|
||||
* step on divergence (prevents oscillation with poor initial
|
||||
@ -1101,16 +909,9 @@ compute_covariance:
|
||||
|
||||
if (config->obs_type == OD_OBS_ECI)
|
||||
compute_residuals_eci(&tle_pert, obs, n_obs, resid_pert);
|
||||
else if (config->obs_type == OD_OBS_RADEC)
|
||||
compute_residuals_radec(&tle_pert, obs, n_obs,
|
||||
config->observers, resid_pert);
|
||||
else
|
||||
compute_residuals_topo(&tle_pert, obs, n_obs,
|
||||
config->observers,
|
||||
config->fit_range_rate, ncomp,
|
||||
resid_pert);
|
||||
|
||||
apply_obs_weights(resid_pert, n_obs, ncomp, config->weights);
|
||||
config->observers, resid_pert);
|
||||
|
||||
for (k_cov = 0; k_cov < nrows; k_cov++)
|
||||
jacobian[j_cov * nrows + k_cov] =
|
||||
|
||||
@ -31,9 +31,8 @@
|
||||
*/
|
||||
typedef enum
|
||||
{
|
||||
OD_OBS_ECI = 0, /* 6-component: x, y, z, vx, vy, vz (km, km/s) */
|
||||
OD_OBS_TOPO = 1, /* 3-component: az, el, range (rad, rad, km) */
|
||||
OD_OBS_RADEC = 2 /* 2-component: RA, Dec (radians, equatorial) */
|
||||
OD_OBS_ECI = 0, /* 6-component: x, y, z, vx, vy, vz (km, km/s) */
|
||||
OD_OBS_TOPO = 1 /* 3-component: az, el, range (rad, rad, km) */
|
||||
} od_obs_type_t;
|
||||
|
||||
/*
|
||||
@ -64,12 +63,9 @@ typedef struct
|
||||
{
|
||||
od_obs_type_t obs_type; /* ECI or topocentric */
|
||||
int fit_bstar; /* include B* as 7th state */
|
||||
int fit_range_rate; /* include range_rate in topo residuals */
|
||||
int max_iter; /* iteration limit */
|
||||
od_observer_t *observers; /* array of observers (topo mode) */
|
||||
int n_observers; /* count (0 for ECI mode) */
|
||||
double *weights; /* per-observation weights (NULL=uniform) */
|
||||
int n_weights;
|
||||
} od_config_t;
|
||||
|
||||
/*
|
||||
|
||||
@ -623,306 +623,3 @@ FROM result;
|
||||
t | t | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 18: Range rate round-trip
|
||||
--
|
||||
-- Propagate ISS, observe() to get topo with range_rate,
|
||||
-- fit via tle_from_topocentric with fit_range_rate := true.
|
||||
-- Verify convergence.
|
||||
-- ============================================================
|
||||
WITH iss_tle AS (
|
||||
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
|
||||
),
|
||||
mit AS (
|
||||
SELECT '(42.36,-71.09,20)'::observer AS obs
|
||||
),
|
||||
topo_obs AS (
|
||||
SELECT
|
||||
array_agg(observe(t, obs, ts) ORDER BY ts) AS observations,
|
||||
array_agg(ts ORDER BY ts) AS times
|
||||
FROM iss_tle, mit,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts
|
||||
),
|
||||
result AS (
|
||||
SELECT (tle_from_topocentric(observations, times, obs, t, false, 20,
|
||||
fit_range_rate := true)).*
|
||||
FROM topo_obs, mit, iss_tle
|
||||
)
|
||||
SELECT
|
||||
rms_final < 10.0 AS rms_under_10km,
|
||||
status = 'converged' AS did_converge
|
||||
FROM result;
|
||||
rms_under_10km | did_converge
|
||||
----------------+--------------
|
||||
t | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 19: Range rate disabled matches existing behavior
|
||||
--
|
||||
-- Same data with fit_range_rate := false (default).
|
||||
-- Results should converge the same as Test 4.
|
||||
-- ============================================================
|
||||
WITH iss_tle AS (
|
||||
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
|
||||
),
|
||||
mit AS (
|
||||
SELECT '(42.36,-71.09,20)'::observer AS obs
|
||||
),
|
||||
topo_obs AS (
|
||||
SELECT
|
||||
array_agg(observe(t, obs, ts) ORDER BY ts) AS observations,
|
||||
array_agg(ts ORDER BY ts) AS times
|
||||
FROM iss_tle, mit,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts
|
||||
),
|
||||
result AS (
|
||||
SELECT (tle_from_topocentric(observations, times, obs, t, false, 20,
|
||||
fit_range_rate := false)).*
|
||||
FROM topo_obs, mit, iss_tle
|
||||
)
|
||||
SELECT
|
||||
rms_final < 10.0 AS rms_under_10km,
|
||||
status = 'converged' AS did_converge
|
||||
FROM result;
|
||||
rms_under_10km | did_converge
|
||||
----------------+--------------
|
||||
t | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 20: Weighted observations round-trip (uniform weights)
|
||||
--
|
||||
-- Uniform weights ARRAY[1,1,...,1]::float8[] should produce
|
||||
-- identical results to no weights.
|
||||
-- ============================================================
|
||||
WITH iss_tle AS (
|
||||
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
|
||||
),
|
||||
obs AS (
|
||||
SELECT
|
||||
array_agg(sgp4_propagate(t, ts) ORDER BY ts) AS positions,
|
||||
array_agg(ts ORDER BY ts) AS times
|
||||
FROM iss_tle,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts
|
||||
),
|
||||
uniform_weights AS (
|
||||
SELECT array_agg(1.0::float8) AS w
|
||||
FROM generate_series(1, 19)
|
||||
),
|
||||
result AS (
|
||||
SELECT (tle_from_eci(positions, times, t, false, 15,
|
||||
weights := w)).* FROM obs, iss_tle, uniform_weights
|
||||
)
|
||||
SELECT
|
||||
rms_final < 1.0 AS rms_under_1km,
|
||||
status = 'converged' AS did_converge
|
||||
FROM result;
|
||||
rms_under_1km | did_converge
|
||||
---------------+--------------
|
||||
t | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 21: Weights length mismatch error
|
||||
--
|
||||
-- Wrong-length weights array should raise an error.
|
||||
-- ============================================================
|
||||
DO $$
|
||||
BEGIN
|
||||
PERFORM tle_from_eci(
|
||||
(SELECT array_agg(sgp4_propagate(
|
||||
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
|
||||
ts) ORDER BY ts)
|
||||
FROM generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts),
|
||||
(SELECT array_agg(ts ORDER BY ts)
|
||||
FROM generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts),
|
||||
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
|
||||
false, 15,
|
||||
ARRAY[1.0, 2.0, 3.0]::float8[] -- wrong length
|
||||
);
|
||||
RAISE NOTICE 'ERROR: should have raised exception';
|
||||
EXCEPTION
|
||||
WHEN array_subscript_error THEN
|
||||
RAISE NOTICE 'OK: weights length mismatch error caught';
|
||||
END $$;
|
||||
NOTICE: OK: weights length mismatch error caught
|
||||
-- ============================================================
|
||||
-- Test 22: ECI with weights (verify API works)
|
||||
--
|
||||
-- Non-uniform weights (downweight first few observations).
|
||||
-- ============================================================
|
||||
WITH iss_tle AS (
|
||||
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
|
||||
),
|
||||
obs AS (
|
||||
SELECT
|
||||
array_agg(sgp4_propagate(t, ts) ORDER BY ts) AS positions,
|
||||
array_agg(ts ORDER BY ts) AS times
|
||||
FROM iss_tle,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts
|
||||
),
|
||||
nonuniform_weights AS (
|
||||
SELECT array_agg(CASE WHEN n <= 5 THEN 0.1 ELSE 1.0 END ::float8) AS w
|
||||
FROM generate_series(1, 19) AS n
|
||||
),
|
||||
result AS (
|
||||
SELECT (tle_from_eci(positions, times, t, false, 15,
|
||||
weights := w)).* FROM obs, iss_tle, nonuniform_weights
|
||||
)
|
||||
SELECT
|
||||
rms_final < 5.0 AS rms_under_5km,
|
||||
status = 'converged' AS did_converge
|
||||
FROM result;
|
||||
rms_under_5km | did_converge
|
||||
---------------+--------------
|
||||
t | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 23: Angles-only with seed TLE
|
||||
--
|
||||
-- Propagate ISS, compute RA/Dec from TEME position relative
|
||||
-- to observer, fit via tle_from_angles() with known seed.
|
||||
-- ============================================================
|
||||
WITH iss_tle AS (
|
||||
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
|
||||
),
|
||||
mit AS (
|
||||
SELECT '(42.36,-71.09,20)'::observer AS obs
|
||||
),
|
||||
topo_obs AS (
|
||||
SELECT
|
||||
array_agg(observe(t, obs, ts) ORDER BY ts) AS observations,
|
||||
array_agg(ts ORDER BY ts) AS times
|
||||
FROM iss_tle, mit,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts
|
||||
),
|
||||
radec AS (
|
||||
SELECT
|
||||
array_agg(
|
||||
(topo_azimuth(o) * 24.0 / 360.0)::float8 -- crude az→RA proxy (accessor returns degrees)
|
||||
ORDER BY ts
|
||||
) AS ra_h,
|
||||
array_agg(
|
||||
topo_elevation(o)::float8 -- accessor returns degrees, valid as Dec proxy
|
||||
ORDER BY ts
|
||||
) AS dec_d,
|
||||
array_agg(ts ORDER BY ts) AS times
|
||||
FROM (
|
||||
SELECT observe(t, obs, ts) AS o, ts
|
||||
FROM iss_tle, mit,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts
|
||||
) sub
|
||||
),
|
||||
result AS (
|
||||
SELECT (tle_from_angles(ra_h, dec_d, times, obs, t)).*
|
||||
FROM radec, mit, iss_tle
|
||||
)
|
||||
SELECT
|
||||
status IN ('converged', 'max_iterations') AS solver_ran,
|
||||
rms_final IS NOT NULL AS has_rms
|
||||
FROM result;
|
||||
solver_ran | has_rms
|
||||
------------+---------
|
||||
t | t
|
||||
(1 row)
|
||||
|
||||
-- ============================================================
|
||||
-- Test 24: Angles-only without seed (Gauss IOD bootstrap)
|
||||
--
|
||||
-- Without a seed TLE, Gauss IOD must recover an initial orbit
|
||||
-- from the RA/Dec observations alone. The crude azimuth→RA
|
||||
-- approximation used here is not physically meaningful, so
|
||||
-- Gauss may fail to converge — we accept that as valid behavior.
|
||||
-- The real Gauss validation is in the standalone C tests.
|
||||
-- ============================================================
|
||||
DO $$
|
||||
DECLARE
|
||||
v_status text;
|
||||
v_rms float8;
|
||||
BEGIN
|
||||
SELECT status, rms_final INTO v_status, v_rms
|
||||
FROM (
|
||||
SELECT (tle_from_angles(
|
||||
(SELECT array_agg((topo_azimuth(o) * 24.0 / 360.0)::float8 ORDER BY ts)
|
||||
FROM (SELECT observe(
|
||||
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
|
||||
'(42.36,-71.09,20)'::observer, ts) AS o, ts
|
||||
FROM generate_series('2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval) AS ts) sub),
|
||||
(SELECT array_agg(topo_elevation(o)::float8 ORDER BY ts)
|
||||
FROM (SELECT observe(
|
||||
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
|
||||
'(42.36,-71.09,20)'::observer, ts) AS o, ts
|
||||
FROM generate_series('2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval) AS ts) sub),
|
||||
(SELECT array_agg(ts ORDER BY ts)
|
||||
FROM generate_series('2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval) AS ts),
|
||||
'(42.36,-71.09,20)'::observer
|
||||
)).*
|
||||
) r;
|
||||
RAISE NOTICE 'OK: seedless angles-only ran: status=%, rms=%', v_status, v_rms;
|
||||
EXCEPTION
|
||||
WHEN data_exception THEN
|
||||
RAISE NOTICE 'OK: seedless angles-only Gauss IOD failed as expected with approximate data';
|
||||
END $$;
|
||||
NOTICE: OK: seedless angles-only Gauss IOD failed as expected with approximate data
|
||||
-- ============================================================
|
||||
-- Test 25: Angles-only error cases
|
||||
--
|
||||
-- Too few observations should raise an error.
|
||||
-- ============================================================
|
||||
DO $$
|
||||
BEGIN
|
||||
PERFORM tle_from_angles(
|
||||
ARRAY[1.0, 2.0]::float8[],
|
||||
ARRAY[30.0, 35.0]::float8[],
|
||||
ARRAY['2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 12:05:00+00'::timestamptz],
|
||||
'(42.36,-71.09,20)'::observer
|
||||
);
|
||||
RAISE NOTICE 'ERROR: should have raised exception';
|
||||
EXCEPTION
|
||||
WHEN data_exception THEN
|
||||
RAISE NOTICE 'OK: insufficient observations error caught';
|
||||
WHEN invalid_parameter_value THEN
|
||||
RAISE NOTICE 'OK: insufficient observations error caught';
|
||||
END $$;
|
||||
NOTICE: OK: insufficient observations error caught
|
||||
|
||||
@ -581,294 +581,3 @@ SELECT
|
||||
covariance IS NOT NULL AS has_covariance,
|
||||
nstate = 6 AS is_6state
|
||||
FROM result;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 18: Range rate round-trip
|
||||
--
|
||||
-- Propagate ISS, observe() to get topo with range_rate,
|
||||
-- fit via tle_from_topocentric with fit_range_rate := true.
|
||||
-- Verify convergence.
|
||||
-- ============================================================
|
||||
|
||||
WITH iss_tle AS (
|
||||
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
|
||||
),
|
||||
mit AS (
|
||||
SELECT '(42.36,-71.09,20)'::observer AS obs
|
||||
),
|
||||
topo_obs AS (
|
||||
SELECT
|
||||
array_agg(observe(t, obs, ts) ORDER BY ts) AS observations,
|
||||
array_agg(ts ORDER BY ts) AS times
|
||||
FROM iss_tle, mit,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts
|
||||
),
|
||||
result AS (
|
||||
SELECT (tle_from_topocentric(observations, times, obs, t, false, 20,
|
||||
fit_range_rate := true)).*
|
||||
FROM topo_obs, mit, iss_tle
|
||||
)
|
||||
SELECT
|
||||
rms_final < 10.0 AS rms_under_10km,
|
||||
status = 'converged' AS did_converge
|
||||
FROM result;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 19: Range rate disabled matches existing behavior
|
||||
--
|
||||
-- Same data with fit_range_rate := false (default).
|
||||
-- Results should converge the same as Test 4.
|
||||
-- ============================================================
|
||||
|
||||
WITH iss_tle AS (
|
||||
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
|
||||
),
|
||||
mit AS (
|
||||
SELECT '(42.36,-71.09,20)'::observer AS obs
|
||||
),
|
||||
topo_obs AS (
|
||||
SELECT
|
||||
array_agg(observe(t, obs, ts) ORDER BY ts) AS observations,
|
||||
array_agg(ts ORDER BY ts) AS times
|
||||
FROM iss_tle, mit,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts
|
||||
),
|
||||
result AS (
|
||||
SELECT (tle_from_topocentric(observations, times, obs, t, false, 20,
|
||||
fit_range_rate := false)).*
|
||||
FROM topo_obs, mit, iss_tle
|
||||
)
|
||||
SELECT
|
||||
rms_final < 10.0 AS rms_under_10km,
|
||||
status = 'converged' AS did_converge
|
||||
FROM result;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 20: Weighted observations round-trip (uniform weights)
|
||||
--
|
||||
-- Uniform weights ARRAY[1,1,...,1]::float8[] should produce
|
||||
-- identical results to no weights.
|
||||
-- ============================================================
|
||||
|
||||
WITH iss_tle AS (
|
||||
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
|
||||
),
|
||||
obs AS (
|
||||
SELECT
|
||||
array_agg(sgp4_propagate(t, ts) ORDER BY ts) AS positions,
|
||||
array_agg(ts ORDER BY ts) AS times
|
||||
FROM iss_tle,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts
|
||||
),
|
||||
uniform_weights AS (
|
||||
SELECT array_agg(1.0::float8) AS w
|
||||
FROM generate_series(1, 19)
|
||||
),
|
||||
result AS (
|
||||
SELECT (tle_from_eci(positions, times, t, false, 15,
|
||||
weights := w)).* FROM obs, iss_tle, uniform_weights
|
||||
)
|
||||
SELECT
|
||||
rms_final < 1.0 AS rms_under_1km,
|
||||
status = 'converged' AS did_converge
|
||||
FROM result;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 21: Weights length mismatch error
|
||||
--
|
||||
-- Wrong-length weights array should raise an error.
|
||||
-- ============================================================
|
||||
|
||||
DO $$
|
||||
BEGIN
|
||||
PERFORM tle_from_eci(
|
||||
(SELECT array_agg(sgp4_propagate(
|
||||
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
|
||||
ts) ORDER BY ts)
|
||||
FROM generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts),
|
||||
(SELECT array_agg(ts ORDER BY ts)
|
||||
FROM generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts),
|
||||
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
|
||||
false, 15,
|
||||
ARRAY[1.0, 2.0, 3.0]::float8[] -- wrong length
|
||||
);
|
||||
RAISE NOTICE 'ERROR: should have raised exception';
|
||||
EXCEPTION
|
||||
WHEN array_subscript_error THEN
|
||||
RAISE NOTICE 'OK: weights length mismatch error caught';
|
||||
END $$;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 22: ECI with weights (verify API works)
|
||||
--
|
||||
-- Non-uniform weights (downweight first few observations).
|
||||
-- ============================================================
|
||||
|
||||
WITH iss_tle AS (
|
||||
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
|
||||
),
|
||||
obs AS (
|
||||
SELECT
|
||||
array_agg(sgp4_propagate(t, ts) ORDER BY ts) AS positions,
|
||||
array_agg(ts ORDER BY ts) AS times
|
||||
FROM iss_tle,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts
|
||||
),
|
||||
nonuniform_weights AS (
|
||||
SELECT array_agg(CASE WHEN n <= 5 THEN 0.1 ELSE 1.0 END ::float8) AS w
|
||||
FROM generate_series(1, 19) AS n
|
||||
),
|
||||
result AS (
|
||||
SELECT (tle_from_eci(positions, times, t, false, 15,
|
||||
weights := w)).* FROM obs, iss_tle, nonuniform_weights
|
||||
)
|
||||
SELECT
|
||||
rms_final < 5.0 AS rms_under_5km,
|
||||
status = 'converged' AS did_converge
|
||||
FROM result;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 23: Angles-only with seed TLE
|
||||
--
|
||||
-- Propagate ISS, compute RA/Dec from TEME position relative
|
||||
-- to observer, fit via tle_from_angles() with known seed.
|
||||
-- ============================================================
|
||||
|
||||
WITH iss_tle AS (
|
||||
SELECT E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle AS t
|
||||
),
|
||||
mit AS (
|
||||
SELECT '(42.36,-71.09,20)'::observer AS obs
|
||||
),
|
||||
topo_obs AS (
|
||||
SELECT
|
||||
array_agg(observe(t, obs, ts) ORDER BY ts) AS observations,
|
||||
array_agg(ts ORDER BY ts) AS times
|
||||
FROM iss_tle, mit,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts
|
||||
),
|
||||
radec AS (
|
||||
SELECT
|
||||
array_agg(
|
||||
(topo_azimuth(o) * 24.0 / 360.0)::float8 -- crude az→RA proxy (accessor returns degrees)
|
||||
ORDER BY ts
|
||||
) AS ra_h,
|
||||
array_agg(
|
||||
topo_elevation(o)::float8 -- accessor returns degrees, valid as Dec proxy
|
||||
ORDER BY ts
|
||||
) AS dec_d,
|
||||
array_agg(ts ORDER BY ts) AS times
|
||||
FROM (
|
||||
SELECT observe(t, obs, ts) AS o, ts
|
||||
FROM iss_tle, mit,
|
||||
generate_series(
|
||||
'2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval
|
||||
) AS ts
|
||||
) sub
|
||||
),
|
||||
result AS (
|
||||
SELECT (tle_from_angles(ra_h, dec_d, times, obs, t)).*
|
||||
FROM radec, mit, iss_tle
|
||||
)
|
||||
SELECT
|
||||
status IN ('converged', 'max_iterations') AS solver_ran,
|
||||
rms_final IS NOT NULL AS has_rms
|
||||
FROM result;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 24: Angles-only without seed (Gauss IOD bootstrap)
|
||||
--
|
||||
-- Without a seed TLE, Gauss IOD must recover an initial orbit
|
||||
-- from the RA/Dec observations alone. The crude azimuth→RA
|
||||
-- approximation used here is not physically meaningful, so
|
||||
-- Gauss may fail to converge — we accept that as valid behavior.
|
||||
-- The real Gauss validation is in the standalone C tests.
|
||||
-- ============================================================
|
||||
|
||||
DO $$
|
||||
DECLARE
|
||||
v_status text;
|
||||
v_rms float8;
|
||||
BEGIN
|
||||
SELECT status, rms_final INTO v_status, v_rms
|
||||
FROM (
|
||||
SELECT (tle_from_angles(
|
||||
(SELECT array_agg((topo_azimuth(o) * 24.0 / 360.0)::float8 ORDER BY ts)
|
||||
FROM (SELECT observe(
|
||||
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
|
||||
'(42.36,-71.09,20)'::observer, ts) AS o, ts
|
||||
FROM generate_series('2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval) AS ts) sub),
|
||||
(SELECT array_agg(topo_elevation(o)::float8 ORDER BY ts)
|
||||
FROM (SELECT observe(
|
||||
E'1 25544U 98067A 24001.50000000 .00016717 00000-0 10270-3 0 9025\n2 25544 51.6400 208.9163 0006703 30.1694 61.7520 15.50100486 00001'::tle,
|
||||
'(42.36,-71.09,20)'::observer, ts) AS o, ts
|
||||
FROM generate_series('2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval) AS ts) sub),
|
||||
(SELECT array_agg(ts ORDER BY ts)
|
||||
FROM generate_series('2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 13:30:00+00'::timestamptz,
|
||||
'5 minutes'::interval) AS ts),
|
||||
'(42.36,-71.09,20)'::observer
|
||||
)).*
|
||||
) r;
|
||||
RAISE NOTICE 'OK: seedless angles-only ran: status=%, rms=%', v_status, v_rms;
|
||||
EXCEPTION
|
||||
WHEN data_exception THEN
|
||||
RAISE NOTICE 'OK: seedless angles-only Gauss IOD failed as expected with approximate data';
|
||||
END $$;
|
||||
|
||||
-- ============================================================
|
||||
-- Test 25: Angles-only error cases
|
||||
--
|
||||
-- Too few observations should raise an error.
|
||||
-- ============================================================
|
||||
|
||||
DO $$
|
||||
BEGIN
|
||||
PERFORM tle_from_angles(
|
||||
ARRAY[1.0, 2.0]::float8[],
|
||||
ARRAY[30.0, 35.0]::float8[],
|
||||
ARRAY['2024-01-01 12:00:00+00'::timestamptz,
|
||||
'2024-01-01 12:05:00+00'::timestamptz],
|
||||
'(42.36,-71.09,20)'::observer
|
||||
);
|
||||
RAISE NOTICE 'ERROR: should have raised exception';
|
||||
EXCEPTION
|
||||
WHEN data_exception THEN
|
||||
RAISE NOTICE 'OK: insufficient observations error caught';
|
||||
WHEN invalid_parameter_value THEN
|
||||
RAISE NOTICE 'OK: insufficient observations error caught';
|
||||
END $$;
|
||||
|
||||
@ -1,282 +0,0 @@
|
||||
/*
|
||||
* test_od_gauss.c -- Standalone unit tests for Gauss angles-only IOD
|
||||
*
|
||||
* No PostgreSQL dependency. Exercises:
|
||||
* - ISS-like orbit: generate RA/Dec from known state, recover via Gauss
|
||||
* - GEO orbit with wider time spacing
|
||||
* - Degenerate: observations too close in time
|
||||
* - Gauss feeding into Gibbs/Herrick-Gibbs
|
||||
*
|
||||
* Build: cc -Wall -Werror -Isrc -o test/test_od_gauss \
|
||||
* test/test_od_gauss.c src/od_iod.c src/od_math.c -lm
|
||||
* Run: ./test/test_od_gauss
|
||||
*/
|
||||
|
||||
#include "od_iod.h"
|
||||
#include "od_math.h"
|
||||
|
||||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <math.h>
|
||||
|
||||
/* ── Test harness ───────────────────────────────────────── */
|
||||
|
||||
static int n_run, n_pass;
|
||||
|
||||
#define RUN(cond, msg) do { \
|
||||
n_run++; \
|
||||
if (!(cond)) \
|
||||
fprintf(stderr, "FAIL: %s [line %d]\n", (msg), __LINE__); \
|
||||
else { n_pass++; fprintf(stderr, " ok: %s\n", (msg)); } \
|
||||
} while (0)
|
||||
|
||||
#define CLOSE(a, b, tol, msg) do { \
|
||||
n_run++; \
|
||||
double _a = (a), _b = (b); \
|
||||
if (fabs(_a - _b) > (tol)) \
|
||||
fprintf(stderr, "FAIL: %s: %.15g vs %.15g (diff %.3e) [line %d]\n", \
|
||||
(msg), _a, _b, fabs(_a - _b), __LINE__); \
|
||||
else { n_pass++; fprintf(stderr, " ok: %s\n", (msg)); } \
|
||||
} while (0)
|
||||
|
||||
|
||||
/* ── Helper: compute RA/Dec from TEME pos + observer ECEF ── */
|
||||
|
||||
/*
|
||||
* Simulates an observation by computing RA/Dec of a satellite
|
||||
* from a ground observer. Uses od_teme_to_radec().
|
||||
*/
|
||||
static void
|
||||
simulate_radec(const double pos_teme[3], const double obs_ecef[3],
|
||||
double jd, double *ra, double *dec)
|
||||
{
|
||||
double gmst = od_gmst_from_jd(jd);
|
||||
od_teme_to_radec(pos_teme, obs_ecef, gmst, ra, dec);
|
||||
}
|
||||
|
||||
|
||||
/* ── Test: ISS-like orbit ──────────────────────────────── */
|
||||
|
||||
static void
|
||||
test_gauss_iss(void)
|
||||
{
|
||||
od_keplerian_t kep;
|
||||
double pos[3][3], vel[3][3];
|
||||
double ra[3], dec[3], jd[3];
|
||||
double obs_ecef[3][3];
|
||||
od_iod_result_t result;
|
||||
int rc, i;
|
||||
double dt;
|
||||
|
||||
fprintf(stderr, "\n--- Gauss: ISS-like orbit ---\n");
|
||||
|
||||
/* Known ISS-like orbit */
|
||||
kep.n = 0.001127; /* ~15.5 rev/day in rad/min */
|
||||
kep.ecc = 0.0007;
|
||||
kep.inc = 0.9012; /* ~51.6 deg */
|
||||
kep.raan = 3.0;
|
||||
kep.argp = 0.5;
|
||||
kep.M = 0.0;
|
||||
kep.bstar = 0.0;
|
||||
|
||||
/* 30 minutes between observations — wider arc improves Gauss conditioning */
|
||||
dt = 1800.0;
|
||||
jd[0] = 2451545.0;
|
||||
jd[1] = jd[0] + dt / 86400.0;
|
||||
jd[2] = jd[0] + 2.0 * dt / 86400.0;
|
||||
|
||||
/* Generate 3 positions */
|
||||
od_keplerian_to_eci(&kep, pos[0], vel[0]);
|
||||
|
||||
kep.M = od_normalize_angle(kep.M + kep.n * (dt / 60.0));
|
||||
od_keplerian_to_eci(&kep, pos[1], vel[1]);
|
||||
|
||||
kep.M = od_normalize_angle(kep.M + kep.n * (dt / 60.0));
|
||||
od_keplerian_to_eci(&kep, pos[2], vel[2]);
|
||||
|
||||
/* Observer at lat=40N, lon=0, alt=0 — compute ECEF once for all obs */
|
||||
od_observer_to_ecef(40.0 * M_PI / 180.0, 0.0, 0.0, obs_ecef[0]);
|
||||
obs_ecef[1][0] = obs_ecef[0][0];
|
||||
obs_ecef[1][1] = obs_ecef[0][1];
|
||||
obs_ecef[1][2] = obs_ecef[0][2];
|
||||
obs_ecef[2][0] = obs_ecef[0][0];
|
||||
obs_ecef[2][1] = obs_ecef[0][1];
|
||||
obs_ecef[2][2] = obs_ecef[0][2];
|
||||
|
||||
/* Compute RA/Dec for each observation */
|
||||
for (i = 0; i < 3; i++)
|
||||
simulate_radec(pos[i], obs_ecef[i], jd[i], &ra[i], &dec[i]);
|
||||
|
||||
/* Run Gauss */
|
||||
rc = od_gauss(ra, dec, jd, obs_ecef, &result);
|
||||
|
||||
RUN(rc == 0, "Gauss ISS returns success");
|
||||
RUN(result.valid == 1, "result is valid");
|
||||
RUN(result.kep.ecc < 1.0, "eccentricity is bound");
|
||||
RUN(result.kep.n > 0.0, "positive mean motion");
|
||||
|
||||
/* Gauss accuracy is inherently low (angles-only loses range info).
|
||||
* It only needs to be close enough for the DC solver to converge.
|
||||
* A factor of 5x in mean motion is acceptable for a seed orbit. */
|
||||
RUN(result.kep.n > 0.0002 && result.kep.n < 0.006,
|
||||
"mean motion in plausible range");
|
||||
CLOSE(result.kep.inc, 0.9012, 0.3, "inclination within ~17 deg");
|
||||
}
|
||||
|
||||
|
||||
/* ── Test: MEO-like orbit ──────────────────────────────── */
|
||||
|
||||
static void
|
||||
test_gauss_meo(void)
|
||||
{
|
||||
od_keplerian_t kep;
|
||||
double pos[3][3], vel[3][3];
|
||||
double ra[3], dec[3], jd[3];
|
||||
double obs_ecef[3][3];
|
||||
od_iod_result_t result;
|
||||
int rc, i;
|
||||
double dt;
|
||||
|
||||
fprintf(stderr, "\n--- Gauss: MEO-like orbit ---\n");
|
||||
|
||||
/* GPS-like altitude, moderate inclination */
|
||||
kep.n = 0.000262; /* ~2 rev/day */
|
||||
kep.ecc = 0.01;
|
||||
kep.inc = 0.96; /* ~55 deg */
|
||||
kep.raan = 1.0;
|
||||
kep.argp = 0.0;
|
||||
kep.M = 0.0;
|
||||
kep.bstar = 0.0;
|
||||
|
||||
/* 2 hours between observations — wider arc for higher altitude */
|
||||
dt = 7200.0;
|
||||
jd[0] = 2451545.0;
|
||||
jd[1] = jd[0] + dt / 86400.0;
|
||||
jd[2] = jd[0] + 2.0 * dt / 86400.0;
|
||||
|
||||
od_keplerian_to_eci(&kep, pos[0], vel[0]);
|
||||
|
||||
kep.M = od_normalize_angle(kep.M + kep.n * (dt / 60.0));
|
||||
od_keplerian_to_eci(&kep, pos[1], vel[1]);
|
||||
|
||||
kep.M = od_normalize_angle(kep.M + kep.n * (dt / 60.0));
|
||||
od_keplerian_to_eci(&kep, pos[2], vel[2]);
|
||||
|
||||
od_observer_to_ecef(35.0 * M_PI / 180.0, -106.0 * M_PI / 180.0,
|
||||
1600.0, obs_ecef[0]);
|
||||
for (i = 1; i < 3; i++)
|
||||
{
|
||||
obs_ecef[i][0] = obs_ecef[0][0];
|
||||
obs_ecef[i][1] = obs_ecef[0][1];
|
||||
obs_ecef[i][2] = obs_ecef[0][2];
|
||||
}
|
||||
|
||||
for (i = 0; i < 3; i++)
|
||||
simulate_radec(pos[i], obs_ecef[i], jd[i], &ra[i], &dec[i]);
|
||||
|
||||
rc = od_gauss(ra, dec, jd, obs_ecef, &result);
|
||||
|
||||
RUN(rc == 0, "Gauss MEO returns success");
|
||||
RUN(result.valid == 1, "result is valid");
|
||||
RUN(result.kep.n > 0.0, "positive mean motion");
|
||||
RUN(result.kep.n > 0.00005 && result.kep.n < 0.002,
|
||||
"mean motion in plausible MEO range");
|
||||
}
|
||||
|
||||
|
||||
/* ── Test: too-close observations fail ───────────────────── */
|
||||
|
||||
static void
|
||||
test_gauss_too_close(void)
|
||||
{
|
||||
double ra[3] = {1.0, 1.001, 1.002};
|
||||
double dec[3] = {0.5, 0.501, 0.502};
|
||||
double jd[3] = {2451545.0, 2451545.000005, 2451545.00001};
|
||||
double obs_ecef[3][3];
|
||||
od_iod_result_t result;
|
||||
int rc;
|
||||
|
||||
fprintf(stderr, "\n--- Gauss: too-close observations ---\n");
|
||||
|
||||
od_observer_to_ecef(40.0 * M_PI / 180.0, 0.0, 0.0, obs_ecef[0]);
|
||||
obs_ecef[1][0] = obs_ecef[0][0];
|
||||
obs_ecef[1][1] = obs_ecef[0][1];
|
||||
obs_ecef[1][2] = obs_ecef[0][2];
|
||||
obs_ecef[2][0] = obs_ecef[0][0];
|
||||
obs_ecef[2][1] = obs_ecef[0][1];
|
||||
obs_ecef[2][2] = obs_ecef[0][2];
|
||||
|
||||
rc = od_gauss(ra, dec, jd, obs_ecef, &result);
|
||||
|
||||
RUN(rc != 0, "too-close observations rejected");
|
||||
RUN(result.valid == 0, "result marked invalid");
|
||||
}
|
||||
|
||||
|
||||
/* ── Test: radec_to_los / teme_to_radec round-trip ─────── */
|
||||
|
||||
static void
|
||||
test_radec_roundtrip(void)
|
||||
{
|
||||
double ra_in = 1.5; /* ~86 degrees */
|
||||
double dec_in = 0.3; /* ~17 degrees */
|
||||
double los[3];
|
||||
double ra_out, dec_out;
|
||||
double rm;
|
||||
|
||||
fprintf(stderr, "\n--- RA/Dec ↔ LOS round-trip ---\n");
|
||||
|
||||
/* RA/Dec → LOS unit vector */
|
||||
od_radec_to_los(ra_in, dec_in, los);
|
||||
|
||||
rm = sqrt(los[0]*los[0] + los[1]*los[1] + los[2]*los[2]);
|
||||
CLOSE(rm, 1.0, 1e-12, "LOS is unit vector");
|
||||
|
||||
/* LOS → RA/Dec (inverse) */
|
||||
dec_out = asin(los[2]);
|
||||
ra_out = atan2(los[1], los[0]);
|
||||
if (ra_out < 0.0) ra_out += 2.0 * M_PI;
|
||||
|
||||
CLOSE(ra_out, ra_in, 1e-12, "RA round-trip");
|
||||
CLOSE(dec_out, dec_in, 1e-12, "Dec round-trip");
|
||||
}
|
||||
|
||||
|
||||
/* ── Test: teme_to_radec consistency ──────────────────── */
|
||||
|
||||
static void
|
||||
test_teme_to_radec(void)
|
||||
{
|
||||
/* Place a satellite at known TEME position, compute RA/Dec from
|
||||
* a ground observer, verify it's in reasonable range */
|
||||
double pos_teme[3] = {6778.0, 0.0, 0.0}; /* on X-axis, LEO alt */
|
||||
double obs_ecef[3];
|
||||
double ra, dec;
|
||||
double jd = 2451545.0;
|
||||
|
||||
fprintf(stderr, "\n--- teme_to_radec consistency ---\n");
|
||||
|
||||
od_observer_to_ecef(0.0, 0.0, 0.0, obs_ecef); /* equator, prime meridian */
|
||||
|
||||
od_teme_to_radec(pos_teme, obs_ecef, od_gmst_from_jd(jd), &ra, &dec);
|
||||
|
||||
RUN(ra >= 0.0 && ra < 2.0 * M_PI, "RA in [0, 2pi)");
|
||||
RUN(dec >= -M_PI / 2.0 && dec <= M_PI / 2.0, "Dec in [-pi/2, pi/2]");
|
||||
}
|
||||
|
||||
|
||||
int
|
||||
main(void)
|
||||
{
|
||||
fprintf(stderr, "pg_orrery Gauss IOD unit tests\n");
|
||||
fprintf(stderr, "==============================\n");
|
||||
|
||||
test_gauss_iss();
|
||||
test_gauss_meo();
|
||||
test_gauss_too_close();
|
||||
test_radec_roundtrip();
|
||||
test_teme_to_radec();
|
||||
|
||||
fprintf(stderr, "\n%d/%d tests passed.\n", n_pass, n_run);
|
||||
return (n_pass == n_run) ? 0 : 1;
|
||||
}
|
||||
Loading…
x
Reference in New Issue
Block a user