tetra3 0.5.1

Rust implementation of Tetra3: Fast and robust star plate solver
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122

use rkyv::{Archive, Deserialize, Serialize};

/// A "Generic" star type that will be used for star matching
/// The star RA & Dec assume proper motion has already been applied to the observation epoch.
/// The magnitude is a generic "brightness" value that can be used for filtering, but the exact meaning is catalog-dependent.
#[derive(Debug, Clone, PartialEq, Archive, Serialize, Deserialize)]
pub struct Star {
    pub id: i64,
    pub ra_rad: f32,
    pub dec_rad: f32,
    pub mag: f32,
}

impl Star {
    /// Unit vector pointing to the star's position on the celestial sphere.
    pub fn uvec(&self) -> numeris::Vector3<f32> {
        let ra = self.ra_rad;
        let dec = self.dec_rad;
        // fast cosine, sine at once:
        let (rasin, racos) = ra.sin_cos();
        let (decsin, deccos) = dec.sin_cos();
        numeris::Vector3::from_array([deccos * racos, deccos * rasin, decsin])
    }
}

/// Convert a Hipparcos star to a generic Star, optionally propagating proper motion.
///
/// `epoch_year`: Target year for proper motion propagation (e.g. 2025.0).
/// If None, the catalog position at the Hipparcos reference epoch (J1991.25) is used.
///
/// Proper motion near the celestial poles (|dec| > ~87°) is ignored because
/// the cos(dec) divisor becomes numerically unstable, following the same
/// convention as tetra3/cedar-solve.
#[cfg(feature = "hipparcos")]
pub fn star_from_hipparcos(
    star: &crate::catalogs::hipparcos::HipparcosStar,
    epoch_year: Option<f64>,
) -> Star {
    // Hipparcos reference epoch is J1991.25
    const HIPPARCOS_EPOCH_YEAR: f64 = 1991.25;
    // Convert milliarcseconds/year to radians/year
    const MAS_PER_YR_TO_RAD_PER_YR: f64 =
        2.0 * std::f64::consts::PI / (3600.0 * 1000.0 * 360.0);

    let (ra, dec) = if let Some(target_year) = epoch_year {
        let dt_years = target_year - HIPPARCOS_EPOCH_YEAR;
        let cos_dec = star.dec_rad.cos();

        // Near the poles, cos(dec) -> 0 makes the RA correction blow up.
        // Following tetra3: skip proper motion for |dec| > ~87 degrees.
        let (mu_ra, mu_dec) = if cos_dec.abs() > 0.05 {
            let mu_alpha_cos_delta = star.pm_ra * MAS_PER_YR_TO_RAD_PER_YR;
            let mu_delta = star.pm_dec * MAS_PER_YR_TO_RAD_PER_YR;
            // pm_ra from Hipparcos is mu_alpha*cos(delta), so divide by cos(dec)
            (mu_alpha_cos_delta / cos_dec, mu_delta)
        } else {
            (0.0, 0.0)
        };

        let ra = star.ra_rad + mu_ra * dt_years;
        let dec = star.dec_rad + mu_dec * dt_years;
        (ra as f32, dec as f32)
    } else {
        (star.ra_rad as f32, star.dec_rad as f32)
    };

    Star {
        id: star.hip as i64,
        ra_rad: ra,
        dec_rad: dec,
        mag: star.hp_to_v(),
    }
}

/// Convert a Gaia star to a generic Star, optionally propagating proper motion.
///
/// `epoch_year`: Target year for proper motion propagation (e.g. 2025.0).
/// If None, the catalog position at the Gaia DR3 reference epoch (J2016.0) is used.
///
/// The source_id is stored directly as `i64` — negative values indicate Hipparcos
/// gap-fill stars from the merged catalog.
pub fn star_from_gaia(
    star: &crate::catalogs::gaia::GaiaStar,
    epoch_year: Option<f64>,
) -> Star {
    // Gaia DR3 reference epoch is J2016.0
    const GAIA_EPOCH_YEAR: f64 = 2016.0;
    // Gaia proper motions are in mas/yr, same conversion factor
    const MAS_PER_YR_TO_RAD_PER_YR: f64 =
        2.0 * std::f64::consts::PI / (3600.0 * 1000.0 * 360.0);

    let ra_rad = (star.ra_deg as f64).to_radians();
    let dec_rad = (star.dec_deg as f64).to_radians();

    let (ra, dec) = if let Some(target_year) = epoch_year {
        let dt_years = target_year - GAIA_EPOCH_YEAR;
        let cos_dec = dec_rad.cos();

        let (mu_ra, mu_dec) = match (star.pmra, star.pmdec) {
            (Some(pmra), Some(pmdec)) if cos_dec.abs() > 0.05 => {
                // Gaia pmra is mu_alpha*cos(delta), same convention as Hipparcos
                let mu_alpha_cos_delta = pmra as f64 * MAS_PER_YR_TO_RAD_PER_YR;
                let mu_delta = pmdec as f64 * MAS_PER_YR_TO_RAD_PER_YR;
                (mu_alpha_cos_delta / cos_dec, mu_delta)
            }
            _ => (0.0, 0.0),
        };

        let ra = ra_rad + mu_ra * dt_years;
        let dec = dec_rad + mu_dec * dt_years;
        (ra as f32, dec as f32)
    } else {
        (ra_rad as f32, dec_rad as f32)
    };

    Star {
        id: star.source_id,
        ra_rad: ra,
        dec_rad: dec,
        mag: star.phot_g_mean_mag,
    }
}