05_target_tracking.cpp

Target hierarchy, ephemeris snapshots, Kepler propagation, and transforms.

// SPDX-License-Identifier: AGPL-3.0-or-later
// Copyright (C) 2026 Vallés Puig, Ramon
 
 
#include <siderust/siderust.hpp>
 
#include <cmath>
#include <iomanip>
#include <iostream>
#include <string>
 
using namespace siderust;
using namespace siderust::frames;
using namespace siderust::centers;
using namespace qtty::literals;
 
// ─── Helper: simple coordinate snapshot (mirrors Rust's Target<T>) ──────────
 
template <typename P> struct Snapshot {
  P position;
  Time<TT, JD> time;
 
  void update(P new_pos, Time<TT, JD> new_time) {
    position = new_pos;
    time = new_time;
  }
};
 
// ─── Halley's comet orbit (hardcoded from the Rust `HALLEY` constant) ───────
 
inline KeplerianOrbit halley_orbit() {
  // a = 17.834 AU, e = 0.96714, i = 162.26°, Ω = 58.42°, ω = 111.33°,
  // M = 38.38° at epoch JD 2446467.4 (≈1986 Feb 9).
  return {17.834_au, Eccentricity{0.96714},  162.26_deg, 58.42_deg, 111.33_deg,
          38.38_deg, Time<TT, JD>(2446467.4)};
}
 
// ─── Section 1: Trackable objects ───────────────────────────────────────────
 
void section_trackable_objects(const Time<TT, JD> &jd, const Time<TT, JD> &jd_next) {
  std::cout << "1) Trackable objects (ICRS, star, Sun, planet, Moon)\n";
 
  // ICRS direction — time-invariant target
  spherical::direction::ICRS fixed_icrs(120.0_deg, 22.5_deg);
  ICRSTarget icrs_target(fixed_icrs, jd, "FixedICRS");
 
  // Verify time-invariance: the ICRS direction coordinates are constant.
  std::cout << "  ICRS direction is time-invariant: " << std::fixed << std::setprecision(3)
            << fixed_icrs << std::endl;
 
  // Sirius via the catalog StarTarget
  StarTarget sirius_target(SIRIUS());
  std::cout << "  Sirius via StarTarget: name = " << sirius_target.name() << '\n';
 
  // Sun, Mars, Moon via BodyTarget
  BodyTarget sun_target(Body::Sun);
  BodyTarget mars_target(Body::Mars);
  BodyTarget moon_target(Body::Moon);
 
  // Show Sun barycentric distance at J2000
  auto sun_bary = ephemeris::sun_barycentric(jd);
  std::cout << "  Sun barycentric distance: " << std::fixed << std::setprecision(6)
            << sun_bary.distance() << std::endl;
 
  // Mars heliocentric distance
  auto mars_helio = ephemeris::mars_heliocentric(jd);
  std::cout << "  Mars heliocentric distance: " << std::fixed << std::setprecision(6)
            << mars_helio.distance() << std::endl;
 
  // Moon geocentric distance
  auto moon_geo = ephemeris::moon_geocentric(jd);
  std::cout << "  Moon geocentric distance: " << std::fixed << std::setprecision(1)
            << moon_geo.distance() << "\n"
            << std::endl;
}
 
// ─── Section 2: Target snapshots ────────────────────────────────────────────
 
void section_target_snapshots(const Time<TT, JD> &jd, const Time<TT, JD> &jd_next) {
  std::cout << "2) Target snapshots for arbitrary sky objects\n";
 
  // Mars heliocentric snapshot (VSOP87 ephemeris)
  Snapshot<cartesian::position::EclipticMeanJ2000<qtty::AstronomicalUnit>> mars_snap{
      ephemeris::mars_heliocentric(jd), jd};
  std::cout << "  Mars target at JD " << std::fixed << std::setprecision(1) << mars_snap.time
            << ": r = " << std::setprecision(6) << mars_snap.position.distance() << std::endl;
 
  // Update with next-day ephemeris
  mars_snap.update(ephemeris::mars_heliocentric(jd_next), jd_next);
  std::cout << "  Mars target updated to JD " << std::fixed << std::setprecision(1)
            << mars_snap.time << ": r = " << std::setprecision(6) << mars_snap.position.distance()
            << std::endl;
 
  // Halley's comet — Kepler-propagated snapshot
  auto halley_pos = kepler_position(halley_orbit(), jd);
  Snapshot<cartesian::position::EclipticMeanJ2000<qtty::AstronomicalUnit>> halley_snap{halley_pos,
                                                                                       jd};
  std::cout << "  Halley target at JD " << std::fixed << std::setprecision(1) << halley_snap.time
            << ": r = " << std::setprecision(6) << halley_snap.position.distance() << std::endl;
 
  // DemoSat — satellite-like custom object with a geocentric orbit
  KeplerianOrbit demosat_orbit{
      1.0002_au, Eccentricity{0.001}, 0.1_deg, 35.0_deg, 80.0_deg, 10.0_deg, jd};
  auto demosat_pos = kepler_position(demosat_orbit, jd);
  Snapshot<cartesian::position::EclipticMeanJ2000<qtty::AstronomicalUnit>> demosat_snap{demosat_pos,
                                                                                        jd};
  std::cout << "  DemoSat target at JD " << std::fixed << std::setprecision(1) << demosat_snap.time
            << ": r = " << std::setprecision(6) << demosat_snap.position.distance() << "\n"
            << std::endl;
}
 
// ─── Section 3: Proper motion ───────────────────────────────────────────────
 
inline spherical::direction::ICRS apply_proper_motion(const spherical::direction::ICRS &pos,
                                                      const ProperMotion &pm,
                                                      const Time<TT, JD> &epoch,
                                                      const Time<TT, JD> &target_epoch) {
  constexpr double JULIAN_YEAR = 365.25; // days
  double dt_years = (target_epoch.value() - epoch.value()) / JULIAN_YEAR;
 
  double ra_deg = pos.ra().value();
  double dec_deg = pos.dec().value();
 
  // ProperMotion rates in deg/yr.
  // MuAlphaStar convention: pm_ra is already μα* = μα·cos(δ), so divide by
  // cos(δ)
  double cos_dec = std::cos(dec_deg * constants::pi / 180.0);
  double dra = (cos_dec > 1e-12) ? pm.ra.deg_per_day() * 365.25 * dt_years / cos_dec : 0.0;
  double ddec = pm.dec.deg_per_day() * 365.25 * dt_years;
 
  return spherical::direction::ICRS(qtty::Degree(ra_deg + dra), qtty::Degree(dec_deg + ddec));
}
 
void section_target_with_proper_motion(const Time<TT, JD> &jd) {
  std::cout << "3) Target with proper motion (stellar-style target)\n";
 
  // Betelgeuse approximate ICRS coordinates at J2000
  // (RA ≈ 88.7929°, Dec ≈ +7.4071°)
  spherical::direction::ICRS betelgeuse_pos(88.7929_deg, 7.4071_deg);
 
  // Proper motion: µα* = 27.54 mas/yr, µδ = 10.86 mas/yr
  // Convert mas/yr → deg/yr
  constexpr double MAS_TO_DEG = 1.0 / 3600000.0;
  ProperMotion pm{AngularRate{qtty::Degree(27.54 * MAS_TO_DEG), qtty::Day(365.25)},
                  AngularRate{qtty::Degree(10.86 * MAS_TO_DEG), qtty::Day(365.25)},
                  RaConvention::MuAlphaStar};
 
  std::cout << "  Betelgeuse-like target at J2000: " << betelgeuse_pos << '\n';
 
  // Propagate 25 years
  constexpr double JULIAN_YEAR = 365.25;
  auto jd_future = jd + qtty::Day(25.0 * JULIAN_YEAR);
  auto moved = apply_proper_motion(betelgeuse_pos, pm, jd, jd_future);
 
  std::cout << "  After 25 years: " << moved << "\n\n";
}
 
// ─── Section 4: Frame + center transforms ───────────────────────────────────
 
void section_target_transform(const Time<TT, JD> &jd) {
  std::cout << "4) Target conversion across frame + center\n";
 
  // Mars heliocentric ecliptic → geocentric equatorial
  auto mars_helio = ephemeris::mars_heliocentric(jd);
  auto mars_geoeq = mars_helio.template transform<Geocentric, EquatorialMeanJ2000>(jd);
 
  std::cout << "  Mars heliocentric ecliptic target: r = " << std::fixed << std::setprecision(6)
            << mars_helio.distance() << std::endl;
  std::cout << "  Mars geocentric equatorial target: r = " << std::fixed << std::setprecision(6)
            << mars_geoeq.distance() << std::endl;
}
 
// ─── main
// ─────────────────────────────────────────────────────────────────────
 
int main() {
  auto jd = Time<TT, JD>::J2000();
  auto jd_next = jd + qtty::Day(1.0);
 
  std::cout << "Target + Trackable examples\n"
            << "===========================\n\n";
 
  section_trackable_objects(jd, jd_next);
  section_target_snapshots(jd, jd_next);
  section_target_with_proper_motion(jd);
  section_target_transform(jd);
 
  return 0;
}