outfit 2.1.0

/// Internal representation of the JPL Horizon ephemeris binary file.
///
/// The Horizons ephemeris files (e.g. DE440) are structured as a sequence of blocks.
/// Each block covers a fixed time span and contains Chebyshev polynomial coefficients
/// used to interpolate the position and velocity of Solar System bodies.
///
/// This module provides:
/// - [`HorizonHeader`] describing the global properties of the file (time coverage, record size,
///   Earth–Moon mass ratio, etc.),
/// - [`HorizonRecords`], a collection of parsed data blocks indexed by body identifier,
/// - [`HorizonData`], the main container holding both header and records.
///
/// Bodies are indexed with IDs from 0 to 14 (consistent with the JPL Horizon conventions).
/// For each body, the ephemeris data is represented as a sequence of [`HorizonRecord`] entries,
/// each of which contains:
/// - the validity interval (start/end Julian Date),
/// - the Chebyshev coefficients for interpolation.
///
/// Typical usage flow:
/// 1. Load a binary ephemeris file (see [`EphemFilePath`]),
/// 2. Parse its header to extract metadata and IPT table,
/// 3. Read and store all blocks into [`HorizonRecords`],
/// 4. Perform interpolation using the stored coefficients via [`InterpResult`].
///
/// See also
/// --------
/// * [`JPLHorizonVersion`] – versioning information for Horizons files.
/// * [`HorizonID`] – mapping between integer body IDs and their symbolic identifiers.
/// * [`HorizonRecord`] – single record containing Chebyshev coefficients.
use nom::{
    bytes::complete::take,
    multi::count,
    number::complete::{le_f64, le_i32, le_u32},
    IResult, Parser,
};

use crate::{constants::MJD, jpl_ephem::download_jpl_file::EphemFilePath};

use super::{
    horizon_ids::HorizonID, horizon_records::HorizonRecord, horizon_version::JPLHorizonVersion,
    interpolation_result::InterpResult,
};
use std::{
    collections::HashMap,
    fs::File,
    io::{BufReader, Read, Seek},
};

/// Collection of parsed Horizons ephemeris records.
///
/// The binary ephemeris is organized as a sequence of time-ordered **blocks**.
/// For each block (outer `Vec` index), the data for every supported body is
/// grouped in a `HashMap` keyed by its Horizons integer ID (`0..=14`).
///
/// Concretely, the structure is:
/// - `Vec< … >` — chronological list of blocks,
/// - `HashMap<u8, Vec<HorizonRecord>>` — within each block, ephemeris series per body,
/// - `Vec<HorizonRecord>` — one or more Chebyshev segments covering sub-intervals
///   inside the block.
///
/// Notes
/// -----
/// * This layout preserves fast sequential access by time while keeping per‑body
///   lookups ergonomic.
/// * `HorizonRecord` holds the validity interval and Chebyshev coefficients needed
///   for interpolation.
///
/// See also
/// -----------
/// * [`HorizonRecord`] – single time segment with Chebyshev coefficients.
/// * [`HorizonHeader`] – global file metadata (time coverage, record size, IPT).
type HorizonRecords = Vec<HashMap<u8, Vec<HorizonRecord>>>;

/// IPT table describing how coefficients for each body are laid out inside a block.
///
/// `IPT[b] == [offset, n_coeff, n_sub]` where:
/// - `offset` — byte (or word) offset to the body’s coefficients within a block,
/// - `n_coeff` — number of Chebyshev coefficients per component,
/// - `n_sub` — number of sub‑intervals the block is split into for that body.
///
/// The table has a fixed size of 15 entries, one per Horizons body ID (`0..=14`).
///
/// See also
/// --------
/// * [`HorizonRecord`] – per‑sub‑interval data constructed using IPT metadata.
pub type IPT = [[u32; 3]; 15];

/// Header metadata extracted from a JPL Horizons binary ephemeris.
///
/// This structure summarizes the global properties of the ephemeris file:
/// the version tag (e.g. `DE440`), the IPT table that describes how each
/// body's coefficients are laid out inside a block, the overall time
/// coverage of the file, the fixed block duration used for interpolation,
/// the raw record size on disk, and the Earth–Moon mass ratio (EMRAT).
///
/// Semantics
/// ---------
/// * Times are expressed in Julian Date on the TDB time scale.
/// * The ephemeris is partitioned into fixed‑length blocks (`period_lenght`),
///   and each block may be subdivided into `n_sub` sub‑intervals per body
///   according to the IPT entry.
/// * `recsize` indicates how many bytes compose one raw record as stored
///   in the binary; it is needed to seek precisely to a given block.
///
/// Invariants
/// ----------
/// * `start_period < end_period`
/// * `period_lenght > 0`
/// * `recsize > 0`
///
/// Typical usage
/// -------------
/// The parser fills this header first, then uses `ipt` to locate, decode,
/// and assemble per‑body [`HorizonRecord`] segments across all blocks.
///
/// See also
/// ------------
/// * [`IPT`] – per‑body layout within a block: `[offset, n_coeff, n_sub]`.
/// * [`HorizonRecord`] – atomic segment with Chebyshev coefficients.
/// * [`HorizonData`] – full ephemeris container (header + records).
/// * [`JPLHorizonVersion`] – helpers for Horizons version handling.
#[derive(Debug, PartialEq, Clone)]
pub struct HorizonHeader {
    /// JPL/Horizons version label (e.g., `DE440`, `DE441`).
    jpl_version: String,

    /// IPT metadata table that indexes where and how to read coefficients per body.
    ipt: IPT,

    /// Start of the ephemeris coverage (Julian Date, TDB).
    start_period: f64,

    /// End of the ephemeris coverage (Julian Date, TDB).
    end_period: f64,

    /// Duration of a single block/record (days).
    ///
    /// Interpolation is performed within these fixed-length blocks, potentially
    /// split in `n_sub` sub‑intervals according to the IPT entry for each body.
    period_lenght: f64, // NOTE: consider renaming to `period_length`.

    /// Size of a raw record (in bytes) as stored in the binary file.
    recsize: usize,

    /// Earth–Moon mass ratio (EMRAT) used by the ephemeris.
    earth_moon_mass_ratio: f64,
}

/// Parsed Horizons ephemeris: header + time‑ordered records.
///
/// This is the main container returned by the parser. It bundles the global
/// metadata (file coverage, IPT layout, record sizing) together with all
/// the time‑ordered data blocks for quick access and interpolation.
///
/// Fields
/// ------
/// * `header` — Global metadata extracted from the file header.
/// * `records` — Time-ordered sequence of per‑body Chebyshev segments.
///
/// Typical usage
/// -------------
/// 1. Open an ephemeris file via [`EphemFilePath`],
/// 2. Parse the header to obtain coverage and IPT,
/// 3. Stream/parse the blocks to populate `records`,
/// 4. Interpolate state vectors for a given `JD` using the right block and sub‑interval,
///    returning an [`InterpResult`].
///
/// See also
/// -----------
/// * [`HorizonHeader`] – global properties and IPT table.
/// * [`HorizonRecord`] – atomic interpolation segment (Chebyshev coefficients).
/// * [`InterpResult`] – output of an interpolation at a given epoch.
/// * [`JPLHorizonVersion`] – helper for version/compatibility checks.
/// * [`HorizonID`] – mapping between numeric IDs and named bodies.
#[derive(Debug, Clone)]
pub struct HorizonData {
    header: HorizonHeader,
    records: HorizonRecords,
}

/// Return the dimensionality of the data vector for a given Horizons body.
///
/// Each entry in the ephemeris has a fixed number of components
/// (e.g. 3 for Cartesian position, 1 for time offset, etc.).
///
/// Arguments
/// -----------------
/// * `index` — Horizons integer ID (`0..=14`).
///
/// Return
/// ----------
/// * Number of components (`usize`) for the given body:
///   - `0..=10` → 3 (planets + Sun, Cartesian x/y/z),
///   - `11` → 2 (Earth–Moon barycenter, 2D offset),
///   - `12` → 3 (lunar libration angles),
///   - `13` → 3 (lunar Euler angle rates),
///   - `14` → 1 (TT–TDB offset),
///   - otherwise `0` (invalid / unused).
///
/// See also
/// ------------
/// * [`IPT`] – table that, together with this dimensionality,
///   defines the structure of each record.
fn dimension(index: usize) -> usize {
    match index {
        0..=10 => 3, // Planets + Sun
        11 => 2,     // Earth–Moon barycenter
        12 => 3,     // Lunar librations
        13 => 3,     // Lunar Euler angle rates
        14 => 1,     // TT–TDB offsets
        _ => 0,      // Safety fallback
    }
}

/// Compute the total size (in bytes) of a binary record block.
///
/// This scans the IPT table and sums the contributions of each body:
/// number of sub-intervals × number of coefficients × dimensionality,
/// doubled (for position and velocity), then scaled to bytes.
///
/// Arguments
/// -----------------
/// * `ipt` — IPT table describing the layout of coefficients for each body.
///
/// Return
/// ----------
/// * Total size of one binary record (in bytes).
///
/// Notes
/// ----------
/// * The factor `4` corresponds to 4-byte words in the original
///   Fortran binary representation.
/// * This value should match the `recsize` field in [`HorizonHeader`].
///
/// See also
/// ------------
/// * [`dimension`] – per-body dimensionality.
/// * [`HorizonHeader::recsize`] – stored record size.
fn compute_recsize(ipt: IPT) -> usize {
    let kernel_size: usize = 4
        + (0..15)
            .map(|i| {
                let n_subintervals = ipt[i][1] as usize;
                let n_coeffs = ipt[i][2] as usize;
                let dim = dimension(i);
                2 * n_subintervals * n_coeffs * dim
            })
            .sum::<usize>();

    kernel_size * 4
}

/// Parse the IPT (Index Pointer Table) section from a Horizons binary file.
///
/// The IPT table is a 15×3 matrix of unsigned 32-bit integers, giving
/// the layout of Chebyshev coefficients for each body:
/// - column 0 → offset inside the block,
/// - column 1 → number of coefficients,
/// - column 2 → number of sub-intervals.
///
/// In addition, this function extracts the `numde` identifier and
/// handles a special case for row 12 (libration data).
///
/// Arguments
/// -----------------
/// * `input` — Raw byte slice containing the IPT table to decode.
///
/// Return
/// ----------
/// * An [`IResult`] containing:
///   - the remaining unparsed input slice,
///   - a tuple `(IPT, numde)` where:
///     - `IPT` is the 15×3 coefficient layout table,
///     - `numde` is the ephemeris number (e.g. 440 for DE440).
///
/// Panics
/// ----------
/// * If the libration row (`lpt`) does not contain exactly 3 elements.
///
/// See also
/// ------------
/// * [`IPT`] – per-body coefficient layout.
/// * [`HorizonHeader`] – stores the parsed IPT table for later use.
fn parse_ipt(input: &[u8]) -> IResult<&[u8], (IPT, u32)> {
    let mut ipt: IPT = [[0; 3]; 15];

    let mut input_remaining = input;
    for i in 0..36 {
        let (rest, val) = le_u32(input_remaining)?;
        input_remaining = rest;
        let row = i / 3;
        let col = i % 3;
        ipt[row][col] = val;
    }

    let (remaining_input, numde) = le_u32(input_remaining)?;
    let (input_remaining, lpt) = count(le_u32, 3).parse(remaining_input)?;

    ipt[12] = lpt
        .try_into()
        .expect("Expected lpt to have exactly 3 elements");

    Ok((input_remaining, (ipt, numde)))
}

/// Decode the trailing IPT rows for entries 13 and 14 (libration & TT–TDB).
///
/// This helper reads two consecutive 3×`u32` triplets from the given byte slice,
/// corresponding to `IPT[13]` and `IPT[14]`. Each triplet follows the usual IPT
/// convention: `[offset, n_coeff, n_sub]`.
///
/// Arguments
/// -----------------
/// * `input` — Raw byte slice containing the remaining IPT bytes to decode
///   (must hold exactly 24 bytes: 2 × 3 × 4).
///
/// Return
/// ----------
/// * An [`IResult`] with:
///   - the remaining unparsed input slice,
///   - a `[[u32; 3]; 2]` array where:
///     - index `0` is `IPT[13]` (lunar Euler angle rates / related set),
///     - index `1` is `IPT[14]` (TT–TDB offset series).
///
/// Errors
/// ----------
/// * Returns a nom error if the input does not contain enough bytes.
///
/// See also
/// ------------
/// * [`IPT`] – main 15×3 IPT table used by the ephemeris layout.
/// * [`parse_ipt`] – parses the leading IPT rows and `numde`.
fn parse_ipt_13_14(input: &[u8]) -> IResult<&[u8], [[u32; 3]; 2]> {
    fn parse_three(mut input: &[u8]) -> IResult<&[u8], [u32; 3]> {
        let mut result = [0u32; 3];
        for slot in &mut result {
            let (rest, val) = le_u32(input)?;
            *slot = val;
            input = rest;
        }
        Ok((input, result))
    }

    let (input, ipt_13) = parse_three(input)?;
    let (input, ipt_14) = parse_three(input)?;
    Ok((input, [ipt_13, ipt_14]))
}

/// Read `IPT[13]` and `IPT[14]` from the constants area of a DE binary (DE440+).
///
/// For DE versions ≥ DE440 and when the constants count `ncon` exceeds 400,
/// the two extra IPT rows (13, 14) are stored in the constants section after
/// the 400‑th constant name. This function seeks to that location and decodes
/// the two triplets.
///
/// Arguments
/// -----------------
/// * `file` — Buffered file reader positioned anywhere (seek is performed).
/// * `ncon` — Number of constants present in the ephemeris file.
/// * `de_version` — Parsed DE version (e.g., `DE440`).
///
/// Return
/// ----------
/// * `Ok(Some([[u32; 3]; 2]))` containing the decoded `IPT[13]` and `IPT[14]`
///   triplets when available; `Ok(None)` if the DE version is older than DE440
///   or `ncon <= 400`.
///
/// Errors
/// ----------
/// * I/O errors on seek/read,
/// * `InvalidData` if the trailing IPT rows cannot be parsed.
///
/// Notes
/// ----------
/// * The read offset is computed from the start of the 400‑th constant name,
///   then advanced by `(ncon - 400) * 6` bytes to reach the IPT rows area.
///
/// See also
/// ------------
/// * [`parse_ipt_13_14`] – decoder for the 24‑byte trailing IPT payload.
/// * [`IPT`] – ephemeris layout table used across the module.
fn read_ipt_13_14(
    file: &mut BufReader<File>,
    ncon: u32,
    de_version: &JPLHorizonVersion,
) -> std::io::Result<Option<[[u32; 3]; 2]>> {
    if *de_version < JPLHorizonVersion::DE440 || ncon <= 400 {
        return Ok(None); // Not present prior to DE440 or when ncon <= 400.
    }

    const START_400TH_CONSTANT_NAME: u64 = 2856;
    let offset = START_400TH_CONSTANT_NAME + (ncon as u64 - 400) * 6;

    file.seek(std::io::SeekFrom::Start(offset))?;
    let mut buffer = [0u8; 24];
    file.read_exact(&mut buffer)?;

    let (_, ipt_extra) = parse_ipt_13_14(&buffer).map_err(|_| {
        std::io::Error::new(
            std::io::ErrorKind::InvalidData,
            "Failed to parse ipt[13][14]",
        )
    })?;
    Ok(Some(ipt_extra))
}

/// Parse a block of `ncoeff` f64 values from the input (Chebyshev kernel slice).
///
/// This utility reads a contiguous sequence of `ncoeff` little‑endian `f64`
/// values, which typically represent the Chebyshev coefficients for a given
/// body, component, and sub‑interval inside a DATA RECORD.
///
/// Arguments
/// -----------------
/// * `input` — Raw byte slice containing at least `8 * ncoeff` bytes.
/// * `ncoeff` — Number of `f64` values to decode.
///
/// Return
/// ----------
/// * An [`IResult`] with:
///   - the remaining unparsed input slice,
///   - a `Vec<f64>` of length `ncoeff` containing the decoded coefficients.
///
/// Errors
/// ----------
/// * Returns a nom error if fewer than `ncoeff` doubles are available.
///
/// See also
/// ------------
/// * [`HorizonRecord`] – holds the per‑segment Chebyshev coefficients.
/// * [`IPT`] – provides `n_coeff` and `n_sub` used to size such blocks.
fn parse_block(input: &[u8], ncoeff: usize) -> IResult<&[u8], Vec<f64>> {
    count(le_f64, ncoeff).parse(input)
}

/// Extract all sub-records for a given body from one DATA RECORD block.
///
/// Each Horizons DATA RECORD contains Chebyshev coefficients for all bodies,
/// arranged according to the IPT entry. This function slices out the coefficients
/// corresponding to a single body and produces a list of [`HorizonRecord`]s,
/// one per sub-interval.
///
/// Arguments
/// -----------------
/// * `block` — Full decoded block (vector of `f64`), including leading JD start/end
///   followed by concatenated coefficients for all bodies.
/// * `ipt` — IPT entry `[offset, n_coeff, n_sub]` for this body:
///   - `offset`: starting index (in 4-byte words) to the body’s coefficients,
///   - `n_coeff`: number of Chebyshev coefficients per component,
///   - `n_sub`: number of sub-intervals within the block.
///
/// Return
/// ----------
/// * A `Vec<HorizonRecord>` with length equal to `n_sub`. Each record holds:
///   - block start/end epoch (`jd_start`, `jd_end`),
///   - the slice of Chebyshev coefficients for that sub-interval.
///
/// Notes
/// ----------
/// * The first two entries of `block` are reserved for `[jd_start, jd_end]`
///   and are skipped.
/// * The returned [`HorizonRecord`]s share references into the same `block` slice.
///
/// See also
/// ------------
/// * [`parse_all_blocks`] – drives this extraction for all bodies and all records.
/// * [`HorizonRecord`] – container for one sub-interval’s coefficients.
fn extract_body_records(block: &[f64], ipt: [u32; 3]) -> Vec<HorizonRecord> {
    let offset = ipt[0] as usize;
    let n_coeffs = ipt[1] as usize;
    let n_subs = ipt[2] as usize;

    let jd_start = block[0];
    let jd_end = block[1];

    let coeffs = &block[2..]; // skip JD start/end
    let mut records = Vec::with_capacity(n_subs);

    for subinterval_number in 0..n_subs {
        let record = HorizonRecord::new(
            jd_start,
            jd_end,
            coeffs,
            offset,
            subinterval_number,
            n_coeffs,
        );
        records.push(record);
    }
    records
}

/// Parse and decode all DATA RECORD blocks from a Horizons binary file.
///
/// This routine seeks past the first two header records, then reads the
/// remaining file contents into memory. It partitions the file into fixed-size
/// DATA RECORD blocks and decodes each block into a set of body records
/// using [`extract_body_records`].
///
/// Arguments
/// -----------------
/// * `file` — Opened `BufReader<File>` pointing to the Horizons binary.
/// * `recsize` — Size of each DATA RECORD in **bytes** (see `HorizonHeader::recsize`).
/// * `ncoeff` — Number of `f64` coefficients to parse per block (sum over bodies).
/// * `ipt` — Full IPT table describing the offset/size/sub-intervals per body.
///
/// Return
/// ----------
/// * `Ok(HorizonRecords)` – Vector of hashmaps, one per block:
///   - outer index = block index (time-ordered),
///   - inner key = `u8` body ID (`0..=14`),
///   - value = `Vec<HorizonRecord>` (all sub-intervals for that body in this block).
/// * `Err(std::io::Error)` if a block cannot be read or parsed.
///
/// Notes
/// ----------
/// * Skips the first two header records (`2 * recsize` bytes) before reading.
/// * Currently extracts bodies `0..11` (planets + Sun + EMB); extension may
///   be needed for libration, nutation, TT–TDB.
/// * Uses [`parse_block`] to decode raw `f64` values from each slice.
///
/// See also
/// ------------
/// * [`extract_body_records`] – builds the per-body records within a block.
/// * [`HorizonRecord`] – container for one sub-interval’s coefficients.
/// * [`IPT`] – index pointer table giving offsets and sizes.
fn parse_all_blocks(
    file: &mut BufReader<File>,
    recsize: usize,
    ncoeff: usize,
    ipt: IPT,
) -> std::io::Result<HorizonRecords> {
    let offset = 2 * recsize;
    file.seek(std::io::SeekFrom::Start(offset as u64))?;

    let mut buffer = Vec::new();
    file.read_to_end(&mut buffer)?;

    let n_blocks = buffer.len() / recsize;

    (0..n_blocks)
        .map(|i| {
            let start = i * recsize;
            let end = start + recsize;
            let slice = &buffer[start..end];

            let (_, block) = parse_block(slice, ncoeff).map_err(|_| {
                std::io::Error::new(
                    std::io::ErrorKind::InvalidData,
                    format!("Failed to parse block {i}"),
                )
            })?;

            let records = (0..12)
                .map(|body| {
                    let ipt_body = ipt[body];
                    let body_rec = extract_body_records(&block, ipt_body);
                    (body as u8, body_rec)
                })
                .collect::<HashMap<u8, _>>();

            Ok(records)
        })
        .collect()
}

impl HorizonData {
    /// Read and parse a JPL Horizons legacy ephemeris binary file.
    ///
    /// This function:
    /// 1. Opens the ephemeris file at the given path,
    /// 2. Parses the header sections (TTL banner, CNAM constants, SS array),
    /// 3. Reads the IPT (Index Pointer Table) and, if available, the extended
    ///    IPT rows (13 and 14 for librations and TT–TDB),
    /// 4. Computes the record size (`recsize`) and number of coefficients,
    /// 5. Parses all DATA RECORD blocks into [`HorizonRecord`]s.
    ///
    /// Arguments
    /// -----------------
    /// * `ephem_path` — Path to the Horizons DE binary file
    ///   (e.g. DE440, DE441) wrapped in [`EphemFilePath`].
    ///
    /// Return
    /// ----------
    /// * A fully constructed [`HorizonData`] containing:
    ///   - [`HorizonHeader`] metadata (version, IPT, coverage, EMRAT, etc.),
    ///   - All time-ordered [`HorizonRecord`]s grouped by body ID.
    ///
    /// Panics
    /// ----------
    /// * If the file cannot be opened or read,
    /// * If any parsing step (TTL, CNAM, SS, IPT, blocks) fails.
    ///
    /// See also
    /// ------------
    /// * [`HorizonHeader`] – global header extracted from the file.
    /// * [`HorizonRecord`] – container for one sub-interval’s coefficients.
    pub fn read_horizon_file(ephem_path: &EphemFilePath) -> Self {
        let version = match ephem_path {
            EphemFilePath::JPLHorizon(_, version) => version,
            _ => panic!("Invalid ephemeris file path"),
        };

        const OLD_MAX: usize = 400;

        /// Parse une chaîne Fortran CHAR*6 (6 bytes)
        fn parse_char6(input: &[u8]) -> IResult<&[u8], String> {
            let (rest, raw) = take(6usize)(input)?;
            Ok((rest, String::from_utf8_lossy(raw).trim_end().to_string()))
        }

        /// Parse les en-têtes TTL (14×3 CHAR*6)
        fn parse_ttl(input: &[u8]) -> IResult<&[u8], Vec<String>> {
            count(parse_char6, 14 * 3).parse(input)
        }

        /// Parse le champ CNAM (400 chaînes CHAR*6)
        fn parse_cnam(input: &[u8]) -> IResult<&[u8], Vec<String>> {
            count(parse_char6, OLD_MAX).parse(input)
        }

        let mut file = BufReader::new(
            File::open(ephem_path.path())
                .unwrap_or_else(|_| panic!("Failed to open the JPL ephemeris file: {ephem_path}")),
        );

        let mut file_data = vec![0u8; 1 << 12];
        file.read_exact(&mut file_data)
            .expect("Failed to read the JPL ephemeris file.");

        let (input, _) = parse_ttl(&file_data).expect("failed to parse TTL");

        let (input, _) = parse_cnam(input).expect("failed to parse CNAM");

        let (input, ss) = count(le_f64::<_, nom::error::Error<_>>, 3)
            .parse(input)
            .expect("failed to parse SS");

        let (input, ncon) = le_i32::<_, nom::error::Error<_>>(input).expect("failed to parse NCON");

        let (input, _) =
            take::<_, _, nom::error::Error<_>>(8usize)(input).expect("failed to parse 4 bytes");
        let (_, earth_moon_mass_ratio) = le_f64::<_, nom::error::Error<_>>(input)
            .expect("failed to parse Earth-Moon mass ratio");

        file.seek(std::io::SeekFrom::Start(2696))
            .expect("Failed to seek to the JPL ephemeris file IPT.");
        let mut buffer = [0u8; 204]; // JPL_HEADER_SIZE
        file.read_exact(&mut buffer)
            .expect("Failed to read the JPL ephemeris file IPT.");

        let (_, (mut ipt, numde)) = parse_ipt(&buffer).expect("failed to parse IPT header");

        let ipt_extra = read_ipt_13_14(&mut file, ncon as u32, version)
            .expect("Failed to read IPT[13] and IPT[14]");

        ipt[13] = ipt_extra
            .as_ref()
            .map(|ipt_extra| ipt_extra[0])
            .unwrap_or([0; 3]);

        ipt[14] = ipt_extra
            .as_ref()
            .map(|ipt_extra| ipt_extra[1])
            .unwrap_or([0; 3]);

        let recsize = compute_recsize(ipt);

        let ncoeff = recsize / 8;
        let blocks = parse_all_blocks(&mut file, recsize, ncoeff, ipt).unwrap();

        HorizonData {
            header: HorizonHeader {
                jpl_version: format!("DE{numde}"),
                ipt,
                start_period: ss[0],
                end_period: ss[1],
                period_lenght: ss[2],
                recsize,
                earth_moon_mass_ratio,
            },
            records: blocks,
        }
    }

    /// Locate the DATA RECORD index and normalized time fraction (`tau`)
    /// for a given ephemeris time.
    ///
    /// The Horizons binary is divided into fixed-length blocks of duration
    /// `header.period_length`. This method finds which block contains the
    /// requested Modified Julian Date (MJD) and computes the normalized
    /// Chebyshev parameter `tau ∈ [0,1]` inside that block.
    ///
    /// Arguments
    /// -----------------
    /// * `et` — Ephemeris time (Modified Julian Date, TDB).
    ///
    /// Return
    /// ----------
    /// * `(nr, tau)` where:
    ///   - `nr`: zero-based index of the block inside `records`,
    ///   - `tau`: normalized fractional time within that block.
    ///
    /// Panics
    /// ----------
    /// * If `et` lies outside the coverage of the ephemeris file.
    ///
    /// See also
    /// ------------
    /// * [`HorizonData::get_record_horizon`] – retrieves the actual record for a body.
    fn get_record_index(&self, et: MJD) -> (usize, f64) {
        // ephem_start and ephem_end are in JD
        let (ephem_start, ephem_end, ephem_step) = (
            self.header.start_period,
            self.header.end_period,
            self.header.period_lenght,
        );

        let jd_base = 2_400_000.5;
        let et_jd = jd_base + et.trunc();

        if et_jd < ephem_start || et_jd > ephem_end {
            panic!("Time outside ephemeris range");
        }

        let mut nr = ((et_jd - ephem_start) / ephem_step).floor() as usize;

        if (et_jd - ephem_end).abs() < 1e-10 {
            nr -= 1;
        }

        let interval_start = (nr as f64) * ephem_step + ephem_start;
        let tau = ((et_jd - interval_start) + et.fract()) / ephem_step;
        (nr, tau)
    }

    /// Retrieve the [`HorizonRecord`] and time fraction for a given body and epoch.
    ///
    /// This function:
    /// 1. Determines the block index and normalized time fraction with
    ///    \[`get_record_index`\],
    /// 2. Looks up the list of [`HorizonRecord`]s for the requested body,
    /// 3. Chooses the correct sub-interval based on `tau`.
    ///
    /// Arguments
    /// -----------------
    /// * `body` — Horizons body ID (`0..=14`).
    /// * `et` — Ephemeris time (Modified Julian Date, TDB).
    ///
    /// Return
    /// ----------
    /// * `Some((&HorizonRecord, tau))` where:
    ///   - `&HorizonRecord` — record covering the corresponding sub-interval,
    ///   - `tau` — normalized time fraction inside the parent block.
    ///
    /// Panics
    /// ----------
    /// * If the requested body or block is not present in `records`.
    ///
    /// See also
    /// ------------
    /// * [`HorizonRecord::interpolate`] – evaluate Chebyshev polynomials.
    fn get_record_horizon(&self, body: u8, et: MJD) -> Option<(&HorizonRecord, f64)> {
        let (nr, tau) = self.get_record_index(et);

        let records = &self.records[nr];
        let record_body = records
            .get(&body)
            .unwrap_or_else(|| panic!("Failed to get record for body {body} in block {nr}"));

        let ipt_body = self.header.ipt[body as usize];
        let n_subs = ipt_body[2] as usize;

        let sub_index = (tau * n_subs as f64).floor().min(n_subs as f64 - 1.0) as usize;

        Some((&record_body[sub_index], tau))
    }

    /// Interpolate the state vector of a target body relative to a center body.
    ///
    /// This is the main high-level query function. It retrieves the records
    /// for both `target` and `center`, evaluates their Chebyshev polynomials
    /// at the requested time, and returns the difference.
    ///
    /// Special handling:
    /// * If `target == Earth`, the Moon’s contribution is removed using the
    ///   Earth–Moon mass ratio (EMRAT).
    ///
    /// Arguments
    /// -----------------
    /// * `target` — Target body as [`HorizonID`].
    /// * `center` — Center body as [`HorizonID`].
    /// * `et` — Ephemeris time (Modified Julian Date, TDB).
    /// * `compute_velocity` — If `true`, also compute velocity.
    /// * `compute_acceleration` — If `true`, also compute acceleration.
    ///
    /// Return
    /// ----------
    /// * An [`InterpResult`] containing:
    ///   - `position` \[km\],
    ///   - `velocity` \[km/day\] (if requested),
    ///   - `acceleration` \[km/day²\] (if requested),
    ///     all expressed relative to the center body.
    ///
    /// See also
    /// ------------
    /// * [`InterpResult::to_au`] – convert to AU / AU/day.
    /// * [`HorizonRecord::interpolate`] – polynomial evaluation.
    pub fn ephemeris(
        &self,
        target: HorizonID,
        center: HorizonID,
        et: MJD,
        compute_velocity: bool,
        compute_acceleration: bool,
    ) -> InterpResult {
        let ipt_target = self.header.ipt[target as usize];
        let ipt_center = self.header.ipt[center as usize];
        let (record, tau) = self.get_record_horizon(target.into(), et).unwrap();
        let mut interp_target = record.interpolate(
            tau,
            compute_velocity,
            compute_acceleration,
            ipt_target[2] as usize,
        );

        if target == HorizonID::Earth {
            let ipt_moon = self.header.ipt[HorizonID::Moon as usize];
            let (record, tau) = self.get_record_horizon(HorizonID::Moon.into(), et).unwrap();
            let interp_moon = record.interpolate(
                tau,
                compute_velocity,
                compute_acceleration,
                ipt_moon[2] as usize,
            );
            interp_target = interp_target - interp_moon / (1. + self.header.earth_moon_mass_ratio);
        }

        let (record, tau) = self.get_record_horizon(center.into(), et).unwrap();
        let interp_center = record.interpolate(
            tau,
            compute_velocity,
            compute_acceleration,
            ipt_center[2] as usize,
        );

        interp_target - interp_center
    }
}

#[cfg(test)]
mod test_horizon_reader {
    #[cfg(feature = "jpl-download")]
    use super::*;

    #[cfg(feature = "jpl-download")]
    use crate::unit_test_global::JPL_EPHEM_HORIZON;

    #[test]
    #[cfg(feature = "jpl-download")]
    fn test_jpl_reader_from_horizon() {
        assert_eq!(
            JPL_EPHEM_HORIZON.header,
            HorizonHeader {
                jpl_version: "DE440".to_string(),
                ipt: [
                    [3, 14, 4,],
                    [171, 10, 2,],
                    [231, 13, 2,],
                    [309, 11, 1,],
                    [342, 8, 1,],
                    [366, 7, 1,],
                    [387, 6, 1,],
                    [405, 6, 1,],
                    [423, 6, 1,],
                    [441, 13, 8,],
                    [753, 11, 2,],
                    [819, 10, 4,],
                    [899, 10, 4,],
                    [1019, 0, 0,],
                    [1019, 0, 0,],
                ],
                start_period: 2287184.5,
                end_period: 2688976.5,
                period_lenght: 32.0,
                recsize: 8144,
                earth_moon_mass_ratio: 81.30056822149722,
            }
        );

        assert_eq!(JPL_EPHEM_HORIZON.records.len(), 12556);
        assert_eq!(
            JPL_EPHEM_HORIZON
                .records
                .iter()
                .fold(0, |acc, x| acc + x.len()),
            150672
        );

        assert_eq!(
            &JPL_EPHEM_HORIZON.records[0].get(&0).unwrap()[0],
            &HorizonRecord {
                start_jd: 2287184.5,
                end_jd: 2287216.5,
                x: vec![
                    -45337704.29199142,
                    -11420952.2182182,
                    1231640.71525489,
                    13474.74253284046,
                    -4516.491868368539,
                    255.9277522073883,
                    -0.6852443171165454,
                    -1.186291661304585,
                    0.1035908763014567,
                    -0.002445022374643364,
                    -0.0003798510942101873,
                    4.866105543177588e-5,
                    -2.092814231764934e-6,
                    -1.118766077984058e-7
                ],
                y: vec![
                    19369902.32537584,
                    -12637978.7311934,
                    -560491.8087798061,
                    75269.29929452346,
                    -1778.157215889909,
                    -139.1111450981129,
                    18.58942667472901,
                    -0.8470743363574242,
                    -0.02513657259988719,
                    0.006832859866311067,
                    -0.0004590542406706136,
                    5.317330201500909e-7,
                    2.764140697666422e-6,
                    -2.579709214481493e-7
                ],
                z: vec![
                    15085546.93080799,
                    -5541484.32081567,
                    -428281.1733889182,
                    38734.17578050081,
                    -474.2794185515511,
                    -101.0728095721337,
                    9.98761508428064,
                    -0.3272813126206938,
                    -0.02428417941794788,
                    0.003901383011452982,
                    -0.000204980181894197,
                    -4.825494149048871e-6,
                    1.694428982754255e-6,
                    -1.260377702474755e-7
                ]
            }
        );
    }

    #[test]
    #[cfg(feature = "jpl-download")]
    fn test_get_record_from_horizon() {
        let (record, tau) = JPL_EPHEM_HORIZON
            .get_record_horizon(4, 57028.479297592596)
            .unwrap();

        assert_eq!(tau, 0.6399780497686152);

        assert_eq!(
            record,
            &HorizonRecord {
                start_jd: 2457008.5,
                end_jd: 2457040.5,
                x: vec![
                    -558258575.175456,
                    -13081137.99303625,
                    70236.79190208673,
                    240.5982017483768,
                    -0.8633795271054904,
                    -0.0002571950240619499,
                    4.544021170810181e-6,
                    -1.839536840631949e-8
                ],
                y: vec![
                    515866454.9405554,
                    -10984765.90508825,
                    -64874.14676298688,
                    261.5078409285777,
                    0.4491808227783525,
                    -0.001923759413667307,
                    2.250631406524948e-6,
                    -6.23548505589978e-8
                ],
                z: vec![
                    234693317.6472925,
                    -4389890.166800962,
                    -29516.72247187927,
                    106.2324198779676,
                    0.2135457077378355,
                    -0.000818018169817293,
                    9.216197575953068e-7,
                    -2.484668250013335e-8
                ]
            },
        );
    }

    #[test]
    #[cfg(feature = "jpl-download")]
    fn test_get_record_index() {
        let (index, tau) = JPL_EPHEM_HORIZON.get_record_index(57028.479297592596);

        assert_eq!(index, 5307);
        assert_eq!(tau, 0.6399780497686152);
    }

    #[test]
    #[cfg(feature = "jpl-download")]
    fn test_interpolation_from_horizon() {
        let (record, tau) = JPL_EPHEM_HORIZON
            .get_record_horizon(10, 57028.479297592596)
            .unwrap();

        let res = record.interpolate(tau, true, true, 2);

        assert_eq!(
            res,
            InterpResult {
                position: [[428149.04652929713, -105270.2354841389, -68083.3417807072]].into(),
                velocity: Some([[589.5451313546943, 729.3492107653134, 300.3651374864013]].into()),
                acceleration: Some(
                    [[-0.919215789187584, 0.8829808730528121, 0.38984144066801635]].into()
                )
            }
        );

        let (record, tau) = JPL_EPHEM_HORIZON
            .get_record_horizon(10, 57_049.231_857_592_59)
            .unwrap();

        let res = record.interpolate(tau, true, true, 2);

        assert_eq!(
            res,
            InterpResult {
                position: [[440183.15997455275, -89933.41046658876, -61760.61450611751]].into(),
                velocity: Some([[569.7900066838628, 749.1773000798205, 309.2398409126585]].into()),
                acceleration: Some(
                    [[-1.038549415842939, 0.9990469757389817, 0.4553928281168678]].into()
                )
            }
        );

        let (record, tau) = JPL_EPHEM_HORIZON
            .get_record_horizon(10, 60781.51949044435)
            .unwrap();
        let res = record.interpolate(tau, true, true, 2);

        assert_eq!(
            res,
            InterpResult {
                position: [[-742814.3000341875, -727671.3536844663, -288321.53733077314]].into(),
                velocity: Some(
                    [[1085.6256324761375, -327.2648113240611, -162.13166222548898]].into()
                ),
                acceleration: Some(
                    [[-0.06525161089395584, 1.261197328255844, 0.5376907387973575]].into()
                )
            }
        )
    }

    #[test]
    #[cfg(feature = "jpl-download")]
    fn test_target_center_interpolation() {
        let interp = JPL_EPHEM_HORIZON.ephemeris(
            HorizonID::Earth,
            HorizonID::Sun,
            60781.51949044435,
            true,
            true,
        );

        assert_eq!(
            interp.to_au(),
            InterpResult {
                position: [[
                    -0.8988034555610618,
                    -0.4096429669081428,
                    -0.17756458835186803
                ]]
                .into(),
                velocity: Some(
                    [[
                        0.007368756100393145,
                        -0.014193450423040196,
                        -0.006152419296313352
                    ]]
                    .into()
                ),
                acceleration: Some(
                    [[
                        0.0002624924486094171,
                        0.00011873474463671509,
                        5.13298993717986e-5
                    ]]
                    .into()
                )
            }
        );

        let interp = JPL_EPHEM_HORIZON.ephemeris(
            HorizonID::Mars,
            HorizonID::Sun,
            60781.51949044435,
            true,
            true,
        );

        assert_eq!(
            interp.to_au(),
            InterpResult {
                position: [[-1.5211490583088887, 0.601249019351933, 0.31680918257233887]].into(),
                velocity: Some(
                    [[
                        -0.005170287251815057,
                        -0.010584792921766313,
                        -0.0047155422859661905
                    ]]
                    .into()
                ),
                acceleration: Some(
                    [[
                        9.733944178212551e-5,
                        -3.84704427625744e-5,
                        -2.0271140253173187e-5
                    ]]
                    .into()
                )
            }
        );

        let interp = JPL_EPHEM_HORIZON.ephemeris(
            HorizonID::Earth,
            HorizonID::Sun,
            52550.18467592593,
            true,
            true,
        );

        assert_eq!(
            interp.to_au(),
            InterpResult {
                position: [[0.9861809593158523, 0.1557130752386633, 0.06750574448131498]].into(),
                velocity: Some(
                    [[
                        -0.003193635816506268,
                        0.015502851910183767,
                        0.006721021432383671
                    ]]
                    .into()
                ),
                acceleration: Some(
                    [[
                        -0.0002926932324959577,
                        -4.506599368500431e-5,
                        -1.9370434372826522e-5
                    ]]
                    .into()
                )
            }
        );
    }
}