indicatrix 0.9.0

//! Acceptance tests of chromophore-color-plan.md §11, criteria 1-6 and 10-14, each asserted at
//! the stated threshold.

#![expect(
    clippy::doc_markdown,
    clippy::suboptimal_flops,
    clippy::cast_precision_loss,
    reason = "test code quoting the plan's symbols and thresholds"
)]

use std::time::Instant;

use super::*;
use crate::{
    color::{
        body_color::{BodyColors, Illuminant, body_color, body_colors, delta_e_2000, xyz_to_lab},
        cie1931::cie_1931_cmf,
    },
    geometry::cuts::StandardGemCuts,
    optics::{
        absorption::{AbsorptionBand, AbsorptionTensor},
        materials::GemMaterial,
    },
    render_setup::{
        materials::{MaterialOverrides, apply_material_overrides, effective_stone_width_mm},
        measure_model_width,
    },
};

use super::primary_data_tests::{PPMA_CM3, hue_deg as gia_hue_deg, sigma_lab, table};

const COLORLESS: [f64; 3] = [100.0, 0.0, 0.0];
/// The 5 mm reference path of the plan's "blue recipe".
const REF_MM: f64 = 5.0;
/// [`REF_MM`] as the `f32` path length the solver API takes.
#[expect(
    clippy::cast_possible_truncation,
    reason = "5.0 is exactly representable in f32"
)]
const REF_MM_F32: f32 = REF_MM as f32;

fn cat() -> &'static ChromophoreCatalogue {
    ChromophoreCatalogue::global()
}

fn recipe_of(host: &str, entries: &[(&str, f64)]) -> ColorRecipe {
    let mut r = ColorRecipe::new(host, cat().data_version);
    for (id, amount) in entries {
        assert!(r.set_amount(id, *amount));
    }
    r
}

fn colors_of(recipe: &ColorRecipe, path_mm: f64, ill: Illuminant) -> BodyColors {
    let (tensor, _) = resolve(recipe, cat()).expect("resolves cleanly");
    body_colors(&tensor, path_mm, ill)
}

fn hue_deg(lab: [f64; 3]) -> f64 {
    lab[2].atan2(lab[1]).to_degrees().rem_euclid(360.0)
}

fn chroma(lab: [f64; 3]) -> f64 {
    lab[1].hypot(lab[2])
}

/// splitmix64: a reproducible seeded generator.
struct Rng(u64);

impl Rng {
    fn next_u64(&mut self) -> u64 {
        self.0 = self.0.wrapping_add(0x9E37_79B9_7F4A_7C15);
        let mut z = self.0;
        z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
        z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
        z ^ (z >> 31)
    }

    /// Uniform in `[0, 1)`.
    fn unit(&mut self) -> f64 {
        (self.next_u64() >> 11) as f64 / (1u64 << 53) as f64
    }

    fn below(&mut self, n: usize) -> usize {
        usize::try_from(self.next_u64() % n as u64).expect("a remainder below n fits usize")
    }
}

/// A random recipe of 1..=`max_items` distinct selectable items with amounts up to
/// `max_frac * conc_max`; end-member fractions are scaled to sum <= 1. `None` for a host without
/// selectable items.
fn random_recipe(
    host_id: &str,
    rng: &mut Rng,
    max_items: usize,
    max_frac: f64,
) -> Option<ColorRecipe> {
    let host = cat().host(host_id)?;
    let mut ids = cat().selectable_elements(host_id);
    if ids.is_empty() {
        return None;
    }
    let k = 1 + rng.below(max_items.min(ids.len()));
    let mut recipe = ColorRecipe::new(host_id, cat().data_version);
    for _ in 0..k {
        let id = ids.remove(rng.below(ids.len()));
        let amount = host.element_conc_max(&id) * max_frac * (0.02 + 0.98 * rng.unit());
        recipe.set_amount(&id, amount);
    }
    let em_sum: f64 = host.end_members.iter().map(|m| recipe.amount(&m.id)).sum();
    if em_sum > 0.95 {
        for m in &host.end_members {
            let a = recipe.amount(&m.id);
            if a > 0.0 {
                recipe.set_amount(&m.id, a * 0.95 / em_sum);
            }
        }
    }
    Some(recipe)
}

/// Integrated o-ray absorbance (mm^-1 nm) of a recipe over 380..780 nm.
fn integrated_alpha(recipe: &ColorRecipe) -> f64 {
    let (t, _) = resolve(recipe, cat()).expect("resolves");
    (380..=780)
        .map(|l| {
            t.o_ray
                .iter()
                .map(|b| f64::from(b.evaluate(l as f32)))
                .sum::<f64>()
        })
        .sum()
}

/// Share of the integrated absorbance contributed by `id` alone.
fn alpha_share(recipe: &ColorRecipe, id: &str) -> f64 {
    let alone = recipe_of(&recipe.host, &[(id, recipe.amount(id))]);
    integrated_alpha(&alone) / integrated_alpha(recipe).max(1e-30)
}

// ---------------------------------------------------------------------------------------------
// Criterion 1
// ---------------------------------------------------------------------------------------------

/// Cr amounts of the five criterion-1 steps (0.05 -> 0.5 wt% Cr2O3).
fn cr_steps() -> [f64; 5] {
    std::array::from_fn(|i| 0.05 + (0.5 - 0.05) * i as f64 / 4.0)
}

/// Criterion 1: the Cr hue lies in [0, 30] U [330, 360) at every step (5 mm reference path).
///
/// `data_version` 4: the band is derived from the GIA 2020 cross-section. GIA's calculated
/// unpolarised colors (Appendix 2, D65, 100-3000 ppma*cm) have hues 339, 340, 343, 352, 8 and
/// 26 degrees, i.e. inside the plan's band, and the five steps (67-670 ppma*cm) sit at
/// 336-354 degrees (checked against the cross-section in
/// `criterion_1_chroma_and_hue_follow_the_gia_cross_section`).
#[test]
fn criterion_1_cr_hue() {
    for cr in cr_steps() {
        let lab = colors_of(
            &recipe_of("corundum", &[("Cr", cr)]),
            REF_MM,
            Illuminant::D65,
        )
        .unpolarised
        .lab;
        let h = hue_deg(lab);
        assert!(
            (0.0..=30.0).contains(&h) || (330.0..360.0).contains(&h),
            "Cr {cr} wt%: hue {h:.1} deg outside [0, 30] U [330, 360), Lab {lab:?}"
        );
    }
}

/// Criterion 1: chroma C* increases monotonically over the 5 Cr steps at the 5 mm reference path.
///
/// Un-ignored for `data_version` 4. The old ignore rested on the research's stated CIELAB (C*
/// 32.5 / 73.8 / 40.2 at 0.05 / 0.25 / 0.72 wt%), which disagrees with its own spectra by dE00
/// 16-26 and is no target (`reference_validation`). From the GIA cross-section the C* of the
/// unpolarised color rises monotonically over 67-670 ppma*cm (the five steps) and saturates
/// near 700-3000 ppma*cm (GIA Appendix 2 "both", D65: 27.8, 47.6, 65.0, 70.8, 69.5, 70.0 at 100,
/// 200, 375, 750, 1500, 3000 ppma*cm); the recipe gives 22.5, 47.1, 57.1, 61.6, 63.7 (the cross-section itself 22.7, 48.3, 58.2, 62.5, 64.3).
#[test]
fn criterion_1_cr_chroma_is_monotonic() {
    let mut last = 0.0;
    for cr in cr_steps() {
        let lab = colors_of(
            &recipe_of("corundum", &[("Cr", cr)]),
            REF_MM,
            Illuminant::D65,
        )
        .unpolarised
        .lab;
        let c = chroma(lab);
        assert!(
            c > last,
            "Cr {cr} wt%: C* {c:.2} does not exceed the previous step's {last:.2}"
        );
        last = c;
    }
}

/// Criterion 1: the blue recipe (Fe 1000 + Ti 100 ppm at 5 mm) is > 20 from colorless, its o/e
/// pleochroism is > 3, and Ti alone stays < 1 for any amount up to `conc_max`.
#[test]
fn criterion_1_blue_recipe_ti_alone_and_pleochroism() {
    let blue = colors_of(
        &recipe_of("corundum", &[("Fe", 1000.0), ("Ti", 100.0)]),
        REF_MM,
        Illuminant::D65,
    );
    let de = delta_e_2000(blue.unpolarised.lab, COLORLESS);
    assert!(
        de > 20.0,
        "blue recipe dE from colorless {de:.2} (Lab {:?})",
        blue.unpolarised.lab
    );
    let de_oe = delta_e_2000(blue.o_ray.lab, blue.e_ray.lab);
    assert!(de_oe > 3.0, "blue recipe dE(o, e) {de_oe:.2}");

    let host = cat().host("corundum").expect("corundum");
    for frac in [0.1, 0.5, 1.0] {
        let ti = host.element_conc_max("Ti") * frac;
        let lab = colors_of(
            &recipe_of("corundum", &[("Ti", ti)]),
            REF_MM,
            Illuminant::D65,
        )
        .unpolarised
        .lab;
        let de_ti = delta_e_2000(lab, COLORLESS);
        assert!(de_ti < 1.0, "Ti {ti} alone: dE from colorless {de_ti:.3}");
    }
}

/// Criterion 1: the blue recipe's hue h_ab.
///
/// `data_version` 4: the plan's [220, 270] was set by the research's S-Met color (b* -41, hue
/// 261). The GIA Fe2+-Ti4+ cross-section gives a violet-blue instead: the unpolarised color
/// (2 E_perp + E_par) / 3 of the recipe's 100 pairs over 5 mm (20 ppma*cm) has hue 291 degrees
/// (E_perp_c 289-313, E_par_c 245-270 over 50-200 pairs and 2-5 mm). The recipe's own Fe3+ and
/// Fe2+-Fe3+ bands (Fe 1000 ppm) shift the hue by about -25 degrees to 266. The band is the
/// GIA hue -50 / +20 degrees, i.e. [241, 311], inside the blue range [220, 310].
#[test]
fn criterion_1_blue_recipe_hue() {
    let lab = colors_of(
        &recipe_of("corundum", &[("Fe", 1000.0), ("Ti", 100.0)]),
        REF_MM,
        Illuminant::D65,
    )
    .unpolarised
    .lab;
    let h = hue_deg(lab);
    // Pairs = min(Fe2+ = 850, Ti 100) = 100 ppm of Al sites = 40 ppma of all atoms; over 5 mm.
    let g = table("gia_v3plus_fe2plus_ti4plus_sigma_5nm");
    let cols = (
        g.col("sigma_FeTi_E_perp_c_cm2"),
        g.col("sigma_FeTi_E_par_c_cm2"),
    );
    let gia = gia_hue_deg(sigma_lab(&g, cols, "both", 40.0 * 0.5, Illuminant::D65));
    eprintln!("blue recipe hue {h:.1}, GIA Fe-Ti hue {gia:.1}");
    assert!(
        (gia - 50.0..=gia + 20.0).contains(&h) && (220.0..=310.0).contains(&h),
        "blue recipe hue {h:.1} deg outside [GIA {gia:.1} - 50, + 20]"
    );
}

/// Criterion 1 (primary data): chroma and hue of the Cr recipe at the five criterion-1 steps
/// follow the color of the GIA cross-section (2 E_perp + E_par) / 3 within C* 2.5 and 5 degrees (measured <= 1.2 and <= 1),
/// and the GIA chroma itself rises monotonically over the steps.
#[test]
fn criterion_1_chroma_and_hue_follow_the_gia_cross_section() {
    let g = table("gia_cr3plus_sample1110_sigma");
    let cols = (g.col("sigma_E_perp_c_cm2"), g.col("sigma_E_par_c_cm2"));
    let n_wt = cat()
        .host("corundum")
        .expect("corundum")
        .n_site_for_unit("wt_pct_oxide:Cr2O3");
    let mut last = 0.0;
    for cr in cr_steps() {
        let ours = colors_of(
            &recipe_of("corundum", &[("Cr", cr)]),
            REF_MM,
            Illuminant::D65,
        )
        .unpolarised
        .lab;
        // wt% x n_wt = Cr per cm^3; over 0.5 cm in GIA's ppma*cm.
        let gia = sigma_lab(
            &g,
            cols,
            "both",
            cr * n_wt * 0.5 / PPMA_CM3,
            Illuminant::D65,
        );
        eprintln!(
            "Cr {cr:.4}: recipe C* {:.1} h {:.0}, GIA sigma C* {:.1} h {:.0}",
            chroma(ours),
            hue_deg(ours),
            chroma(gia),
            hue_deg(gia)
        );
        assert!((chroma(ours) - chroma(gia)).abs() <= 2.5);
        assert!((hue_deg(ours) - hue_deg(gia)).abs() <= 5.0);
        assert!(chroma(gia) > last, "the GIA chroma rises over the steps");
        last = chroma(gia);
    }
}

// ---------------------------------------------------------------------------------------------
// Criterion 2
// ---------------------------------------------------------------------------------------------

/// Peak of the E-perp-c Fe2+-Ti4+ band (cm^-1) of the blue recipe under the catalogue `cat`.
fn blue_pair_peak_cm1(cat: &ChromophoreCatalogue) -> f64 {
    let r = recipe_of("corundum", &[("Fe", 1000.0), ("Ti", 100.0)]);
    let (t, _) = resolve(&r, cat).expect("resolves");
    // The E-perp-c Fe2+-Ti4+ band (the strongest o band between 570 and 600 nm).
    let band = t
        .o_ray
        .iter()
        .filter(|b| (570.0..=600.0).contains(&b.center_nm))
        .map(|b| f64::from(b.peak))
        .fold(0.0, f64::max);
    band * 10.0
}

/// Criterion 2: random pairing would give < 0.1 cm^-1, the clustered law >= 2 and within 5 % of
/// the GIA value; the pure host transmits Y >= 0.98 at 10 mm.
///
/// `data_version` 4: the "empirical 5 cm^-1" of the research (0.05 per ppm of Al sites, fitted to
/// its sapphire tables) is replaced by the GIA 2020 cross-section per pair, sigma = 1.98e-18 cm^2
/// at 580 nm: 100 pairs (ppm of Al sites) = 4.70e18 cm^-3 give 9.31 cm^-1 (E_perp_c), 0.093 per
/// pair-ppm; the data_version 3 value was 5.0.
#[test]
fn criterion_2_pairing_regression_and_pure_host() {
    let clustered = blue_pair_peak_cm1(cat());
    let g = table("gia_v3plus_fe2plus_ti4plus_sigma_5nm");
    let n_site = cat()
        .host("corundum")
        .expect("corundum")
        .n_site_for_unit("ppm_site");
    let gia = g.at(g.col("sigma_FeTi_E_perp_c_cm2"), 580.0) * 100.0 * n_site;
    assert!(
        clustered >= 2.0 && (clustered / gia - 1.0).abs() <= 0.05,
        "clustered pair peak {clustered} cm^-1 (GIA {gia:.2})"
    );

    let mut random = cat().clone();
    let pair = random
        .hosts
        .iter_mut()
        .find(|h| h.id == "corundum")
        .expect("corundum")
        .chromophores
        .iter_mut()
        .find(|c| c.kind == "ivct_pair" && c.partners[1] == "Ti4+")
        .expect("Fe2+-Ti4+ pair");
    pair.pairing = Some("random".to_string());
    pair.pair_z = Some(6.0);
    let random_peak = blue_pair_peak_cm1(&random);
    assert!(
        random_peak < 0.1,
        "random pairing would give {random_peak} cm^-1 (must be < 0.1)"
    );

    for host in &cat().hosts {
        if cat().selectable_elements(&host.id).is_empty() {
            continue;
        }
        let pure = colors_of(
            &ColorRecipe::new(&host.id, cat().data_version),
            10.0,
            Illuminant::D65,
        );
        assert!(
            pure.unpolarised.xyz[1] >= 0.98,
            "pure {} at 10 mm: Y = {:.4}",
            host.id,
            pure.unpolarised.xyz[1]
        );
    }
}

// ---------------------------------------------------------------------------------------------
// Criterion 3
// ---------------------------------------------------------------------------------------------

/// Criterion 3: 200 reproducibly seeded recipes of <= 3 items (amounts <= 0.5 c_max) in corundum,
/// beryl, garnet and tanzanite: dE <= 1 for >= 95 %, <= 2 for 100 %, `reachable` iff dE <= 1.
#[test]
fn criterion_3_round_trip_of_200_seeded_recipes() {
    let mut rng = Rng(0x00C0_1042_5EED);
    let hosts = ["corundum", "beryl", "garnet_pyralspite", "tanzanite"];
    let (mut total, mut within_1) = (0usize, 0usize);
    let mut worst: Vec<String> = Vec::new();
    for host in hosts {
        for _ in 0..50 {
            let recipe = random_recipe(host, &mut rng, 3, 0.5).expect("host has items");
            let target = colors_of(&recipe, REF_MM, Illuminant::D65).unpolarised.lab;
            let res = solve_physics(cat(), host, target, None, REF_MM_F32, &[]);
            total += 1;
            within_1 += usize::from(res.delta_e <= 1.0);
            assert_eq!(
                res.reachable,
                res.delta_e <= 1.0,
                "{host}: reachable flag disagrees with dE {}",
                res.delta_e
            );
            if res.delta_e > 1.0 {
                worst.push(format!(
                    "{host} {:?} -> dE {:.2}",
                    recipe.entries, res.delta_e
                ));
            }
            assert!(
                res.delta_e <= 2.0,
                "{host}: round trip of {:?} left dE {:.2} > 2.0",
                recipe.entries,
                res.delta_e
            );
        }
    }
    assert_eq!(total, 200);
    assert!(
        within_1 * 100 >= total * 95,
        "only {within_1}/{total} round trips reach dE <= 1.0; misses: {worst:#?}"
    );
}

/// Criterion 3: a ruby pick (L* 40, a* 55, b* 30 +- 5) yields a Cr-dominated recipe
/// (Cr share of alpha >= 70 %), and so do the GIA ruby colors.
///
/// `data_version` 4: a ruby with b* +30 (orange-red) is not reachable by the GIA Cr3+ cross-section
/// (b* turns positive only at 1500+ ppma*cm; the solver leaves dE00 12-16 and the closest recipe
/// is still Cr 96-99 %); the GIA Appendix 2 unpolarised D65 colors at 100-375 ppma*cm (purplish
/// red, b* -10 to -19) are reached with Cr alone to dE00 0.7, 1.1 and 2.7 (the GIA "both" column is
/// itself 2.5 off its cross-section at 375 ppma*cm), so the limit is 3.
#[test]
fn criterion_3_ruby_pick_is_cr_dominated() {
    for (da, db) in [
        (0.0, 0.0),
        (5.0, 5.0),
        (5.0, -5.0),
        (-5.0, 5.0),
        (-5.0, -5.0),
    ] {
        let target = [40.0, 55.0 + da, 30.0 + db];
        let res = solve_physics(cat(), "corundum", target, None, REF_MM_F32, &[]);
        let share = alpha_share(&res.recipe, "Cr");
        assert!(
            share >= 0.7,
            "ruby pick {target:?}: Cr share of alpha {share:.2} (recipe {:?}, dE {:.2})",
            res.recipe.entries,
            res.delta_e
        );
    }
    let g = table("gia_cr3plus_sample1110_sigma");
    for a in g.appendix2().iter().filter(|a| {
        a.spec == "both" && matches!(a.illuminant, Illuminant::D65) && a.ppma_cm <= 375.0
    }) {
        let res = solve_physics(cat(), "corundum", a.lab, None, REF_MM_F32, &[]);
        let share = alpha_share(&res.recipe, "Cr");
        assert!(
            share >= 0.7 && res.delta_e <= 3.0,
            "GIA ruby {:?}: Cr share {share:.2}, dE {:.2} (recipe {:?})",
            a.lab,
            res.delta_e,
            res.recipe.entries
        );
    }
}

/// Criterion 3: an unreachable green reports its dE and `reachable = false`; locks are honoured
/// to 1e-12.
#[test]
fn criterion_3_unreachable_green_and_locks() {
    let green = [95.0, -90.0, 90.0];
    let res = solve_physics(cat(), "corundum", green, None, REF_MM_F32, &[]);
    assert!(!res.reachable);
    assert!(
        res.delta_e > 1.0,
        "unreachable pick reports dE {:.2}",
        res.delta_e
    );

    let target = colors_of(
        &recipe_of("corundum", &[("Cr", 0.3)]),
        REF_MM,
        Illuminant::D65,
    )
    .unpolarised
    .lab;
    for locked in [("Cr".to_string(), 0.22), ("Fe".to_string(), 123.456)] {
        let res = solve_physics(
            cat(),
            "corundum",
            target,
            None,
            REF_MM_F32,
            std::slice::from_ref(&locked),
        );
        assert!(
            (res.recipe.amount(&locked.0) - locked.1).abs() < 1e-12,
            "lock {} changed to {}",
            locked.0,
            res.recipe.amount(&locked.0)
        );
    }
}

/// Criterion 3: color-change target pair (D65 + 3200 K) for alexandrite: both dE <= 3.
#[test]
fn criterion_3_color_change_dual_target() {
    let source = recipe_of("chrysoberyl", &[("Cr", 0.12), ("Fe", 0.05)]);
    let (tensor, _) = resolve(&source, cat()).expect("resolves");
    let d65 = body_colors(&tensor, REF_MM, Illuminant::D65)
        .unpolarised
        .lab;
    let a = body_colors(&tensor, REF_MM, Illuminant::Planckian(3200.0))
        .unpolarised
        .lab;
    let req = SolveRequest {
        target_a_lab: Some(a),
        ..SolveRequest::new("chrysoberyl", d65, REF_MM_F32)
    };
    let never = std::sync::atomic::AtomicBool::new(false);
    let res = solve_physics_with(cat(), &req, &never).expect("not cancelled");
    let de_a = res.delta_e_a.expect("dual target reports dE_A");
    assert!(res.delta_e <= 3.0, "alexandrite dE_D65 {:.2}", res.delta_e);
    assert!(de_a <= 3.0, "alexandrite dE_A {de_a:.2}");
}

// ---------------------------------------------------------------------------------------------
// Criterion 4
// ---------------------------------------------------------------------------------------------

/// The host with the most selectable items.
fn largest_host() -> &'static str {
    cat()
        .hosts
        .iter()
        .max_by_key(|h| {
            (
                cat().selectable_elements(&h.id).len(),
                std::cmp::Reverse(h.id.clone()),
            )
        })
        .map(|h| h.id.as_str())
        .expect("hosts")
}

/// Criterion 4: evaluation budget (<= 40 000), wall time (<= 250 ms, release only), stage-1 exit
/// (<= 400 evaluations), |S| = 1 for a pure-Cr target, and determinism.
#[test]
fn criterion_4_budget_time_stage_exit_and_determinism() {
    let host = largest_host();
    let ids = cat().selectable_elements(host);
    assert!(ids.len() >= 6, "{host} offers {} items", ids.len());

    // Worst case: an unreachable target forces every stage over all subsets.
    let hard = [60.0, -70.0, 70.0];
    let start = Instant::now();
    let worst = solve_physics(cat(), host, hard, None, REF_MM_F32, &[]);
    let elapsed = start.elapsed();
    assert!(
        worst.evals <= MAX_EVALS,
        "{host}: {} evaluations",
        worst.evals
    );
    eprintln!(
        "{host}: worst case {} evals, {elapsed:?}{}",
        worst.evals,
        if cfg!(debug_assertions) {
            " (timing asserted in release only)"
        } else {
            ""
        }
    );
    if !cfg!(debug_assertions) {
        assert!(
            elapsed.as_millis() <= 250,
            "{host}: worst-case solve took {elapsed:?} ({} evals)",
            worst.evals
        );
    }

    // Stage 1 exits after <= 400 evaluations for a size-1 target (corundum and the largest host).
    for (host, entry) in [("corundum", "Cr"), (host, ids[0].as_str())] {
        let host_data = cat().host(host).expect("host");
        let amount = 0.3 * host_data.element_conc_max(entry);
        let target = colors_of(
            &recipe_of(host, &[(entry, amount)]),
            REF_MM,
            Illuminant::D65,
        )
        .unpolarised
        .lab;
        let res = solve_physics(cat(), host, target, None, REF_MM_F32, &[]);
        assert!(res.delta_e <= 1.0, "{host}/{entry}: dE {:.3}", res.delta_e);
        assert!(
            res.evals <= 400,
            "{host}/{entry}: {} evaluations before exiting stage 1",
            res.evals
        );
        assert_eq!(
            res.recipe.entries.len(),
            1,
            "{host}/{entry}: {:?}",
            res.recipe.entries
        );
    }

    // Pure Cr in corundum: |S| = 1 and it is Cr.
    let target = colors_of(
        &recipe_of("corundum", &[("Cr", 0.25)]),
        REF_MM,
        Illuminant::D65,
    )
    .unpolarised
    .lab;
    let res = solve_physics(cat(), "corundum", target, None, REF_MM_F32, &[]);
    let solved: Vec<&str> = res.recipe.entries.iter().map(|e| e.id.as_str()).collect();
    assert_eq!(
        solved,
        ["Cr"],
        "a pure-Cr target must never select more items"
    );

    // Determinism: identical input, identical bits.
    let again = solve_physics(cat(), "corundum", target, None, REF_MM_F32, &[]);
    assert_eq!(res.delta_e.to_bits(), again.delta_e.to_bits());
    assert_eq!(res.recipe, again.recipe);
    for (a, b) in res.recipe.entries.iter().zip(&again.recipe.entries) {
        assert_eq!(a.amount.to_bits(), b.amount.to_bits());
    }
    assert_eq!(
        res.achieved.lab.map(f64::to_bits),
        again.achieved.lab.map(f64::to_bits)
    );
}

// ---------------------------------------------------------------------------------------------
// Criterion 5
// ---------------------------------------------------------------------------------------------

/// Criterion 5: for 50 random recipes per host the swatch from `body_colors(resolved)` equals the
/// swatch recomputed from the stored `resolved_bands` bitwise, and the solver's reported
/// `achieved` is within 0.01 of the swatch of its stored recipe.
#[test]
fn criterion_5_wysiwyg_for_50_random_recipes_per_host() {
    let mut rng = Rng(0x57A7_1C05);
    for host in &cat().hosts {
        if cat().selectable_elements(&host.id).is_empty() {
            continue;
        }
        for _ in 0..50 {
            let recipe = random_recipe(&host.id, &mut rng, 3, 1.0).expect("items");
            let (tensor, _) = resolve(&recipe, cat()).expect("resolves cleanly");
            let swatch = body_colors(&tensor, REF_MM, Illuminant::D65);
            let stored = ResolvedBands::from_tensor(&tensor).to_tensor();
            let from_bands = body_colors(&stored, REF_MM, Illuminant::D65);
            assert_eq!(
                swatch.unpolarised.lab.map(f64::to_bits),
                from_bands.unpolarised.lab.map(f64::to_bits),
                "{}: swatch of resolved vs stored bands differ bitwise ({:?})",
                host.id,
                recipe.entries
            );
            assert_eq!(
                delta_e_2000(swatch.unpolarised.lab, from_bands.unpolarised.lab),
                0.0
            );

            let res = solve_physics(
                cat(),
                &host.id,
                swatch.unpolarised.lab,
                None,
                REF_MM_F32,
                &[],
            );
            let solved = colors_of(&res.recipe, REF_MM, Illuminant::D65).unpolarised;
            let de = delta_e_2000(res.achieved.lab, solved.lab);
            assert!(
                de < 0.01,
                "{}: solver-reported achieved vs swatch dE {de} ({:?})",
                host.id,
                res.recipe.entries
            );
            let stored_swatch = body_colors(
                &res.recipe.resolved_bands.to_tensor(),
                REF_MM,
                Illuminant::D65,
            );
            assert!(
                delta_e_2000(res.achieved.lab, stored_swatch.unpolarised.lab) < 0.01,
                "{}: achieved vs stored resolved_bands",
                host.id
            );
        }
    }
}

/// Hosts allowed to drop or merge bands at `conc_typical` (explicit allow-list, §3.3) with the
/// dE00 their pinned `BandsDropped`/`BandsMerged` must stay below. The last four host ids
/// overlap 8-20 chromophores at `conc_typical`; their merges (cheapest color error first, see
/// `resolve::budget_bands`) measure dE00 0.09-0.46 after the research §5 pruning. Beryl's
/// all-elements recipe gained the two Mn chromophores when they got a `conc_typical` on
/// 2026-10-12 (red beryl 0.1-0.35 wt% MnO is a strong absorber at 560 nm): its rays then carry
/// 14 (ordinary) and 13 (extraordinary) bands, and the merges (six and five) measure up to
/// dE00 1.74 for that kitchen-sink stone (Cr + V + Fe + Mn at once, nearly opaque at 5 mm), so
/// its bound is 2.
const BAND_DROP_ALLOW_LIST: &[(&str, f64)] = &[
    ("beryl", 2.0),
    ("corundum", 1.0),
    ("yag", 1.0),
    ("cubic_zirconia", 1.0),
    ("garnet_pyralspite", 1.0),
    ("glass_silicate", 1.0),
    ("spinel", 1.0),
];

/// The allow-listed dE00 bound of `host_id`, `None` when the host may not budget at all.
fn band_drop_allowance(host_id: &str) -> Option<f64> {
    BAND_DROP_ALLOW_LIST
        .iter()
        .find(|(id, _)| *id == host_id)
        .map(|(_, limit)| *limit)
}

/// All elements of `host` at the middle of their `conc_typical` range, resolved: the band
/// budget warnings and the tensor.
fn conc_typical_resolve(host: &HostData) -> (AbsorptionTensor, Vec<ResolveWarning>) {
    let mut recipe = ColorRecipe::new(&host.id, cat().data_version);
    for id in cat().selectable_elements(&host.id) {
        let Some(unit) = host.element_unit(&id) else {
            // End member: the middle of its typical range.
            let typical = host
                .chromophores
                .iter()
                .filter(|c| c.end_member.as_deref() == Some(id.as_str()))
                .filter_map(|c| c.conc_typical)
                .map(|t| f64::midpoint(t[0], t[1]))
                .fold(0.0, f64::max);
            recipe.set_amount(&id, typical.min(0.3));
            continue;
        };
        let n_in = host.n_site_for_unit(&unit);
        let typical = host
            .chromophores
            .iter()
            .filter(|c| c.is_offered() && c.elements.iter().any(|e| e == &id))
            .filter_map(|c| {
                c.conc_typical
                    .map(|t| f64::midpoint(t[0], t[1]) * host.n_site_for_unit(&c.conc_unit) / n_in)
            })
            .fold(0.0, f64::max);
        recipe.set_amount(&id, typical.min(host.element_conc_max(&id)));
    }
    resolve(&recipe, cat()).expect("resolves cleanly")
}

/// The over-budget findings of `hosts`: bands per mode above 8, or dropped or merged bands
/// whose dE00 reaches the host's allowance (or a host on no allow-list).
fn band_budget_violations<'a>(hosts: impl Iterator<Item = &'a HostData>) -> Vec<String> {
    let mut violations = Vec::new();
    for host in hosts {
        let (tensor, warnings) = conc_typical_resolve(host);
        for (mode, bands) in [
            ("o", Some(&tensor.o_ray)),
            ("e", Some(&tensor.e_ray)),
            ("beta", tensor.beta_ray.as_ref()),
        ] {
            if bands.is_some_and(|b| b.len() > 8) {
                violations.push(format!("{}: {mode} has more than 8 bands", host.id));
            }
        }
        for w in &warnings {
            let (what, mode, count, delta_e) = match w {
                ResolveWarning::BandsDropped {
                    mode,
                    count,
                    delta_e,
                } => ("drops", mode, count, delta_e),
                ResolveWarning::BandsMerged {
                    mode,
                    count,
                    delta_e,
                } => ("merges", mode, count, delta_e),
                _ => continue,
            };
            if !band_drop_allowance(&host.id).is_some_and(|limit| *delta_e < limit) {
                violations.push(format!(
                    "{}: {what} {count} {mode} band(s) (dE00 {delta_e:.2}), not allow-listed with that dE",
                    host.id
                ));
            }
        }
    }
    violations
}

/// Criterion 5: <= 8 bands per mode for `conc_typical` recipes (allow-list with a pinned
/// per-host dE00 bound) for every host that needs no budgeting beyond the data (no warning at
/// all).
#[test]
fn criterion_5_conc_typical_fits_the_band_budget() {
    let within: Vec<&HostData> = cat()
        .hosts
        .iter()
        .filter(|h| band_drop_allowance(&h.id).is_none())
        .collect();
    let v = band_budget_violations(within.iter().copied());
    assert!(v.is_empty(), "{v:#?}");
    for host in within {
        let (_, warnings) = conc_typical_resolve(host);
        assert!(
            !warnings.iter().any(|w| matches!(
                w,
                ResolveWarning::BandsDropped { .. } | ResolveWarning::BandsMerged { .. }
            )),
            "{} is not allow-listed but needs budgeting: {warnings:?}",
            host.id
        );
    }
}

/// Criterion 5: the same for ALL hosts.
#[test]
fn criterion_5_conc_typical_fits_the_band_budget_for_all_hosts() {
    let v = band_budget_violations(cat().hosts.iter());
    assert!(v.is_empty(), "{v:#?}");
}

// ---------------------------------------------------------------------------------------------
// Criterion 6
// ---------------------------------------------------------------------------------------------

#[test]
fn criterion_6_stone_size_default_and_resolved_scale() {
    assert_eq!(effective_stone_width_mm(0.0), 7.0);
    assert_eq!(effective_stone_width_mm(5.0), 5.0);

    let model_width = measure_model_width(&StandardGemCuts::standard_round_brilliant())
        .expect("round brilliant measures");
    let base = GemMaterial::all_materials()
        .into_iter()
        .find(|m| m.name == "Ruby")
        .expect("Ruby built-in");
    let overrides = MaterialOverrides {
        inclusion_sigma_s: 0.0,
        c_axis_override: None,
        edge_rounding_radius: 0.0,
        stone_width_mm: 0.0,
    };
    // A physics recipe is per millimetre: the 7 mm default applies while no size is set.
    let physics = apply_material_overrides(
        base.clone()
            .with_chromophore_absorption(base.absorption.clone()),
        &overrides,
        Some(model_width),
    );
    let expected = 7.0 / model_width;
    let got = f64::from(physics.absorption_path_scale);
    assert!(
        ((got - expected) / expected).abs() < 1e-6,
        "scale {got} vs 7/model_width {expected}"
    );
    // A per-model-unit material (the built-in table) with no size set gets exactly the face-up
    // calibration `1 / K` (its swatch colour over `K` units of path), whatever the model width.
    let plain = apply_material_overrides(base, &overrides, Some(model_width));
    assert_eq!(
        plain.absorption_path_scale.to_bits(),
        (1.0 / crate::render_setup::MODEL_UNIT_FACE_UP_PATH).to_bits()
    );
}

/// The built-in materials' own unpolarised color at the 5 mm reference path (bands scaled from
/// model units to mm with W = 7 mm and the measured model width) and the best physics recipe the
/// solver finds for it, with the treatments of the host allowed: `(name, target Lab, solve)`.
fn builtin_vs_recipe() -> Vec<(&'static str, [f64; 3], SolveResult)> {
    let model_width = measure_model_width(&StandardGemCuts::standard_round_brilliant())
        .expect("round brilliant measures");
    let per_mm = model_width / f64::from(effective_stone_width_mm(0.0));
    #[expect(
        clippy::cast_possible_truncation,
        reason = "band peaks are stored as f32; the scaled wavelength is in the normal float range"
    )]
    let scaled = |bands: &[AbsorptionBand]| -> Vec<AbsorptionBand> {
        bands
            .iter()
            .map(|b| AbsorptionBand {
                peak: (f64::from(b.peak) * per_mm) as f32,
                ..*b
            })
            .collect()
    };
    let names = [
        "Ruby",
        "Sapphire",
        "Emerald",
        "Tanzanite",
        "Peridot",
        "Amethyst",
        "Aquamarine",
        "Citrine",
    ];
    let materials = GemMaterial::all_materials();
    names
        .into_iter()
        .map(|name| {
            let m = materials
                .iter()
                .find(|m| m.name == name)
                .expect("built-in exists");
            let host = cat().host_for_material(name).expect("host for built-in");
            let t = &m.absorption;
            let tensor = if t.is_pleochroic {
                AbsorptionTensor::uniaxial(scaled(&t.o_ray), scaled(&t.e_ray))
            } else {
                AbsorptionTensor::isotropic(scaled(&t.o_ray))
            };
            let tensor = t.beta_ray.as_ref().map_or(tensor, |beta| {
                AbsorptionTensor::biaxial(scaled(&t.o_ray), scaled(beta), scaled(&t.e_ray))
            });
            let target = body_colors(&tensor, REF_MM, Illuminant::D65)
                .unpolarised
                .lab;
            let all: Vec<String> = host.treatments.iter().map(|t| t.id.clone()).collect();
            let req = SolveRequest {
                treatments: TreatmentPolicy::Allow(&all),
                ..SolveRequest::new(&host.id, target, REF_MM_F32)
            };
            let res = solve_physics_with(cat(), &req, &std::sync::atomic::AtomicBool::new(false))
                .expect("not cancelled");
            eprintln!(
                "{name}: target {target:.1?} -> dE {:.2} {:?} {:?} lab {:.1?}",
                res.delta_e, res.recipe.entries, res.recipe.treatments, res.achieved.lab
            );
            (name, target, res)
        })
        .collect()
}

/// Criterion 6 / §2.4 item 4 against the PRIMARY-data recipes: the built-in colors (unchanged,
/// bit for bit) against the best physics recipe of `data_version` 4, pinned per built-in.
///
/// Re-pinned on 2026-10-19 after the measured color refits (owner decision; values from the
/// owner's Windows run). Agreement (dE00 <= 3.5): Peridot 0.15, Sapphire 0.22, Aquamarine
/// 0.75, Ruby 1.00 (a Cr + Fe recipe), Emerald 1.01, Citrine 2.41 (2.69 until `data_version`
/// 5 removed the fixed 10 % Fe2+ share that tinted every iron quartz). The remaining two
/// really differ from the primary-data recipes: Tanzanite 7.27 (the heated built-in is still
/// more saturated than the V3+ recipe at conc_max) and Amethyst 11.25 (built-in C* 82, h 321
/// against the Caltech-shaped [FeO4]0 recipe C* 44 at conc_max).
#[test]
fn criterion_6_built_in_colors_vs_primary_data_recipes_are_pinned() {
    let pinned = [
        ("Ruby", 1.00),
        ("Sapphire", 0.22),
        ("Emerald", 1.01),
        ("Tanzanite", 7.27),
        ("Peridot", 0.15),
        ("Amethyst", 11.25),
        ("Aquamarine", 0.75),
        ("Citrine", 2.41),
    ];
    let results = builtin_vs_recipe();
    for ((name, _, res), (pin_name, pin)) in results.iter().zip(pinned) {
        assert_eq!(*name, pin_name);
        assert!(
            (res.delta_e - pin).abs() < 0.25,
            "{name}: built-in vs recipe dE00 {:.2}, pinned {pin}",
            res.delta_e
        );
    }
    let agree = results.iter().filter(|(_, _, r)| r.delta_e <= 3.5).count();
    assert_eq!(
        agree, 6,
        "Ruby, Sapphire, Emerald, Peridot, Aquamarine and Citrine agree to 3.5"
    );
    let ruby = &results[0].2;
    assert!(ruby.delta_e <= 1.25 && ruby.recipe.entries.iter().any(|e| e.id == "Cr"));
}

/// Criterion 6 / §2.4 item 4: a plausible recipe reproduces the built-in color (bands scaled
/// from model units to mm with W = 7 mm and the measured model width) within dE00 <= 3 for at
/// least 6 of the 8 built-ins. Ignored until the measured color refits of 2026-10-19; since
/// then six reproduce within 3 (Ruby, Sapphire, Emerald, Peridot, Aquamarine, Citrine), and
/// Tanzanite and Amethyst are the two that differ.
#[test]
fn criterion_6_built_in_sanity_at_least_6_of_8() {
    let results = builtin_vs_recipe();
    let passed = results.iter().filter(|(_, _, r)| r.delta_e <= 3.0).count();
    let report: Vec<String> = results
        .iter()
        .map(|(n, _, r)| format!("{n}: dE {:.2}", r.delta_e))
        .collect();
    assert!(
        passed >= 6,
        "only {passed}/8 built-ins reproduced within dE00 3: {report:?}"
    );
}

// ---------------------------------------------------------------------------------------------
// Criterion 10
// ---------------------------------------------------------------------------------------------

/// Pinned band fixture of corundum + 0.20 wt% Cr2O3 at 5 mm, `data_version` 5 (the corundum
/// bands are those of version 4; version 5 changed the glass and quartz valence states only):
/// `(centre nm, fwhm cm^-1, peak mm^-1)` per mode, sorted by centre.
const PINNED_DATA_VERSION: u32 = 5;
const PINNED_O: [(f32, f32, f32); 5] = [
    (345.0, 3687.4, 0.032_346),
    (404.36, 3286.7002, 1.02598),
    (538.39, 2874.5999, 0.41608),
    (559.87, 1583.9, 0.7246),
    (693.3, 55.0, 0.056_016),
];
const PINNED_E: [(f32, f32, f32); 4] = [
    (345.0, 12000.0, 0.021_274),
    (400.44, 3257.5, 1.7583),
    (537.81, 2416.0, 0.34726),
    (693.3, 55.0, 0.007_282_08),
];

#[test]
fn criterion_10_pinned_snapshot() {
    assert_eq!(
        cat().data_version,
        PINNED_DATA_VERSION,
        "changing a catalogue number requires bumping data_version and this fixture together"
    );
    let recipe = recipe_of("corundum", &[("Cr", 0.20)]);
    let (tensor, _) = resolve(&recipe, cat()).expect("resolves cleanly");
    let rel = |a: f32, b: f32| (f64::from(a) - f64::from(b)).abs() / f64::from(b).abs();
    for (mode, bands, pinned) in [
        ("o", &tensor.o_ray, &PINNED_O[..]),
        ("e", &tensor.e_ray, &PINNED_E[..]),
    ] {
        assert_eq!(bands.len(), pinned.len(), "{mode}: band count");
        for (b, (centre, fwhm, peak)) in bands.iter().zip(pinned.iter().copied()) {
            assert!(
                rel(b.center_nm, centre) < 1e-6,
                "{mode} centre {} vs {centre}",
                b.center_nm
            );
            assert!(
                rel(b.width_nm * 2.354_82, fwhm) < 1e-6,
                "{mode} width {} vs {fwhm}",
                b.width_nm * 2.354_82
            );
            assert!(rel(b.peak, peak) < 1e-6, "{mode} peak {} vs {peak}", b.peak);
        }
    }
}

// ---------------------------------------------------------------------------------------------
// Criterion 11
// ---------------------------------------------------------------------------------------------

#[test]
fn criterion_11_offer_rule_and_band_limits() {
    let offered = |host: &str| cat().selectable_elements(host);
    let cor = offered("corundum");
    for el in ["Cr", "Fe", "Ti", "V"] {
        assert!(
            cor.iter().any(|e| e == el),
            "corundum must offer {el}: {cor:?}"
        );
    }
    assert!(offered("chrysoberyl").iter().any(|e| e == "Cr"));
    // Diamond offers its visible N centres: N drives an offered chromophore with a visible band.
    let diamond = cat().host("diamond").expect("diamond");
    assert!(offered("diamond").iter().any(|e| e == "N"));
    assert!(diamond.chromophores.iter().any(|c| {
        c.is_offered()
            && c.elements.iter().any(|e| e == "N")
            && c.usable_bands()
                .any(|b| (380.0..=780.0).contains(&b.centre_nm))
    }));

    for host in &cat().hosts {
        for chromo in &host.chromophores {
            for band in &chromo.bands {
                assert!(
                    (250.0..=1200.0).contains(&band.centre_nm),
                    "{}/{}: band at {} nm outside 250-1200 nm",
                    host.id,
                    chromo.id,
                    band.centre_nm
                );
            }
        }
    }
}

// ---------------------------------------------------------------------------------------------
// Criterion 12
// ---------------------------------------------------------------------------------------------

/// A corundum-like host whose Cr3+ rows carry the Report's 25.9 and the correct 51 cm^-1/wt%.
const CR_FIXTURE: &str = r#"
data_version = 1

[[host]]
id = "ruby_fixture"
name = "Ruby fixture"
formula = "Al2O3"
optical = "uniaxial"
materials = ["Fixture"]
density_g_cm3 = 3.98
formula_weight_g_mol = 101.96
sites_per_fu = 2.0
cation_site_density_cm3 = 4.701475e+22
total_atom_density_cm3 = 1.175369e+23

  [[host.chromophore]]
  id = "Cr3+ (Report 25.9)"
  kind = "ion"
  elements = ["Cr"]
  conc_unit = "wt_pct_oxide:Cr2O3"
    [[host.chromophore.band]]
    centre_nm = 556.0
    fwhm_cm1 = 2800.0
    shape = "gaussian_energy"
    peak_coeff = 25.9
    pol = { o = 1.0, e = 0.75 }
    sigma_cm2 = 1.62e-19

  [[host.chromophore]]
  id = "Cr3+ (51)"
  kind = "ion"
  elements = ["Cr"]
  conc_unit = "wt_pct_oxide:Cr2O3"
    [[host.chromophore.band]]
    centre_nm = 556.0
    fwhm_cm1 = 2800.0
    shape = "gaussian_energy"
    peak_coeff = 51.0
    pol = { o = 1.0, e = 0.75 }
    sigma_cm2 = 1.62e-19
"#;

#[test]
fn criterion_12_consistency_checks() {
    // The Report's 25.9 is flagged by the loader, 51 passes (loaded through from_toml).
    let fixture = ChromophoreCatalogue::from_toml(CR_FIXTURE).expect("fixture loads");
    let host = fixture.host("ruby_fixture").expect("fixture host");
    let band_of = |id: &str| {
        &host
            .chromophores
            .iter()
            .find(|c| c.id == id)
            .expect(id)
            .bands[0]
    };
    assert!(
        band_of("Cr3+ (Report 25.9)").suspect.is_some(),
        "25.9 must be flagged suspect"
    );
    assert!(
        band_of("Cr3+ (51)").suspect.is_none(),
        "51 must pass the 2x check"
    );
    let offered = |id: &str| {
        host.chromophores
            .iter()
            .find(|c| c.id == id)
            .expect(id)
            .is_offered()
    };
    assert!(
        !offered("Cr3+ (Report 25.9)") && offered("Cr3+ (51)"),
        "the flagged row is hidden, the 51 row offered"
    );

    // The shipped catalogue carries the 51 row and flags GR1 and [h.-Fe3+].
    let corundum = cat().host("corundum").expect("corundum");
    let cr3 = corundum
        .chromophores
        .iter()
        .find(|c| c.id == "Cr3+")
        .expect("Cr3+");
    assert!(cr3.suspect.is_none() && cr3.bands.iter().all(|b| b.suspect.is_none()));
    // data_version 4: the Cr3+ bands are the GIA cross-section (each carries its sigma, which the
    // loader checks against peak_coeff within 2x). The U band region (E-perp-c, 560 nm) carries
    // the Dodd/Dubinsky anchor sigma 1.62e-19 cm^2 = 51.1 cm^-1/wt% within 8 %.
    let n_wt = corundum.n_site_for_unit("wt_pct_oxide:Cr2O3");
    let u_560: f64 = cr3
        .usable_bands()
        .filter(|b| b.pol.get("o").copied().unwrap_or(0.0) > 0.0)
        .map(|b| {
            let d = 1e7 / 560.0 - 1e7 / b.centre_nm;
            b.peak_coeff.expect("peak")
                * (-4.0 * std::f64::consts::LN_2 * (d / b.fwhm_cm1.expect("fwhm")).powi(2)).exp()
        })
        .sum();
    assert!(
        (u_560 / (1.62e-19 * n_wt) - 1.0).abs() < 0.08,
        "the U band carries ~51 cm^-1/wt% at 560 nm: {u_560}"
    );
    assert!(cr3.bands.iter().all(|b| b.sigma_cm2.is_some()));

    let flagged = |host: &str, needle: &str| {
        let c = cat()
            .host(host)
            .expect("host")
            .chromophores
            .iter()
            .find(|c| c.id.contains(needle))
            .unwrap_or_else(|| panic!("{host}/{needle} exists"));
        c.suspect.is_some() || c.bands.iter().any(|b| b.suspect.is_some())
    };
    assert!(flagged("diamond", "GR1"), "GR1 must be flagged");
    assert!(
        flagged("corundum", "h\u{2022}"),
        "[h.-Fe3+] must be flagged"
    );
}

// ---------------------------------------------------------------------------------------------
// Criterion 13
// ---------------------------------------------------------------------------------------------

/// Criterion 13: ruby R1 (0.7 nm FWHM) at 1 nm differs from a 0.1 nm reference integral by
/// dE00 < 0.3.
#[test]
fn criterion_13_narrow_lines_ruby_r1() {
    let (centre, sigma, peak) = (694.3, 0.7 / 2.354_82, 3.0);
    let alpha = |lam: f64| peak * (-0.5 * ((lam - centre) / sigma).powi(2)).exp();
    let col_1nm = body_color(alpha, REF_MM, Illuminant::D65);

    let step = 0.1;
    let (mut xyz, mut white) = ([0.0; 3], [0.0; 3]);
    #[expect(
        clippy::cast_possible_truncation,
        clippy::cast_sign_loss,
        reason = "the rounded step count (4000) is a small positive number"
    )]
    let steps = ((780.0_f64 - 380.0) / step).round() as usize;
    for i in 0..=steps {
        let lambda = 380.0 + i as f64 * step;
        #[expect(
            clippy::cast_possible_truncation,
            reason = "the CMF takes f32 nanometres; lambda is 380..=780"
        )]
        let cmf = cie_1931_cmf(lambda as f32);
        let s = Illuminant::D65.spectral_power(lambda);
        let t = (-alpha(lambda) * REF_MM).exp();
        for k in 0..3 {
            xyz[k] += t * s * f64::from(cmf[k]) * step;
            white[k] += s * f64::from(cmf[k]) * step;
        }
    }
    let norm = white[1];
    let lab_ref = xyz_to_lab(xyz.map(|v| v / norm), white.map(|v| v / norm));
    let de = delta_e_2000(col_1nm.lab, lab_ref);
    assert!(de < 0.3, "1 nm swatch vs 0.1 nm reference: dE00 {de:.4}");
}

// ---------------------------------------------------------------------------------------------
// Criterion 14
// ---------------------------------------------------------------------------------------------

fn color_change_at(recipe: &ColorRecipe, path_mm: f64) -> f64 {
    let (tensor, _) = resolve(recipe, cat()).expect("resolves");
    let d65 = body_colors(&tensor, path_mm, Illuminant::D65);
    let a = body_colors(&tensor, path_mm, Illuminant::Planckian(3200.0));
    delta_e_2000(d65.unpolarised.lab, a.unpolarised.lab)
}

fn color_change(recipe: &ColorRecipe) -> f64 {
    color_change_at(recipe, REF_MM)
}

/// Criterion 14: an alexandrite-like V + Cr corundum recipe has dE_cc (D65 vs 3200 K) > 8.
///
/// `data_version` 4: the V3+ cross-section is GIA's, which gives a D65 vs A color change of
/// 8-25 dE00 (peak in E_par_c at a few hundred ppma*cm; the report-1 table's 2.9-5.8 is wrong).
/// V 2000 ppm of Al sites (800 ppma, conc_max) over 5 mm is 400 ppma*cm: dE_cc 17.1 alone; the
/// color-change stone is V-dominated, so Cr 0.02 wt% (54 ppma) leaves 13.0. (The data_version
/// 3 recipe Cr 0.1 wt% + V 2000 ppm is now a ruby-like 5.65: Cr dominates.)
#[test]
fn criterion_14_alexandrite_like_color_change() {
    let alexandrite_like = recipe_of("corundum", &[("Cr", 0.02), ("V", 2000.0)]);
    let cc = color_change(&alexandrite_like);
    assert!(cc > 8.0, "V+Cr corundum dE_cc {cc:.2} (needs > 8)");
    let v_alone = color_change(&recipe_of("corundum", &[("V", 2000.0)]));
    assert!(
        (8.0..=25.0).contains(&v_alone),
        "V corundum dE_cc {v_alone:.2}"
    );
    // Cr takes the color change away again.
    assert!(color_change(&recipe_of("corundum", &[("Cr", 0.1), ("V", 2000.0)])) < cc);
}

/// Criterion 14: a Cr-only ruby recipe has dE_cc < 7 (D65 vs 3200 K).
///
/// `data_version` 4: the plan's "< 5" is replaced by the GIA-derived value. The color of the
/// GIA Cr3+ cross-section (2 E_perp + E_par) / 3 changes by 2.4 / 2.8 / 5.0 / 5.9 dE00 at
/// 0.05 / 0.1 / 0.3 / 0.5 wt% over 5 mm (ruby is darker and redder under incandescent light), not
/// "no change": the recipe is within 0.4 of it and the limit 7 keeps a Cr-only ruby apart from
/// the color-change recipes above (> 8).
#[test]
fn criterion_14_ruby_has_no_color_change() {
    let g = table("gia_cr3plus_sample1110_sigma");
    let cols = (g.col("sigma_E_perp_c_cm2"), g.col("sigma_E_par_c_cm2"));
    let n_wt = cat()
        .host("corundum")
        .expect("corundum")
        .n_site_for_unit("wt_pct_oxide:Cr2O3");
    let planck = Illuminant::Planckian(3200.0);
    for cr in [0.05, 0.1, 0.3, 0.5] {
        let cc = color_change(&recipe_of("corundum", &[("Cr", cr)]));
        let d = cr * n_wt * 0.5 / PPMA_CM3;
        // GIA sigma colors under the same two illuminants, over the same 5 mm.
        let n = d * PPMA_CM3;
        let t = |l: f64, col: usize| (-g.at(col, l) * n).exp();
        let mix = |l: f64| -((2.0 * t(l, cols.0) + t(l, cols.1)) / 3.0).ln() / 5.0;
        let gia_cc = delta_e_2000(
            body_color(mix, 5.0, Illuminant::D65).lab,
            body_color(mix, 5.0, planck).lab,
        );
        eprintln!("Cr {cr} wt%: recipe dE_cc {cc:.2}, GIA sigma {gia_cc:.2}");
        assert!(
            (cc - gia_cc).abs() < 0.5,
            "Cr {cr}: {cc:.2} vs GIA {gia_cc:.2}"
        );
        assert!(
            cc < 7.0,
            "Cr-only ruby ({cr} wt%) dE_cc {cc:.2} (needs < 7)"
        );
    }
}