use core::cmp::Ordering;
use skymath::{position_angle, precess, separation, Angle, Epoch, Equatorial};
use crate::optics::{Field, RadiusPolicy};
pub trait SkyObject {
fn position(&self) -> Equatorial;
}
#[derive(Debug, Clone, Copy, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub enum Membership {
Circular {
radius: Angle,
},
Rectangle {
fov: (Angle, Angle),
},
Rotated {
fov: (Angle, Angle),
position_angle: Angle,
},
}
#[derive(Debug, Clone, Copy, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub enum Query {
AllWithinField,
NearestOne,
NearestN {
n: usize,
max_radius: Option<Angle>,
},
}
#[derive(Debug, Clone, Copy, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct Constraint {
pub membership: Membership,
pub query: Query,
pub pixel_scale: Option<(f64, f64)>,
}
impl Constraint {
#[must_use]
pub fn within(field: &Field, policy: RadiusPolicy) -> Self {
Self {
membership: Membership::Circular {
radius: field.radius(policy),
},
query: Query::AllWithinField,
pixel_scale: field.pixel_scale(),
}
}
#[must_use]
pub fn circular(radius: Angle) -> Self {
Self {
membership: Membership::Circular { radius },
query: Query::AllWithinField,
pixel_scale: None,
}
}
#[must_use]
pub fn frame(field: &Field) -> Self {
Self {
membership: Membership::Rectangle {
fov: (field.width(), field.height()),
},
query: Query::AllWithinField,
pixel_scale: field.pixel_scale(),
}
}
#[must_use]
pub fn frame_rotated(field: &Field, position_angle: Angle) -> Self {
Self {
membership: Membership::Rotated {
fov: (field.width(), field.height()),
position_angle,
},
query: Query::AllWithinField,
pixel_scale: field.pixel_scale(),
}
}
#[must_use]
pub fn all(mut self) -> Self {
self.query = Query::AllWithinField;
self
}
#[must_use]
pub fn nearest_one(mut self) -> Self {
self.query = Query::NearestOne;
self
}
#[must_use]
pub fn nearest_n(mut self, n: usize) -> Self {
self.query = Query::NearestN {
n,
max_radius: None,
};
self
}
#[must_use]
pub fn nearest_n_within(mut self, n: usize, max_radius: Angle) -> Self {
self.query = Query::NearestN {
n,
max_radius: Some(max_radius),
};
self
}
}
#[derive(Debug, Clone, Copy, PartialEq)]
#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
pub struct Offset {
pub sky: (Angle, Angle),
pub frame: Option<(Angle, Angle)>,
pub pixels: Option<(f64, f64)>,
}
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct Match<'a, T> {
pub object: &'a T,
pub separation: Angle,
pub in_frame: bool,
pub offset: Offset,
pub position_angle: Angle,
}
fn tangent_components(sep: Angle, pa: Angle) -> (f64, f64) {
let (s, r) = (pa.radians(), sep.radians());
(r * s.sin(), r * s.cos())
}
fn rotate(east: f64, north: f64, pa_rad: f64) -> (f64, f64) {
let (s, c) = pa_rad.sin_cos();
(east * c - north * s, east * s + north * c)
}
fn bound_radius(m: Membership) -> f64 {
match m {
Membership::Circular { radius } => radius.radians(),
Membership::Rectangle { fov } | Membership::Rotated { fov, .. } => {
(fov.0.radians() / 2.0).hypot(fov.1.radians() / 2.0)
}
}
}
fn contains(sep: Angle, east: f64, north: f64, m: Membership) -> bool {
match m {
Membership::Circular { radius } => {
let r = radius.radians();
r.is_finite() && r >= 0.0 && sep.radians() <= r
}
Membership::Rectangle { fov } => within_rect(sep, east, north, fov, 0.0),
Membership::Rotated {
fov,
position_angle,
} => within_rect(sep, east, north, fov, position_angle.radians()),
}
}
fn within_rect(sep: Angle, east: f64, north: f64, fov: (Angle, Angle), pa: f64) -> bool {
let (hx, hy) = (fov.0.radians() / 2.0, fov.1.radians() / 2.0);
let circum = hx.hypot(hy);
if sep.radians() > circum || sep.radians().is_nan() {
return false;
}
let (x, y) = rotate(east, north, pa);
x.abs() <= hx && y.abs() <= hy
}
fn build_offset(east: f64, north: f64, m: Membership, scale: Option<(f64, f64)>) -> Offset {
let sky = (Angle::from_radians(east), Angle::from_radians(north));
match m {
Membership::Circular { .. } => Offset {
sky,
frame: None,
pixels: None,
},
Membership::Rectangle { .. } | Membership::Rotated { .. } => {
let pa = match m {
Membership::Rotated { position_angle, .. } => position_angle.radians(),
_ => 0.0,
};
let (x, y) = rotate(east, north, pa);
let frame = Some((Angle::from_radians(x), Angle::from_radians(y)));
let pixels = scale.map(|(sx, sy)| {
(
Angle::from_radians(x).arcseconds() / sx,
Angle::from_radians(y).arcseconds() / sy,
)
});
Offset { sky, frame, pixels }
}
}
}
fn evaluate<'a, T: SkyObject>(
pointing: Equatorial,
obj: &'a T,
m: Membership,
scale: Option<(f64, f64)>,
) -> Match<'a, T> {
let pos = obj.position();
let sep = separation(pointing, pos);
let pa = position_angle(pointing, pos);
let (east, north) = tangent_components(sep, pa);
Match {
object: obj,
separation: sep,
in_frame: contains(sep, east, north, m),
offset: build_offset(east, north, m, scale),
position_angle: pa,
}
}
fn keep<T>(m: &Match<'_, T>, query: Query) -> bool {
match query {
Query::AllWithinField | Query::NearestOne => m.in_frame,
Query::NearestN { max_radius, .. } => {
max_radius.map_or(true, |r| m.separation.radians() <= r.radians())
}
}
}
fn rank_candidates<'a, T: SkyObject, I>(
pointing: Equatorial,
candidates: I,
c: &Constraint,
) -> Vec<Match<'a, T>>
where
I: Iterator<Item = (usize, &'a T)>,
{
let mut scored: Vec<(usize, Match<'a, T>)> = candidates
.map(|(i, o)| (i, evaluate(pointing, o, c.membership, c.pixel_scale)))
.filter(|(_, m)| keep(m, c.query))
.collect();
scored.sort_by(|a, b| {
a.1.separation
.radians()
.partial_cmp(&b.1.separation.radians())
.unwrap_or(Ordering::Equal)
.then(a.0.cmp(&b.0))
});
let mut out: Vec<Match<'a, T>> = scored.into_iter().map(|(_, m)| m).collect();
match c.query {
Query::NearestOne => out.truncate(1),
Query::NearestN { n, .. } => out.truncate(n),
Query::AllWithinField => {}
}
out
}
#[must_use]
pub fn rank<T: SkyObject>(pointing: Equatorial, objects: &[T], c: Constraint) -> Vec<Match<'_, T>> {
let p = precess(pointing, Epoch::J2000);
rank_candidates(p, objects.iter().enumerate(), &c)
}
#[must_use]
pub fn is_framed<T: SkyObject>(
pointing: Equatorial,
object: &T,
membership: Membership,
) -> Match<'_, T> {
let p = precess(pointing, Epoch::J2000);
evaluate(p, object, membership, None)
}
pub struct Matcher<T> {
storage: Vec<T>,
sorted: Vec<(f64, usize)>,
}
impl<T: SkyObject> Matcher<T> {
#[must_use]
pub fn from_objects(objects: Vec<T>) -> Self {
let mut sorted: Vec<(f64, usize)> = objects
.iter()
.enumerate()
.map(|(i, o)| (o.position().dec().degrees(), i))
.collect();
sorted.sort_by(|a, b| {
a.0.partial_cmp(&b.0)
.unwrap_or(Ordering::Equal)
.then(a.1.cmp(&b.1))
});
Self {
storage: objects,
sorted,
}
}
#[must_use]
pub fn objects(&self) -> &[T] {
&self.storage
}
#[must_use]
pub fn query(&self, pointing: Equatorial, c: Constraint) -> Vec<Match<'_, T>> {
let p = precess(pointing, Epoch::J2000);
let r = match c.query {
Query::NearestN { max_radius, .. } => max_radius.map_or(f64::INFINITY, |a| a.radians()),
_ => bound_radius(c.membership),
};
let idxs = self.band(p.dec().degrees(), r);
rank_candidates(p, idxs.into_iter().map(|i| (i, &self.storage[i])), &c)
}
#[must_use]
pub fn is_framed<'a>(
&self,
pointing: Equatorial,
object: &'a T,
m: Membership,
) -> Match<'a, T> {
is_framed(pointing, object, m)
}
fn band(&self, dec0_deg: f64, r_rad: f64) -> Vec<usize> {
if r_rad.is_infinite() && r_rad > 0.0 {
return self.sorted.iter().map(|&(_, i)| i).collect();
}
if !r_rad.is_finite() || r_rad < 0.0 {
return Vec::new();
}
let r_deg = r_rad.to_degrees();
let (lo, hi) = (dec0_deg - r_deg, dec0_deg + r_deg);
let start = self.sorted.partition_point(|&(d, _)| d < lo);
let end = self.sorted.partition_point(|&(d, _)| d <= hi);
self.sorted[start..end].iter().map(|&(_, i)| i).collect()
}
}
#[cfg(test)]
mod tests {
use super::*;
#[derive(Clone)]
struct Obj {
name: &'static str,
ra: f64,
dec: f64,
}
impl SkyObject for Obj {
fn position(&self) -> Equatorial {
Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
}
}
fn catalog() -> Vec<Obj> {
vec![
Obj {
name: "M 31",
ra: 10.6847,
dec: 41.2688,
},
Obj {
name: "M 110",
ra: 10.0921,
dec: 41.6853,
},
Obj {
name: "M 33",
ra: 23.4621,
dec: 30.6599,
},
Obj {
name: "M 42",
ra: 83.8221,
dec: -5.3911,
},
]
}
fn m31() -> Equatorial {
Equatorial::j2000(Angle::from_degrees(10.6847), Angle::from_degrees(41.2688)).unwrap()
}
#[test]
fn nearest_one_circular() {
let cat = catalog();
let c = Constraint::circular(Angle::from_degrees(2.0)).nearest_one();
let hits = rank(m31(), &cat, c);
assert_eq!(hits.len(), 1);
assert_eq!(hits[0].object.name, "M 31");
assert!(hits[0].separation.arcseconds() < 1.0);
}
#[test]
fn all_within_field_ranked() {
let cat = catalog();
let c = Constraint::circular(Angle::from_degrees(1.0)).all();
let hits = rank(m31(), &cat, c);
assert_eq!(
hits.len(),
2,
"{hits:?}",
hits = hits.iter().map(|h| h.object.name).collect::<Vec<_>>()
);
assert_eq!(hits[0].object.name, "M 31");
assert_eq!(hits[1].object.name, "M 110");
assert!(hits[0].separation.radians() <= hits[1].separation.radians());
}
#[test]
fn nearest_n_bounds_and_counts() {
let cat = catalog();
let c = Constraint::circular(Angle::from_degrees(1.0)).nearest_n(3);
let hits = rank(m31(), &cat, c);
assert_eq!(hits.len(), 3, "top-3 by separation regardless of frame");
assert_eq!(hits[0].object.name, "M 31");
let c2 = Constraint::circular(Angle::from_degrees(1.0))
.nearest_n_within(3, Angle::from_degrees(1.0));
assert_eq!(rank(m31(), &cat, c2).len(), 2, "only M31 + M110 within 1°");
}
#[test]
fn coordinates_only_never_name() {
let cat = vec![
Obj {
name: "M 31",
ra: 200.0,
dec: -40.0,
},
Obj {
name: "Some Galaxy",
ra: 10.6847,
dec: 41.2688,
},
];
let c = Constraint::circular(Angle::from_degrees(2.0)).nearest_one();
let hits = rank(m31(), &cat, c);
assert_eq!(hits.len(), 1);
assert_eq!(hits[0].object.name, "Some Galaxy");
}
#[test]
fn rectangle_excludes_circle_only_corner() {
let field = Field::from_fov(Angle::from_degrees(2.0), Angle::from_degrees(1.0)).unwrap();
let north_obj = Obj {
name: "N",
ra: 10.6847,
dec: 41.2688 + 0.4,
}; let high_obj = Obj {
name: "H",
ra: 10.6847,
dec: 41.2688 + 0.9,
}; let cat = vec![north_obj, high_obj];
let c = Constraint::frame(&field).all();
let hits = rank(m31(), &cat, c);
assert_eq!(hits.len(), 1);
assert_eq!(hits[0].object.name, "N");
}
#[test]
fn matcher_matches_rank_exactly() {
let cat = catalog();
let c = Constraint::circular(Angle::from_degrees(5.0)).all();
let via_rank: Vec<_> = rank(m31(), &cat, c).iter().map(|m| m.object.name).collect();
let matcher = Matcher::from_objects(cat.clone());
let via_index: Vec<_> = matcher
.query(m31(), c)
.iter()
.map(|m| m.object.name)
.collect();
assert_eq!(via_rank, via_index);
assert_eq!(matcher.objects().len(), 4);
}
#[test]
fn is_framed_reports_geometry() {
let m110 = catalog()[1].clone();
let m = is_framed(
m31(),
&m110,
Membership::Circular {
radius: Angle::from_degrees(1.0),
},
);
assert!(m.in_frame);
assert!((0.4..0.9).contains(&m.separation.degrees()));
}
#[test]
fn empty_catalog_and_zero_radius() {
let cat = catalog();
assert!(rank(
m31(),
&[] as &[Obj],
Constraint::circular(Angle::from_degrees(1.0))
)
.is_empty());
let c = Constraint::circular(Angle::from_degrees(-1.0)).all();
assert!(
rank(m31(), &cat, c).is_empty(),
"negative radius matches nothing"
);
}
}