indicatrix 0.9.0

//! Dispersion regression tests: Abbe numbers, `n_d` at the sodium D line, physical
//! bounds at the sampled band edges, and the wavelength-dependent extraordinary
//! index.

use crate::optics::materials::{GemMaterial, OpticalCharacter};

/// Verifies every built-in material's computed Abbe number
/// `V_d = (n_d - 1) / (n_F - n_C)` against a trusted-source reference value, so a
/// future mistyped coefficient (Sellmeier or Cauchy) gets caught instead of
/// silently shipping -- see each material's own doc comment in `all_materials` for
/// the primary/reputable source and the derivation of its target Abbe number.
///
/// Source priority for these targets (re-verified against trusted online sources,
/// which now OUTRANK `GEMSTONE_RENDERING_BLUEPRINT.md` -- the document is used only
/// as a tiebreaker/corroboration, never as the primary authority, because it is a
/// research compilation whose own "Fire Delta n(F-C)" column demonstrably mixes
/// Fraunhofer B-G and F-C conventions row-to-row, e.g. Quartz lists 0.013, the
/// well-known B-G figure, in a column labelled F-C, contradicting that same row's
/// own `V_d`):
///   1. Primary optical measurements (Sellmeier/Cauchy fits from refractiveindex.info
///      and the papers it cites).
///   2. Reputable gemological references (GIA, Mindat, Handbook of Mineralogy,
///      International Gem Society), corroborated by 2+ independent sources where
///      possible, with their "dispersion" figures explicitly treated as Fraunhofer
///      B-G (never plugged directly into an F-C slot).
///   3. The document, as a tiebreaker only.
///
/// Two confidence tiers, reflected in the tolerance:
///   - TIGHT (0.5): Diamond, Sapphire, Ruby, Quartz, Spinel, Cubic Zirconia,
///     Alexandrite, Synthetic Moissanite -- dispersion coefficients transcribed
///     directly from a primary optical-constants paper (Peter 1923; Malitson &
///     Dodge 1972; Ghosh 1999; Tropf & Thomas 1991; Wood & Nassau 1982; Walling
///     1980; Wang 2013 corroborated by Singh 1971), target computed from that exact
///     published formula. Alexandrite and Moissanite CONTRADICT the document here
///     (68.0 -> 73.90 and 21.5 -> 25.94 respectively) -- per the source priority
///     above the primary measurement wins; see each entry's comment in
///     `all_materials` for the two-source corroboration.
///   - LOWER CONFIDENCE (2.5-4.0, roughly 5-12% relative): Zircon, Topaz,
///     Tourmaline, Tanzanite, Emerald -- no primary Sellmeier/Cauchy fit exists in
///     the optics literature for these species (confirmed absent from
///     refractiveindex.info's database), so `n_d` is taken from 2+ corroborating
///     reputable gemological sources and the dispersion shape is a 2-parameter
///     Cauchy fit solved from that `n_d` and a Delta n(F-C) obtained by converting
///     each species' well-corroborated gemological B-G "dispersion" figure via a
///     B-G->F-C ratio (0.579, range 0.569-0.587) computed directly from the 8
///     genuine primary Sellmeier curves above -- not from cross-comparing
///     gemological tables against each other, which is how the earlier 0.591
///     factor was derived. All five CONTRADICT the document (Zircon 28.0 -> 40.96, Topaz 64.0 ->
///     76.83, Tourmaline 55.0 -> 64.58, Tanzanite 45.0 -> 40.21, Emerald's estimate
///     shifts modestly from 69.99 to 70.93 with the refined ratio). Zircon's `n_d`
///     also carried an unrelated, unsourced error, corrected to 1.925 from
///     1.956878 (see its entry's comment) -- critical-angle-relevant since `n_d`
///     drives it directly.
#[test]
fn builtin_material_abbe_numbers_match_published_values() {
    // (name, sourced Abbe number V_d, tolerance)
    let expected: &[(&str, f32, f32)] = &[
        // -- tight tier: primary Sellmeier source for both n_d and V_d --
        ("Diamond", 55.27, 0.5),
        ("Sapphire", 72.27, 0.5),
        ("Ruby", 72.27, 0.5),
        ("Quartz", 69.65, 0.5),
        ("Spinel", 60.63, 0.5),
        ("Cubic Zirconia", 33.52, 0.5),
        ("Alexandrite", 73.90, 0.5),
        ("Synthetic Moissanite", 25.94, 0.5),
        // -- lower-confidence tier: no primary fit exists; n_d from corroborated
        // gemological sources, dispersion shape from a gemological B-G figure
        // converted via the physically-derived 0.579 ratio -- see each entry's
        // comment in `all_materials` --
        ("Zircon", 40.96, 3.0),
        ("Topaz", 76.83, 4.0),
        ("Tourmaline", 64.58, 3.5),
        ("Tanzanite", 40.21, 3.0),
        ("Emerald", 70.93, 3.5),
    ];

    for &(name, published_abbe, tol) in expected {
        let material = GemMaterial::by_name(name)
            .unwrap_or_else(|| panic!("{name} must be a built-in material"));
        let n_f = material.dispersion.evaluate(486.1);
        let n_d = material.dispersion.evaluate(589.3);
        let n_c = material.dispersion.evaluate(656.3);
        let dn = n_f - n_c;
        assert!(
            dn > 0.0,
            "{name}: Delta n(F-C) must be positive (normal dispersion), got {dn}"
        );
        let computed_abbe = (n_d - 1.0) / dn;
        assert!(
            (computed_abbe - published_abbe).abs() <= tol,
            "{name}: computed Abbe V_d={computed_abbe:.2} (n_d={n_d:.5}, n_F={n_f:.5}, n_C={n_c:.5}) \
                 does not match published V_d={published_abbe:.2} within tolerance {tol} -- check the \
                 dispersion coefficients against the source cited in all_materials()"
        );
    }
}
/// Companion to `builtin_material_abbe_numbers_match_published_values`: pins down
/// `n_d` at the sodium D line (589.3nm) for every re-sourced material below, so a
/// future coefficient edit that silently shifts `n_d` gets caught. `n_d` drives
/// the critical angle and therefore brilliance/windowing/extinction, so any
/// intentional change here is significant downstream and must be deliberate, not
/// an accident of refitting the dispersion shape.
///
/// For five of the six re-sourced materials (Synthetic Moissanite, Topaz,
/// Tourmaline, Tanzanite, Alexandrite), `n_d` is unchanged from before -- only the
/// dispersion shape (`V_d`) moved, per
/// `builtin_material_abbe_numbers_match_published_values`.
///
/// Zircon is the exception, and the one material in this file whose `n_d`
/// DID change materially: from 1.956878 to 1.925, a -0.0319 shift, far above
/// the ~0.002 threshold at which critical-angle-driven behaviour (brilliance,
/// windowing, extinction) is expected to move. The old value was never actually
/// sourced from anywhere (not the document, which gives `n_o=1.9250`, nor any
/// external reference found); 1.925 is corroborated by the Handbook of Mineralogy,
/// International Gem Society, and the document itself. See Zircon's comment in
/// `all_materials` for the full derivation. Any test elsewhere in the workspace
/// that hardcodes zircon's brilliance/windowing/pose behaviour will need
/// re-stabilising against this new `n_d`.
#[test]
fn builtin_material_n_d_matches_sourced_values() {
    // (name, sourced n_d, tolerance)
    let expected_n_d: &[(&str, f32, f32)] = &[
        ("Zircon", 1.925, 1e-4),
        ("Synthetic Moissanite", 2.647_434, 1e-4),
        ("Topaz", 1.627_178, 1e-4),
        ("Tourmaline", 1.639_405, 1e-4),
        ("Tanzanite", 1.700_858, 1e-4),
        ("Alexandrite", 1.742_73, 1e-4),
    ];

    for &(name, expected, tol) in expected_n_d {
        let material = GemMaterial::by_name(name)
            .unwrap_or_else(|| panic!("{name} must be a built-in material"));
        let n_d = material.dispersion.evaluate(589.3);
        assert!(
            (n_d - expected).abs() < tol,
            "{name}: n_d={n_d:.6} does not match the sourced value {expected:.6} -- check the \
                 dispersion coefficients against the source cited in all_materials(); if this is a \
                 deliberate re-sourcing, update this test's expected value and flag the n_d change \
                 in the task report if it moves by more than ~0.002 (critical-angle significance)"
        );
    }
}
/// Every built-in material's dispersion curve (Sellmeier or Cauchy) must
/// evaluate to a finite, physically-sane (`n >= 1.0`) index at BOTH ends of this
/// renderer's actual sampled visible band (380nm violet, 780nm red) -- not just
/// near the sodium D line every other dispersion test in this file checks. Most
/// Cauchy fits here are only solved/verified near D and the F/C Fraunhofer lines
/// (486-656nm); this pins the `.max(1.0)` floor added to
/// `DispersionModel::Cauchy::evaluate` (see that variant's own doc comment) so an
/// out-of-fit-range extrapolation can never silently produce `n < 1` or a NaN.
#[test]
fn every_builtin_dispersion_curve_stays_physical_at_the_sampled_band_edges() {
    for material in GemMaterial::all_materials() {
        for lambda_nm in [380.0f32, 780.0] {
            let n = material.dispersion.evaluate(lambda_nm);
            assert!(
                n.is_finite(),
                "{}: dispersion.evaluate({lambda_nm}) must be finite, got {n}",
                material.name
            );
            assert!(
                n >= 1.0,
                "{}: dispersion.evaluate({lambda_nm}) = {n} must be >= 1.0 \
                     (physically-impossible index)",
                material.name
            );
            // The extraordinary ray (when a genuine independent curve is
            // present, or the constant-offset fallback otherwise) must be
            // equally well-behaved.
            if material.optical_character == OpticalCharacter::UniaxialPositive
                || material.optical_character == OpticalCharacter::UniaxialNegative
            {
                let n_e = material.extraordinary_index_at(lambda_nm, n);
                assert!(
                    n_e.is_finite() && n_e >= 1.0,
                    "{}: extraordinary_index_at({lambda_nm}, {n}) = {n_e} must be \
                         finite and >= 1.0",
                    material.name
                );
            }
        }
    }
}
/// Every new built-in must resolve by its own exact name
/// (the same `by_name` round-trip `by_name_round_trips_every_builtin_material`
/// already covers exhaustively -- this test's real job is the tolerance check
/// below) and its `n(589.3nm)` must match the target table below to within
/// 0.002.
#[test]
fn m4_new_species_match_their_target_n_d_within_tolerance() {
    // (name, target n_D, tolerance)
    let expected: &[(&str, f32, f32)] = &[
        ("Aquamarine", 1.577, 0.002),
        ("Morganite", 1.577, 0.002),
        ("Chrysoberyl (Yellow)", 1.746, 0.002),
        ("Amethyst", 1.544, 0.001),
        ("Citrine", 1.544, 0.001),
        ("Pyrope Garnet", 1.714, 0.002),
        ("Almandine Garnet", 1.790, 0.002),
        ("Spessartine Garnet", 1.800, 0.002),
        ("Color-Change Garnet (Pyrope-Spessartine)", 1.760, 0.001),
        ("Grossular Garnet (Tsavorite)", 1.734, 0.002),
        ("Andradite Garnet (Demantoid)", 1.887, 0.002),
        ("Peridot", 1.654, 0.002),
        ("YAG", 1.833, 0.002),
        ("Benitoite", 1.757, 0.002),
        ("Andalusite", 1.634, 0.002),
        ("Opal", 1.45, 0.002),
        ("Glass (N-BK7)", 1.5168, 0.002),
        ("Glass (F2)", 1.620, 0.002),
    ];

    for &(name, target_n_d, tol) in expected {
        let material =
            GemMaterial::by_name(name).unwrap_or_else(|| panic!("{name} must resolve via by_name"));
        assert_eq!(material.name, name, "{name} must resolve to itself exactly");
        let n_d = material.dispersion.evaluate(589.3);
        assert!(
            (n_d - target_n_d).abs() <= tol,
            "{name}: n_d={n_d:.5} does not match target {target_n_d} within {tol}"
        );
    }

    // GGG is checked separately: its target n_D (1.970) is only reproduced to
    // within ~0.002 by construction from this entry's own Cauchy fit (see that
    // entry's comment for why no primary Sellmeier fit was available), so this
    // uses the same tolerance but is called out on its own for clarity.
    let ggg = GemMaterial::by_name("GGG").expect("GGG must resolve via by_name");
    let ggg_n_d = ggg.dispersion.evaluate(589.3);
    assert!(
        (ggg_n_d - 1.970).abs() <= 0.002,
        "GGG: n_d={ggg_n_d:.5} does not match target 1.970 within 0.002"
    );
}

/// GGG's stored dispersion is the F-C convention: the gemological B-G figure 0.045 times
/// the file-wide 0.579 B-G->F-C ratio gives `Delta n(F-C) = 0.026055`, which the Cauchy
/// fit reproduces between the Fraunhofer F (486.1nm) and C (656.3nm) lines. Feeding 0.045
/// straight into the F-C slot would give 0.045 here, 1.7x too dispersive.
#[test]
fn ggg_dispersion_uses_the_converted_f_minus_c_value() {
    let ggg = GemMaterial::by_name("GGG").expect("GGG must resolve via by_name");
    let delta_f_c = ggg.dispersion.evaluate(486.1) - ggg.dispersion.evaluate(656.3);
    assert!(
        (delta_f_c - 0.026_055).abs() <= 1e-3,
        "GGG: n(F)-n(C)={delta_f_c:.5} should be 0.045 * 0.579 = 0.026055 within 1e-3"
    );
}
/// Quartz's genuine per-axis (Ghosh o/e) dispersion must make its
/// extraordinary-ray index actually VARY with wavelength in a way the old
/// constant-offset approximation could not (the offset `n_o(lambda) +
/// birefringence_delta` tracks `n_o`'s own curvature exactly, so `n_e(lambda) -
/// n_o(lambda)` is constant under the old model but need not be under a genuine
/// independent curve); a material with no `uniaxial_extraordinary_dispersion` set
/// (e.g. Sapphire) must keep the exact old constant-delta behaviour, bit-identical.
#[test]
fn quartz_extraordinary_index_has_wavelength_dependent_birefringence_while_legacy_entries_stay_constant()
 {
    let quartz = GemMaterial::by_name("Quartz").expect("Quartz must resolve");
    assert!(
        quartz.uniaxial_extraordinary_dispersion.is_some(),
        "test premise: Quartz must carry a genuine per-axis e-ray dispersion curve"
    );

    let delta_at = |lambda_nm: f32| {
        let n_o = quartz.dispersion.evaluate(lambda_nm);
        let n_e = quartz.extraordinary_index_at(lambda_nm, n_o);
        n_e - n_o
    };
    let delta_380 = delta_at(380.0);
    let delta_780 = delta_at(780.0);
    assert!(
        (delta_380 - delta_780).abs() > 1e-4,
        "Quartz's n_e - n_o must genuinely vary across the visible band: \
             delta(380nm)={delta_380:.6}, delta(780nm)={delta_780:.6}"
    );

    // A legacy entry (no per-axis curve) must keep the OLD constant-delta
    // behaviour exactly: n_e - n_o == birefringence_delta at every wavelength.
    let sapphire = GemMaterial::by_name("Sapphire").expect("Sapphire must resolve");
    assert!(
        sapphire.uniaxial_extraordinary_dispersion.is_none(),
        "test premise: Sapphire must NOT carry a per-axis e-ray dispersion curve"
    );
    for lambda_nm in [380.0f32, 589.3, 780.0] {
        let n_o = sapphire.dispersion.evaluate(lambda_nm);
        let n_e = sapphire.extraordinary_index_at(lambda_nm, n_o);
        assert!(
            (n_e - n_o - sapphire.birefringence_delta).abs() < 1e-6,
            "Sapphire (legacy entry) at {lambda_nm}nm: n_e - n_o must equal the \
                 constant birefringence_delta exactly"
        );
    }
}

/// Rutile's principal indices follow `DeVore` (1951): `n_o(D) = 2.6129`, `n_e(D) =
/// 2.9086`, with Fraunhofer `Delta n(F-C)` of 0.1636 (o) and 0.2072 (e) -- far below
/// the gemological B-G figure (~0.30) the earlier Cauchy fit mistook for F-C -- and the
/// stored `birefringence_delta` equals `n_e(D) - n_o(D)`.
#[test]
fn rutile_matches_devore_indices_and_dispersion() {
    let rutile = GemMaterial::by_name("Rutile").expect("Rutile must be a built-in material");
    let e_ray = rutile
        .uniaxial_extraordinary_dispersion
        .expect("Rutile must carry its own e-ray dispersion curve");

    let n_o = rutile.dispersion.evaluate(589.3);
    let n_e = e_ray.evaluate(589.3);
    assert!((n_o - 2.613).abs() <= 0.002, "n_o(D) = {n_o}");
    assert!((n_e - 2.909).abs() <= 0.002, "n_e(D) = {n_e}");

    let fc_o = rutile.dispersion.evaluate(486.1) - rutile.dispersion.evaluate(656.3);
    let fc_e = e_ray.evaluate(486.1) - e_ray.evaluate(656.3);
    assert!((fc_o - 0.164).abs() <= 0.003, "o-ray F-C = {fc_o}");
    assert!((fc_e - 0.207).abs() <= 0.003, "e-ray F-C = {fc_e}");

    assert!(
        (rutile.birefringence_delta - (n_e - n_o)).abs() <= 1e-3,
        "birefringence_delta {} must equal n_e(D) - n_o(D) = {}",
        rutile.birefringence_delta,
        n_e - n_o
    );
}