warehack.ing/pg_orrery
2
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Phase 3: Planetary moons and Jupiter radio burst prediction

Browse Source 
Add observation functions for 19 planetary moons across four systems:
- Galilean moons (Io, Europa, Ganymede, Callisto) via clean-room L1.2 theory
- Saturn moons (Mimas through Hyperion) via TASS 1.7
- Uranus moons (Miranda through Oberon) via GUST86
- Mars moons (Phobos, Deimos) via MarsSat

Add Jupiter decametric radio burst prediction for Radio JOVE operators:
- io_phase_angle() — Io orbital phase from superior conjunction
- jupiter_cml() — System III Central Meridian Longitude with light-time correction
- jupiter_burst_probability() — Carr et al. (1983) source regions A, B, C, D

L1.2 Galilean theory is a clean-room MIT implementation from the published
IMCCE FORTRAN coefficients. All other ephemeris libraries are MIT-licensed
extractions from Stellarium with static caching removed for PARALLEL SAFE.

All 10 regression tests pass. Extension .so grows from 2.4MB to 2.5MB.
This commit is contained in:
Ryan Malloy 2026-02-16 01:55:13 -07:00
parent 0544a78276
commit ad7209d0db
15 changed files with 6224 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Makefile
View File

@ -7,7 +7,9 @@ OBJS = src/pg_orbit.o src/tle_type.o src/eci_type.o src/observer_type.o \
src/sgp4_funcs.o src/coord_funcs.o src/pass_funcs.o src/gist_tle.o \ src/sgp4_funcs.o src/coord_funcs.o src/pass_funcs.o src/gist_tle.o \
src/star_funcs.o src/kepler_funcs.o \ src/star_funcs.o src/kepler_funcs.o \
src/vsop87.o src/elp82b.o src/elliptic_to_rectangular.o \ src/vsop87.o src/elp82b.o src/elliptic_to_rectangular.o \
src/precession.o src/sidereal_time.o src/planet_funcs.o src/precession.o src/sidereal_time.o src/planet_funcs.o \
src/tass17.o src/gust86.o src/marssat.o src/l12.o \
src/moon_funcs.o src/radio_funcs.o
# sat_code C++ sources (compiled with g++, linked with extern "C" symbols) # sat_code C++ sources (compiled with g++, linked with extern "C" symbols)
SAT_CODE_DIR = lib/sat_code SAT_CODE_DIR = lib/sat_code
@ -21,7 +23,7 @@ OBJS += $(SAT_CODE_OBJS)
# Regression tests # Regression tests
REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index convenience \ REGRESS = tle_parse sgp4_propagate coord_transforms pass_prediction gist_index convenience \
star_observe kepler_comet planet_observe star_observe kepler_comet planet_observe moon_observe
REGRESS_OPTS = --inputdir=test REGRESS_OPTS = --inputdir=test
# Need C++ runtime for sat_code # Need C++ runtime for sat_code

45
sql/pg_orbit--0.1.0--0.2.0.sql
View File

@ -120,3 +120,48 @@ CREATE FUNCTION moon_observe(observer, timestamptz) RETURNS topocentric
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
COMMENT ON FUNCTION moon_observe(observer, timestamptz) IS 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.'; 'Observe the Moon via ELP2000-82B from (observer, time). Returns topocentric az/el with geocentric range in km.';
-- ============================================================
-- Phase 3: 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.';
-- ============================================================
-- Phase 3: 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.';

45
sql/pg_orbit--0.2.0.sql
View File

@ -636,3 +636,48 @@ CREATE FUNCTION moon_observe(observer, timestamptz) RETURNS topocentric
AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE; AS 'MODULE_PATHNAME' LANGUAGE C IMMUTABLE STRICT PARALLEL SAFE;
COMMENT ON FUNCTION moon_observe(observer, timestamptz) IS 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.'; 'Observe the Moon via ELP2000-82B from (observer, time). Returns topocentric az/el with geocentric range in km.';
-- ============================================================
-- Phase 3: 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.';
-- ============================================================
-- Phase 3: 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.';

440
src/gust86.c Normal file
View File

@ -0,0 +1,440 @@
/************************************************************************
COMPUTATION OF THE COORDINATES OF THE URANIAN SATELLITES (GUST86),
version 0.1 (1988,1995) by LASKAR J. and JACOBSON, R. can be found at
ftp://ftp.imcce.fr/pub/ephem/satel/gust86
I (Johannes Gajdosik) have just taken the Fortran code and data
obtained from above and rearranged it into this piece of software.
I can neither allow nor forbid the usage of the GUST86 theory.
The copyright notice below covers not the works of LASKAR J. and JACOBSON, R.,
but just my work, that is the compilation of the GUST86 theory
into the software supplied in this file.
Copyright (c) 2005 Johannes Gajdosik
Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense,
and/or sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
Derived from Stellarium's GUST86 implementation.
Modified for pg_orbit: removed all static mutable state for thread safety
(PostgreSQL PARALLEL SAFE). The original used static caching arrays and
CalcInterpolatedElements for performance; this version computes fresh on
each call which is acceptable for SQL query workloads.
1) Rotate results to "dynamical equinox and ecliptic J2000",
the reference frame of VSOP87 and VSOP87A.
2) units: julian day, AU, rad
3) use EllipticToRectangularN from elliptic_to_rectangular.h
****************************************************************/
#include "gust86.h"
#include "elliptic_to_rectangular.h"
#include <math.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
static
const double fqn[5] = {4.44519055,
2.492952519,
1.516148111,
0.721718509,
0.46669212};
static
const double fqe[5] = {20.082*M_PI/(180*365.25),
6.217*M_PI/(180*365.25),
2.865*M_PI/(180*365.25),
2.078*M_PI/(180*365.25),
0.386*M_PI/(180*365.25)};
static
const double fqi[5] = {-20.309*M_PI/(180*365.25),
-6.288*M_PI/(180*365.25),
-2.836*M_PI/(180*365.25),
-1.843*M_PI/(180*365.25),
-0.259*M_PI/(180*365.25)};
static
const double phn[5] = {-0.238051,
3.098046,
2.285402,
0.856359,
-0.915592};
static
const double phe[5] = {0.611392,
2.408974,
2.067774,
0.735131,
0.426767};
static
const double phi[5] = {5.702313,
0.395757,
0.589326,
1.746237,
4.206896};
static
void CalcGust86Elem(double t, double elem[30]) {
double an[5], ae[5], ai[5];
int i;
for (i = 0; i < 5; i++) {
an[i] = fmod(fqn[i] * t + phn[i], 2*M_PI);
ae[i] = fmod(fqe[i] * t + phe[i], 2*M_PI);
ai[i] = fmod(fqi[i] * t + phi[i], 2*M_PI);
}
elem[0*6+0] = 4.44352267
- cos(an[0] - an[1] * 3. + an[2] * 2.) * 3.492e-5
+ cos(an[0] * 2. - an[1] * 6. + an[2] * 4.) * 8.47e-6
+ cos(an[0] * 3. - an[1] * 9. + an[2] * 6.) * 1.31e-6
- cos(an[0] - an[1] ) * 5.228e-5
- cos(an[0] * 2. - an[1] * 2. ) * 1.3665e-4;
elem[0*6+1] =
sin(an[0] - an[1] * 3. + an[2] * 2.) * .02547217
- sin(an[0] * 2. - an[1] * 6. + an[2] * 4.) * .00308831
- sin(an[0] * 3. - an[1] * 9. + an[2] * 6.) * 3.181e-4
- sin(an[0] * 4. - an[1] * 12 + an[2] * 8.) * 3.749e-5
- sin(an[0] - an[1] ) * 5.785e-5
- sin(an[0] * 2. - an[1] * 2. ) * 6.232e-5
- sin(an[0] * 3. - an[1] * 3. ) * 2.795e-5
+ t * 4.44519055 - .23805158;
elem[0*6+2] = cos(ae[0]) * .00131238
+ cos(ae[1]) * 7.181e-5
+ cos(ae[2]) * 6.977e-5
+ cos(ae[3]) * 6.75e-6
+ cos(ae[4]) * 6.27e-6
+ cos(an[0]) * 1.941e-4
- cos(-an[0] + an[1] * 2.) * 1.2331e-4
+ cos(an[0] * -2. + an[1] * 3.) * 3.952e-5;
elem[0*6+3] = sin(ae[0]) * .00131238
+ sin(ae[1]) * 7.181e-5
+ sin(ae[2]) * 6.977e-5
+ sin(ae[3]) * 6.75e-6
+ sin(ae[4]) * 6.27e-6
+ sin(an[0]) * 1.941e-4
- sin(-an[0] + an[1] * 2.) * 1.2331e-4
+ sin(an[0] * -2. + an[1] * 3.) * 3.952e-5;
elem[0*6+4] = cos(ai[0]) * .03787171
+ cos(ai[1]) * 2.701e-5
+ cos(ai[2]) * 3.076e-5
+ cos(ai[3]) * 1.218e-5
+ cos(ai[4]) * 5.37e-6;
elem[0*6+5] = sin(ai[0]) * .03787171
+ sin(ai[1]) * 2.701e-5
+ sin(ai[2]) * 3.076e-5
+ sin(ai[3]) * 1.218e-5
+ sin(ai[4]) * 5.37e-6;
elem[1*6+0] = 2.49254257
+ cos(an[0] - an[1] * 3. + an[2] * 2.) * 2.55e-6
- cos( an[1] - an[2] ) * 4.216e-5
- cos( an[1] * 2. - an[2] * 2.) * 1.0256e-4;
elem[1*6+1] =
- sin(an[0] - an[1] * 3. + an[2] * 2.) * .0018605
+ sin(an[0] * 2. - an[1] * 6. + an[2] * 4.) * 2.1999e-4
+ sin(an[0] * 3. - an[1] * 9. + an[2] * 6.) * 2.31e-5
+ sin(an[0] * 4. - an[1] * 12 + an[2] * 8.) * 4.3e-6
- sin( an[1] - an[2] ) * 9.011e-5
- sin( an[1] * 2. - an[2] * 2.) * 9.107e-5
- sin( an[1] * 3. - an[2] * 3.) * 4.275e-5
- sin( an[1] * 2. - an[3] * 2.) * 1.649e-5
+ t * 2.49295252 + 3.09804641;
elem[1*6+2] = cos(ae[0]) * -3.35e-6
+ cos(ae[1]) * .00118763
+ cos(ae[2]) * 8.6159e-4
+ cos(ae[3]) * 7.15e-5
+ cos(ae[4]) * 5.559e-5
- cos(-an[1] + an[2] * 2.) * 8.46e-5
+ cos(an[1] * -2. + an[2] * 3.) * 9.181e-5
+ cos(-an[1] + an[3] * 2.) * 2.003e-5
+ cos(an[1]) * 8.977e-5;
elem[1*6+3] = sin(ae[0]) * -3.35e-6
+ sin(ae[1]) * .00118763
+ sin(ae[2]) * 8.6159e-4
+ sin(ae[3]) * 7.15e-5
+ sin(ae[4]) * 5.559e-5
- sin(-an[1] + an[2] * 2.) * 8.46e-5
+ sin(an[1] * -2. + an[2] * 3.) * 9.181e-5
+ sin(-an[1] + an[3] * 2.) * 2.003e-5
+ sin(an[1]) * 8.977e-5;
elem[1*6+4] = cos(ai[0]) * -1.2175e-4
+ cos(ai[1]) * 3.5825e-4
+ cos(ai[2]) * 2.9008e-4
+ cos(ai[3]) * 9.778e-5
+ cos(ai[4]) * 3.397e-5;
elem[1*6+5] = sin(ai[0]) * -1.2175e-4
+ sin(ai[1]) * 3.5825e-4
+ sin(ai[2]) * 2.9008e-4
+ sin(ai[3]) * 9.778e-5
+ sin(ai[4]) * 3.397e-5;
elem[2*6+0] = 1.5159549
+ cos(an[2] - an[3] * 2. + ae[2]) * 9.74e-6
- cos(an[1] - an[2]) * 1.06e-4
+ cos(an[1] * 2. - an[2] * 2.) * 5.416e-5
- cos(an[2] - an[3]) * 2.359e-5
- cos(an[2] * 2. - an[3] * 2.) * 7.07e-5
- cos(an[2] * 3. - an[3] * 3.) * 3.628e-5;
elem[2*6+1] =
sin(an[0] - an[1] * 3. + an[2] * 2.) * 6.6057e-4
- sin(an[0] * 2. - an[1] * 6. + an[2] * 4.) * 7.651e-5
- sin(an[0] * 3. - an[1] * 9. + an[2] * 6.) * 8.96e-6
- sin(an[0] * 4. - an[1] * 12. + an[2] * 8.) * 2.53e-6
- sin(an[2] - an[3] * 4. + an[4] * 3.) * 5.291e-5
- sin(an[2] - an[3] * 2. + ae[4]) * 7.34e-6
- sin(an[2] - an[3] * 2. + ae[3]) * 1.83e-6
+ sin(an[2] - an[3] * 2. + ae[2]) * 1.4791e-4
+ sin(an[2] - an[3] * 2. + ae[1]) * -7.77e-6
+ sin(an[1] - an[2]) * 9.776e-5
+ sin(an[1] * 2. - an[2] * 2.) * 7.313e-5
+ sin(an[1] * 3. - an[2] * 3.) * 3.471e-5
+ sin(an[1] * 4. - an[2] * 4.) * 1.889e-5
- sin(an[2] - an[3]) * 6.789e-5
- sin(an[2] * 2. - an[3] * 2.) * 8.286e-5
+ sin(an[2] * 3. - an[3] * 3.) * -3.381e-5
- sin(an[2] * 4. - an[3] * 4.) * 1.579e-5
- sin(an[2] - an[4]) * 1.021e-5
- sin(an[2] * 2. - an[4] * 2.) * 1.708e-5
+ t * 1.51614811 + 2.28540169;
elem[2*6+2] = cos(ae[0]) * -2.1e-7
- cos(ae[1]) * 2.2795e-4
+ cos(ae[2]) * .00390469
+ cos(ae[3]) * 3.0917e-4
+ cos(ae[4]) * 2.2192e-4
+ cos(an[1]) * 2.934e-5
+ cos(an[2]) * 2.62e-5
+ cos(-an[1] + an[2] * 2.) * 5.119e-5
- cos(an[1] * -2. + an[2] * 3.) * 1.0386e-4
- cos(an[1] * -3. + an[2] * 4.) * 2.716e-5
+ cos(an[3]) * -1.622e-5
+ cos(-an[2] + an[3] * 2.) * 5.4923e-4
+ cos(an[2] * -2. + an[3] * 3.) * 3.47e-5
+ cos(an[2] * -3. + an[3] * 4.) * 1.281e-5
+ cos(-an[2] + an[4] * 2.) * 2.181e-5
+ cos(an[2]) * 4.625e-5;
elem[2*6+3] = sin(ae[0]) * -2.1e-7
- sin(ae[1]) * 2.2795e-4
+ sin(ae[2]) * .00390469
+ sin(ae[3]) * 3.0917e-4
+ sin(ae[4]) * 2.2192e-4
+ sin(an[1]) * 2.934e-5
+ sin(an[2]) * 2.62e-5
+ sin(-an[1] + an[2] * 2.) * 5.119e-5
- sin(an[1] * -2. + an[2] * 3.) * 1.0386e-4
- sin(an[1] * -3. + an[2] * 4.) * 2.716e-5
+ sin(an[3]) * -1.622e-5
+ sin(-an[2] + an[3] * 2.) * 5.4923e-4
+ sin(an[2] * -2. + an[3] * 3.) * 3.47e-5
+ sin(an[2] * -3. + an[3] * 4.) * 1.281e-5
+ sin(-an[2] + an[4] * 2.) * 2.181e-5
+ sin(an[2]) * 4.625e-5;
elem[2*6+4] = cos(ai[0]) * -1.086e-5
- cos(ai[1]) * 8.151e-5
+ cos(ai[2]) * .00111336
+ cos(ai[3]) * 3.5014e-4
+ cos(ai[4]) * 1.065e-4;
elem[2*6+5] = sin(ai[0]) * -1.086e-5
- sin(ai[1]) * 8.151e-5
+ sin(ai[2]) * .00111336
+ sin(ai[3]) * 3.5014e-4
+ sin(ai[4]) * 1.065e-4;
elem[3*6+0] = .72166316
- cos(an[2] - an[3] * 2. + ae[2]) * 2.64e-6
- cos(an[3] * 2. - an[4] * 3. + ae[4]) * 2.16e-6
+ cos(an[3] * 2. - an[4] * 3. + ae[3]) * 6.45e-6
- cos(an[3] * 2. - an[4] * 3. + ae[2]) * 1.11e-6
+ cos(an[1] - an[3]) * -6.223e-5
- cos(an[2] - an[3]) * 5.613e-5
- cos(an[3] - an[4]) * 3.994e-5
- cos(an[3] * 2. - an[4] * 2.) * 9.185e-5
- cos(an[3] * 3. - an[4] * 3.) * 5.831e-5
- cos(an[3] * 4. - an[4] * 4.) * 3.86e-5
- cos(an[3] * 5. - an[4] * 5.) * 2.618e-5
- cos(an[3] * 6. - an[4] * 6.) * 1.806e-5;
elem[3*6+1] =
sin(an[2] - an[3] * 4. + an[4] * 3.) * 2.061e-5
- sin(an[2] - an[3] * 2. + ae[4]) * 2.07e-6
- sin(an[2] - an[3] * 2. + ae[3]) * 2.88e-6
- sin(an[2] - an[3] * 2. + ae[2]) * 4.079e-5
+ sin(an[2] - an[3] * 2. + ae[1]) * 2.11e-6
- sin(an[3] * 2. - an[4] * 3. + ae[4]) * 5.183e-5
+ sin(an[3] * 2. - an[4] * 3. + ae[3]) * 1.5987e-4
+ sin(an[3] * 2. - an[4] * 3. + ae[2]) * -3.505e-5
- sin(an[3] * 3. - an[4] * 4. + ae[4]) * 1.56e-6
+ sin(an[1] - an[3]) * 4.054e-5
+ sin(an[2] - an[3]) * 4.617e-5
- sin(an[3] - an[4]) * 3.1776e-4
- sin(an[3] * 2. - an[4] * 2.) * 3.0559e-4
- sin(an[3] * 3. - an[4] * 3.) * 1.4836e-4
- sin(an[3] * 4. - an[4] * 4.) * 8.292e-5
+ sin(an[3] * 5. - an[4] * 5.) * -4.998e-5
- sin(an[3] * 6. - an[4] * 6.) * 3.156e-5
- sin(an[3] * 7. - an[4] * 7.) * 2.056e-5
- sin(an[3] * 8. - an[4] * 8.) * 1.369e-5
+ t * .72171851 + .85635879;
elem[3*6+2] = cos(ae[0]) * -2e-8
- cos(ae[1]) * 1.29e-6
- cos(ae[2]) * 3.2451e-4
+ cos(ae[3]) * 9.3281e-4
+ cos(ae[4]) * .00112089
+ cos(an[1]) * 3.386e-5
+ cos(an[3]) * 1.746e-5
+ cos(-an[1] + an[3] * 2.) * 1.658e-5
+ cos(an[2]) * 2.889e-5
- cos(-an[2] + an[3] * 2.) * 3.586e-5
+ cos(an[3]) * -1.786e-5
- cos(an[4]) * 3.21e-5
- cos(-an[3] + an[4] * 2.) * 1.7783e-4
+ cos(an[3] * -2. + an[4] * 3.) * 7.9343e-4
+ cos(an[3] * -3. + an[4] * 4.) * 9.948e-5
+ cos(an[3] * -4. + an[4] * 5.) * 4.483e-5
+ cos(an[3] * -5. + an[4] * 6.) * 2.513e-5
+ cos(an[3] * -6. + an[4] * 7.) * 1.543e-5;
elem[3*6+3] = sin(ae[0]) * -2e-8
- sin(ae[1]) * 1.29e-6
- sin(ae[2]) * 3.2451e-4
+ sin(ae[3]) * 9.3281e-4
+ sin(ae[4]) * .00112089
+ sin(an[1]) * 3.386e-5
+ sin(an[3]) * 1.746e-5
+ sin(-an[1] + an[3] * 2.) * 1.658e-5
+ sin(an[2]) * 2.889e-5
- sin(-an[2] + an[3] * 2.) * 3.586e-5
+ sin(an[3]) * -1.786e-5
- sin(an[4]) * 3.21e-5
- sin(-an[3] + an[4] * 2.) * 1.7783e-4
+ sin(an[3] * -2. + an[4] * 3.) * 7.9343e-4
+ sin(an[3] * -3. + an[4] * 4.) * 9.948e-5
+ sin(an[3] * -4. + an[4] * 5.) * 4.483e-5
+ sin(an[3] * -5. + an[4] * 6.) * 2.513e-5
+ sin(an[3] * -6. + an[4] * 7.) * 1.543e-5;
elem[3*6+4] = cos(ai[0]) * -1.43e-6
- cos(ai[1]) * 1.06e-6
- cos(ai[2]) * 1.4013e-4
+ cos(ai[3]) * 6.8572e-4
+ cos(ai[4]) * 3.7832e-4;
elem[3*6+5] = sin(ai[0]) * -1.43e-6
- sin(ai[1]) * 1.06e-6
- sin(ai[2]) * 1.4013e-4
+ sin(ai[3]) * 6.8572e-4
+ sin(ai[4]) * 3.7832e-4;
elem[4*6+0] = .46658054
+ cos(an[3] * 2. - an[4] * 3. + ae[4]) * 2.08e-6
- cos(an[3] * 2. - an[4] * 3. + ae[3]) * 6.22e-6
+ cos(an[3] * 2. - an[4] * 3. + ae[2]) * 1.07e-6
- cos(an[1] - an[4]) * 4.31e-5
+ cos(an[2] - an[4]) * -3.894e-5
- cos(an[3] - an[4]) * 8.011e-5
+ cos(an[3] * 2. - an[4] * 2.) * 5.906e-5
+ cos(an[3] * 3. - an[4] * 3.) * 3.749e-5
+ cos(an[3] * 4. - an[4] * 4.) * 2.482e-5
+ cos(an[3] * 5. - an[4] * 5.) * 1.684e-5;
elem[4*6+1] =
- sin(an[2] - an[3] * 4. + an[4] * 3.) * 7.82e-6
+ sin(an[3] * 2. - an[4] * 3. + ae[4]) * 5.129e-5
- sin(an[3] * 2. - an[4] * 3. + ae[3]) * 1.5824e-4
+ sin(an[3] * 2. - an[4] * 3. + ae[2]) * 3.451e-5
+ sin(an[1] - an[4]) * 4.751e-5
+ sin(an[2] - an[4]) * 3.896e-5
+ sin(an[3] - an[4]) * 3.5973e-4
+ sin(an[3] * 2. - an[4] * 2.) * 2.8278e-4
+ sin(an[3] * 3. - an[4] * 3.) * 1.386e-4
+ sin(an[3] * 4. - an[4] * 4.) * 7.803e-5
+ sin(an[3] * 5. - an[4] * 5.) * 4.729e-5
+ sin(an[3] * 6. - an[4] * 6.) * 3e-5
+ sin(an[3] * 7. - an[4] * 7.) * 1.962e-5
+ sin(an[3] * 8. - an[4] * 8.) * 1.311e-5
+ t * .46669212 - .9155918;
elem[4*6+2] = cos(ae[1]) * -3.5e-7
+ cos(ae[2]) * 7.453e-5
- cos(ae[3]) * 7.5868e-4
+ cos(ae[4]) * .00139734
+ cos(an[1]) * 3.9e-5
+ cos(-an[1] + an[4] * 2.) * 1.766e-5
+ cos(an[2]) * 3.242e-5
+ cos(an[3]) * 7.975e-5
+ cos(an[4]) * 7.566e-5
+ cos(-an[3] + an[4] * 2.) * 1.3404e-4
- cos(an[3] * -2. + an[4] * 3.) * 9.8726e-4
- cos(an[3] * -3. + an[4] * 4.) * 1.2609e-4
- cos(an[3] * -4. + an[4] * 5.) * 5.742e-5
- cos(an[3] * -5. + an[4] * 6.) * 3.241e-5
- cos(an[3] * -6. + an[4] * 7.) * 1.999e-5
- cos(an[3] * -7. + an[4] * 8.) * 1.294e-5;
elem[4*6+3] = sin(ae[1]) * -3.5e-7
+ sin(ae[2]) * 7.453e-5
- sin(ae[3]) * 7.5868e-4
+ sin(ae[4]) * .00139734
+ sin(an[1]) * 3.9e-5
+ sin(-an[1] + an[4] * 2.) * 1.766e-5
+ sin(an[2]) * 3.242e-5
+ sin(an[3]) * 7.975e-5
+ sin(an[4]) * 7.566e-5
+ sin(-an[3] + an[4] * 2.) * 1.3404e-4
- sin(an[3] * -2. + an[4] * 3.) * 9.8726e-4
- sin(an[3] * -3. + an[4] * 4.) * 1.2609e-4
- sin(an[3] * -4. + an[4] * 5.) * 5.742e-5
- sin(an[3] * -5. + an[4] * 6.) * 3.241e-5
- sin(an[3] * -6. + an[4] * 7.) * 1.999e-5
- sin(an[3] * -7. + an[4] * 8.) * 1.294e-5;
elem[4*6+4] = cos(ai[0]) * -4.4e-7
- cos(ai[1]) * 3.1e-7
+ cos(ai[2]) * 3.689e-5
- cos(ai[3]) * 5.9633e-4
+ cos(ai[4]) * 4.5169e-4;
elem[4*6+5] = sin(ai[0]) * -4.4e-7
- sin(ai[1]) * 3.1e-7
+ sin(ai[2]) * 3.689e-5
- sin(ai[3]) * 5.9633e-4
+ sin(ai[4]) * 4.5169e-4;
}
static
const double gust86_rmu[5] = {
1.291892353675174e-08,
1.291910570526396e-08,
1.291910102284198e-08,
1.291942656265575e-08,
1.291935967091320e-08};
static const double GUST86toVsop87[9] = {
9.753206632086812015e-01, 6.194425668001473004e-02, 2.119257251551559653e-01,
-2.006444610981783542e-01,-1.519328516640849367e-01, 9.678110398294910731e-01,
9.214881523275189928e-02,-9.864478281437795399e-01,-1.357544776485127136e-01
};
void GetGust86Coor(double jd, int body, double *xyz, double *xyzdot) {
double elem[5*6];
double x[6];
const double t = jd - 2444239.5;
CalcGust86Elem(t, elem);
EllipticToRectangularN(gust86_rmu[body], elem + (body * 6), 0.0, x);
xyz[0] = GUST86toVsop87[0]*x[0] + GUST86toVsop87[1]*x[1] + GUST86toVsop87[2]*x[2];
xyz[1] = GUST86toVsop87[3]*x[0] + GUST86toVsop87[4]*x[1] + GUST86toVsop87[5]*x[2];
xyz[2] = GUST86toVsop87[6]*x[0] + GUST86toVsop87[7]*x[1] + GUST86toVsop87[8]*x[2];
if (xyzdot) {
xyzdot[0] = GUST86toVsop87[0]*x[3] + GUST86toVsop87[1]*x[4] + GUST86toVsop87[2]*x[5];
xyzdot[1] = GUST86toVsop87[3]*x[3] + GUST86toVsop87[4]*x[4] + GUST86toVsop87[5]*x[5];
xyzdot[2] = GUST86toVsop87[6]*x[3] + GUST86toVsop87[7]*x[4] + GUST86toVsop87[8]*x[5];
}
}

78
src/gust86.h Normal file
View File

@ -0,0 +1,78 @@
/************************************************************************
COMPUTATION OF THE COORDINATES OF THE URANIAN SATELLITES (GUST86),
version 0.1 (1988,1995) by LASKAR J. and JACOBSON, R. can be found at
ftp://ftp.imcce.fr/pub/ephem/satel/gust86
I (Johannes Gajdosik) have just taken the Fortran code and data
obtained from above and rearranged it into this piece of software.
I can neither allow nor forbid the usage of the GUST86 theory.
The copyright notice below covers not the works of LASKAR J. and JACOBSON, R.,
but just my work, that is the compilation of the GUST86 theory
into the software supplied in this file.
Copyright (c) 2005 Johannes Gajdosik
Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense,
and/or sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
Derived from Stellarium's GUST86 implementation.
Modified for pg_orbit: removed static mutable state and
CalcInterpolatedElements caching for thread safety
(PostgreSQL PARALLEL SAFE). Elements are computed fresh on each call.
1) Rotate results to "dynamical equinox and ecliptic J2000",
the reference frame of VSOP87 and VSOP87A.
2) units: julian day, AU, rad
****************************************************************/
#ifndef PG_ORBIT_GUST86_H
#define PG_ORBIT_GUST86_H
#ifdef __cplusplus
extern "C" {
#endif
#define GUST86_MIRANDA 0
#define GUST86_ARIEL 1
#define GUST86_UMBRIEL 2
#define GUST86_TITANIA 3
#define GUST86_OBERON 4
void GetGust86Coor(double jd, int body, double *xyz, double *xyzdot);
/* Return the rectangular coordinates of the given satellite
and the given julian date jd expressed in dynamical time (TAI+32.184s).
The origin of the xyz-coordinates is the center of Uranus.
The reference frame is "dynamical equinox and ecliptic J2000",
which is the reference frame in VSOP87 and VSOP87A.
body: 0=Miranda, 1=Ariel, 2=Umbriel, 3=Titania, 4=Oberon
xyz[3]: position (AU), relative to Uranus center
xyzdot[3]: velocity (AU/day), may be NULL
*/
#ifdef __cplusplus
}
#endif
#endif

918
src/l12.c Normal file
View File

@ -0,0 +1,918 @@
/************************************************************************
L1.2 Galilean satellite theory -- Lainey, Duriez & Vienne
Clean-room implementation for pg_orbit.
Positions and velocities of Io, Europa, Ganymede, and Callisto
relative to Jupiter's center, in VSOP87 ecliptic J2000 coordinates.
Reference:
Lainey V., Duriez L., Vienne A.
"New accurate ephemerides for the Galilean satellites of Jupiter"
Astronomy & Astrophysics, 2004
ftp://ftp.imcce.fr/pub/ephem/satel/galilean/L1/L1.2/
The theory expresses each moon's orbit using modified Delaunay
variables {a, L, K, H, Q, P}, where:
a = semi-major axis (AU)
L = mean longitude (radians)
K = e * cos(varpi) (eccentricity vector, real part)
H = e * sin(varpi) (eccentricity vector, imaginary part)
Q = sin(i/2) * cos(Omega) (inclination vector, real part)
P = sin(i/2) * sin(Omega) (inclination vector, imaginary part)
Each variable is computed as a Fourier series in time since the
theory epoch (JD 2433282.5 = 1950 Jan 1.0 TT).
Copyright (c) 2026 Ryan Malloy <ryan@supported.systems>
Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense,
and/or sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
Thread-safe: all functions are reentrant with no static mutable state.
****************************************************************/
#include "l12.h"
#include <math.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846264338327950288
#endif
#define TWO_PI (2.0 * M_PI)
/* L1.2 theory epoch: JD 2433282.5 (1950 Jan 1.0 TT) */
#define L12_EPOCH_JD 2433282.5
/* Kepler equation convergence threshold */
#define KEPLER_TOL 1.0e-12
/* Maximum series lengths across all four moons */
#define MAX_TERMS_A 38
#define MAX_TERMS_L 36
#define MAX_TERMS_Z 50
#define MAX_TERMS_ZETA 25
/*
* A single Fourier term: amplitude * trig(phase + frequency * t)
*
* For semi-major axis 'a': amplitude * cos(phase + freq * t)
* For mean longitude 'L': amplitude * sin(phase + freq * t)
* For complex eccentricity 'z': amplitude * cos/sin(phase + freq * t)
* For complex inclination 'zeta': amplitude * cos/sin(phase + freq * t)
*/
typedef struct {
double amp;
double phi;
double nu;
} l12_fourier_term;
/*
* All theory data for one Galilean moon.
*/
typedef struct {
double grav_param; /* mu: GM_Jupiter in AU^3/day^2 */
double lon0; /* mean longitude at epoch (rad) */
double lon_rate; /* mean longitude rate (rad/day) */
int n_a; /* number of semi-major axis terms */
int n_l; /* number of mean longitude terms */
int n_z; /* number of eccentricity terms */
int n_zeta; /* number of inclination terms */
l12_fourier_term fa[MAX_TERMS_A];
l12_fourier_term fl[MAX_TERMS_L];
l12_fourier_term fz[MAX_TERMS_Z];
l12_fourier_term fzeta[MAX_TERMS_ZETA];
} l12_moon_data;
/*
* Rotation matrix: L1.2 reference frame -> VSOP87 ecliptic J2000
*
* This is the product of the equatorial-to-ecliptic rotation
* evaluated at J2000, representing the physical orientation of
* Jupiter's Laplacian plane relative to the ecliptic. These nine
* values are an astronomical coordinate transform (physical constant).
*
* Stored row-major: row i, col j = rot_l12_to_vsop87[3*i + j]
*/
static const double rot_l12_to_vsop87[9] = {
9.994327815023905713e-01,
3.039550993390781261e-02,
-1.449924943755843383e-02,
-3.089770442223671880e-02,
9.988822846893227815e-01,
-3.577028369016394015e-02,
1.339578739122566807e-02,
3.619798764705610479e-02,
9.992548516622136737e-01
};
/*
* Convert modified Delaunay elements to position and velocity.
*
* Elements: {a, L, K, H, Q, P} where
* a = semi-major axis
* L = mean longitude
* K = e*cos(varpi), H = e*sin(varpi)
* Q = sin(i/2)*cos(Omega), P = sin(i/2)*sin(Omega)
*
* The mean anomaly in modified Delaunay variables is simply L
* (since the perturbation series already accounts for the
* longitude of periapse through K and H).
*
* Kepler's equation in K,H form:
* E = L + K*sin(E) - H*cos(E)
* where E is the eccentric longitude (not anomaly).
*/
static void
delaunay_to_cartesian(double grav_param, const double elems[6],
double pos[3], double vel[3])
{
double sma, mean_lon, kk, hh, qq, pp;
double mean_mot, ecc_lon, cos_e, sin_e;
double delta, dle, rsm1, inv_r_over_a;
double phi_ecc, psi_fac;
double xp, yp, vxp, vyp;
double f2, one_m_2pp, one_m_2qq, two_pq;
int iter;
sma = elems[0];
mean_lon = elems[1];
kk = elems[2];
hh = elems[3];
qq = elems[4];
pp = elems[5];
/* mean motion from Kepler's third law */
mean_mot = sqrt(grav_param / (sma * sma * sma));
/* solve generalized Kepler equation:
* E = L + K*sin(E) - H*cos(E)
* using Newton-Raphson iteration */
ecc_lon = mean_lon + kk * sin(mean_lon) - hh * cos(mean_lon);
for (iter = 0; iter < 20; iter++) {
cos_e = cos(ecc_lon);
sin_e = sin(ecc_lon);
delta = (mean_lon - ecc_lon + kk * sin_e - hh * cos_e)
/ (1.0 - kk * cos_e - hh * sin_e);
ecc_lon += delta;
if (fabs(delta) <= KEPLER_TOL)
break;
}
cos_e = cos(ecc_lon);
sin_e = sin(ecc_lon);
/* auxiliary quantities */
dle = hh * cos_e - kk * sin_e;
rsm1 = -kk * cos_e - hh * sin_e; /* (r/a - 1) */
inv_r_over_a = 1.0 / (1.0 + rsm1); /* a/r */
/* eccentricity-related factor */
phi_ecc = sqrt(1.0 - kk * kk - hh * hh);
psi_fac = 1.0 / (1.0 + phi_ecc);
/* position in orbital plane */
xp = sma * (cos_e - kk - psi_fac * hh * dle);
yp = sma * (sin_e - hh + psi_fac * kk * dle);
/* velocity in orbital plane */
vxp = mean_mot * inv_r_over_a * sma * (-sin_e - psi_fac * hh * rsm1);
vyp = mean_mot * inv_r_over_a * sma * ( cos_e + psi_fac * kk * rsm1);
/* rotate from orbital plane to 3D using Q, P (inclination vector) */
f2 = 2.0 * sqrt(1.0 - qq * qq - pp * pp);
one_m_2pp = 1.0 - 2.0 * pp * pp;
one_m_2qq = 1.0 - 2.0 * qq * qq;
two_pq = 2.0 * pp * qq;
pos[0] = xp * one_m_2pp + yp * two_pq;
pos[1] = xp * two_pq + yp * one_m_2qq;
pos[2] = (qq * yp - xp * pp) * f2;
vel[0] = vxp * one_m_2pp + vyp * two_pq;
vel[1] = vxp * two_pq + vyp * one_m_2qq;
vel[2] = (qq * vyp - vxp * pp) * f2;
}
/*
* Apply 3x3 rotation matrix to a 3-vector.
*/
static void
rotate_vec(const double mat[9], const double in[3], double out[3])
{
out[0] = mat[0] * in[0] + mat[1] * in[1] + mat[2] * in[2];
out[1] = mat[3] * in[0] + mat[4] * in[1] + mat[5] * in[2];
out[2] = mat[6] * in[0] + mat[7] * in[1] + mat[8] * in[2];
}
/* forward declaration -- data defined below */
static const l12_moon_data galilean_moons[4];
/*
* Evaluate the L1.2 Fourier series and compute Cartesian
* position/velocity for one Galilean moon.
*/
void
GetL12Coor(double jd, int body, double *xyz, double *xyzdot)
{
const l12_moon_data *md;
double dt, angle, accum, re, im;
double elems[6];
double body_pos[3], body_vel[3];
int j;
/* select moon data (bounds check) */
if (body < 0 || body > 3)
return;
md = &galilean_moons[body];
/* time since L1.2 epoch in days */
dt = jd - L12_EPOCH_JD;
/*
* 1. Semi-major axis: sum of cosine terms
* a = sum_j amp_j * cos(phi_j + nu_j * dt)
*/
accum = 0.0;
for (j = 0; j < md->n_a; j++) {
angle = md->fa[j].phi + md->fa[j].nu * dt;
accum += md->fa[j].amp * cos(angle);
}
elems[0] = accum;
/*
* 2. Mean longitude: linear trend + sine series
* L = lon0 + lon_rate * dt + sum_j amp_j * sin(phi_j + nu_j * dt)
*/
accum = md->lon0 + md->lon_rate * dt;
for (j = 0; j < md->n_l; j++) {
angle = md->fl[j].phi + md->fl[j].nu * dt;
accum += md->fl[j].amp * sin(angle);
}
accum = fmod(accum, TWO_PI);
if (accum < 0.0)
accum += TWO_PI;
elems[1] = accum;
/*
* 3. Complex eccentricity z = K + iH = sum_j amp_j * exp(i*(phi_j + nu_j*dt))
* K = sum amp * cos(angle), H = sum amp * sin(angle)
*/
re = 0.0;
im = 0.0;
for (j = 0; j < md->n_z; j++) {
angle = md->fz[j].phi + md->fz[j].nu * dt;
re += md->fz[j].amp * cos(angle);
im += md->fz[j].amp * sin(angle);
}
elems[2] = re;
elems[3] = im;
/*
* 4. Complex inclination zeta = Q + iP = sum_j amp_j * exp(i*(phi_j + nu_j*dt))
* Q = sum amp * cos(angle), P = sum amp * sin(angle)
*/
re = 0.0;
im = 0.0;
for (j = 0; j < md->n_zeta; j++) {
angle = md->fzeta[j].phi + md->fzeta[j].nu * dt;
re += md->fzeta[j].amp * cos(angle);
im += md->fzeta[j].amp * sin(angle);
}
elems[4] = re;
elems[5] = im;
/* convert elements to position/velocity in L1.2 frame */
delaunay_to_cartesian(md->grav_param, elems, body_pos, body_vel);
/* rotate from L1.2 frame to VSOP87 ecliptic J2000 */
rotate_vec(rot_l12_to_vsop87, body_pos, xyz);
if (xyzdot)
rotate_vec(rot_l12_to_vsop87, body_vel, xyzdot);
}
/* ================================================================
* Theory coefficients for the four Galilean moons.
*
* These are astronomical constants derived from fitting
* observations of the Galilean satellite system. They represent
* physical measurements of orbital parameters and are scientific
* data, not copyrightable expression.
*
* Source: L1.2 theory data files
* ftp://ftp.imcce.fr/pub/ephem/satel/galilean/L1/L1.2/
* ================================================================
*/
static const l12_moon_data galilean_moons[4] = {
/* ---- Io (J-I) ---- */
{
/* grav_param */ 0.2824894284338140e-06,
/* lon0 */ 0.1446213296021224e+01,
/* lon_rate */ 0.3551552286182400e+01,
/* n_a */ 38, /* n_l */ 32, /* n_z */ 23, /* n_zeta */ 15,
/* semi-major axis cosine series (38 terms) */
{
{ 0.0028210960212903, 0.00000000000000e+00, 0.00000000000000e+00},
{ 0.0000000762024588, 0.36392902322306e+01, 0.35644591656241e+01},
{ 0.0000000180900324, 0.99554707056522e+00, 0.71289183312483e+01},
{ 0.0000000172337652, 0.18196487820921e+01, 0.17822295777568e+01},
{ 0.0000000101726080, 0.28150559763861e+01, 0.89111478635073e+01},
{ 0.0000000094794086, 0.34760224933239e+01, 0.80200331112799e+01},
{ 0.0000000092196266, 0.46347004953370e+01, 0.10693377436209e+02},
{ 0.0000000058581604, 0.11586746335276e+01, 0.26733443704266e+01},
{-0.0000000036218148, 0.23173675289588e+01, 0.53466887181044e+01},
{ 0.0000000034892754, 0.17122470613669e+00, 0.12475607079684e+02},
{ 0.0000000030842852, 0.36170311370435e+01, 0.63501320826717e+01},
{ 0.0000000020794650, 0.19906655633153e+01, 0.14257836656755e+02},
{ 0.0000000013655244, 0.49369712857369e+01, 0.13584836518140e-01},
{ 0.0000000011682572, 0.57934065580556e+01, 0.13366721796637e+02},
{-0.0000000008031976, 0.66879731833041e+00, 0.16040066232595e+02},
{ 0.0000000007309510, 0.56300556878949e+01, 0.17822295806244e+02},
{ 0.0000000007014118, 0.43297377080515e+01, 0.71002044886497e+01},
{ 0.0000000006561624, 0.43188797534991e+01, 0.13034138433510e-01},
{ 0.0000000005753088, 0.54252179509841e+01, 0.95251981240076e+01},
{ 0.0000000004359548, 0.11670110887440e+01, 0.19604525331797e+02},
{ 0.0000000003711992, 0.14936154077537e+01, 0.12938928912340e-01},
{-0.0000000003412576, 0.26346374300664e+01, 0.15571117257960e-01},
{ 0.0000000003432980, 0.17994723387341e+01, 0.31750663461810e+01},
{ 0.0000000003228344, 0.29861854159944e+01, 0.21386754933987e+02},
{ 0.0000000003014418, 0.19871924348983e+00, 0.24675315510310e-01},
{ 0.0000000001707670, 0.50718778620273e+01, 0.35514255456604e+01},
{ 0.0000000001655832, 0.29783205832994e+01, 0.44555739317536e+01},
{ 0.0000000001612910, 0.48058392680935e+01, 0.23168984521460e+02},
{ 0.0000000001527992, 0.18275651107267e+01, 0.18713410599600e+02},
{ 0.0000000001523312, 0.46323297275220e+01, 0.44686092108182e+01},
{ 0.0000000001449720, 0.19079860214667e+01, 0.30506511533200e-02},
{ 0.0000000001188688, 0.53321680658912e+01, 0.70987549449082e+01},
{ 0.0000000001129258, 0.95031497804420e+00, 0.12700264165343e+02},
{ 0.0000000000986086, 0.34190944178580e+00, 0.24951214111224e+02},
{-0.0000000000877720, 0.36228267942948e+01, 0.17958145576535e+01},
{ 0.0000000000857194, 0.33682834215727e+01, 0.59736711266730e+01},
{-0.0000000000545492, 0.19473964103154e+01, 0.22929425718040e-01},
{ 0.0000000000326102, 0.24880420823571e+01, 0.25119610718870e-01}
},
/* mean longitude sine series (32 terms) */
{
{-0.0001925258348666, 0.49369589722645e+01, 0.13584836583050e-01},
{-0.0000970803596076, 0.43188796477322e+01, 0.13034138432430e-01},
{-0.0000898817416500, 0.19080016428617e+01, 0.30506486715800e-02},
{-0.0000553101050262, 0.14936156681569e+01, 0.12938928911550e-01},
{-0.0000503584426150, 0.36410196089987e+01, 0.35644591049605e+01},
{-0.0000444412770116, 0.18196478828985e+01, 0.17822295777568e+01},
{ 0.0000418078870490, 0.26346334480977e+01, 0.15571117221300e-01},
{ 0.0000372356597388, 0.21402440902650e+01, 0.14500977488900e-02},
{-0.0000234440533016, 0.19871945729267e+00, 0.24675315507400e-01},
{-0.0000160313164240, 0.28203470990931e+01, 0.95196190000000e-04},
{-0.0000119049755698, 0.99521552502799e+00, 0.71289183312483e+01},
{-0.0000109014269320, 0.11586742711973e+01, 0.26733443704266e+01},
{ 0.0000087217118104, 0.22995085327344e+01, 0.44456805185000e-03},
{ 0.0000082229455492, 0.84723690387904e+00, 0.54980078903000e-03},
{ 0.0000075365481720, 0.30644603245150e+01, 0.64826749624000e-03},
{-0.0000061452803962, 0.28150499448772e+01, 0.89111478635073e+01},
{-0.0000057575824778, 0.34760236756099e+01, 0.80200331112799e+01},
{-0.0000053196302672, 0.14952058549171e+01, 0.29001290992300e-02},
{-0.0000051181206936, 0.46347077042449e+01, 0.10693377436209e+02},
{-0.0000047817413326, 0.49236512419835e+01, 0.30554833877800e-02},
{ 0.0000045554015322, 0.19585097634352e+01, 0.22928625941210e-01},
{ 0.0000043204134698, 0.15842888614383e+01, 0.15677112478190e-01},
{ 0.0000037684098282, 0.23173652780077e+01, 0.53466887181044e+01},
{-0.0000031403738248, 0.22184076281042e+01, 0.25155489165510e-01},
{ 0.0000024336535428, 0.85320650238886e+00, 0.25426968834400e-02},
{-0.0000020289901692, 0.36168998565188e+01, 0.63501320826717e+01},
{ 0.0000018665438704, 0.48458061649481e+01, 0.13589674898130e-01},
{-0.0000018552431038, 0.17086811529922e+00, 0.12475607079684e+02},
{-0.0000016229875536, 0.62803206871082e+01, 0.60414171604000e-03},
{-0.0000013160987604, 0.14718125754925e+01, 0.14358460012320e-01},
{ 0.0000008070729808, 0.38735416148641e+00, 0.37658680379100e-02},
{ 0.0000002602397658, 0.14337589305551e+01, 0.45692429208000e-02}
},
/* complex eccentricity cos/sin series (23 terms) */
{
{ 0.0041510849668155, 0.40899396355450e+01, -0.12906864146660e-01},
{ 0.0006260521444113, 0.14461888986270e+01, 0.35515522949802e+01},
{ 0.0000352747346169, 0.21256287034578e+01, 0.12727416567000e-03},
{ 0.0000198194483636, 0.55835619926762e+01, 0.32065751140000e-04},
{ 0.0000146399842989, 0.44137212696837e+00, 0.26642533547700e-02},
{ 0.0000098749504021, 0.45076118781320e+00, -0.35773660260022e+01},
{-0.0000096819265753, 0.59097266550442e+01, 0.17693227079462e+01},
{-0.0000083063168209, 0.28751474873012e+00, 0.87820791951527e+00},
{ 0.0000059689735869, 0.50740752477871e+01, 0.71160118048918e+01},
{-0.0000052220588690, 0.27460731023666e+01, 0.67796100936000e-03},
{ 0.0000046538995236, 0.49143203339385e+01, -0.53595956184347e+01},
{ 0.0000045951340101, 0.42533513770304e+01, -0.44684808200578e+01},
{-0.0000037711061757, 0.54120093562773e+01, -0.17951364445825e+01},
{ 0.0000037405126681, 0.30946737347297e+01, -0.71418251640916e+01},
{ 0.0000022044764663, 0.54360702580001e+01, -0.26491700241240e-01},
{ 0.0000018698303790, 0.41124042914226e+01, -0.27985797954589e+01},
{-0.0000015410375360, 0.27141931505529e+01, 0.27731236679900e-02},
{ 0.0000013214613496, 0.12750177723530e+01, -0.89240547799787e+01},
{-0.0000012707585609, 0.51141075152507e+01, 0.37654227378982e+00},
{ 0.0000012193607962, 0.59977053365953e+01, -0.98566169956900e-02},
{-0.0000011886104747, 0.32658350285168e+01, 0.53337818460633e+01},
{ 0.0000008742035177, 0.23903528311144e+01, 0.25194921818800e-02},
{-0.0000007689215742, 0.38308837306225e+01, -0.27293225500100e-02}
},
/* complex inclination cos/sin series (15 terms) */
{
{ 0.0003142172466014, 0.27964219722923e+01, -0.23150960980000e-02},
{ 0.0000904169207946, 0.10477061879627e+01, -0.56920638196000e-03},
{ 0.0000175695395780, 0.24150809680215e+01, 0.00000000000000e+00},
{ 0.0000164452324013, 0.33368861773902e+01, -0.12491307197000e-03},
{ 0.0000055424829091, 0.59720202381027e+01, -0.30561164720000e-04},
{ 0.0000035856270353, 0.84898736841329e+00, -0.25244521900630e-01},
{ 0.0000024180760140, 0.55603770950923e+01, 0.29003681445800e-02},
{-0.0000008673084930, 0.28496686106299e+00, -0.14500593353200e-02},
{-0.0000003176227277, 0.53834633036029e+01, -0.23498632298700e-01},
{ 0.0000003152816608, 0.45569499027478e+01, 0.43504654304000e-02},
{ 0.0000002338676726, 0.17633292120047e+01, 0.14501339138600e-02},
{ 0.0000001754553689, 0.48429319984493e+01, -0.25688816532440e-01},
{ 0.0000001286319583, 0.57543347143871e+01, -0.25813660979740e-01},
{ 0.0000000967213304, 0.11503592426900e+01, -0.29001471397800e-02},
{ 0.0000000000692310, 0.40745966852008e+01, -0.32506757319070e-01}
}
},
/* ---- Europa (J-II) ---- */
{
/* grav_param */ 0.2824832743928930e-06,
/* lon0 */ -0.3735263437471362e+00,
/* lon_rate */ 0.1769322711123470e+01,
/* n_a */ 38, /* n_l */ 36, /* n_z */ 41, /* n_zeta */ 25,
/* semi-major axis cosine series (38 terms) */
{
{ 0.0044871037804314, 0.00000000000000e+00, 0.00000000000000e+00},
{ 0.0000004324367498, 0.18196456062910e+01, 0.17822295777568e+01},
{ 0.0000001603614750, 0.43002726529577e+01, 0.26733443704266e+01},
{-0.0000001019146786, 0.54589480865442e+01, 0.53466887181044e+01},
{ 0.0000000924734786, 0.56222139048906e+01, 0.89111478887838e+00},
{-0.0000000523665800, 0.36392846323417e+01, 0.35644591656241e+01},
{ 0.0000000511509000, 0.29783307371014e+01, 0.44555739317536e+01},
{-0.0000000311907780, 0.99466557754027e+00, 0.71289183312483e+01},
{-0.0000000272859938, 0.28144480309092e+01, 0.89111478635073e+01},
{ 0.0000000232225828, 0.62608434364366e+01, 0.27856729211550e+01},
{-0.0000000181310770, 0.43188692380649e+01, 0.13034138308860e-01},
{ 0.0000000174960544, 0.16563941638726e+01, 0.62378035398422e+01},
{-0.0000000122874072, 0.46421290370833e+01, 0.10693377254218e+02},
{-0.0000000095367130, 0.14936536615312e+01, 0.12938928820690e-01},
{-0.0000000084863836, 0.17146854643555e+00, 0.12475607079684e+02},
{ 0.0000000071939342, 0.49376739095661e+01, 0.13584833017030e-01},
{ 0.0000000069122354, 0.62488746138492e+01, 0.41785094280464e+01},
{ 0.0000000061377568, 0.33434976298081e+00, 0.80200331112799e+01},
{-0.0000000045343054, 0.19892156959655e+01, 0.14257836656755e+02},
{ 0.0000000044574684, 0.34597804303324e+00, 0.35357692227539e+01},
{ 0.0000000042350072, 0.62719655202169e+01, 0.13928364636651e+01},
{-0.0000000028783772, 0.38108811302610e+01, 0.16040066232595e+02},
{ 0.0000000024354662, 0.99587190880214e+00, 0.90414989297380e+00},
{ 0.0000000022532940, 0.52958965893939e+01, 0.98022627054664e+01},
{ 0.0000000021573570, 0.62379050559630e+01, 0.55713458670107e+01},
{-0.0000000016530062, 0.56456686036734e+01, 0.17822294694543e+02},
{ 0.0000000016464798, 0.26346435392424e+01, 0.15571117126760e-01},
{ 0.0000000011589838, 0.32732388195745e+01, 0.17691951440716e+01},
{-0.0000000010251826, 0.19079858535660e+01, 0.30506497838200e-02},
{-0.0000000010203510, 0.11692020351116e+01, 0.19604525194172e+02},
{ 0.0000000007614982, 0.16862812414995e+01, 0.35342961443230e+01},
{ 0.0000000007104494, 0.59112717191092e+01, 0.24092214574831e+01},
{-0.0000000006957184, 0.24879412197796e+01, 0.25119609730670e-01},
{-0.0000000005817914, 0.19872303312324e+00, 0.24675315511270e-01},
{-0.0000000003792178, 0.15765189821595e+01, 0.25244441830200e-01},
{ 0.0000000003397378, 0.58126953372535e+01, 0.25973067138760e-01},
{ 0.0000000003159492, 0.23545476741301e+01, 0.26068277099550e-01},
{ 0.0000000002538154, 0.19471441186087e+01, 0.22929424919760e-01}
},
/* mean longitude sine series (36 terms) */
{
{ 0.0008576433172936, 0.43188693178264e+01, 0.13034138308050e-01},
{ 0.0004549582875086, 0.14936531751079e+01, 0.12938928819620e-01},
{ 0.0003248939825174, 0.18196494533458e+01, 0.17822295777568e+01},
{-0.0003074250079334, 0.49377037005911e+01, 0.13584832867240e-01},
{ 0.0001982386144784, 0.19079869054760e+01, 0.30510121286900e-02},
{ 0.0001834063551804, 0.21402853388529e+01, 0.14500978933800e-02},
{-0.0001434383188452, 0.56222140366630e+01, 0.89111478887838e+00},
{-0.0000771939140944, 0.43002724372350e+01, 0.26733443704266e+01},
{-0.0000632289777196, 0.26346392822098e+01, 0.15571117084700e-01},
{ 0.0000446766477388, 0.54589448561143e+01, 0.53466887181044e+01},
{ 0.0000436574731410, 0.36392908617709e+01, 0.35644591656241e+01},
{ 0.0000349172750296, 0.28289867162553e+01, 0.29885749150000e-04},
{-0.0000325709094646, 0.53721409780230e+01, 0.12495233774000e-03},
{ 0.0000205826473860, 0.15258464215508e+01, 0.29001315522200e-02},
{-0.0000192706087556, 0.29783311531879e+01, 0.44555739317536e+01},
{ 0.0000168028316254, 0.24879414119403e+01, 0.25119609725650e-01},
{-0.0000141628733606, 0.29183576504413e+01, 0.64930403718000e-03},
{ 0.0000140713155600, 0.19872319369353e+00, 0.24675315510310e-01},
{ 0.0000131946915760, 0.99584744364935e+00, 0.71289183312483e+01},
{ 0.0000106598617620, 0.53356907396678e+01, 0.30233219231900e-02},
{-0.0000104011727738, 0.62608296198866e+01, 0.27856729211550e+01},
{ 0.0000100746080234, 0.44288900030073e+01, 0.55297871931000e-03},
{ 0.0000097414019416, 0.27312462188296e+01, 0.89111510230745e+01},
{-0.0000094651366640, 0.25010358163865e+01, 0.93478322470000e-04},
{ 0.0000091108073324, 0.15765182522628e+01, 0.25244441682120e-01},
{-0.0000087720567668, 0.15376962386886e+01, 0.15676393315070e-01},
{-0.0000078429703340, 0.58128473756772e+01, 0.25973069246350e-01},
{-0.0000075566039418, 0.30586251688920e+01, 0.43252872559000e-03},
{-0.0000066580990752, 0.19591270593390e+01, 0.22928567412490e-01},
{-0.0000065854142774, 0.18617673337640e+01, 0.26093058384670e-01},
{-0.0000058131135230, 0.16563893807978e+01, 0.62378035398422e+01},
{ 0.0000055720865276, 0.39565695752204e+01, 0.25481216339900e-02},
{-0.0000048198508906, 0.62720230965345e+01, 0.13928364605775e+01},
{ 0.0000042728431266, 0.46220912946918e+01, 0.10693377982182e+02},
{ 0.0000042175545304, 0.13509343368359e+01, 0.30164435787800e-02},
{ 0.0000037707624520, 0.51034507119889e+01, 0.25219658202250e-01}
},
/* complex eccentricity cos/sin series (41 terms) */
{
{-0.0093589104136341, 0.40899396509039e+01, -0.12906864146660e-01},
{ 0.0002988994545555, 0.59097265185595e+01, 0.17693227079462e+01},
{ 0.0002139036390350, 0.21256289300016e+01, 0.12727418407000e-03},
{ 0.0001980963564781, 0.27435168292650e+01, 0.67797343009000e-03},
{ 0.0001210388158965, 0.55839943711203e+01, 0.32056614900000e-04},
{ 0.0000837042048393, 0.16094538368039e+01, -0.90402165808846e+00},
{ 0.0000823525166369, 0.14461887708689e+01, 0.35515522949802e+01},
{-0.0000315906532820, 0.28751224400811e+00, 0.87820791951527e+00},
{-0.0000294503681314, 0.45078002968967e+00, -0.35773660260022e+01},
{-0.0000278946698536, 0.22704374310903e+01, -0.17951364497113e+01},
{ 0.0000144958688621, 0.29313956641719e+01, -0.26862512422390e+01},
{ 0.0000139052321679, 0.60542576187622e+01, -0.25941002404500e-01},
{ 0.0000108374431350, 0.59320761116863e+01, -0.10163502160128e+01},
{-0.0000082175838585, 0.49144730088838e+01, -0.53595956184347e+01},
{ 0.0000073925894084, 0.25962855881215e+01, -0.25845792927870e-01},
{ 0.0000062618381566, 0.62252936384007e+01, -0.71418248393794e+01},
{-0.0000051968296512, 0.54353355159239e+01, -0.26491696725040e-01},
{-0.0000043507065743, 0.51150292346242e+01, 0.37654221060604e+00},
{ 0.0000042081682285, 0.31202613836361e+01, 0.44427230757158e+01},
{ 0.0000041298266970, 0.42533371370636e+01, -0.44684808200578e+01},
{-0.0000036991221930, 0.52487564172390e+01, 0.26604375002057e+01},
{-0.0000027357551003, 0.12734806685602e+01, -0.89240546532297e+01},
{-0.0000026854901206, 0.75596663258784e+00, 0.15806953180460e-01},
{ 0.0000023074479953, 0.19438998534712e+01, 0.71160114825227e+01},
{ 0.0000020163445050, 0.58484195254467e+01, -0.24091843401454e+01},
{-0.0000018506530067, 0.26838225102582e+01, 0.68236609075000e-03},
{ 0.0000018159137522, 0.26048690461733e+01, 0.62248966818788e+01},
{-0.0000017894118824, 0.57385537790777e+01, -0.10706284328884e+02},
{ 0.0000016518864520, 0.32658492478888e+01, 0.53337818460633e+01},
{-0.0000015660692561, 0.61789350505156e+01, -0.11453588847960e-01},
{ 0.0000014426949422, 0.60014075911383e+01, -0.17664769149811e+01},
{ 0.0000013196935928, 0.55753025652974e+01, -0.62507103665413e+01},
{-0.0000011726743714, 0.50242932747650e+01, -0.14351210704140e-01},
{-0.0000009550285338, 0.28409403047363e+01, 0.14257595044900e-02},
{-0.0000007569857746, 0.38098760906143e+01, -0.27271627793800e-02},
{-0.0000007495662748, 0.29896372346394e+01, 0.20097553243900e-02},
{ 0.0000007091149133, 0.27139331814919e+01, -0.19924932783300e-02},
{ 0.0000005646670312, 0.21683602575236e+01, -0.12871803232940e-01},
{-0.0000002004455524, 0.12893849410519e+01, -0.32724923477800e-02},
{-0.0000001623489363, 0.24189454629613e+00, 0.44609337678800e-02},
{ 0.0000001058862562, 0.45356953407129e+01, 0.39269908172100e-02}
},
/* complex inclination cos/sin series (25 terms) */
{
{ 0.0040404917832303, 0.10477063169425e+01, -0.56920640540000e-03},
{ 0.0002200421034564, 0.33368857864364e+01, -0.12491307307000e-03},
{ 0.0001662544744719, 0.24134862374711e+01, 0.00000000000000e+00},
{ 0.0000590282470983, 0.59719930968366e+01, -0.30561602250000e-04},
{-0.0000105030331400, 0.27964978379152e+01, -0.23150966123800e-02},
{-0.0000102943248250, 0.84898796322150e+00, -0.25244521901650e-01},
{ 0.0000072600013020, 0.55603730312676e+01, 0.29003676713100e-02},
{ 0.0000018391258758, 0.28480515491153e+00, -0.14500579196900e-02},
{ 0.0000014880605763, 0.48429974929766e+01, -0.25688815138710e-01},
{-0.0000008828196274, 0.65011185407635e+00, 0.34696170683100e-02},
{ 0.0000008714042768, 0.17639430319108e+01, 0.14501352157600e-02},
{ 0.0000008536188044, 0.45568506666427e+01, 0.43504641410100e-02},
{ 0.0000006846214331, 0.57542117253981e+01, -0.25813768702630e-01},
{ 0.0000004471826348, 0.53834694321520e+01, -0.23498632366370e-01},
{ 0.0000003034392168, 0.22078201315180e+01, -0.25783170906020e-01},
{ 0.0000001799083735, 0.31858868501531e+01, 0.88086056517000e-03},
{-0.0000001792048645, 0.51949494917342e+01, -0.20193236931900e-02},
{-0.0000001098546626, 0.59286821904995e+01, 0.49197316579700e-02},
{-0.0000001083128732, 0.45808061408794e+01, -0.59959459406000e-03},
{ 0.0000001062153531, 0.38387102863271e+01, -0.53795085847000e-03},
{ 0.0000000768496749, 0.35553768729770e+01, 0.58005587812700e-02},
{-0.0000000692273841, 0.46440611341931e+01, 0.30253219029200e-02},
{ 0.0000000676969224, 0.13621319661456e+00, -0.44430413602000e-03},
{-0.0000000621559952, 0.30093497179950e+01, -0.13603287200690e-01},
{ 0.0000000000608298, 0.40529569532600e+01, -0.32510869900940e-01}
}
},
/* ---- Ganymede (J-III) ---- */
{
/* grav_param */ 0.2824981841847230e-06,
/* lon0 */ 0.2874089391143348e+00,
/* lon_rate */ 0.8782079235893280e+00,
/* n_a */ 38, /* n_l */ 31, /* n_z */ 50, /* n_zeta */ 18,
/* semi-major axis cosine series (38 terms) */
{
{ 0.0071566594572575, 0.00000000000000e+00, 0.00000000000000e+00},
{ 0.0000013930299110, 0.11586745884981e+01, 0.26733443704266e+01},
{ 0.0000006449829346, 0.56222145702102e+01, 0.89111478887838e+00},
{ 0.0000002298059520, 0.12995924044108e+01, 0.10034433456729e+01},
{-0.0000001221434370, 0.49612436330515e+01, 0.17822295777568e+01},
{ 0.0000001095798176, 0.19486708778350e+01, 0.15051650461529e+01},
{ 0.0000000701435616, 0.64978508114196e+00, 0.50172166963138e+00},
{ 0.0000000547868566, 0.25992050672074e+01, 0.20068866945508e+01},
{-0.0000000394635858, 0.23173535605652e+01, 0.53466887181044e+01},
{-0.0000000363221428, 0.36393008632056e+01, 0.35644591656241e+01},
{ 0.0000000290949364, 0.20123392230605e+01, 0.17535157350384e+01},
{ 0.0000000281244968, 0.32490010762048e+01, 0.25086083721948e+01},
{-0.0000000207924698, 0.29783308899923e+01, 0.44555739317536e+01},
{ 0.0000000146896774, 0.38988244013504e+01, 0.30103300418262e+01},
{-0.0000000119930042, 0.16563968316083e+01, 0.62378035398422e+01},
{ 0.0000000112067460, 0.43188665692819e+01, 0.13034138285340e-01},
{-0.0000000109535132, 0.49372826282154e+01, 0.13584834937940e-01},
{ 0.0000000099867772, 0.96700720263958e+00, 0.62699633776977e+00},
{ 0.0000000077668260, 0.45486373016444e+01, 0.35120517182683e+01},
{ 0.0000000074143972, 0.16140449852661e+00, 0.12531361661566e+00},
{ 0.0000000066346638, 0.33441073536010e+00, 0.80200331112799e+01},
{ 0.0000000057842684, 0.14936630646671e+01, 0.12938928799370e-01},
{-0.0000000055768352, 0.44651777597613e+01, 0.11287386144710e+01},
{-0.0000000049395106, 0.61563894598809e+01, 0.17520662093793e+01},
{ 0.0000000041439704, 0.51984558307998e+01, 0.40137734147421e+01},
{-0.0000000040765630, 0.99543742426922e+00, 0.71289183312483e+01},
{-0.0000000036862062, 0.46386836178626e+01, 0.10693377254218e+02},
{ 0.0000000033617538, 0.37493658441448e+01, 0.87808669180168e+00},
{ 0.0000000033348284, 0.22668196818990e+01, 0.16304394485673e+01},
{-0.0000000025754698, 0.33293196902303e+00, 0.17952648120307e+01},
{ 0.0000000024363084, 0.19604838407749e+01, 0.11232854513197e+00},
{ 0.0000000022265432, 0.58482745704418e+01, 0.45154950951905e+01},
{ 0.0000000020032676, 0.29166648062069e+01, 0.21321610765333e+01},
{-0.0000000018115706, 0.99782757414001e+00, 0.90414978368384e+00},
{ 0.0000000014535006, 0.18748212041600e+01, 0.89112137093506e+01},
{-0.0000000006819260, 0.19871670124324e+00, 0.24675315493830e-01},
{ 0.0000000004433776, 0.24880003003965e+01, 0.25119610196650e-01},
{-0.0000000002836658, 0.58126277034761e+01, 0.25973068607520e-01}
},
/* mean longitude sine series (31 terms) */
{
{ 0.0002310797886226, 0.21402987195942e+01, 0.14500978438400e-02},
{-0.0001828635964118, 0.43188672736968e+01, 0.13034138282630e-01},
{ 0.0001512378778204, 0.49373102372298e+01, 0.13584834812520e-01},
{-0.0001163720969778, 0.43002659861490e+01, 0.26733443704266e+01},
{-0.0000955478069846, 0.14936612842567e+01, 0.12938928798570e-01},
{ 0.0000815246854464, 0.56222137132535e+01, 0.89111478887838e+00},
{-0.0000801219679602, 0.12995922951532e+01, 0.10034433456729e+01},
{-0.0000607017260182, 0.64978769669238e+00, 0.50172167043264e+00},
{ 0.0000543922473002, 0.27927547440639e+01, 0.29880873700000e-04},
{-0.0000489253646474, 0.53711728089803e+01, 0.12495278292000e-03},
{-0.0000427574981536, 0.18196513407448e+01, 0.17822295777568e+01},
{-0.0000307360417826, 0.19498372703786e+01, 0.15051650064903e+01},
{-0.0000169767346458, 0.19078637281659e+01, 0.30507678226700e-02},
{ 0.0000154725890508, 0.56912713028984e+01, 0.65164073556000e-03},
{-0.0000145268863648, 0.18863875475387e+00, 0.12530827181195e+00},
{-0.0000135654458738, 0.27930238268852e+01, 0.55663681407000e-03},
{-0.0000134648621904, 0.25991972928128e+01, 0.20068866945508e+01},
{ 0.0000095524017320, 0.23173520454449e+01, 0.53466887181044e+01},
{ 0.0000087955125170, 0.36393024031096e+01, 0.35644591656241e+01},
{ 0.0000075462003630, 0.53560617584395e+01, 0.92426977490000e-04},
{-0.0000071146195958, 0.20120561622463e+01, 0.17535157644008e+01},
{ 0.0000064153141218, 0.15526366820734e+01, 0.29001309732400e-02},
{-0.0000063221625942, 0.32490122452649e+01, 0.25086083721948e+01},
{-0.0000056564973024, 0.24862139082596e+01, 0.44834622386000e-03},
{ 0.0000052570245720, 0.19871532348033e+00, 0.24675315501580e-01},
{ 0.0000047020767994, 0.29783317790630e+01, 0.44555739317536e+01},
{-0.0000047004229470, 0.96617050453708e+00, 0.62699712737505e+00},
{-0.0000046565198820, 0.36125113449716e+01, 0.43633231340000e-03},
{-0.0000042349322008, 0.19604744669606e+01, 0.11232854282257e+00},
{-0.0000038755741918, 0.22619624763183e+01, 0.25146663939730e-01},
{-0.0000032577733688, 0.56861827246039e+01, 0.17074576501600e-02}
},
/* complex eccentricity cos/sin series (50 terms) */
{
{ 0.0014289811307319, 0.21256295942739e+01, 0.12727413029000e-03},
{ 0.0007710931226760, 0.55836330003496e+01, 0.32064341100000e-04},
{ 0.0005925911780766, 0.40899396636448e+01, -0.12906864146660e-01},
{ 0.0002045597496146, 0.52713683670372e+01, -0.12523544076106e+00},
{ 0.0001785118648258, 0.28743156721063e+00, 0.87820792442520e+00},
{ 0.0001131999784893, 0.14462127277818e+01, 0.35515522949802e+01},
{-0.0000658778169210, 0.22702423990985e+01, -0.17951364394537e+01},
{ 0.0000497058888328, 0.59096792204858e+01, 0.17693227129285e+01},
{-0.0000316384926978, 0.16093054939404e+01, -0.90402165028424e+00},
{ 0.0000287801237327, 0.46217321268757e+01, -0.62695712341840e+00},
{-0.0000181744317896, 0.59210641379360e+01, 0.37648623991673e+00},
{ 0.0000105558175161, 0.39720191398746e+01, -0.11286788041058e+01},
{-0.0000070808673396, 0.60542548894164e+01, -0.25941002415210e-01},
{-0.0000070804404020, 0.27978433776854e+01, 0.67774258703000e-03},
{-0.0000061046181888, 0.14151685760988e+01, -0.87530769416913e+00},
{-0.0000057610853129, 0.42530537622646e+01, -0.44684807882788e+01},
{-0.0000057310334964, 0.29311803223072e+01, -0.26862512192699e+01},
{ 0.0000048299146941, 0.27138294508149e+01, 0.27731329671900e-02},
{ 0.0000046610005483, 0.33222980229554e+01, -0.16304004832039e+01},
{-0.0000038142769361, 0.25962943627643e+01, -0.25845792955510e-01},
{ 0.0000034982417330, 0.15866568011217e+01, 0.18816512920593e+01},
{-0.0000030091617315, 0.35921173988567e+01, -0.35773660056343e+01},
{-0.0000024732926446, 0.53461730094807e+01, 0.25122576835111e+00},
{ 0.0000024416432533, 0.47049477027963e+01, -0.25049613834712e+00},
{ 0.0000024171568015, 0.34508032389167e+01, 0.00000000000000e+00},
{ 0.0000023143850535, 0.55385759257808e+01, 0.28683339028800e-02},
{ 0.0000022651772374, 0.55608006706192e+01, 0.14501892967800e-02},
{ 0.0000022247695560, 0.26725424635341e+01, -0.21321221654766e+01},
{ 0.0000020947921969, 0.22350374116258e+01, 0.23833730192673e+01},
{-0.0000014042712722, 0.93718044411525e+00, 0.13799296041822e+01},
{ 0.0000011932531874, 0.28861941414418e+01, 0.28850946416201e+01},
{-0.0000011180389240, 0.49139919849718e+01, -0.53595955727170e+01},
{ 0.0000011076384510, 0.20227538540345e+01, -0.26338438454332e+01},
{-0.0000010371714944, 0.40722739402948e+00, -0.87385759089274e+00},
{-0.0000008993128501, 0.30942691883530e+01, -0.71418251640916e+01},
{ 0.0000007268381420, 0.54334774230433e+01, -0.26491687896550e-01},
{-0.0000007178049665, 0.52487423493616e+01, 0.26604375002057e+01},
{ 0.0000006908412319, 0.40596134184175e+01, -0.75221793556997e+00},
{-0.0000006784151570, 0.38846818226669e+01, 0.42496535534400e-02},
{ 0.0000006772314920, 0.23013479896873e+01, 0.26317235358158e+01},
{ 0.0000006659820028, 0.35359530295550e+01, 0.33868163258510e+01},
{-0.0000006339665249, 0.39268665697903e+01, 0.44426670658559e+01},
{-0.0000006286307601, 0.19440608894162e+01, 0.71160114019304e+01},
{-0.0000006128705113, 0.25027415074658e+01, 0.62249001971556e+01},
{ 0.0000005660807396, 0.13729316457251e+01, -0.31355655165873e+01},
{-0.0000005206551413, 0.55749300982469e+01, -0.62507103665413e+01},
{-0.0000004718481418, 0.45366605084874e+01, 0.16786677353000e-03},
{-0.0000004583970422, 0.19351070248496e+01, -0.98151695574500e+01},
{-0.0000004577854173, 0.62350780976534e+01, 0.17563373058434e+01},
{ 0.0000003466029660, 0.75412427489767e+00, 0.15807097495700e-01}
},
/* complex inclination cos/sin series (18 terms) */
{
{ 0.0015932721570848, 0.33368862796665e+01, -0.12491307058000e-03},
{ 0.0008533093128905, 0.24133881688166e+01, 0.00000000000000e+00},
{ 0.0003513347911037, 0.59720789850127e+01, -0.30561017710000e-04},
{-0.0001441929255483, 0.10477061764435e+01, -0.56920632124000e-03},
{ 0.0000157303527750, 0.55604041197704e+01, 0.29003665011200e-02},
{ 0.0000025161319881, 0.28477653709685e+00, -0.14500554486800e-02},
{ 0.0000020438305183, 0.17652628559998e+01, 0.14501383926500e-02},
{ 0.0000017939612784, 0.45568977341583e+01, 0.43504621590400e-02},
{ 0.0000013614276895, 0.84898872627945e+00, -0.25244521900630e-01},
{-0.0000008996109017, 0.46441156003340e+01, 0.30253214588300e-02},
{-0.0000008702078430, 0.27972000093551e+01, -0.23150965645100e-02},
{-0.0000004371144064, 0.48429530385679e+01, -0.25688816011500e-01},
{-0.0000002174259374, 0.57543785603741e+01, -0.25813642993310e-01},
{-0.0000001926397869, 0.20118539705648e+01, 0.29330596864500e-02},
{ 0.0000001589279656, 0.35554727018503e+01, 0.58005577768400e-02},
{-0.0000001432228753, 0.11966574544002e+01, -0.15750124983800e-02},
{-0.0000000926213408, 0.22052538606469e+01, -0.25782797426020e-01},
{ 0.0000000000106902, 0.45764213311755e+01, -0.32611614716800e-01}
}
},
/* ---- Callisto (J-IV) ---- */
{
/* grav_param */ 0.2824921448899090e-06,
/* lon0 */ -0.3620341291375704e+00,
/* lon_rate */ 0.3764862334338280e+00,
/* n_a */ 22, /* n_l */ 19, /* n_z */ 46, /* n_zeta */ 18,
/* semi-major axis cosine series (22 terms) */
{
{ 0.0125879701715314, 0.00000000000000e+00, 0.00000000000000e+00},
{ 0.0000035952049470, 0.64965776007116e+00, 0.50172168165034e+00},
{ 0.0000027580210652, 0.18084235781510e+01, 0.31750660413359e+01},
{ 0.0000012874896172, 0.62718908285025e+01, 0.13928364698403e+01},
{-0.0000004173729106, 0.12990650292663e+01, 0.10034433697108e+01},
{ 0.0000002790757718, 0.71428870045577e+00, 0.75007225869130e+00},
{-0.0000001998252258, 0.19489881012004e+01, 0.15051650461529e+01},
{-0.0000001001149838, 0.25987168731338e+01, 0.20068867266014e+01},
{-0.0000000513967092, 0.32484798706247e+01, 0.25086084022422e+01},
{-0.0000000475687992, 0.48635521917696e+01, 0.74862216593606e+00},
{ 0.0000000348242240, 0.15082713497295e+00, 0.37645917070525e+00},
{ 0.0000000283840630, 0.51672973364888e+01, 0.12530678073049e+00},
{-0.0000000263234638, 0.33499822210495e+01, 0.30103491232578e+01},
{ 0.0000000239106346, 0.43573519442736e+01, 0.62698238798737e+00},
{ 0.0000000219977422, 0.15075404808879e+01, 0.27986109086768e+01},
{-0.0000000171144478, 0.62607361864777e+01, 0.27856729335053e+01},
{-0.0000000141956834, 0.45481077718910e+01, 0.35120517575303e+01},
{-0.0000000120003630, 0.18583887479127e+01, 0.11287042579152e+01},
{ 0.0000000108418904, 0.54873138800427e+01, 0.67395238593127e+01},
{ 0.0000000108218254, 0.59772630516669e+01, 0.10163811590412e+01},
{ 0.0000000002477642, 0.56894071957878e+01, 0.65165021654000e-03},
{-0.0000000001874576, 0.28598333265121e+01, 0.55639542661000e-03}
},
/* mean longitude sine series (19 terms) */
{
{ 0.0005586040123824, 0.21404207189815e+01, 0.14500979323100e-02},
{-0.0003805813868176, 0.27358844897853e+01, 0.29729650620000e-04},
{ 0.0002205152863262, 0.64979652596400e+00, 0.50172167243580e+00},
{ 0.0001877895151158, 0.18084787604005e+01, 0.31750660413359e+01},
{ 0.0000766916975242, 0.62720114319755e+01, 0.13928364636651e+01},
{ 0.0000747056855106, 0.12995916202344e+01, 0.10034433456729e+01},
{-0.0000388323297366, 0.71289234751879e+00, 0.75007236972328e+00},
{ 0.0000335036484314, 0.53712641184981e+01, 0.12494011725000e-03},
{ 0.0000293032677938, 0.19493939340593e+01, 0.15051650209131e+01},
{ 0.0000185940935472, 0.14630998372377e+01, 0.29001339405200e-02},
{-0.0000170438022886, 0.56893382353856e+01, 0.65165044781000e-03},
{ 0.0000151393833114, 0.28749516044614e+01, 0.55646069067000e-03},
{-0.0000148825637256, 0.33321074618840e+01, 0.12530790075011e+00},
{ 0.0000129927896682, 0.25991973549465e+01, 0.20068866945508e+01},
{-0.0000116117398772, 0.56192268627131e+01, 0.93166256720000e-04},
{ 0.0000066211702894, 0.48564958193206e+01, 0.74862286166569e+00},
{ 0.0000065387442486, 0.35580120361824e+01, 0.16550513741900e-02},
{ 0.0000061580798140, 0.32490037889701e+01, 0.25086083721948e+01},
{ 0.0000046797140778, 0.96612169919707e+00, 0.62699716616712e+00}
},
/* complex eccentricity cos/sin series (46 terms) */
{
{ 0.0073755808467977, 0.55836071576084e+01, 0.32065099140000e-04},
{ 0.0002065924169942, 0.59209831565786e+01, 0.37648624194703e+00},
{ 0.0001589869764021, 0.28744006242623e+00, 0.87820792442520e+00},
{-0.0001561131605348, 0.21257397865089e+01, 0.12727441285000e-03},
{ 0.0001486043380971, 0.14462134301023e+01, 0.35515522949802e+01},
{ 0.0000635073108731, 0.59096803285954e+01, 0.17693227129285e+01},
{ 0.0000599351698525, 0.41125517584798e+01, -0.27985797954589e+01},
{ 0.0000540660842731, 0.55390350845569e+01, 0.28683408228300e-02},
{-0.0000489596900866, 0.46218149483338e+01, -0.62695712529519e+00},
{ 0.0000333682283528, 0.52066975238880e+01, -0.37358601734497e+00},
{ 0.0000295832427279, 0.59322697896516e+01, -0.10163502275209e+01},
{ 0.0000292325461337, 0.52707623402008e+01, -0.12523542448602e+00},
{ 0.0000197588369441, 0.33317768022759e+01, 0.00000000000000e+00},
{-0.0000183551029746, 0.39720443249757e+01, -0.11286788041058e+01},
{ 0.0000090411191759, 0.55606719963947e+01, 0.14501837490800e-02},
{-0.0000081987970452, 0.33223313720086e+01, -0.16304004832039e+01},
{-0.0000060406575087, 0.13970265485562e+01, 0.43191832032300e-02},
{ 0.0000056895636122, 0.41990956668120e+01, -0.37213592656720e+00},
{-0.0000040434854859, 0.47008406172134e+01, -0.25049602889288e+00},
{-0.0000039403527376, 0.26725832255243e+01, -0.21321221654766e+01},
{ 0.0000036901291978, 0.35207772267753e+00, 0.11265585018525e+01},
{-0.0000028551622596, 0.55601265129356e+01, -0.31584886140000e-04},
{-0.0000026588026505, 0.25969882784477e+00, 0.14182025553300e-02},
{-0.0000019711212463, 0.20228019680496e+01, -0.26338438454332e+01},
{ 0.0000019322089806, 0.51418595457408e+01, 0.25123117908847e+00},
{-0.0000018673159813, 0.93674892088247e+00, 0.13799296163047e+01},
{ 0.0000016838424078, 0.60796033426941e+01, 0.75294520843775e+00},
{-0.0000016695689644, 0.15867810488422e+01, 0.18816512864243e+01},
{ 0.0000016317841395, 0.45789534393209e+01, 0.14822429153000e-02},
{-0.0000016159095087, 0.30157253757329e+00, -0.14180284447900e-02},
{-0.0000014034621874, 0.59433512039442e+01, -0.24091866865037e+01},
{-0.0000012029942283, 0.27137754880270e+01, 0.27731373092900e-02},
{-0.0000011758260607, 0.40581098970285e+01, -0.75221789504525e+00},
{-0.0000010798624964, 0.22364861319452e+01, 0.23833729650229e+01},
{-0.0000010108880552, 0.13729872033949e+01, -0.31355655165873e+01},
{-0.0000008876681807, 0.50534107615010e+01, -0.50169281575682e+00},
{ 0.0000008869382117, 0.50147420853991e+01, 0.68353864231000e-03},
{-0.0000008194699011, 0.62190878357566e+01, -0.37503615934219e+00},
{ 0.0000007093782158, 0.44118312641559e+01, -0.24221246131672e+01},
{ 0.0000006728737059, 0.31910016062920e+01, -0.37068584881726e+00},
{ 0.0000006297345982, 0.13595719733984e+01, 0.11251084135564e+01},
{ 0.0000006128899757, 0.51402161299290e+01, 0.71160095483087e+01},
{-0.0000005580987049, 0.34117733109010e+01, -0.12539396666771e+01},
{ 0.0000005321318002, 0.35377046967957e+01, 0.57685340271300e-02},
{-0.0000004739086661, 0.21645217929478e+01, 0.58845474482000e-03},
{ 0.0000004518928658, 0.44963664372727e+01, 0.29023325111200e-02}
},
/* complex inclination cos/sin series (18 terms) */
{
{ 0.0038422977898495, 0.24133922085557e+01, 0.00000000000000e+00},
{ 0.0022453891791894, 0.59721736773277e+01, -0.30561255250000e-04},
{-0.0002604479450559, 0.33368746306409e+01, -0.12491309972000e-03},
{ 0.0000332112143230, 0.55604137742337e+01, 0.29003768850700e-02},
{ 0.0000049727136261, 0.28488229706820e+00, -0.14500571761900e-02},
{-0.0000049416729114, 0.10476908456459e+01, -0.56920298857000e-03},
{ 0.0000043945193428, 0.17684273746003e+01, 0.14501344524700e-02},
{ 0.0000037630501589, 0.45567680530533e+01, 0.43504645407000e-02},
{-0.0000030823418750, 0.20094360655956e+01, 0.29313051376700e-02},
{ 0.0000004719790711, 0.18055417618741e+01, 0.14195445432000e-02},
{-0.0000004637177865, 0.38277528822158e+01, -0.14808731001600e-02},
{ 0.0000003497224175, 0.46444360330108e+01, 0.30253130162300e-02},
{-0.0000003467132626, 0.10120757927163e+01, 0.43816126822900e-02},
{ 0.0000003324412570, 0.35549391686606e+01, 0.58005379032100e-02},
{ 0.0000001945374351, 0.61251687150860e+01, 0.28808264872800e-02},
{ 0.0000001727743329, 0.11900773236610e+01, -0.29001068524700e-02},
{-0.0000001485176585, 0.62335834706368e+01, 0.14807679092700e-02},
{ 0.0000000000666922, 0.40616225761771e+01, -0.32724923474890e-01}
}
}
};

73
src/l12.h Normal file
View File

@ -0,0 +1,73 @@
/************************************************************************
L1.2 Galilean satellite theory -- Lainey, Duriez & Vienne
Clean-room implementation for pg_orbit.
The L1.2 theory provides positions and velocities of Jupiter's four
Galilean moons (Io, Europa, Ganymede, Callisto) relative to Jupiter's
center, in the VSOP87 ecliptic J2000 reference frame.
Reference:
Lainey V., Duriez L., Vienne A.
"New accurate ephemerides for the Galilean satellites of Jupiter"
Astronomy & Astrophysics, 2004
ftp://ftp.imcce.fr/pub/ephem/satel/galilean/L1/L1.2/
The theory coefficients are astronomical constants derived from
fitting observations of the Galilean system. They represent
physical measurements and are not copyrightable.
Copyright (c) 2026 Ryan Malloy <ryan@supported.systems>
Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense,
and/or sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
Thread-safe: all functions are reentrant with no static mutable state.
****************************************************************/
#ifndef L12_H
#define L12_H
#ifdef __cplusplus
extern "C" {
#endif
#define L12_IO 0
#define L12_EUROPA 1
#define L12_GANYMEDE 2
#define L12_CALLISTO 3
/*
* Compute the position (and optionally velocity) of a Galilean moon
* using the L1.2 semi-analytic theory.
*
* jd: Julian date in Terrestrial Time (TT).
* body: Moon index: 0=Io, 1=Europa, 2=Ganymede, 3=Callisto.
* xyz: Output position in AU, VSOP87 ecliptic J2000 frame,
* relative to Jupiter's center. Must point to double[3].
* xyzdot: Output velocity in AU/day, same frame. May be NULL
* to skip velocity computation.
*/
void GetL12Coor(double jd, int body, double *xyz, double *xyzdot);
#ifdef __cplusplus
}
#endif
#endif /* L12_H */

425
src/marssat.c Normal file
View File

@ -0,0 +1,425 @@
/************************************************************************
The Ephemerides of the Martian satellites
(adjustement from 1877 to 2005, Version 1.0)
by Valery Lainey can be obtained from Valery Lainey:
V.Lainey (Lainey@oma.be)
ROB- 3, Avenue Circulaire, B-1180 Bruxelles (Belgium)
IMCCE - 77, Avenue Denfert-Rochereau 75014 Paris (France)
-----------------------------------------------------------------------
I (Johannes Gajdosik) have just taken Valery Laineys Fortran code,
MarsSatV1-0.f, which he kindly supplied, and rearranged it into
this piece of software.
I can neither allow nor forbid the usage of Valery Laineys
Ephemerides of the Martian satellites.
The copyright notice below covers not the work of Valery Lainey
but just my work, that is the compilation of Valery Laineys
Ephemerides of the Martian satellites into the software supplied in this file.
Copyright (c) 2006 Johannes Gajdosik
Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense,
and/or sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
Derived from Stellarium's MARSSAT implementation.
Modified for pg_orbit: removed all static mutable state for thread safety
(PostgreSQL PARALLEL SAFE). The original used static caching arrays and
CalcInterpolatedElements for performance; this version computes fresh on
each call which is acceptable for SQL query workloads.
1) do not calculate constant terms at runtime but beforehand
2) unite terms with the same frequencies
****************************************************************/
#include "marssat.h"
#include "elliptic_to_rectangular.h"
#include <math.h>
#include <string.h>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
struct MarsSatTerm {
double phase, frequency, amplitude;
};
struct MarsSatTermList {
const struct MarsSatTerm *terms;
int size;
};
struct MarsSatBody {
double mu, l, acc;
double constants[6];
const struct MarsSatTermList lists[4];
};
static const struct MarsSatTerm mars_sat_phobos_0[16] = {
{ 5.013490350126586e+00, 2.715567382195733e+01, 4.537539999999999e-09},
{ 6.839910590780000e-01, 1.969446151585829e+01, 7.312000000000000e-10},
{ 4.514031245592586e+00, 4.073350897245015e+01, 7.785599999999993e-10},
{ 5.531351678889586e+00, 1.357783661756435e+01, 3.529399999999985e-10},
{ 5.700307292822586e+00, 4.685013393001369e+01, 2.094599999999994e-10},
{ 4.228062220664587e+00, 5.431134764391466e+01, 9.139999999999999e-11},
{ 1.167907583017500e+00, 7.461211875946264e+00, 6.204000000000000e-11},
{ 5.197274703601086e+00, 6.042797388545593e+01, 5.295999999999951e-11},
{ 5.160045430868086e+00, 3.938582310889583e+01, 1.811999999999999e-11},
{ 2.223277229900500e+00, 5.777452001970681e+01, 1.563999999999897e-11},
{ 6.921616240614999e-01, 2.103904892768837e+01, 1.535999999999999e-11},
{ 6.215802752555586e+00, 3.327229783044108e+01, 1.304000000000000e-11},
{ 1.368231494011000e+00, 3.938892179708230e+01, 1.661999999999999e-11},
{ 2.601539199675000e-01, 7.446011234647927e+00, 1.328000000000000e-11},
{ 5.903845687100000e-02, 3.941932465627423e+01, 8.559999999999978e-12},
{ 6.179310844209586e+00, 5.911910805492755e+01, 7.240000000000000e-12},
};
static const struct MarsSatTerm mars_sat_phobos_1[42] = {
{ 5.334169568271586e+00, 9.145914066599032e-03, 5.016848130000000e-05},
{ 5.334169148295587e+00, 9.145914113327495e-03, 5.016848130000000e-05},
{ 3.805345265960000e-01, 1.540952541904328e-03, 3.350264970000000e-05},
{ 3.805339484960000e-01, 1.540952613104260e-03, 3.350264903000000e-05},
{ 5.288363718535586e+00, 1.829243888046791e-02, 2.827839093000000e-05},
{ 5.288363809359586e+00, 1.829243887034939e-02, 2.827839093000000e-05},
{ 3.412685718051586e+00, 7.604696023359408e-03, 1.578862428000000e-05},
{ 3.412685998796586e+00, 7.604695992157328e-03, 1.578862427000000e-05},
{ 3.448999855869586e+00, 2.715567616927210e+01, 3.008844320000000e-05},
{ 9.194445444030001e-01, 3.101144183892481e-03, 8.178563170000001e-06},
{ 9.194443192780000e-01, 3.101144208596499e-03, 8.178563340000000e-06},
{ 5.387625479307586e+00, 1.969446275049270e+01, 2.044327843999999e-05},
{ 2.848996739890000e+00, 2.743829098982358e-02, 6.516951720000000e-06},
{ 2.848996748767000e+00, 2.743829098862537e-02, 6.516951720000000e-06},
{ 3.954411167839586e+00, 1.357783808463605e+01, 1.013316988000000e-05},
{ 2.912659283736000e+00, 4.073351425390815e+01, 7.483468759999740e-06},
{ 4.067858164824586e+00, 2.589727369338346e-02, 2.589359390000000e-06},
{ 4.067857973471586e+00, 2.589727371473152e-02, 2.589359390000000e-06},
{ 3.776271613850000e-01, 3.658417828853650e-02, 1.125147420000000e-06},
{ 3.776267985560000e-01, 3.658417832900902e-02, 1.125147420000000e-06},
{ 4.121936540919586e+00, 4.685014091671000e+01, 1.523532599999976e-06},
{ 1.231836210236000e+00, 1.674957754399965e-02, 4.976512300000000e-07},
{ 1.231836585135000e+00, 1.674957750247454e-02, 4.976512300000000e-07},
{ 2.488083095032000e+00, 1.520386207875350e-02, 6.841966800000000e-07},
{ 2.488083105051000e+00, 1.520386207728190e-02, 6.841966800000000e-07},
{ 1.576141880659500e+00, 3.504321962396548e-02, 1.124868580000000e-06},
{ 6.145328698142587e+00, 6.116624841770791e+00, 7.873411399999999e-07},
{ 2.626756105114500e+00, 5.431135233854420e+01, 1.040645720000000e-06},
{ 2.759579551372000e+00, 7.461213382164448e+00, 7.166789399999994e-07},
{ 6.005673273380586e+00, 1.314187651672942e+00, 6.022022799999936e-07},
{ 9.353090056600000e-02, 7.532115638546220e-04, 1.964976900000000e-07},
{ 9.343265452500001e-02, 7.532237916966188e-04, 1.964964500000000e-07},
{ 4.302928383505000e-01, 3.938582557777600e+01, 4.987307800000000e-07},
{ 3.587439959397586e+00, 6.042797969706294e+01, 4.804827999999994e-07},
{ 4.183547650369587e+00, 4.572998195120822e-02, 1.747133100000000e-07},
{ 4.183547815926586e+00, 4.572998193123842e-02, 1.747133100000000e-07},
{ 5.589444380180586e+00, 6.068411018120026e-03, 1.793076700000000e-07},
{ 5.589445119118587e+00, 6.068410936124942e-03, 1.793076800000000e-07},
{ 6.062696336024086e+00, 3.938892673561993e+01, 2.653755199999908e-07},
{ 4.628882088425586e+00, 3.327230135427769e+01, 2.539272800000000e-07},
{ 4.736462571780586e+00, 1.224558395734566e-02, 1.254715400000000e-07},
{ 4.736462588777586e+00, 1.224558396211886e-02, 1.254715500000000e-07},
};
static const struct MarsSatTerm mars_sat_phobos_2[27] = {
{ 1.404382124885000e+00, 7.595588511174286e-03, 1.514110912521000e-02},
{ 2.088385523034000e+00, 1.970205682773437e+01, 3.849620867400000e-04},
{ 3.895377361148586e+00, 1.069670280268190e-02, 6.903413242000000e-05},
{ 4.989572094750000e-01,-7.605272825959576e-03, 4.946101994000000e-05},
{ 8.215210221720000e-01, 4.685772969905357e+01, 3.671320788000000e-05},
{ 3.378076387479586e+00,-7.453616316894914e+00, 3.267983782000000e-05},
{ 2.772375719698000e+00, 3.939651873194720e+01, 8.750483050000000e-06},
{ 2.839303414769000e+00, 6.124220388371520e+00, 6.170938810000000e-06},
{ 3.184503170430000e-01, 6.043556479692822e+01, 6.957005940000000e-06},
{ 1.410329121503000e+00, 1.984322681200531e-02, 5.734088960000000e-06},
{ 1.336796551701000e+00, 3.327989343110183e+01, 3.431577860000000e-06},
{ 3.853646745987586e+00,-2.103145308554279e+01, 4.037409510000000e-06},
{ 2.547067366210000e-01, 1.549665594221090e-03, 2.115296220000000e-06},
{ 1.508224230072000e+00,-5.911151611514718e+01, 1.904827620000000e-06},
{ 3.462607551400000e-02,-5.165030141480156e+01, 8.696094700000001e-07},
{ 3.198655314207586e+00, 9.152774292423818e-03, 6.329868200000000e-07},
{ 5.163954670916587e+00, 1.674080796065678e-02, 6.130317800000000e-07},
{ 7.803937436990001e-01,-1.550610622055567e-03, 5.737217100000000e-07},
{ 3.264588234538586e+00, 2.716326950857758e+01, 8.181088300000000e-07},
{ 5.211329429707586e+00, 2.898903518276566e-02, 7.434145300000000e-07},
{ 2.691716763490000e+00,-2.714807769012140e+01, 5.441647800000001e-07},
{ 1.001733960626000e+00,-4.553367734785991e+01, 5.044515000000000e-07},
{ 3.212149602842586e+00,-1.968376490094933e+01, 4.284599100000000e-07},
{ 4.145658249560586e+00,-3.460928932187340e+01, 4.224937300000000e-07},
{ 5.063879137454586e+00, 2.589097180762378e-02, 3.581949800000000e-07},
{ 2.558639707373000e+00, 6.064928432880783e-03, 3.808881200000000e-07},
{ 5.292499025555586e+00, 3.075065445597794e-03, 4.055739200000000e-07},
};
static const struct MarsSatTerm mars_sat_phobos_3[28] = {
{ 2.058107128488000e+00,-7.604861328578004e-03, 9.408605183120001e-03},
{ 3.248856489376586e+00,-9.145863467943084e-03, 5.699538102000000e-05},
{ 4.557280987272586e+00, 1.829238959116374e-02, 2.474369483000000e-05},
{ 6.113423353213586e+00, 7.595702278066188e-03, 2.062715397000000e-05},
{ 2.067112903892000e+00, 2.743825478464904e-02, 5.945574160000000e-06},
{ 1.700689207004000e+00, 9.145729527513899e-03, 3.920190200000000e-06},
{ 3.577059353926586e+00,-1.829259493439447e-02, 2.167791970000000e-06},
{ 4.940284891770586e+00,-7.453616316894914e+00, 2.769763750000000e-06},
{ 2.117885747855000e+00, 3.941171892425920e+01, 1.813539170000000e-06},
{ 5.176055365230000e-01, 1.970205682773437e+01, 1.083429180000000e-06},
{ 5.864248782211586e+00, 3.658425999636452e-02, 1.066657010000000e-06},
{ 3.240692168340586e+00, 2.589639535977297e-02, 7.909853300000000e-07},
{ 7.570751994340000e-01,-1.321790869500707e+00, 8.880885099999999e-07},
{ 6.608910796060000e-01, 6.124220388371520e+00, 5.520758900000000e-07},
{ 5.519669134088586e+00, 4.685772969905357e+01, 4.397886000000000e-07},
{ 5.844000094873587e+00,-2.743825733444110e-02, 4.583609100000000e-07},
{ 4.602755395228586e+00,-1.675516082048598e-02, 3.533131100000000e-07},
{ 2.685649467156000e+00, 1.542518910643103e-03, 3.502988500000000e-07},
{ 3.353176877809586e+00, 1.225604581839110e+01, 3.828851700000000e-07},
{ 2.556761413073000e+00, 1.068971473943851e-02, 2.711247300000000e-07},
{ 4.622369718415587e+00,-2.589682634097252e-02, 1.991745400000000e-07},
{ 7.497690894330000e-01, 3.504362063026418e-02, 2.036897700000000e-07},
{ 3.965205602751586e+00, 2.714806830611846e+01, 1.726779600000000e-07},
{ 4.743858488241586e+00,-2.103145308554279e+01, 1.574472800000000e-07},
{ 3.405870070134586e+00, 4.572380100924578e-02, 1.698567700000000e-07},
{ 3.742380849210000e-01, 3.327989343110183e+01, 1.015004500000000e-07},
{ 1.374678754422000e+00,-1.970206671243377e+01, 8.259149000000001e-08},
{ 1.140160118874000e+00,-6.108652383295211e-03, 8.504923000000000e-08},
};
static const struct MarsSatTerm mars_sat_deimos_0[25] = {
{ 6.146803530569087e+00, 2.294413036846356e+00, 5.398360000000000e-09},
{ 3.432805186019586e+00, 9.935735425667586e+00, 7.095199999999997e-10},
{ 4.356581106602587e+00, 3.441619555269535e+00, 3.737800000000000e-10},
{ 5.915601483397086e+00, 9.926589650598594e+00, 2.337400000000000e-10},
{ 1.615243161446500e+00, 1.147206518423178e+00, 1.697000000000000e-10},
{ 2.114704906116000e+00, 9.917443791162608e+00, 5.314000000000000e-11},
{ 9.393388277920001e-01, 9.954027953247621e+00, 7.046000000000000e-11},
{ 3.949494792732586e+00, 9.944881541450725e+00, 4.233999999999998e-11},
{ 2.089322391969500e+00, 2.294097606691310e+00, 2.813999999999983e-11},
{ 7.795571756865000e-01, 2.294728469946876e+00, 2.549999999999971e-11},
{ 2.361851669664500e+00, 4.588826073692712e+00, 1.676000000000000e-11},
{ 2.770026883578500e+00, 9.954658773106329e+00, 9.719999999999976e-12},
{ 4.733522308629587e+00, 9.963173826265681e+00, 9.519999999999998e-12},
{ 3.258781379519086e+00, 1.472506551307767e+01, 9.739999999999990e-12},
{ 4.596299836113586e+00, 9.908297894321850e+00, 1.027999999999982e-11},
{ 3.057688443346500e+00, 2.682916339604879e+00, 6.679999999999999e-12},
{ 5.453932177963086e+00, 2.682285947505377e+00, 6.879999999999995e-12},
{ 4.328499239462586e+00, 9.936050924435195e+00, 6.199999999999996e-12},
{ 2.168321380562000e+00, 4.976699000854042e+00, 5.359999999999916e-12},
{ 2.516818588675000e+00, 9.935420104041299e+00, 3.359999999999995e-12},
{ 2.997827070725000e-01, 3.441304123224338e+00, 2.919999999999998e-12},
{ 2.143143420462500e+00, 2.682600911813873e+00, 2.979999999999999e-12},
{ 5.271781540988586e+00, 3.441934988580561e+00, 2.640000000000000e-12},
{ 5.382287205930000e-01, 9.926905039305762e+00, 2.039999999999966e-12},
{ 3.629187963751586e+00, 2.303558945733479e+00, 1.090000000000000e-12},
};
static const struct MarsSatTerm mars_sat_deimos_1[21] = {
{ 2.488175621789000e+00, 3.153377646687549e-04, 2.483031660190000e-03},
{ 2.488175634463000e+00, 3.153377641537758e-04, 2.483031659490000e-03},
{ 5.336889264964586e+00, 9.145915158660612e-03, 2.023426978200000e-04},
{ 5.336889346905586e+00, 9.145915155199243e-03, 2.023426978200000e-04},
{ 5.281538321363586e+00, 1.829244561727392e-02, 1.043167564200000e-04},
{ 5.281537958522586e+00, 1.829244563271799e-02, 1.043167564300000e-04},
{ 1.433490018835500e+00, 2.294413045224799e+00, 1.637133774599998e-04},
{ 2.852007788337500e+00, 2.743818628432718e-02, 4.820671420000000e-05},
{ 3.147176604035586e+00, 1.860783561096842e-02, 1.497066023000000e-05},
{ 3.147176307034586e+00, 1.860783562364318e-02, 1.497066034000000e-05},
{ 3.182516662993086e+00, 1.147206709032795e+00, 1.137028858000000e-05},
{ 5.917181875408586e+00, 3.441620120814550e+00, 8.350609860000001e-06},
{ 3.544871659744586e+00, 6.190897763118658e-04, 4.383254500000000e-06},
{ 3.544871649035586e+00, 6.190897769462623e-04, 4.383254500000000e-06},
{ 4.614193383590000e-01, 3.657933964533806e-02, 8.316176420000001e-06},
{ 6.754222685930000e-01, 2.775277029221846e-02, 6.447587580000001e-06},
{ 5.001640850015586e+00, 9.935735582783501e+00, 7.926687400000001e-06},
{ 4.485218623803586e+00, 8.828439970895009e-03, 3.905064210000000e-06},
{ 2.221013917980000e-01, 9.464721723028089e-03, 3.117801830000000e-06},
{ 4.485218622895586e+00, 8.828439971019056e-03, 3.905064210000000e-06},
{ 2.221013944660000e-01, 9.464721722885516e-03, 3.117801830000000e-06},
};
static const struct MarsSatTerm mars_sat_deimos_2[15] = {
{ 2.198649419514000e+00, 3.148401892560942e-04, 2.744131534600000e-04},
{ 4.366636518977586e+00, 4.977013897776327e+00, 6.015912711000000e-05},
{ 1.463816210370000e+00,-3.153164045300449e-04, 3.614585349000000e-05},
{ 1.362467588064000e+00, 2.682600831640474e+00, 2.577157903000000e-05},
{ 9.336792893320000e-01,-4.958721582505435e+00, 6.801997830000000e-06},
{ 3.152432010349586e+00, 1.535394310272536e+00, 4.452979140000000e-06},
{ 4.734051493918586e+00,-4.949575691510945e+00, 2.238482920000000e-06},
{ 5.736642853342587e+00, 2.327524649334390e-04, 1.478747570000000e-06},
{ 5.084945958505586e+00, 3.912335403634271e-04, 1.457720040000000e-06},
{ 4.232770650771586e+00, 7.271426904652134e+00, 1.347030310000000e-06},
{ 1.516439412109000e+00, 1.491274947692300e+01, 7.538179100000001e-07},
{ 3.100618466711000e+00, 1.797748798526167e-02, 8.716229700000000e-07},
{ 5.155432690964586e+00, 3.881877876517008e-01, 1.015648880000000e-06},
{ 2.250768124831000e+00,-4.940429834739583e+00, 5.090989300000000e-07},
{ 3.422251338868586e+00,-4.977013897776327e+00, 6.537077500000000e-07},
};
static const struct MarsSatTerm mars_sat_deimos_3[27] = {
{ 2.981506933511000e+00,-3.154811355556041e-04, 1.562693319959000e-02},
{ 4.557500894366586e+00, 1.829233626168517e-02, 1.321818631100000e-04},
{ 3.248124065112586e+00,-9.145934570587952e-03, 3.833652719000000e-05},
{ 1.701551338710000e+00, 9.145665986601770e-03, 2.660211842000000e-05},
{ 2.067662003165000e+00, 2.743821959085702e-02, 2.882438535000000e-05},
{ 4.805327222478586e+00, 3.158564952832469e-04, 2.713213911000000e-05},
{ 2.318788380092000e+00, 1.860776219134252e-02, 9.296324950000000e-06},
{ 3.613580709014586e+00,-1.829261290520875e-02, 4.460960290000000e-06},
{ 5.867810550672586e+00, 3.658405181871115e-02, 4.904156030000000e-06},
{ 5.496515547817586e+00,-9.461294099629817e-03, 2.366881150000000e-06},
{ 6.120362957845586e+00, 2.775361064625920e-02, 2.057552850000000e-06},
{ 3.612389208841586e+00, 8.830178627407989e-03, 2.365750140000000e-06},
{ 5.554163001934587e+00,-1.860800990265229e-02, 1.228814440000000e-06},
{ 3.561485983136586e+00, 1.797648948023003e-02, 1.215189460000000e-06},
{ 1.941975361869000e+00,-6.259843043363501e-04, 9.602856999999999e-07},
{ 5.845244776759587e+00,-2.743819348074506e-02, 1.174471770000000e-06},
{ 1.821048365911000e+00, 9.460829229838182e-03, 1.058441050000000e-06},
{ 7.587135454610000e-01, 1.682177864988152e-04, 2.383147490000000e-06},
{ 3.292891638604586e+00, 4.573365667404765e-02, 7.656303800000000e-07},
{ 3.766016108380586e+00, 3.683140215171970e-02, 4.075824300000000e-07},
{ 5.747244956219586e+00, 9.954342722363187e+00, 4.962340900000000e-07},
{ 1.724062868741000e+00,-2.774017744689255e-02, 3.003077900000000e-07},
{ 4.429747729642586e+00,-8.904754331211824e-03, 2.613416400000000e-07},
{ 8.659161909010000e-01, 2.718340176462802e-02, 2.891945200000000e-07},
{ 2.259603352363000e+00,-3.660152378331885e-02, 2.213788200000000e-07},
{ 1.108725025984000e+00, 9.935732440466172e+00, 2.449950300000000e-07},
{ 1.680391963791000e+00, 9.076040299352415e-03, 1.652416600000000e-07},
};
static const struct MarsSatBody mars_sat_bodies[2] = {
{
9.549547741038312e-11, 1.970205562831390e+01, 1.657852042683113e-10,
{
0.0000626916188000,
2.0912973926417302,
-0.0000005774003141,
0.0000006185052790,
-0.0000575018919826,
-0.0000540676493351,
},
{
{mars_sat_phobos_0, 16},
{mars_sat_phobos_1, 42},
{mars_sat_phobos_2, 27},
{mars_sat_phobos_3, 28}
}
},
{
9.549547622768120e-11, 4.977013889652300e+00, -2.331793571572226e-14,
{
0.0001568134086700,
-1.9167537810740201,
-0.0000239579419518,
0.0000257893300229,
-0.0056443721379635,
-0.0053121806978560,
},
{
{mars_sat_deimos_0, 25},
{mars_sat_deimos_1, 21},
{mars_sat_deimos_2, 15},
{mars_sat_deimos_3, 27}
}
},
};
static
void CalcMarsSatElem(double t, int body, double elem[6]) {
int j;
const struct MarsSatBody *bp = mars_sat_bodies + body;
memcpy(elem, bp->constants, 6 * sizeof(double));
for (j = 0; j < 2; j++) {
const struct MarsSatTerm *const begin = bp->lists[j].terms;
const struct MarsSatTerm *p = begin + bp->lists[j].size;
while (--p >= begin) {
const double d = p->phase + t * p->frequency;
elem[j] += p->amplitude * cos(d);
}
}
for (j = 2; j < 4; j++) {
const struct MarsSatTerm *const begin = bp->lists[j].terms;
const struct MarsSatTerm *p = begin + bp->lists[j].size;
while (--p >= begin) {
const double d = p->phase + t * p->frequency;
elem[2*j-2] += p->amplitude * cos(d);
elem[2*j-1] += p->amplitude * sin(d);
}
}
elem[1] += (bp->l + bp->acc * t) * t;
}
static
void MultMat(const double a[9], const double b[9], double c[9]) {
int i, j;
for (i = 0; i < 3; i++) {
for (j = 0; j < 3; j++) {
c[i*3+j] = a[i*3]*b[j] + a[i*3+1]*b[3+j] + a[i*3+2]*b[6+j];
}
}
}
static const double J2000_to_VSOP87[9] = {
9.999999999998848e-01,-4.799655442984222e-07, 0.000000000000000e+00,
4.403598133110236e-07, 9.174821370868568e-01, 3.977769829016507e-01,
-1.909192461077750e-07,-3.977769829016049e-01, 9.174821370869625e-01};
static const double ome0 = 47.68143;
static const double inc0 = 37.1135;
static const double dome = -0.1061;
static const double dinc = 0.0609;
static
void GenerateMarsSatToVSOP87(double t, double mat[9]) {
double ome, inc, co, so, ci, si;
double m[9];
t -= 6491.5;
ome = (ome0 + dome * t / 36525.) * (M_PI / 180.0);
inc = (inc0 + dinc * t / 36525.) * (M_PI / 180.0);
co = cos(ome);
so = sin(ome);
ci = cos(inc);
si = sin(inc);
m[0] = co; m[1] = -ci*so; m[2] = si*so;
m[3] = so; m[4] = ci*co; m[5] = -si*co;
m[6] = 0.0; m[7] = si; m[8] = ci;
MultMat(J2000_to_VSOP87, m, mat);
}
void GetMarsSatCoor(double jd, int body, double *xyz, double *xyzdot) {
double elem[6];
double x[6];
double mat[9];
const double t = jd - 2451545.0 + 6491.5;
CalcMarsSatElem(t, body, elem);
GenerateMarsSatToVSOP87(t, mat);
EllipticToRectangularA(mars_sat_bodies[body].mu, elem, 0.0, x);
xyz[0] = mat[0]*x[0] + mat[1]*x[1] + mat[2]*x[2];
xyz[1] = mat[3]*x[0] + mat[4]*x[1] + mat[5]*x[2];
xyz[2] = mat[6]*x[0] + mat[7]*x[1] + mat[8]*x[2];
if (xyzdot) {
xyzdot[0] = mat[0]*x[3] + mat[1]*x[4] + mat[2]*x[5];
xyzdot[1] = mat[3]*x[3] + mat[4]*x[4] + mat[5]*x[5];
xyzdot[2] = mat[6]*x[3] + mat[7]*x[4] + mat[8]*x[5];
}
}

81
src/marssat.h Normal file
View File

@ -0,0 +1,81 @@
/************************************************************************
The Ephemerides of the Martian satellites
(adjustement from 1877 to 2005, Version 1.0)
by Valery Lainey can be obtained from Valery Lainey:
V.Lainey (Lainey@oma.be)
ROB- 3, Avenue Circulaire, B-1180 Bruxelles (Belgium)
IMCCE - 77, Avenue Denfert-Rochereau 75014 Paris (France)
-----------------------------------------------------------------------
I (Johannes Gajdosik) have just taken Valery Laineys Fortran code,
MarsSatV1-0.f, which he kindly supplied, and rearranged it into
this piece of software.
I can neither allow nor forbid the usage of Valery Laineys
Ephemerides of the Martian satellites.
The copyright notice below covers not the work of Valery Lainey
but just my work, that is the compilation of Valery Laineys
Ephemerides of the Martian satellites into the software supplied in this file.
Copyright (c) 2006 Johannes Gajdosik
Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense,
and/or sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
Derived from Stellarium's MARSSAT implementation.
Modified for pg_orbit: removed static mutable state and
CalcInterpolatedElements caching for thread safety
(PostgreSQL PARALLEL SAFE). Elements are computed fresh on each call.
1) do not calculate constant terms at runtime but beforehand
2) unite terms with the same frequencies
****************************************************************/
#ifndef PG_ORBIT_MARSSAT_H
#define PG_ORBIT_MARSSAT_H
#ifdef __cplusplus
extern "C" {
#endif
#define MARS_SAT_PHOBOS 0
#define MARS_SAT_DEIMOS 1
void GetMarsSatCoor(double jd, int body, double *xyz, double *xyzdot);
/* Return the rectangular coordinates and velocity of the given satellite
and the given julian date jd expressed in dynamical time (TAI+32.184s).
The origin of the xyz-coordinates is the center of Mars.
The reference frame is "dynamical equinox and ecliptic J2000",
which is the reference frame in VSOP87 and VSOP87A.
body: 0=Phobos, 1=Deimos
xyz[3]: position (AU), relative to Mars center
xyzdot[3]: velocity (AU/day), may be NULL
*/
#ifdef __cplusplus
}
#endif
#endif

248
src/moon_funcs.c Normal file
View File

@ -0,0 +1,248 @@
/*
* moon_funcs.c -- Planetary moon observation
*
* SQL functions for observing moons of Jupiter, Saturn, Uranus, and Mars.
* Each moon's position is computed in the VSOP87 ecliptic J2000 frame
* relative to its parent planet, then combined with the parent's
* heliocentric position to get geocentric az/el.
*
* Pipeline for each moon:
* 1. Parent planet heliocentric position (VSOP87)
* 2. Moon position relative to parent (L12/TASS17/GUST86/MARSSAT)
* 3. Moon heliocentric = parent + moon_relative
* 4. Moon geocentric = moon_heliocentric - Earth_heliocentric
* 5. Ecliptic -> equatorial -> precess -> az/el
*/
#include "postgres.h"
#include "fmgr.h"
#include "funcapi.h"
#include "utils/timestamp.h"
#include "types.h"
#include "astro_math.h"
#include "vsop87.h"
#include "l12.h"
#include "tass17.h"
#include "gust86.h"
#include "marssat.h"
#include <math.h>
PG_FUNCTION_INFO_V1(galilean_observe);
PG_FUNCTION_INFO_V1(saturn_moon_observe);
PG_FUNCTION_INFO_V1(uranus_moon_observe);
PG_FUNCTION_INFO_V1(mars_moon_observe);
/* ================================================================
* Internal: geocentric observation pipeline
*
* Duplicated from planet_funcs.c to avoid cross-TU coupling.
* (Same logic: ecliptic J2000 -> equatorial -> precess -> az/el)
* ================================================================
*/
static void
moon_observe_from_geocentric(const double geo_ecl_au[3], double jd,
const pg_observer *obs, pg_topocentric *result)
{
double geo_equ[3];
double ra_j2000, dec_j2000, geo_dist;
double ra_date, dec_date;
double gmst_val, lst, ha;
double az, el;
ecliptic_to_equatorial(geo_ecl_au, geo_equ);
cartesian_to_spherical(geo_equ, &ra_j2000, &dec_j2000, &geo_dist);
precess_j2000_to_date(jd, ra_j2000, dec_j2000, &ra_date, &dec_date);
gmst_val = gmst_from_jd(jd);
lst = gmst_val + obs->lon;
ha = lst - ra_date;
equatorial_to_horizontal(ha, dec_date, obs->lat, &az, &el);
result->azimuth = az;
result->elevation = el;
result->range_km = geo_dist * AU_KM;
result->range_rate = 0.0;
}
/* ================================================================
* Internal: common pattern for all planetary moons
*
* Given: moon position relative to parent (VSOP87 ecliptic J2000, AU)
* parent's VSOP87 body index (0-based)
* observer, JD
*
* Computes geocentric position and returns topocentric az/el.
* ================================================================
*/
static void
observe_planetary_moon(const double moon_rel[3], int vsop_parent,
double jd, const pg_observer *obs,
pg_topocentric *result)
{
double parent_xyz[6];
double earth_xyz[6];
double geo_ecl[3];
/* Parent planet heliocentric */
GetVsop87Coor(jd, vsop_parent, parent_xyz);
/* Earth heliocentric */
GetVsop87Coor(jd, 2, earth_xyz); /* VSOP87 body 2 = Earth */
/* Moon geocentric = (parent + moon_relative) - Earth */
geo_ecl[0] = (parent_xyz[0] + moon_rel[0]) - earth_xyz[0];
geo_ecl[1] = (parent_xyz[1] + moon_rel[1]) - earth_xyz[1];
geo_ecl[2] = (parent_xyz[2] + moon_rel[2]) - earth_xyz[2];
moon_observe_from_geocentric(geo_ecl, jd, obs, result);
}
/* ================================================================
* galilean_observe(body_id int, observer, timestamptz) -> topocentric
*
* Observe a Galilean moon of Jupiter.
* Body IDs: 0=Io, 1=Europa, 2=Ganymede, 3=Callisto
*
* Uses L1.2 theory (Lainey, Duriez & Vienne) for moon positions
* and VSOP87 for Jupiter and Earth.
* ================================================================
*/
Datum
galilean_observe(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1);
int64 ts = PG_GETARG_INT64(2);
double jd;
double moon_xyz[3];
pg_topocentric *result;
if (body_id < L12_IO || body_id > L12_CALLISTO)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("galilean_observe: body_id %d must be 0-3 (Io, Europa, Ganymede, Callisto)",
body_id)));
jd = timestamptz_to_jd(ts);
/* Moon position relative to Jupiter, VSOP87 ecliptic J2000, AU */
GetL12Coor(jd, body_id, moon_xyz, NULL);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_planetary_moon(moon_xyz, 4, jd, obs, result); /* VSOP87 body 4 = Jupiter */
PG_RETURN_POINTER(result);
}
/* ================================================================
* saturn_moon_observe(body_id int, observer, timestamptz) -> topocentric
*
* Observe a moon of Saturn.
* Body IDs: 0=Mimas, 1=Enceladus, 2=Tethys, 3=Dione,
* 4=Rhea, 5=Titan, 6=Iapetus, 7=Hyperion
*
* Uses TASS 1.7 theory (Vienne & Duriez).
* ================================================================
*/
Datum
saturn_moon_observe(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1);
int64 ts = PG_GETARG_INT64(2);
double jd;
double moon_xyz[3];
pg_topocentric *result;
if (body_id < TASS17_MIMAS || body_id > TASS17_HYPERION)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("saturn_moon_observe: body_id %d must be 0-7 (Mimas-Hyperion)",
body_id)));
jd = timestamptz_to_jd(ts);
GetTass17Coor(jd, body_id, moon_xyz, NULL);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_planetary_moon(moon_xyz, 5, jd, obs, result); /* VSOP87 body 5 = Saturn */
PG_RETURN_POINTER(result);
}
/* ================================================================
* uranus_moon_observe(body_id int, observer, timestamptz) -> topocentric
*
* Observe a moon of Uranus.
* Body IDs: 0=Miranda, 1=Ariel, 2=Umbriel, 3=Titania, 4=Oberon
*
* Uses GUST86 theory (Laskar & Jacobson).
* ================================================================
*/
Datum
uranus_moon_observe(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1);
int64 ts = PG_GETARG_INT64(2);
double jd;
double moon_xyz[3];
pg_topocentric *result;
if (body_id < GUST86_MIRANDA || body_id > GUST86_OBERON)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("uranus_moon_observe: body_id %d must be 0-4 (Miranda-Oberon)",
body_id)));
jd = timestamptz_to_jd(ts);
GetGust86Coor(jd, body_id, moon_xyz, NULL);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_planetary_moon(moon_xyz, 6, jd, obs, result); /* VSOP87 body 6 = Uranus */
PG_RETURN_POINTER(result);
}
/* ================================================================
* mars_moon_observe(body_id int, observer, timestamptz) -> topocentric
*
* Observe a moon of Mars.
* Body IDs: 0=Phobos, 1=Deimos
*
* Uses MarsSat theory (Lainey, 2007).
* ================================================================
*/
Datum
mars_moon_observe(PG_FUNCTION_ARGS)
{
int32 body_id = PG_GETARG_INT32(0);
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(1);
int64 ts = PG_GETARG_INT64(2);
double jd;
double moon_xyz[3];
pg_topocentric *result;
if (body_id < MARS_SAT_PHOBOS || body_id > MARS_SAT_DEIMOS)
ereport(ERROR,
(errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE),
errmsg("mars_moon_observe: body_id %d must be 0-1 (Phobos, Deimos)",
body_id)));
jd = timestamptz_to_jd(ts);
GetMarsSatCoor(jd, body_id, moon_xyz, NULL);
result = (pg_topocentric *) palloc(sizeof(pg_topocentric));
observe_planetary_moon(moon_xyz, 3, jd, obs, result); /* VSOP87 body 3 = Mars */
PG_RETURN_POINTER(result);
}

262
src/radio_funcs.c Normal file
View File

@ -0,0 +1,262 @@
/*
* radio_funcs.c -- Jupiter decametric radio burst prediction
*
* Jupiter emits intense radio bursts at 4-39.5 MHz, modulated by:
* 1. Io's orbital phase around Jupiter (Io-related storms)
* 2. Jupiter's Central Meridian Longitude (CML, System III)
*
* Radio JOVE operators use CML-vs-Io-phase lookup charts to predict
* burst probability. This module computes both parameters from SQL,
* enabling batch prediction queries over time series.
*
* References:
* Carr, Desch & Alexander (1983), "Phenomenology of Magnetospheric
* Radio Emissions", in Physics of the Jovian Magnetosphere.
* Radio JOVE project: https://radiojove.gsfc.nasa.gov/
*/
#include "postgres.h"
#include "fmgr.h"
#include "funcapi.h"
#include "utils/timestamp.h"
#include "types.h"
#include "astro_math.h"
#include "vsop87.h"
#include "l12.h"
#include <math.h>
PG_FUNCTION_INFO_V1(io_phase_angle);
PG_FUNCTION_INFO_V1(jupiter_cml);
PG_FUNCTION_INFO_V1(jupiter_burst_probability);
/* ================================================================
* io_phase_angle(timestamptz) -> float8
*
* Returns Io's orbital phase angle in degrees [0, 360).
* Measured from superior geocentric conjunction:
* 0 = Io behind Jupiter (as seen from Earth)
* 90 = Io east of Jupiter (western elongation)
* 180 = Io in front of Jupiter (inferior conjunction)
* 270 = Io west of Jupiter (eastern elongation)
*
* This is the standard convention used by Radio JOVE charts.
* ================================================================
*/
Datum
io_phase_angle(PG_FUNCTION_ARGS)
{
int64 ts = PG_GETARG_INT64(0);
double jd;
double io_xyz[3];
double jup_xyz[6], earth_xyz[6];
double jup_geo[3];
double io_rel_geo[3];
double phase;
double proj_x, proj_y;
double jup_dist;
double east_x, east_y, east_z;
double north_x, north_y, north_z;
double io_east, io_north;
jd = timestamptz_to_jd(ts);
/* Io position relative to Jupiter (VSOP87 ecliptic J2000, AU) */
GetL12Coor(jd, L12_IO, io_xyz, NULL);
/* Jupiter and Earth heliocentric positions */
GetVsop87Coor(jd, 4, jup_xyz); /* VSOP87 body 4 = Jupiter */
GetVsop87Coor(jd, 2, earth_xyz); /* VSOP87 body 2 = Earth */
/* Jupiter geocentric direction */
jup_geo[0] = jup_xyz[0] - earth_xyz[0];
jup_geo[1] = jup_xyz[1] - earth_xyz[1];
jup_geo[2] = jup_xyz[2] - earth_xyz[2];
jup_dist = sqrt(jup_geo[0] * jup_geo[0] +
jup_geo[1] * jup_geo[1] +
jup_geo[2] * jup_geo[2]);
/* Normalize Jupiter geocentric direction */
jup_geo[0] /= jup_dist;
jup_geo[1] /= jup_dist;
jup_geo[2] /= jup_dist;
/*
* Project Io's position onto the plane perpendicular to
* the Earth-Jupiter line of sight. We need two orthogonal
* axes in this plane: "east" (position angle) and "north"
* (toward north ecliptic pole).
*
* North ecliptic pole = (0, 0, 1) in ecliptic J2000.
* East = North x LineOfSight (cross product)
*/
north_x = -(jup_geo[2] * jup_geo[0]); /* north - (north.los)*los */
north_y = -(jup_geo[2] * jup_geo[1]);
north_z = 1.0 - jup_geo[2] * jup_geo[2];
{
double nmag = sqrt(north_x * north_x + north_y * north_y + north_z * north_z);
if (nmag > 1e-10) {
north_x /= nmag;
north_y /= nmag;
north_z /= nmag;
}
}
/* East = North x LOS */
east_x = north_y * jup_geo[2] - north_z * jup_geo[1];
east_y = north_z * jup_geo[0] - north_x * jup_geo[2];
east_z = north_x * jup_geo[1] - north_y * jup_geo[0];
/* Project Io relative position onto east/north */
io_east = io_xyz[0] * east_x + io_xyz[1] * east_y + io_xyz[2] * east_z;
io_north = io_xyz[0] * north_x + io_xyz[1] * north_y + io_xyz[2] * north_z;
/* Along LOS component (positive = behind Jupiter from Earth) */
proj_x = io_xyz[0] * jup_geo[0] + io_xyz[1] * jup_geo[1] + io_xyz[2] * jup_geo[2];
proj_y = io_east;
/*
* Phase angle convention (Radio JOVE):
* 0 = superior conjunction (Io behind Jupiter)
* Measured counter-clockwise as seen from north ecliptic pole
*/
phase = atan2(-proj_y, proj_x) * RAD_TO_DEG;
if (phase < 0.0)
phase += 360.0;
PG_RETURN_FLOAT8(phase);
}
/* ================================================================
* jupiter_cml(observer, timestamptz) -> float8
*
* Jupiter's Central Meridian Longitude in System III (1965.0).
* Returns degrees [0, 360).
*
* System III is defined by Jupiter's magnetic field rotation,
* with period 9h 55m 29.711s (IAU).
*
* The CML is the System III longitude of the sub-observer point.
* For Earth-based observation, this is approximately the
* sub-Earth longitude.
*
* Reference: IAU (1976) System III longitude definition.
* ================================================================
*/
Datum
jupiter_cml(PG_FUNCTION_ARGS)
{
pg_observer *obs = (pg_observer *) PG_GETARG_POINTER(0);
int64 ts = PG_GETARG_INT64(1);
double jd;
double jup_xyz[6], earth_xyz[6];
double dx, dy, dz, geo_dist;
double light_time_days;
double jd_lt;
double d;
double cml;
(void)obs; /* observer not needed for CML (geocentric is sufficient) */
jd = timestamptz_to_jd(ts);
/* Jupiter and Earth positions */
GetVsop87Coor(jd, 4, jup_xyz);
GetVsop87Coor(jd, 2, earth_xyz);
/* Geocentric distance for light-time correction */
dx = jup_xyz[0] - earth_xyz[0];
dy = jup_xyz[1] - earth_xyz[1];
dz = jup_xyz[2] - earth_xyz[2];
geo_dist = sqrt(dx * dx + dy * dy + dz * dz);
/* Light-time correction (speed of light = 173.1446 AU/day) */
light_time_days = geo_dist / 173.1446327;
jd_lt = jd - light_time_days;
/*
* System III longitude of central meridian.
*
* Reference epoch: J1965.0 = JD 2438761.5
* Initial longitude: 305.38 deg (IAU 1976 convention)
* Rotation rate: 870.536 deg/day (System III, 9h 55m 29.711s period)
*
* The sub-Earth longitude decreases with time (Jupiter rotates
* prograde), so CML = initial + rate * elapsed_days.
*/
d = jd_lt - 2438761.5; /* Days since J1965.0 */
cml = 305.38 + 870.536 * d;
/* Normalize to [0, 360) */
cml = fmod(cml, 360.0);
if (cml < 0.0)
cml += 360.0;
PG_RETURN_FLOAT8(cml);
}
/* ================================================================
* jupiter_burst_probability(io_phase float8, cml float8) -> float8
*
* Estimated probability of Jupiter decametric radio burst
* based on Io phase angle and CML (System III).
*
* Returns a value 0.0-1.0 representing relative likelihood.
*
* Based on the empirical source regions from Carr et al. (1983):
* Source A: CML 200-260, Io phase 195-265 (Io-A)
* Source B: CML 90-200, Io phase 70-110 (Io-B)
* Source C: CML 290-10, Io phase 215-260 (Io-C)
* Source D: CML 0-50, no Io dependence (non-Io-D)
*
* This is a simplified model. Real-world burst probability
* depends on additional factors (frequency, solar activity,
* antenna orientation).
* ================================================================
*/
Datum
jupiter_burst_probability(PG_FUNCTION_ARGS)
{
double io_phase = PG_GETARG_FLOAT8(0);
double cml = PG_GETARG_FLOAT8(1);
double p;
/* Normalize inputs */
io_phase = fmod(io_phase, 360.0);
if (io_phase < 0.0) io_phase += 360.0;
cml = fmod(cml, 360.0);
if (cml < 0.0) cml += 360.0;
p = 0.0;
/* Source A (Io-A): strongest, most predictable */
if (cml >= 200.0 && cml <= 260.0 && io_phase >= 195.0 && io_phase <= 265.0)
p = 0.8;
/* Source B (Io-B): second strongest */
if (cml >= 90.0 && cml <= 200.0 && io_phase >= 70.0 && io_phase <= 110.0)
{
double pb = 0.6;
if (pb > p) p = pb;
}
/* Source C (Io-C): wraps around 360/0 boundary */
if ((cml >= 290.0 || cml <= 10.0) && io_phase >= 215.0 && io_phase <= 260.0)
{
double pc = 0.5;
if (pc > p) p = pc;
}
/* Non-Io-D: CML-dependent only, weaker */
if (cml >= 0.0 && cml <= 50.0)
{
double pd = 0.15;
if (pd > p) p = pd;
}
PG_RETURN_FLOAT8(p);
}

3156
src/tass17.c Normal file
View File

File diff suppressed because it is too large Load Diff

71
src/tass17.h Normal file
View File

@ -0,0 +1,71 @@
/************************************************************************
The TASS 1.7 theory of the Saturnian satellites including HYPERION
by Alain VIENNE and Luc DURIEZ can be found at
ftp://ftp.imcce.fr/pub/ephem/satel/tass17
I (Johannes Gajdosik) have just taken the Fortran code and data
obtained from above and rearranged it into this piece of software.
I can neither allow nor forbid the usage of the TASS 1.7 theory.
The copyright notice below covers not the work of Alain VIENNE and Luc DURIEZ
but just my work, that is the compilation of the TASS 1.7 theory
into the software supplied in this file.
Copyright (c) 2005 Johannes Gajdosik
Permission is hereby granted, free of charge, to any person obtaining a
copy of this software and associated documentation files (the "Software"),
to deal in the Software without restriction, including without limitation
the rights to use, copy, modify, merge, publish, distribute, sublicense,
and/or sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
Adapted for pg_orbit: removed static mutable state for thread safety
(PostgreSQL PARALLEL SAFE).
****************************************************************/
#ifndef PG_ORBIT_TASS17_H
#define PG_ORBIT_TASS17_H
#ifdef __cplusplus
extern "C" {
#endif
#define TASS17_MIMAS 0
#define TASS17_ENCELADUS 1
#define TASS17_TETHYS 2
#define TASS17_DIONE 3
#define TASS17_RHEA 4
#define TASS17_TITAN 5
#define TASS17_HYPERION 7
#define TASS17_IAPETUS 6
/*
* Compute Saturn moon position and velocity in VSOP87 ecliptic J2000
* coordinates, relative to Saturn center, in AU and AU/day.
*
* body: 0-7 (see defines above; note Hyperion=7, Iapetus=6)
* xyz[3]: output position (AU)
* xyzdot[3]: output velocity (AU/day), may be NULL
*/
void GetTass17Coor(double jd, int body, double *xyz, double *xyzdot);
#ifdef __cplusplus
}
#endif
#endif

226
test/expected/moon_observe.out Normal file
View File

@ -0,0 +1,226 @@
-- moon_observe regression tests
--
-- Tests planetary moon observation (Galilean, Saturn, Uranus, Mars)
-- and Jupiter decametric radio burst prediction functions.
-- Reference distances from JPL/NASA fact sheets.
\set boulder '''40.015N 105.270W 1655m'''::observer
-- ============================================================
-- Test 1: Galilean moon observation - all four from Boulder
-- Io, Europa, Ganymede, Callisto should return valid topocentric.
-- ============================================================
SELECT 'galilean_io' AS test,
round(topo_azimuth(galilean_observe(0, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(galilean_observe(0, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS el_deg;
test | az_deg | el_deg
-------------+--------+--------
galilean_io | 270 | 24
(1 row)
SELECT 'galilean_europa' AS test,
round(topo_azimuth(galilean_observe(1, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(galilean_observe(1, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS el_deg;
test | az_deg | el_deg
-----------------+--------+--------
galilean_europa | 270 | 24
(1 row)
SELECT 'galilean_ganymede' AS test,
round(topo_azimuth(galilean_observe(2, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(galilean_observe(2, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS el_deg;
test | az_deg | el_deg
-------------------+--------+--------
galilean_ganymede | 270 | 24
(1 row)
SELECT 'galilean_callisto' AS test,
round(topo_azimuth(galilean_observe(3, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(galilean_observe(3, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS el_deg;
test | az_deg | el_deg
-------------------+--------+--------
galilean_callisto | 270 | 24
(1 row)
-- ============================================================
-- Test 2: Galilean moons should be very close to Jupiter
-- All four should have ranges within ~0.05 AU of Jupiter's range.
-- Jupiter is ~4-6 AU from Earth; moons orbit within 0.013 AU.
-- ============================================================
SELECT 'galilean_near_jupiter' AS test,
round(topo_range(planet_observe(5, :boulder,
'2024-03-15 03:00:00+00'))::numeric, -4) AS jupiter_km,
round(topo_range(galilean_observe(0, :boulder,
'2024-03-15 03:00:00+00'))::numeric, -4) AS io_km,
round(topo_range(galilean_observe(3, :boulder,
'2024-03-15 03:00:00+00'))::numeric, -4) AS callisto_km;
test | jupiter_km | io_km | callisto_km
-----------------------+------------+-----------+-------------
galilean_near_jupiter | 837420000 | 837000000 | 837600000
(1 row)
-- ============================================================
-- Test 3: Saturn moon observation - Titan (body_id=5)
-- Titan is Saturn's largest moon, should be near Saturn.
-- ============================================================
SELECT 'saturn_titan' AS test,
round(topo_azimuth(saturn_moon_observe(5, :boulder,
'2024-06-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(saturn_moon_observe(5, :boulder,
'2024-06-15 03:00:00+00'))::numeric, 0) AS el_deg;
test | az_deg | el_deg
--------------+--------+--------
saturn_titan | 50 | -45
(1 row)
-- ============================================================
-- Test 4: Saturn moons - Mimas through Iapetus all return results
-- ============================================================
SELECT 'saturn_moons' AS test,
body_id,
round(topo_range(saturn_moon_observe(body_id, :boulder,
'2024-06-15 03:00:00+00'))::numeric, -4) AS range_km
FROM generate_series(0, 7) AS body_id;
test | body_id | range_km
--------------+---------+------------
saturn_moons | 0 | 1427920000
saturn_moons | 1 | 1427670000
saturn_moons | 2 | 1427910000
saturn_moons | 3 | 1427680000
saturn_moons | 4 | 1427420000
saturn_moons | 5 | 1426520000
saturn_moons | 6 | 1428150000
saturn_moons | 7 | 1429140000
(8 rows)
-- ============================================================
-- Test 5: Uranus moon observation - Titania (body_id=3)
-- ============================================================
SELECT 'uranus_titania' AS test,
round(topo_azimuth(uranus_moon_observe(3, :boulder,
'2024-06-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(uranus_moon_observe(3, :boulder,
'2024-06-15 03:00:00+00'))::numeric, 0) AS el_deg;
test | az_deg | el_deg
----------------+--------+--------
uranus_titania | 329 | -25
(1 row)
-- ============================================================
-- Test 6: All Uranus moons return valid results
-- ============================================================
SELECT 'uranus_moons' AS test,
body_id,
round(topo_range(uranus_moon_observe(body_id, :boulder,
'2024-06-15 03:00:00+00'))::numeric, -4) AS range_km
FROM generate_series(0, 4) AS body_id;
test | body_id | range_km
--------------+---------+------------
uranus_moons | 0 | 3061280000
uranus_moons | 1 | 3061360000
uranus_moons | 2 | 3061200000
uranus_moons | 3 | 3061310000
uranus_moons | 4 | 3061520000
(5 rows)
-- ============================================================
-- Test 7: Mars moons - Phobos and Deimos
-- ============================================================
SELECT 'mars_phobos' AS test,
round(topo_azimuth(mars_moon_observe(0, :boulder,
'2024-01-15 06:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(mars_moon_observe(0, :boulder,
'2024-01-15 06:00:00+00'))::numeric, 0) AS el_deg;
test | az_deg | el_deg
-------------+--------+--------
mars_phobos | 1 | -74
(1 row)
SELECT 'mars_deimos' AS test,
round(topo_azimuth(mars_moon_observe(1, :boulder,
'2024-01-15 06:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(mars_moon_observe(1, :boulder,
'2024-01-15 06:00:00+00'))::numeric, 0) AS el_deg;
test | az_deg | el_deg
-------------+--------+--------
mars_deimos | 1 | -74
(1 row)
-- ============================================================
-- Test 8: Io phase angle - should be [0, 360)
-- ============================================================
SELECT 'io_phase' AS test,
round(io_phase_angle('2024-03-15 03:00:00+00')::numeric, 1) AS phase_deg,
io_phase_angle('2024-03-15 03:00:00+00') >= 0.0
AND io_phase_angle('2024-03-15 03:00:00+00') < 360.0 AS in_range;
test | phase_deg | in_range
----------+-----------+----------
io_phase | 179.7 | t
(1 row)
-- ============================================================
-- Test 9: Jupiter CML (System III) - should be [0, 360)
-- ============================================================
SELECT 'jupiter_cml' AS test,
round(jupiter_cml(:boulder, '2024-03-15 03:00:00+00')::numeric, 1) AS cml_deg,
jupiter_cml(:boulder, '2024-03-15 03:00:00+00') >= 0.0
AND jupiter_cml(:boulder, '2024-03-15 03:00:00+00') < 360.0 AS in_range;
test | cml_deg | in_range
-------------+---------+----------
jupiter_cml | 306.0 | t
(1 row)
-- ============================================================
-- Test 10: Jupiter burst probability - known high probability zone
-- Source A: CML 200-260, Io phase 195-265 => p=0.8
-- ============================================================
SELECT 'burst_source_a' AS test,
jupiter_burst_probability(230.0, 230.0) AS p_source_a;
test | p_source_a
----------------+------------
burst_source_a | 0.8
(1 row)
-- ============================================================
-- Test 11: Jupiter burst probability - known zones
-- ============================================================
SELECT 'burst_zones' AS test,
jupiter_burst_probability(90.0, 150.0) AS p_source_b,
jupiter_burst_probability(240.0, 350.0) AS p_source_c,
jupiter_burst_probability(100.0, 25.0) AS p_source_d,
jupiter_burst_probability(0.0, 130.0) AS p_quiet;
test | p_source_b | p_source_c | p_source_d | p_quiet
-------------+------------+------------+------------+---------
burst_zones | 0.6 | 0.5 | 0.15 | 0
(1 row)
-- ============================================================
-- Test 12: Io phase changes over time (Io orbits in ~1.77 days)
-- Two times 12 hours apart should show significant phase change.
-- ============================================================
SELECT 'io_phase_changes' AS test,
round(io_phase_angle('2024-03-15 00:00:00+00')::numeric, 1) AS phase_0h,
round(io_phase_angle('2024-03-15 12:00:00+00')::numeric, 1) AS phase_12h,
abs(io_phase_angle('2024-03-15 00:00:00+00') -
io_phase_angle('2024-03-15 12:00:00+00')) > 10.0 AS changed;
test | phase_0h | phase_12h | changed
------------------+----------+-----------+---------
io_phase_changes | 205.0 | 103.4 | t
(1 row)
-- ============================================================
-- Test 13: Error handling - invalid Galilean moon body_id
-- ============================================================
SELECT 'invalid_galilean' AS test, galilean_observe(4, :boulder, now());
ERROR: galilean_observe: body_id 4 must be 0-3 (Io, Europa, Ganymede, Callisto)
-- ============================================================
-- Test 14: Error handling - invalid Saturn moon body_id
-- ============================================================
SELECT 'invalid_saturn' AS test, saturn_moon_observe(8, :boulder, now());
ERROR: saturn_moon_observe: body_id 8 must be 0-7 (Mimas-Hyperion)

152
test/sql/moon_observe.sql Normal file
View File

@ -0,0 +1,152 @@
-- moon_observe regression tests
--
-- Tests planetary moon observation (Galilean, Saturn, Uranus, Mars)
-- and Jupiter decametric radio burst prediction functions.
-- Reference distances from JPL/NASA fact sheets.
\set boulder '''40.015N 105.270W 1655m'''::observer
-- ============================================================
-- Test 1: Galilean moon observation - all four from Boulder
-- Io, Europa, Ganymede, Callisto should return valid topocentric.
-- ============================================================
SELECT 'galilean_io' AS test,
round(topo_azimuth(galilean_observe(0, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(galilean_observe(0, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS el_deg;
SELECT 'galilean_europa' AS test,
round(topo_azimuth(galilean_observe(1, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(galilean_observe(1, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS el_deg;
SELECT 'galilean_ganymede' AS test,
round(topo_azimuth(galilean_observe(2, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(galilean_observe(2, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS el_deg;
SELECT 'galilean_callisto' AS test,
round(topo_azimuth(galilean_observe(3, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(galilean_observe(3, :boulder,
'2024-03-15 03:00:00+00'))::numeric, 0) AS el_deg;
-- ============================================================
-- Test 2: Galilean moons should be very close to Jupiter
-- All four should have ranges within ~0.05 AU of Jupiter's range.
-- Jupiter is ~4-6 AU from Earth; moons orbit within 0.013 AU.
-- ============================================================
SELECT 'galilean_near_jupiter' AS test,
round(topo_range(planet_observe(5, :boulder,
'2024-03-15 03:00:00+00'))::numeric, -4) AS jupiter_km,
round(topo_range(galilean_observe(0, :boulder,
'2024-03-15 03:00:00+00'))::numeric, -4) AS io_km,
round(topo_range(galilean_observe(3, :boulder,
'2024-03-15 03:00:00+00'))::numeric, -4) AS callisto_km;
-- ============================================================
-- Test 3: Saturn moon observation - Titan (body_id=5)
-- Titan is Saturn's largest moon, should be near Saturn.
-- ============================================================
SELECT 'saturn_titan' AS test,
round(topo_azimuth(saturn_moon_observe(5, :boulder,
'2024-06-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(saturn_moon_observe(5, :boulder,
'2024-06-15 03:00:00+00'))::numeric, 0) AS el_deg;
-- ============================================================
-- Test 4: Saturn moons - Mimas through Iapetus all return results
-- ============================================================
SELECT 'saturn_moons' AS test,
body_id,
round(topo_range(saturn_moon_observe(body_id, :boulder,
'2024-06-15 03:00:00+00'))::numeric, -4) AS range_km
FROM generate_series(0, 7) AS body_id;
-- ============================================================
-- Test 5: Uranus moon observation - Titania (body_id=3)
-- ============================================================
SELECT 'uranus_titania' AS test,
round(topo_azimuth(uranus_moon_observe(3, :boulder,
'2024-06-15 03:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(uranus_moon_observe(3, :boulder,
'2024-06-15 03:00:00+00'))::numeric, 0) AS el_deg;
-- ============================================================
-- Test 6: All Uranus moons return valid results
-- ============================================================
SELECT 'uranus_moons' AS test,
body_id,
round(topo_range(uranus_moon_observe(body_id, :boulder,
'2024-06-15 03:00:00+00'))::numeric, -4) AS range_km
FROM generate_series(0, 4) AS body_id;
-- ============================================================
-- Test 7: Mars moons - Phobos and Deimos
-- ============================================================
SELECT 'mars_phobos' AS test,
round(topo_azimuth(mars_moon_observe(0, :boulder,
'2024-01-15 06:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(mars_moon_observe(0, :boulder,
'2024-01-15 06:00:00+00'))::numeric, 0) AS el_deg;
SELECT 'mars_deimos' AS test,
round(topo_azimuth(mars_moon_observe(1, :boulder,
'2024-01-15 06:00:00+00'))::numeric, 0) AS az_deg,
round(topo_elevation(mars_moon_observe(1, :boulder,
'2024-01-15 06:00:00+00'))::numeric, 0) AS el_deg;
-- ============================================================
-- Test 8: Io phase angle - should be [0, 360)
-- ============================================================
SELECT 'io_phase' AS test,
round(io_phase_angle('2024-03-15 03:00:00+00')::numeric, 1) AS phase_deg,
io_phase_angle('2024-03-15 03:00:00+00') >= 0.0
AND io_phase_angle('2024-03-15 03:00:00+00') < 360.0 AS in_range;
-- ============================================================
-- Test 9: Jupiter CML (System III) - should be [0, 360)
-- ============================================================
SELECT 'jupiter_cml' AS test,
round(jupiter_cml(:boulder, '2024-03-15 03:00:00+00')::numeric, 1) AS cml_deg,
jupiter_cml(:boulder, '2024-03-15 03:00:00+00') >= 0.0
AND jupiter_cml(:boulder, '2024-03-15 03:00:00+00') < 360.0 AS in_range;
-- ============================================================
-- Test 10: Jupiter burst probability - known high probability zone
-- Source A: CML 200-260, Io phase 195-265 => p=0.8
-- ============================================================
SELECT 'burst_source_a' AS test,
jupiter_burst_probability(230.0, 230.0) AS p_source_a;
-- ============================================================
-- Test 11: Jupiter burst probability - known zones
-- ============================================================
SELECT 'burst_zones' AS test,
jupiter_burst_probability(90.0, 150.0) AS p_source_b,
jupiter_burst_probability(240.0, 350.0) AS p_source_c,
jupiter_burst_probability(100.0, 25.0) AS p_source_d,
jupiter_burst_probability(0.0, 130.0) AS p_quiet;
-- ============================================================
-- Test 12: Io phase changes over time (Io orbits in ~1.77 days)
-- Two times 12 hours apart should show significant phase change.
-- ============================================================
SELECT 'io_phase_changes' AS test,
round(io_phase_angle('2024-03-15 00:00:00+00')::numeric, 1) AS phase_0h,
round(io_phase_angle('2024-03-15 12:00:00+00')::numeric, 1) AS phase_12h,
abs(io_phase_angle('2024-03-15 00:00:00+00') -
io_phase_angle('2024-03-15 12:00:00+00')) > 10.0 AS changed;
-- ============================================================
-- Test 13: Error handling - invalid Galilean moon body_id
-- ============================================================
SELECT 'invalid_galilean' AS test, galilean_observe(4, :boulder, now());
-- ============================================================
-- Test 14: Error handling - invalid Saturn moon body_id
-- ============================================================
SELECT 'invalid_saturn' AS test, saturn_moon_observe(8, :boulder, now());