Skip to main content

target_match/
matcher.rs

1//! The matching engine: input trait, constraints, ranking, and a prebuilt index.
2//!
3//! A consumer's catalogue type implements [`SkyObject`] (position only — never a
4//! name). A [`Constraint`] combines a [`Membership`] shape (circular, rectangle,
5//! or rotated rectangle) with a [`Query`] mode. [`rank`] scans a slice; [`Matcher`]
6//! builds a declination-sorted index once and answers repeated queries, returning
7//! results identical to [`rank`].
8//!
9//! # Geometry
10//!
11//! Matching precesses the pointing to J2000 (via [`skymath::precess`]), then works
12//! in the local tangent frame about the pointing: an object's offset is decomposed
13//! from its great-circle separation and position angle (East of North, both from
14//! `skymath`) into East/North components, which are rotated into the camera frame
15//! for rectangular membership. The circumscribed circle pre-filters both rectangle
16//! tests, so the tangent decomposition is only evaluated for objects near the
17//! frame — never on the far side of the sky.
18
19use core::cmp::Ordering;
20
21use skymath::{position_angle, precess, separation, Angle, Epoch, Equatorial};
22
23use crate::optics::{Field, RadiusPolicy};
24
25/// A catalogue object that can be matched by sky position.
26///
27/// The trait exposes **only** a J2000 position — matching never reads a name or
28/// designation. A caller's own type keeps its identity; a [`Match`] borrows it.
29///
30/// # Example
31///
32/// ```
33/// use skymath::{Angle, Equatorial};
34/// use target_match::SkyObject;
35///
36/// struct Target {
37///     name: &'static str,
38///     ra: f64,
39///     dec: f64,
40/// }
41/// impl SkyObject for Target {
42///     fn position(&self) -> Equatorial {
43///         Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
44///     }
45/// }
46///
47/// let m31 = Target { name: "M 31", ra: 10.6847, dec: 41.2688 };
48/// assert!((m31.position().ra().degrees() - 10.6847).abs() < 1e-9);
49/// ```
50pub trait SkyObject {
51    /// The object's J2000 equatorial position — the only thing [`rank`],
52    /// [`is_framed`], and [`Matcher`] read from an implementer.
53    ///
54    /// # Example
55    ///
56    /// ```
57    /// use skymath::{Angle, Equatorial};
58    /// use target_match::SkyObject;
59    ///
60    /// struct Target { ra_deg: f64, dec_deg: f64 }
61    /// impl SkyObject for Target {
62    ///     fn position(&self) -> Equatorial {
63    ///         Equatorial::j2000(Angle::from_degrees(self.ra_deg), Angle::from_degrees(self.dec_deg)).unwrap()
64    ///     }
65    /// }
66    ///
67    /// let m31 = Target { ra_deg: 10.6847, dec_deg: 41.2688 };
68    /// assert!((m31.position().ra().degrees() - 10.6847).abs() < 1e-9);
69    /// ```
70    fn position(&self) -> Equatorial;
71}
72
73/// The shape that decides whether an object is "in frame".
74///
75/// Usually built for you by a [`Constraint`] constructor; pass one directly to
76/// [`is_framed`] to test a single object.
77///
78/// # Example
79///
80/// ```
81/// use skymath::{Angle, Equatorial, ParseMode};
82/// use target_match::{is_framed, Membership, SkyObject};
83///
84/// struct Target {
85///     ra: f64,
86///     dec: f64,
87/// }
88/// impl SkyObject for Target {
89///     fn position(&self) -> Equatorial {
90///         Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
91///     }
92/// }
93///
94/// let m31 = Target { ra: 10.6847, dec: 41.2688 };
95/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
96///
97/// let shape = Membership::Circular { radius: Angle::from_degrees(1.0) };
98/// assert!(is_framed(pointing, &m31, shape).in_frame);
99/// ```
100#[derive(Debug, Clone, Copy, PartialEq)]
101#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
102pub enum Membership {
103    /// Pure angular distance: in frame iff separation ≤ `radius`.
104    Circular {
105        /// Search radius.
106        radius: Angle,
107    },
108    /// Axis-aligned rectangle of the given width×height field of view.
109    Rectangle {
110        /// `(width, height)` field of view.
111        fov: (Angle, Angle),
112    },
113    /// Rectangle rotated by a camera position angle (degrees East of North).
114    Rotated {
115        /// `(width, height)` field of view.
116        fov: (Angle, Angle),
117        /// Camera position angle, East of North.
118        position_angle: Angle,
119    },
120}
121
122/// What to return from a match.
123///
124/// Set on a [`Constraint`] via [`Constraint::all`], [`Constraint::nearest_one`],
125/// [`Constraint::nearest_n`], or [`Constraint::nearest_n_within`].
126///
127/// # Example
128///
129/// ```
130/// use skymath::{Angle, Equatorial, ParseMode};
131/// use target_match::{rank, Constraint, SkyObject};
132///
133/// struct Target {
134///     ra: f64,
135///     dec: f64,
136/// }
137/// impl SkyObject for Target {
138///     fn position(&self) -> Equatorial {
139///         Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
140///     }
141/// }
142///
143/// let catalog = [Target { ra: 10.6847, dec: 41.2688 }];
144/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
145///
146/// // `nearest_n(1)` sets the query to `Query::NearestN { n: 1, .. }`.
147/// let hits = rank(pointing, &catalog, Constraint::circular(Angle::from_degrees(1.0)).nearest_n(1));
148/// assert_eq!(hits.len(), 1);
149/// ```
150#[derive(Debug, Clone, Copy, PartialEq)]
151#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
152pub enum Query {
153    /// Every object inside the membership shape, ranked nearest-first.
154    AllWithinField,
155    /// The single nearest in-frame object.
156    NearestOne,
157    /// The `n` nearest objects by separation, optionally bounded by a radius.
158    NearestN {
159        /// Maximum number of results.
160        n: usize,
161        /// Optional maximum separation; unbounded when `None`.
162        max_radius: Option<Angle>,
163    },
164}
165
166/// A [`Membership`] shape combined with a [`Query`] mode (and the plate scale,
167/// when known, so pixel offsets can be reported).
168///
169/// Build one with a shape constructor ([`within`](Constraint::within),
170/// [`circular`](Constraint::circular), [`frame`](Constraint::frame),
171/// [`frame_rotated`](Constraint::frame_rotated)), then optionally switch the
172/// query mode ([`all`](Constraint::all), [`nearest_one`](Constraint::nearest_one),
173/// [`nearest_n`](Constraint::nearest_n),
174/// [`nearest_n_within`](Constraint::nearest_n_within)). Pass the result to
175/// [`rank`], [`Matcher::query`], or [`is_framed`].
176///
177/// # Example
178///
179/// ```
180/// use skymath::{Angle, Equatorial, ParseMode};
181/// use target_match::{rank, Constraint, Field, Optics, RadiusPolicy, SkyObject};
182///
183/// struct Target {
184///     ra: f64,
185///     dec: f64,
186/// }
187/// impl SkyObject for Target {
188///     fn position(&self) -> Equatorial {
189///         Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
190///     }
191/// }
192///
193/// let catalog = [Target { ra: 10.6847, dec: 41.2688 }];
194/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
195/// let field = Field::from_optics(Optics {
196///     focal_mm: 800.0,
197///     pixel_um: (3.76, 3.76),
198///     binning: (1, 1),
199///     pixels: (6248, 4176),
200/// })
201/// .unwrap();
202///
203/// let c = Constraint::within(&field, RadiusPolicy::Circumscribed).nearest_one();
204/// assert_eq!(rank(pointing, &catalog, c).len(), 1);
205/// ```
206#[derive(Debug, Clone, Copy, PartialEq)]
207#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
208pub struct Constraint {
209    /// The in-frame shape.
210    pub membership: Membership,
211    /// The query mode.
212    pub query: Query,
213    /// Per-axis plate scale (arcsec/px), used only to fill pixel offsets.
214    pub pixel_scale: Option<(f64, f64)>,
215}
216
217impl Constraint {
218    /// Circular membership with a radius derived from a field under `policy`.
219    ///
220    /// # Example
221    ///
222    /// ```
223    /// use target_match::{Constraint, Field, Optics, RadiusPolicy};
224    ///
225    /// let field = Field::from_optics(Optics {
226    ///     focal_mm: 800.0,
227    ///     pixel_um: (3.76, 3.76),
228    ///     binning: (1, 1),
229    ///     pixels: (6248, 4176),
230    /// })
231    /// .unwrap();
232    /// let c = Constraint::within(&field, RadiusPolicy::Circumscribed);
233    /// assert!(c.pixel_scale.is_some(), "the field's plate scale carries through");
234    /// ```
235    #[must_use]
236    pub fn within(field: &Field, policy: RadiusPolicy) -> Self {
237        Self {
238            membership: Membership::Circular {
239                radius: field.radius(policy),
240            },
241            query: Query::AllWithinField,
242            pixel_scale: field.pixel_scale(),
243        }
244    }
245    /// Circular membership with an explicit radius (no plate scale).
246    ///
247    /// # Example
248    ///
249    /// ```
250    /// use skymath::Angle;
251    /// use target_match::Constraint;
252    ///
253    /// let c = Constraint::circular(Angle::from_degrees(2.0));
254    /// assert!(c.pixel_scale.is_none(), "no `Field`, so no plate scale");
255    /// ```
256    #[must_use]
257    pub fn circular(radius: Angle) -> Self {
258        Self {
259            membership: Membership::Circular { radius },
260            query: Query::AllWithinField,
261            pixel_scale: None,
262        }
263    }
264    /// Axis-aligned rectangular membership from a field's width×height.
265    ///
266    /// # Example
267    ///
268    /// ```
269    /// use target_match::{Constraint, Field, Membership, Optics};
270    ///
271    /// let field = Field::from_optics(Optics {
272    ///     focal_mm: 800.0,
273    ///     pixel_um: (3.76, 3.76),
274    ///     binning: (1, 1),
275    ///     pixels: (6248, 4176),
276    /// })
277    /// .unwrap();
278    /// let c = Constraint::frame(&field);
279    /// assert!(matches!(c.membership, Membership::Rectangle { .. }));
280    /// ```
281    #[must_use]
282    pub fn frame(field: &Field) -> Self {
283        Self {
284            membership: Membership::Rectangle {
285                fov: (field.width(), field.height()),
286            },
287            query: Query::AllWithinField,
288            pixel_scale: field.pixel_scale(),
289        }
290    }
291    /// Rotated rectangular membership from a field and a camera position angle
292    /// (degrees East of North).
293    ///
294    /// # Example
295    ///
296    /// An object due north of the pointing sits on the frame's `+y` axis when
297    /// axis-aligned; a 90° East-of-North camera rotation moves it onto the `-x`
298    /// axis.
299    ///
300    /// ```
301    /// use skymath::{Angle, Equatorial, ParseMode};
302    /// use target_match::{rank, Constraint, Field, SkyObject};
303    ///
304    /// struct Target {
305    ///     ra: f64,
306    ///     dec: f64,
307    /// }
308    /// impl SkyObject for Target {
309    ///     fn position(&self) -> Equatorial {
310    ///         Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
311    ///     }
312    /// }
313    ///
314    /// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
315    /// let field = Field::from_fov(Angle::from_degrees(2.0), Angle::from_degrees(2.0)).unwrap();
316    /// let catalog = [Target { ra: 10.6847, dec: 41.2688 + 0.3 }]; // ~0.3° due north
317    ///
318    /// let (ax, ay) = rank(pointing, &catalog, Constraint::frame(&field))[0].offset.frame.unwrap();
319    /// assert!(ay.degrees() > 0.25 && ax.degrees().abs() < 0.05);
320    ///
321    /// let rotated = Constraint::frame_rotated(&field, Angle::from_degrees(90.0));
322    /// let (rx, ry) = rank(pointing, &catalog, rotated)[0].offset.frame.unwrap();
323    /// assert!(rx.degrees() < -0.25 && ry.degrees().abs() < 0.05);
324    /// ```
325    #[must_use]
326    pub fn frame_rotated(field: &Field, position_angle: Angle) -> Self {
327        Self {
328            membership: Membership::Rotated {
329                fov: (field.width(), field.height()),
330                position_angle,
331            },
332            query: Query::AllWithinField,
333            pixel_scale: field.pixel_scale(),
334        }
335    }
336    /// Set the query to all-within-field.
337    ///
338    /// # Example
339    ///
340    /// ```
341    /// use skymath::Angle;
342    /// use target_match::{Constraint, Query};
343    ///
344    /// let c = Constraint::circular(Angle::from_degrees(1.0)).nearest_one().all();
345    /// assert_eq!(c.query, Query::AllWithinField);
346    /// ```
347    #[must_use]
348    pub fn all(mut self) -> Self {
349        self.query = Query::AllWithinField;
350        self
351    }
352    /// Set the query to nearest-one.
353    ///
354    /// # Example
355    ///
356    /// ```
357    /// use skymath::Angle;
358    /// use target_match::{Constraint, Query};
359    ///
360    /// let c = Constraint::circular(Angle::from_degrees(1.0)).nearest_one();
361    /// assert_eq!(c.query, Query::NearestOne);
362    /// ```
363    #[must_use]
364    pub fn nearest_one(mut self) -> Self {
365        self.query = Query::NearestOne;
366        self
367    }
368    /// Set the query to the `n` nearest (unbounded).
369    ///
370    /// # Example
371    ///
372    /// ```
373    /// use skymath::Angle;
374    /// use target_match::{Constraint, Query};
375    ///
376    /// let c = Constraint::circular(Angle::from_degrees(1.0)).nearest_n(3);
377    /// assert_eq!(c.query, Query::NearestN { n: 3, max_radius: None });
378    /// ```
379    #[must_use]
380    pub fn nearest_n(mut self, n: usize) -> Self {
381        self.query = Query::NearestN {
382            n,
383            max_radius: None,
384        };
385        self
386    }
387    /// Set the query to the `n` nearest within `max_radius`.
388    ///
389    /// # Example
390    ///
391    /// ```
392    /// use skymath::Angle;
393    /// use target_match::{Constraint, Query};
394    ///
395    /// let radius = Angle::from_degrees(1.0);
396    /// let c = Constraint::circular(radius).nearest_n_within(3, radius);
397    /// assert_eq!(c.query, Query::NearestN { n: 3, max_radius: Some(radius) });
398    /// ```
399    #[must_use]
400    pub fn nearest_n_within(mut self, n: usize, max_radius: Angle) -> Self {
401        self.query = Query::NearestN {
402            n,
403            max_radius: Some(max_radius),
404        };
405        self
406    }
407}
408
409/// The offset of a matched object relative to the frame centre.
410///
411/// Carried on every [`Match`]. `frame` and `pixels` are only populated for
412/// rectangular [`Membership`] (see [`Constraint::frame`],
413/// [`Constraint::frame_rotated`]) — circular membership has no frame axes.
414///
415/// # Example
416///
417/// ```
418/// use skymath::{Angle, Equatorial, ParseMode};
419/// use target_match::{rank, Constraint, Field, Optics, SkyObject};
420///
421/// struct Target {
422///     ra: f64,
423///     dec: f64,
424/// }
425/// impl SkyObject for Target {
426///     fn position(&self) -> Equatorial {
427///         Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
428///     }
429/// }
430///
431/// let catalog = [Target { ra: 10.6847, dec: 41.2688 }];
432/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
433/// let field = Field::from_optics(Optics {
434///     focal_mm: 800.0,
435///     pixel_um: (3.76, 3.76),
436///     binning: (1, 1),
437///     pixels: (6248, 4176),
438/// })
439/// .unwrap();
440///
441/// let hits = rank(pointing, &catalog, Constraint::frame(&field));
442/// let offset = hits[0].offset;
443/// assert!(offset.frame.is_some(), "rectangular membership reports a frame-aligned offset");
444/// assert!(offset.pixels.is_some(), "plate scale from `Optics` fills the pixel offset");
445/// ```
446#[derive(Debug, Clone, Copy, PartialEq)]
447#[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
448pub struct Offset {
449    /// Sky-tangent offset `(East, North)` — always present.
450    pub sky: (Angle, Angle),
451    /// Frame-aligned offset `(x, y)`, present for rectangular membership.
452    pub frame: Option<(Angle, Angle)>,
453    /// Frame-aligned offset in pixels `(x, y)`, present when a plate scale is known.
454    pub pixels: Option<(f64, f64)>,
455}
456
457/// A ranked match: a borrowed catalogue object plus its computed geometry.
458///
459/// Returned by [`rank`], [`Matcher::query`], and [`is_framed`].
460///
461/// # Example
462///
463/// ```
464/// use skymath::{Angle, Equatorial, ParseMode};
465/// use target_match::{rank, Constraint, Field, Optics, RadiusPolicy, SkyObject};
466///
467/// struct Target {
468///     name: &'static str,
469///     ra: f64,
470///     dec: f64,
471/// }
472/// impl SkyObject for Target {
473///     fn position(&self) -> Equatorial {
474///         Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
475///     }
476/// }
477///
478/// let catalog = [Target { name: "M 31", ra: 10.6847, dec: 41.2688 }];
479/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
480/// let field = Field::from_optics(Optics {
481///     focal_mm: 800.0,
482///     pixel_um: (3.76, 3.76),
483///     binning: (1, 1),
484///     pixels: (6248, 4176),
485/// })
486/// .unwrap();
487///
488/// let hits = rank(pointing, &catalog, Constraint::within(&field, RadiusPolicy::Circumscribed));
489/// let m = &hits[0];
490/// assert_eq!(m.object.name, "M 31");
491/// assert!(m.in_frame);
492/// assert!(m.separation.arcseconds() < 2.0);
493/// ```
494#[derive(Debug, Clone, Copy, PartialEq)]
495pub struct Match<'a, T> {
496    /// The matched object, borrowed from the caller's slice or the [`Matcher`].
497    pub object: &'a T,
498    /// Great-circle separation from the frame centre.
499    pub separation: Angle,
500    /// Whether the object is inside the active membership shape.
501    pub in_frame: bool,
502    /// The object's offset relative to the frame centre.
503    pub offset: Offset,
504    /// Position angle from frame centre to the object (degrees East of North).
505    pub position_angle: Angle,
506}
507
508/// Tangent-frame `(East, North)` offset in radians, decomposed from an
509/// already-computed separation and position angle (the polar decomposition
510/// `skymath::tangent_offset` performs, reusing this evaluation's sep/PA
511/// instead of recomputing them).
512fn tangent_components(sep: Angle, pa: Angle) -> (f64, f64) {
513    let (s, r) = (pa.radians(), sep.radians());
514    (r * s.sin(), r * s.cos())
515}
516
517fn rotate(east: f64, north: f64, pa_rad: f64) -> (f64, f64) {
518    let (s, c) = pa_rad.sin_cos();
519    (east * c - north * s, east * s + north * c)
520}
521
522/// The circumscribed-circle radius (radians) that bounds a membership shape.
523fn bound_radius(m: Membership) -> f64 {
524    match m {
525        Membership::Circular { radius } => radius.radians(),
526        Membership::Rectangle { fov } | Membership::Rotated { fov, .. } => {
527            (fov.0.radians() / 2.0).hypot(fov.1.radians() / 2.0)
528        }
529    }
530}
531
532fn contains(sep: Angle, east: f64, north: f64, m: Membership) -> bool {
533    match m {
534        Membership::Circular { radius } => {
535            let r = radius.radians();
536            r.is_finite() && r >= 0.0 && sep.radians() <= r
537        }
538        Membership::Rectangle { fov } => within_rect(sep, east, north, fov, 0.0),
539        Membership::Rotated {
540            fov,
541            position_angle,
542        } => within_rect(sep, east, north, fov, position_angle.radians()),
543    }
544}
545
546fn within_rect(sep: Angle, east: f64, north: f64, fov: (Angle, Angle), pa: f64) -> bool {
547    let (hx, hy) = (fov.0.radians() / 2.0, fov.1.radians() / 2.0);
548    let circum = hx.hypot(hy);
549    // Circumscribed pre-filter (also rejects NaN separations).
550    if sep.radians() > circum || sep.radians().is_nan() {
551        return false;
552    }
553    let (x, y) = rotate(east, north, pa);
554    x.abs() <= hx && y.abs() <= hy
555}
556
557fn build_offset(east: f64, north: f64, m: Membership, scale: Option<(f64, f64)>) -> Offset {
558    let sky = (Angle::from_radians(east), Angle::from_radians(north));
559    match m {
560        Membership::Circular { .. } => Offset {
561            sky,
562            frame: None,
563            pixels: None,
564        },
565        Membership::Rectangle { .. } | Membership::Rotated { .. } => {
566            let pa = match m {
567                Membership::Rotated { position_angle, .. } => position_angle.radians(),
568                _ => 0.0,
569            };
570            let (x, y) = rotate(east, north, pa);
571            let frame = Some((Angle::from_radians(x), Angle::from_radians(y)));
572            let pixels = scale.map(|(sx, sy)| {
573                (
574                    Angle::from_radians(x).arcseconds() / sx,
575                    Angle::from_radians(y).arcseconds() / sy,
576                )
577            });
578            Offset { sky, frame, pixels }
579        }
580    }
581}
582
583fn evaluate<'a, T: SkyObject>(
584    pointing: Equatorial,
585    obj: &'a T,
586    m: Membership,
587    scale: Option<(f64, f64)>,
588) -> Match<'a, T> {
589    let pos = obj.position();
590    let sep = separation(pointing, pos);
591    let pa = position_angle(pointing, pos);
592    let (east, north) = tangent_components(sep, pa);
593    Match {
594        object: obj,
595        separation: sep,
596        in_frame: contains(sep, east, north, m),
597        offset: build_offset(east, north, m, scale),
598        position_angle: pa,
599    }
600}
601
602/// Whether `obj` should be kept for `query` given its evaluated match.
603fn keep<T>(m: &Match<'_, T>, query: Query) -> bool {
604    match query {
605        Query::AllWithinField | Query::NearestOne => m.in_frame,
606        Query::NearestN { max_radius, .. } => {
607            max_radius.map_or(true, |r| m.separation.radians() <= r.radians())
608        }
609    }
610}
611
612/// Shared core: evaluate candidates, filter, rank, and truncate per the query.
613fn rank_candidates<'a, T: SkyObject, I>(
614    pointing: Equatorial,
615    candidates: I,
616    c: &Constraint,
617) -> Vec<Match<'a, T>>
618where
619    I: Iterator<Item = (usize, &'a T)>,
620{
621    let mut scored: Vec<(usize, Match<'a, T>)> = candidates
622        .map(|(i, o)| (i, evaluate(pointing, o, c.membership, c.pixel_scale)))
623        .filter(|(_, m)| keep(m, c.query))
624        .collect();
625    scored.sort_by(|a, b| {
626        a.1.separation
627            .radians()
628            .partial_cmp(&b.1.separation.radians())
629            .unwrap_or(Ordering::Equal)
630            .then(a.0.cmp(&b.0))
631    });
632    let mut out: Vec<Match<'a, T>> = scored.into_iter().map(|(_, m)| m).collect();
633    match c.query {
634        Query::NearestOne => out.truncate(1),
635        Query::NearestN { n, .. } => out.truncate(n),
636        Query::AllWithinField => {}
637    }
638    out
639}
640
641/// Rank a slice of objects against a pointing under a constraint (stateless scan).
642///
643/// The pointing is precessed to J2000 first. Results are ascending by separation
644/// with ties broken by input order. For repeated queries against one catalogue,
645/// build a [`Matcher`] instead — it returns identical results faster.
646///
647/// # Example
648///
649/// ```
650/// use skymath::{Angle, Equatorial, ParseMode};
651/// use target_match::{rank, Constraint, Field, Optics, RadiusPolicy, SkyObject};
652///
653/// struct Target {
654///     name: &'static str,
655///     ra: f64,
656///     dec: f64,
657/// }
658/// impl SkyObject for Target {
659///     fn position(&self) -> Equatorial {
660///         Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
661///     }
662/// }
663///
664/// let catalog = [
665///     Target { name: "M 31", ra: 10.6847, dec: 41.2688 },
666///     Target { name: "M 33", ra: 23.4621, dec: 30.6599 },
667/// ];
668/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
669/// let field = Field::from_optics(Optics {
670///     focal_mm: 800.0,
671///     pixel_um: (3.76, 3.76),
672///     binning: (1, 1),
673///     pixels: (6248, 4176),
674/// })
675/// .unwrap();
676///
677/// let hits = rank(pointing, &catalog, Constraint::within(&field, RadiusPolicy::Circumscribed).nearest_one());
678/// assert_eq!(hits[0].object.name, "M 31");
679/// ```
680#[must_use]
681pub fn rank<T: SkyObject>(pointing: Equatorial, objects: &[T], c: Constraint) -> Vec<Match<'_, T>> {
682    let p = precess(pointing, Epoch::J2000);
683    rank_candidates(p, objects.iter().enumerate(), &c)
684}
685
686/// Evaluate a single object against a frame, returning its membership + geometry.
687///
688/// The pointing is precessed to J2000 first. Unlike [`rank`], no filtering or
689/// ranking is applied — the returned [`Match`] always describes `object`.
690///
691/// # Example
692///
693/// ```
694/// use skymath::{Angle, Equatorial, ParseMode};
695/// use target_match::{is_framed, Membership, SkyObject};
696///
697/// struct Target {
698///     ra: f64,
699///     dec: f64,
700/// }
701/// impl SkyObject for Target {
702///     fn position(&self) -> Equatorial {
703///         Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
704///     }
705/// }
706///
707/// let m31 = Target { ra: 10.6847, dec: 41.2688 };
708/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
709///
710/// let m = is_framed(pointing, &m31, Membership::Circular { radius: Angle::from_degrees(1.0) });
711/// assert!(m.in_frame);
712/// ```
713#[must_use]
714pub fn is_framed<T: SkyObject>(
715    pointing: Equatorial,
716    object: &T,
717    membership: Membership,
718) -> Match<'_, T> {
719    let p = precess(pointing, Epoch::J2000);
720    evaluate(p, object, membership, None)
721}
722
723/// A prebuilt, declination-sorted index for repeated queries against one catalogue.
724///
725/// Produces results identical to [`rank`] for the same objects, pointing, and
726/// constraint — the index is a performance optimization only. Build it once, then
727/// [`query`](Matcher::query) many pointings.
728///
729/// # Example
730///
731/// ```
732/// use skymath::{Angle, Equatorial, ParseMode};
733/// use target_match::{Constraint, Field, Matcher, Membership, Optics, RadiusPolicy, SkyObject};
734///
735/// struct Target {
736///     name: &'static str,
737///     ra: f64,
738///     dec: f64,
739/// }
740/// impl SkyObject for Target {
741///     fn position(&self) -> Equatorial {
742///         Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
743///     }
744/// }
745///
746/// let matcher = Matcher::from_objects(vec![
747///     Target { name: "M 31", ra: 10.6847, dec: 41.2688 },
748///     Target { name: "M 33", ra: 23.4621, dec: 30.6599 },
749/// ]);
750/// assert_eq!(matcher.objects().len(), 2);
751///
752/// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
753/// let field = Field::from_optics(Optics {
754///     focal_mm: 800.0,
755///     pixel_um: (3.76, 3.76),
756///     binning: (1, 1),
757///     pixels: (6248, 4176),
758/// })
759/// .unwrap();
760///
761/// let hits = matcher.query(pointing, Constraint::within(&field, RadiusPolicy::Circumscribed).nearest_one());
762/// assert_eq!(hits[0].object.name, "M 31");
763///
764/// let m = matcher.is_framed(
765///     pointing,
766///     &matcher.objects()[0],
767///     Membership::Circular { radius: Angle::from_degrees(1.0) },
768/// );
769/// assert!(m.in_frame);
770/// ```
771pub struct Matcher<T> {
772    storage: Vec<T>,
773    /// `(declination_degrees, original_index)` sorted by declination.
774    sorted: Vec<(f64, usize)>,
775}
776
777impl<T: SkyObject> Matcher<T> {
778    /// Build an index from a set of objects (original order is preserved for
779    /// tie-breaking and [`objects`](Matcher::objects)).
780    ///
781    /// # Example
782    ///
783    /// ```
784    /// use skymath::{Angle, Equatorial};
785    /// use target_match::{Matcher, SkyObject};
786    ///
787    /// struct Target { ra_deg: f64, dec_deg: f64 }
788    /// impl SkyObject for Target {
789    ///     fn position(&self) -> Equatorial {
790    ///         Equatorial::j2000(Angle::from_degrees(self.ra_deg), Angle::from_degrees(self.dec_deg)).unwrap()
791    ///     }
792    /// }
793    ///
794    /// let matcher = Matcher::from_objects(vec![
795    ///     Target { ra_deg: 10.6847, dec_deg: 41.2688 },
796    ///     Target { ra_deg: 23.4621, dec_deg: 30.6599 },
797    /// ]);
798    /// assert_eq!(matcher.objects().len(), 2);
799    /// ```
800    #[must_use]
801    pub fn from_objects(objects: Vec<T>) -> Self {
802        let mut sorted: Vec<(f64, usize)> = objects
803            .iter()
804            .enumerate()
805            .map(|(i, o)| (o.position().dec().degrees(), i))
806            .collect();
807        sorted.sort_by(|a, b| {
808            a.0.partial_cmp(&b.0)
809                .unwrap_or(Ordering::Equal)
810                .then(a.1.cmp(&b.1))
811        });
812        Self {
813            storage: objects,
814            sorted,
815        }
816    }
817
818    /// The stored objects, in their original insertion order.
819    ///
820    /// # Example
821    ///
822    /// ```
823    /// use skymath::{Angle, Equatorial};
824    /// use target_match::{Matcher, SkyObject};
825    ///
826    /// struct Target { name: &'static str, ra_deg: f64, dec_deg: f64 }
827    /// impl SkyObject for Target {
828    ///     fn position(&self) -> Equatorial {
829    ///         Equatorial::j2000(Angle::from_degrees(self.ra_deg), Angle::from_degrees(self.dec_deg)).unwrap()
830    ///     }
831    /// }
832    ///
833    /// let matcher = Matcher::from_objects(vec![Target { name: "M 31", ra_deg: 10.6847, dec_deg: 41.2688 }]);
834    /// assert_eq!(matcher.objects()[0].name, "M 31");
835    /// ```
836    #[must_use]
837    pub fn objects(&self) -> &[T] {
838        &self.storage
839    }
840
841    /// Query the index for a pointing under a constraint. Results are
842    /// identical to calling [`rank`] with [`objects`](Matcher::objects).
843    ///
844    /// # Example
845    ///
846    /// ```
847    /// use skymath::{Angle, Equatorial, ParseMode};
848    /// use target_match::{Constraint, Field, Matcher, Optics, RadiusPolicy, SkyObject};
849    ///
850    /// struct Target { name: &'static str, ra_deg: f64, dec_deg: f64 }
851    /// impl SkyObject for Target {
852    ///     fn position(&self) -> Equatorial {
853    ///         Equatorial::j2000(Angle::from_degrees(self.ra_deg), Angle::from_degrees(self.dec_deg)).unwrap()
854    ///     }
855    /// }
856    ///
857    /// let matcher = Matcher::from_objects(vec![Target { name: "M 31", ra_deg: 10.6847, dec_deg: 41.2688 }]);
858    /// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
859    /// let field = Field::from_optics(Optics {
860    ///     focal_mm: 800.0, pixel_um: (3.76, 3.76), binning: (1, 1), pixels: (6248, 4176),
861    /// })
862    /// .unwrap();
863    ///
864    /// let hits = matcher.query(pointing, Constraint::within(&field, RadiusPolicy::Circumscribed).nearest_one());
865    /// assert_eq!(hits[0].object.name, "M 31");
866    /// ```
867    #[must_use]
868    pub fn query(&self, pointing: Equatorial, c: Constraint) -> Vec<Match<'_, T>> {
869        let p = precess(pointing, Epoch::J2000);
870        let r = match c.query {
871            Query::NearestN { max_radius, .. } => max_radius.map_or(f64::INFINITY, |a| a.radians()),
872            _ => bound_radius(c.membership),
873        };
874        let idxs = self.band(p.dec().degrees(), r);
875        rank_candidates(p, idxs.into_iter().map(|i| (i, &self.storage[i])), &c)
876    }
877
878    /// Evaluate a single stored-or-external object against a frame (see [`is_framed`]).
879    ///
880    /// # Example
881    ///
882    /// ```
883    /// use skymath::{Angle, Equatorial, ParseMode};
884    /// use target_match::{Matcher, Membership, SkyObject};
885    ///
886    /// struct Target { ra_deg: f64, dec_deg: f64 }
887    /// impl SkyObject for Target {
888    ///     fn position(&self) -> Equatorial {
889    ///         Equatorial::j2000(Angle::from_degrees(self.ra_deg), Angle::from_degrees(self.dec_deg)).unwrap()
890    ///     }
891    /// }
892    ///
893    /// let matcher = Matcher::from_objects(vec![Target { ra_deg: 10.6847, dec_deg: 41.2688 }]);
894    /// let pointing = Equatorial::parse_j2000("00:42:44.3", "+41:16:09", ParseMode::Strict).unwrap();
895    ///
896    /// let obj = &matcher.objects()[0];
897    /// let m = matcher.is_framed(pointing, obj, Membership::Circular { radius: Angle::from_degrees(1.0) });
898    /// assert!(m.in_frame);
899    /// ```
900    #[must_use]
901    pub fn is_framed<'a>(
902        &self,
903        pointing: Equatorial,
904        object: &'a T,
905        m: Membership,
906    ) -> Match<'a, T> {
907        is_framed(pointing, object, m)
908    }
909
910    /// Original indices whose declination lies within `r_rad` of `dec0_deg`.
911    fn band(&self, dec0_deg: f64, r_rad: f64) -> Vec<usize> {
912        if r_rad.is_infinite() && r_rad > 0.0 {
913            return self.sorted.iter().map(|&(_, i)| i).collect();
914        }
915        if !r_rad.is_finite() || r_rad < 0.0 {
916            return Vec::new();
917        }
918        let r_deg = r_rad.to_degrees();
919        let (lo, hi) = (dec0_deg - r_deg, dec0_deg + r_deg);
920        let start = self.sorted.partition_point(|&(d, _)| d < lo);
921        let end = self.sorted.partition_point(|&(d, _)| d <= hi);
922        self.sorted[start..end].iter().map(|&(_, i)| i).collect()
923    }
924}
925
926#[cfg(test)]
927mod tests {
928    use super::*;
929
930    #[derive(Clone)]
931    struct Obj {
932        name: &'static str,
933        ra: f64,
934        dec: f64,
935    }
936    impl SkyObject for Obj {
937        fn position(&self) -> Equatorial {
938            Equatorial::j2000(Angle::from_degrees(self.ra), Angle::from_degrees(self.dec)).unwrap()
939        }
940    }
941
942    fn catalog() -> Vec<Obj> {
943        vec![
944            Obj {
945                name: "M 31",
946                ra: 10.6847,
947                dec: 41.2688,
948            },
949            Obj {
950                name: "M 110",
951                ra: 10.0921,
952                dec: 41.6853,
953            },
954            Obj {
955                name: "M 33",
956                ra: 23.4621,
957                dec: 30.6599,
958            },
959            Obj {
960                name: "M 42",
961                ra: 83.8221,
962                dec: -5.3911,
963            },
964        ]
965    }
966
967    fn m31() -> Equatorial {
968        Equatorial::j2000(Angle::from_degrees(10.6847), Angle::from_degrees(41.2688)).unwrap()
969    }
970
971    #[test]
972    fn nearest_one_circular() {
973        let cat = catalog();
974        let c = Constraint::circular(Angle::from_degrees(2.0)).nearest_one();
975        let hits = rank(m31(), &cat, c);
976        assert_eq!(hits.len(), 1);
977        assert_eq!(hits[0].object.name, "M 31");
978        assert!(hits[0].separation.arcseconds() < 1.0);
979    }
980
981    #[test]
982    fn all_within_field_ranked() {
983        let cat = catalog();
984        // ~1° radius keeps M31 (self) + M110 (~0.62°); excludes M33/M42.
985        let c = Constraint::circular(Angle::from_degrees(1.0)).all();
986        let hits = rank(m31(), &cat, c);
987        assert_eq!(
988            hits.len(),
989            2,
990            "{hits:?}",
991            hits = hits.iter().map(|h| h.object.name).collect::<Vec<_>>()
992        );
993        assert_eq!(hits[0].object.name, "M 31");
994        assert_eq!(hits[1].object.name, "M 110");
995        assert!(hits[0].separation.radians() <= hits[1].separation.radians());
996    }
997
998    #[test]
999    fn nearest_n_bounds_and_counts() {
1000        let cat = catalog();
1001        let c = Constraint::circular(Angle::from_degrees(1.0)).nearest_n(3);
1002        let hits = rank(m31(), &cat, c);
1003        assert_eq!(hits.len(), 3, "top-3 by separation regardless of frame");
1004        assert_eq!(hits[0].object.name, "M 31");
1005        // Bounded nearest-N respects the radius.
1006        let c2 = Constraint::circular(Angle::from_degrees(1.0))
1007            .nearest_n_within(3, Angle::from_degrees(1.0));
1008        assert_eq!(rank(m31(), &cat, c2).len(), 2, "only M31 + M110 within 1°");
1009    }
1010
1011    #[test]
1012    fn coordinates_only_never_name() {
1013        // A far object literally named "M 31" must NOT match near M31's pointing.
1014        let cat = vec![
1015            Obj {
1016                name: "M 31",
1017                ra: 200.0,
1018                dec: -40.0,
1019            },
1020            Obj {
1021                name: "Some Galaxy",
1022                ra: 10.6847,
1023                dec: 41.2688,
1024            },
1025        ];
1026        let c = Constraint::circular(Angle::from_degrees(2.0)).nearest_one();
1027        let hits = rank(m31(), &cat, c);
1028        assert_eq!(hits.len(), 1);
1029        assert_eq!(hits[0].object.name, "Some Galaxy");
1030    }
1031
1032    #[test]
1033    fn rectangle_excludes_circle_only_corner() {
1034        // An object just outside the rectangle but inside the circumscribed circle.
1035        // Field 2°×1°: half-width 1°, half-height 0.5°. Put an object 0.9° North
1036        // (in frame) vs 0.9° along the diagonal (out of the axis-aligned rect).
1037        let field = Field::from_fov(Angle::from_degrees(2.0), Angle::from_degrees(1.0)).unwrap();
1038        let north_obj = Obj {
1039            name: "N",
1040            ra: 10.6847,
1041            dec: 41.2688 + 0.4,
1042        }; // within 0.5° height
1043        let high_obj = Obj {
1044            name: "H",
1045            ra: 10.6847,
1046            dec: 41.2688 + 0.9,
1047        }; // beyond height, inside circle
1048        let cat = vec![north_obj, high_obj];
1049        let c = Constraint::frame(&field).all();
1050        let hits = rank(m31(), &cat, c);
1051        assert_eq!(hits.len(), 1);
1052        assert_eq!(hits[0].object.name, "N");
1053    }
1054
1055    #[test]
1056    fn matcher_matches_rank_exactly() {
1057        let cat = catalog();
1058        let c = Constraint::circular(Angle::from_degrees(5.0)).all();
1059        let via_rank: Vec<_> = rank(m31(), &cat, c).iter().map(|m| m.object.name).collect();
1060        let matcher = Matcher::from_objects(cat.clone());
1061        let via_index: Vec<_> = matcher
1062            .query(m31(), c)
1063            .iter()
1064            .map(|m| m.object.name)
1065            .collect();
1066        assert_eq!(via_rank, via_index);
1067        assert_eq!(matcher.objects().len(), 4);
1068    }
1069
1070    #[test]
1071    fn is_framed_reports_geometry() {
1072        let m110 = catalog()[1].clone();
1073        let m = is_framed(
1074            m31(),
1075            &m110,
1076            Membership::Circular {
1077                radius: Angle::from_degrees(1.0),
1078            },
1079        );
1080        assert!(m.in_frame);
1081        assert!((0.4..0.9).contains(&m.separation.degrees()));
1082    }
1083
1084    #[test]
1085    fn empty_catalog_and_zero_radius() {
1086        let cat = catalog();
1087        assert!(rank(
1088            m31(),
1089            &[] as &[Obj],
1090            Constraint::circular(Angle::from_degrees(1.0))
1091        )
1092        .is_empty());
1093        let c = Constraint::circular(Angle::from_degrees(-1.0)).all();
1094        assert!(
1095            rank(m31(), &cat, c).is_empty(),
1096            "negative radius matches nothing"
1097        );
1098    }
1099}