spectro_rs/

colorimetry.rs

use crate::{Illuminant, Observer};

/// CIE 1931 2-degree Standard Observer CMFs (380-780nm, 10nm steps)
pub const X_BAR_2: [f32; 41] = [
    0.0014, 0.0042, 0.0143, 0.0435, 0.1344, 0.2839, 0.3483, 0.3362, 0.2908, 0.1954, 0.0956, 0.0320,
    0.0049, 0.0093, 0.0633, 0.1655, 0.2904, 0.4334, 0.5945, 0.7621, 0.9163, 1.0263, 1.0622, 1.0026,
    0.8524, 0.6424, 0.4479, 0.2835, 0.1649, 0.0874, 0.0468, 0.0227, 0.0114, 0.0058, 0.0029, 0.0014,
    0.00069, 0.00033, 0.00017, 0.00008, 0.00004,
];

pub const Y_BAR_2: [f32; 41] = [
    0.0000, 0.0001, 0.0004, 0.0012, 0.0040, 0.0116, 0.0230, 0.0380, 0.0600, 0.0910, 0.1390, 0.2080,
    0.3230, 0.5030, 0.7100, 0.8620, 0.9540, 0.9950, 0.9950, 0.9520, 0.8700, 0.7570, 0.6310, 0.5030,
    0.3810, 0.2650, 0.1750, 0.1070, 0.0610, 0.0320, 0.0170, 0.0082, 0.0041, 0.0021, 0.0010,
    0.00052, 0.00025, 0.00012, 0.00006, 0.00003, 0.00001,
];

pub const Z_BAR_2: [f32; 41] = [
    0.0065, 0.0201, 0.0679, 0.2074, 0.6456, 1.3856, 1.7471, 1.7721, 1.5794, 1.1143, 0.5701, 0.1970,
    0.0415, 0.0052, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
];

/// CIE 1964 10-degree Standard Observer CMFs (380-780nm, 10nm steps)
#[expect(
    clippy::approx_constant,
    reason = "Standard CIE constant with prescribed precision"
)]
pub const X_BAR_10: [f32; 41] = [
    0.0002, 0.0011, 0.0061, 0.0315, 0.1241, 0.3023, 0.5045, 0.6931, 0.8177, 0.7530, 0.5314, 0.3345,
    0.1570, 0.0538, 0.0331, 0.1117, 0.2230, 0.4243, 0.6627, 0.8690, 1.0107, 1.0743, 1.0257, 0.8724,
    0.6553, 0.4456, 0.2800, 0.1622, 0.0869, 0.0434, 0.0218, 0.0107, 0.0053, 0.0026, 0.0013, 0.0006,
    0.0003, 0.0001, 0.0000, 0.0000, 0.0000,
];
pub const Y_BAR_10: [f32; 41] = [
    0.0000, 0.0000, 0.0002, 0.0010, 0.0041, 0.0105, 0.0207, 0.0407, 0.0702, 0.1120, 0.1852, 0.2904,
    0.4190, 0.5764, 0.7435, 0.8872, 0.9666, 0.9983, 0.9873, 0.9331, 0.8420, 0.7163, 0.5596, 0.4203,
    0.3021, 0.2003, 0.1245, 0.0713, 0.0380, 0.0189, 0.0094, 0.0046, 0.0023, 0.0111, 0.0006, 0.0003,
    0.0001, 0.0000, 0.0000, 0.0000, 0.0000,
];
pub const Z_BAR_10: [f32; 41] = [
    0.0007, 0.0045, 0.0259, 0.1343, 0.5285, 1.3003, 2.1932, 3.0334, 3.5534, 3.2392, 2.2235, 1.3400,
    0.5752, 0.1866, 0.0427, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
];

// X_BAR_10, Y_BAR_10, Z_BAR_10 are already pub const.
// No changes needed to visibility.

/// CIE 2015 Physiologically-based LMS Color Matching Functions (2-degree observer, 10nm)
/// These represent the real cone response of the human eye.
pub const L_BAR_2015: [f32; 41] = [
    0.0001, 0.0004, 0.0019, 0.0084, 0.0292, 0.0544, 0.0652, 0.0660, 0.0536, 0.0336, 0.0253, 0.0435,
    0.0906, 0.1834, 0.3541, 0.5363, 0.7024, 0.8358, 0.9328, 0.9859, 1.0000, 0.9575, 0.8524, 0.7081,
    0.5480, 0.3952, 0.2644, 0.1651, 0.0967, 0.0538, 0.0284, 0.0143, 0.0068, 0.0031, 0.0014, 0.0006,
    0.0003, 0.0001, 0.0001, 0.0000, 0.0000,
];
pub const M_BAR_2015: [f32; 41] = [
    0.0000, 0.0001, 0.0006, 0.0028, 0.0121, 0.0298, 0.0450, 0.0526, 0.0519, 0.0440, 0.0494, 0.0772,
    0.1345, 0.2319, 0.3802, 0.5312, 0.6724, 0.7974, 0.8926, 0.9515, 0.9757, 0.9592, 0.8995, 0.7963,
    0.6621, 0.5134, 0.3698, 0.2486, 0.1557, 0.0917, 0.0511, 0.0270, 0.0135, 0.0064, 0.0030, 0.0013,
    0.0006, 0.0003, 0.0001, 0.0001, 0.0000,
];
pub const S_BAR_2015: [f32; 41] = [
    0.0019, 0.0101, 0.0469, 0.1648, 0.4449, 0.8443, 0.9930, 0.8970, 0.6171, 0.3392, 0.1505, 0.0532,
    0.1042, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
    0.0000, 0.0000, 0.0000, 0.0000, 0.0000,
];

/// CIE Standard Illuminants.
/// IMPORTANT: Use the correct observer angle (2° or 10°) matching your CMFs.
pub mod illuminant {
    use super::XYZ;

    // ==================== 2-DEGREE OBSERVER ====================
    /// D50 (Horizon Light, Print Industry - 2°)
    pub const D50: XYZ = XYZ {
        x: 0.96422,
        y: 1.0,
        z: 0.82521,
    };
    /// D55 (Mid-Morning Daylight - 2°)
    pub const D55: XYZ = XYZ {
        x: 0.95682,
        y: 1.0,
        z: 0.92149,
    };
    /// D65 (Noon Daylight, sRGB Standard - 2°)
    pub const D65: XYZ = XYZ {
        x: 0.95047,
        y: 1.0,
        z: 1.08883,
    };
    /// D75 (North Sky Daylight - 2°)
    pub const D75: XYZ = XYZ {
        x: 0.94972,
        y: 1.0,
        z: 1.22638,
    };
    /// Illuminant A (Tungsten 2856K - 2°)
    pub const A: XYZ = XYZ {
        x: 1.09850,
        y: 1.0,
        z: 0.35585,
    };
    /// F2 (Cool White Fluorescent - 2°)
    pub const F2: XYZ = XYZ {
        x: 0.99186,
        y: 1.0,
        z: 0.67393,
    };
    /// F7 (Daylight Fluorescent - 2°)
    pub const F7: XYZ = XYZ {
        x: 0.95041,
        y: 1.0,
        z: 1.08747,
    };
    /// F11 (TL84 Narrow Band - 2°)
    pub const F11: XYZ = XYZ {
        x: 1.00962,
        y: 1.0,
        z: 0.64350,
    };

    // ==================== 10-DEGREE OBSERVER ====================
    /// D50 (10-degree observer)
    pub const D50_10: XYZ = XYZ {
        x: 0.96720,
        y: 1.0,
        z: 0.81427,
    };
    /// D55 (10-degree observer)
    pub const D55_10: XYZ = XYZ {
        x: 0.95799,
        y: 1.0,
        z: 0.90926,
    };
    /// D65 (10-degree observer)
    pub const D65_10: XYZ = XYZ {
        x: 0.94811,
        y: 1.0,
        z: 1.07304,
    };
    /// D75 (10-degree observer)
    pub const D75_10: XYZ = XYZ {
        x: 0.94416,
        y: 1.0,
        z: 1.20641,
    };
    /// Illuminant A (10-degree observer)
    pub const A_10: XYZ = XYZ {
        x: 1.11144,
        y: 1.0,
        z: 0.35200,
    };

    /// CIE 2018 LED Series Illuminants (2-degree).
    /// These replace old F-series for modern lighting analysis.
    pub mod led {
        use super::super::XYZ;
        /// LED-B1: Typical Blue-pumped white LED (2733K)
        pub const B1: XYZ = XYZ {
            x: 1.0967,
            y: 1.0,
            z: 0.3533,
        };
        /// LED-B3: Standard white LED (4103K)
        pub const B3: XYZ = XYZ {
            x: 1.0031,
            y: 1.0,
            z: 0.5361,
        };
        /// LED-B5: High CRI white LED (6598K)
        pub const B5: XYZ = XYZ {
            x: 0.9482,
            y: 1.0,
            z: 1.0642,
        };
        /// LED-BH1: Hybrid warm LED (2851K)
        pub const BH1: XYZ = XYZ {
            x: 1.0824,
            y: 1.0,
            z: 0.3592,
        };
    }

    // Legacy aliases for backward compatibility
    pub const D55_2: XYZ = D55;
    pub const D65_2: XYZ = D65;

    /// Relative Spectral Power Distributions (380-780nm, 10nm steps)
    pub mod spd {
        /// CIE Standard Illuminant A (Tungsten)
        pub const A: [f32; 41] = [
            9.80, 12.09, 14.71, 17.68, 21.00, 24.67, 28.70, 33.09, 37.81, 42.87, 48.24, 53.91,
            59.86, 66.06, 72.50, 79.13, 85.95, 92.91, 100.00, 107.18, 114.44, 121.73, 129.04,
            136.35, 143.62, 150.84, 157.98, 165.03, 171.96, 178.77, 185.43, 191.93, 198.26, 204.41,
            210.37, 216.12, 221.67, 227.00, 232.12, 237.01, 241.68,
        ];

        /// CIE Standard Illuminant D50 (Horizon Daylight)
        pub const D50: [f32; 41] = [
            24.49, 29.87, 49.31, 56.51, 60.03, 57.82, 74.82, 87.25, 90.61, 91.37, 95.11, 91.96,
            95.72, 96.61, 97.13, 102.10, 100.75, 102.32, 100.00, 97.74, 98.92, 93.50, 97.69, 99.27,
            99.04, 95.72, 98.86, 95.67, 98.19, 103.00, 99.13, 87.38, 91.60, 92.89, 76.85, 86.51,
            92.58, 78.23, 57.69, 82.92, 78.27,
        ];

        /// CIE Standard Illuminant D65 (Noon Daylight)
        pub const D65: [f32; 41] = [
            49.98, 54.65, 82.75, 91.49, 93.43, 86.68, 104.87, 117.01, 117.81, 114.86, 115.92,
            108.81, 109.35, 107.80, 104.79, 107.69, 104.41, 104.05, 100.00, 96.33, 95.79, 88.69,
            90.01, 89.60, 87.70, 83.29, 83.70, 80.03, 80.21, 82.28, 78.28, 69.72, 71.61, 74.35,
            61.60, 69.89, 75.09, 63.59, 46.42, 66.81, 63.38,
        ];

        /// CIE Standard Illuminant F2 (Cool White Fluorescent)
        pub const F2: [f32; 41] = [
            1.87, 2.94, 5.17, 6.13, 7.01, 8.56, 43.67, 16.94, 11.35, 12.37, 13.00, 13.23, 13.13,
            12.52, 11.83, 11.22, 11.03, 11.53, 27.74, 17.05, 14.33, 15.52, 19.55, 14.91, 13.22,
            11.12, 8.95, 7.02, 5.42, 4.15, 3.20, 2.47, 1.93, 1.67, 1.29, 1.08, 0.88, 0.77, 0.73,
            0.69, 0.68,
        ];

        /// CIE Standard Illuminant F7 (Broadband Fluorescent)
        pub const F7: [f32; 41] = [
            1.87, 2.92, 5.10, 6.00, 6.85, 8.31, 40.76, 16.06, 10.91, 11.83, 12.40, 12.58, 12.47,
            11.89, 11.33, 10.96, 11.16, 12.12, 27.78, 17.73, 15.20, 16.10, 19.50, 14.64, 12.69,
            10.45, 8.29, 6.41, 4.90, 3.72, 2.83, 2.19, 1.71, 1.43, 1.13, 0.96, 0.78, 0.68, 0.65,
            0.62, 0.62,
        ];

        /// CIE Standard Illuminant F11 (Narrowband Fluorescent)
        pub const F11: [f32; 41] = [
            1.87, 2.35, 2.92, 3.45, 5.10, 18.91, 6.00, 6.11, 6.85, 7.58, 8.31, 40.76, 16.06, 10.32,
            10.91, 11.40, 11.83, 12.17, 12.40, 12.54, 12.58, 12.52, 12.47, 12.20, 11.89, 11.61,
            11.33, 11.10, 10.96, 10.97, 11.16, 11.54, 12.12, 27.78, 17.73, 14.47, 15.20, 15.77,
            16.10, 18.54, 19.50,
        ];
    }
}

impl Illuminant {
    /// Get the relative spectral power distribution for this illuminant.
    pub fn get_spd(&self) -> &'static [f32; 41] {
        match self {
            Illuminant::A => &illuminant::spd::A,
            Illuminant::D50 => &illuminant::spd::D50,
            Illuminant::D55 => &illuminant::spd::D65, // Fallback to D65 if D55 SPD missing
            Illuminant::D65 => &illuminant::spd::D65,
            Illuminant::D75 => &illuminant::spd::D65, // Fallback
            Illuminant::F2 => &illuminant::spd::F2,
            Illuminant::F7 => &illuminant::spd::F7,
            Illuminant::F11 => &illuminant::spd::F11,
        }
    }

    pub fn get_white_point(&self, observer: Observer) -> XYZ {
        match observer {
            Observer::CIE1931_2 => match self {
                Illuminant::D50 => illuminant::D50,
                Illuminant::D55 => illuminant::D55,
                Illuminant::D65 => illuminant::D65,
                Illuminant::D75 => illuminant::D75,
                Illuminant::A => illuminant::A,
                Illuminant::F2 => illuminant::F2,
                Illuminant::F7 => illuminant::F7,
                Illuminant::F11 => illuminant::F11,
            },
            Observer::CIE1964_10 => match self {
                Illuminant::D50 => illuminant::D50_10,
                Illuminant::D55 => illuminant::D55_10,
                Illuminant::D65 => illuminant::D65_10,
                Illuminant::D75 => illuminant::D75_10,
                Illuminant::A => illuminant::A_10,
                // Fallback to 2-degree if 10-degree constants are missing
                Illuminant::F2 => illuminant::F2,
                Illuminant::F7 => illuminant::F7,
                Illuminant::F11 => illuminant::F11,
            },
        }
    }
}

impl Observer {
    pub fn get_cmfs(&self) -> (&'static [f32; 41], &'static [f32; 41], &'static [f32; 41]) {
        match self {
            Observer::CIE1931_2 => (&X_BAR_2, &Y_BAR_2, &Z_BAR_2),
            Observer::CIE1964_10 => (&X_BAR_10, &Y_BAR_10, &Z_BAR_10),
        }
    }
}

/// ASTM E308 Weighting Factors for D65/2° at 10nm.
/// These factors include spectral bandwidth compensation and are the 
/// industry standard for computing tristimulus values from reflectance.
#[rustfmt::skip]
pub mod weighting {
    /// Tristimulus weighting factors for X (D65, 2-degree, 10nm, 380-780nm)
    pub const WX_D65_2_10: [f32; 41] = [
        0.007, 0.022, 0.112, 0.377, 1.188, 2.330, 3.459, 3.724, 3.243, 2.126, 
        1.049, 0.330, 0.051, 0.095, 0.628, 1.687, 2.870, 4.267, 5.628, 6.948, 
        8.310, 8.618, 9.050, 8.505, 7.077, 5.066, 3.549, 2.147, 1.252, 0.681, 
        0.347, 0.150, 0.077, 0.041, 0.017, 0.009, 0.005, 0.002, 0.001, 0.001, 
        0.000
    ];

    /// Tristimulus weighting factors for Y (D65, 2-degree, 10nm, 380-780nm)
    pub const WY_D65_2_10: [f32; 41] = [
        0.000, 0.001, 0.003, 0.010, 0.035, 0.095, 0.228, 0.421, 0.669, 0.989, 
        1.524, 2.141, 3.344, 5.131, 7.041, 8.785, 9.425, 9.792, 9.416, 8.675, 
        7.887, 6.354, 5.374, 4.265, 3.162, 2.089, 1.386, 0.810, 0.463, 0.249, 
        0.126, 0.054, 0.028, 0.015, 0.006, 0.003, 0.002, 0.001, 0.000, 0.000, 
        0.000
    ];

    /// Tristimulus weighting factors for Z (D65, 2-degree, 10nm, 380-780nm)
    pub const WZ_D65_2_10: [f32; 41] = [
        0.035, 0.119, 0.610, 2.059, 6.541, 13.031, 19.880, 22.491, 20.182, 13.888, 
        7.167, 2.325, 0.492, 0.061, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 
        0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 
        0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 
        0.000
    ];

    /// Sum of WY_D65_2_10 weights. 
    /// Corrected to 100.0 to simplify normalization logic.
    pub const SUM_WY_D65_2_10: f32 = 100.000;

    // ==================== D50 (Print Industry Standard) ====================
    
    /// Tristimulus weighting factors for X (D50, 2-degree, 10nm, 380-780nm)
    /// Target White Point: X=96.42, Y=100.00, Z=82.52 (CIE 15:2004)
    pub const WX_D50_2_10: [f32; 41] = [
        0.003, 0.012, 0.067, 0.235, 0.770, 1.566, 2.483, 2.794, 2.510, 1.700, 
        0.866, 0.280, 0.045, 0.086, 0.585, 1.608, 2.785, 4.221, 5.659, 7.090, 
        8.627, 9.137, 9.882, 9.480, 8.042, 5.858, 4.219, 2.585, 1.543, 0.858, 
        0.442, 0.189, 0.100, 0.051, 0.021, 0.012, 0.006, 0.002, 0.001, 0.001, 
        0.000
    ];

    /// Tristimulus weighting factors for Y (D50, 2-degree, 10nm, 380-780nm)
    pub const WY_D50_2_10: [f32; 41] = [
        0.000, 0.000, 0.002, 0.006, 0.023, 0.064, 0.164, 0.316, 0.518, 0.792, 
        1.259, 1.821, 2.943, 4.626, 6.563, 8.376, 9.148, 9.688, 9.469, 8.854, 
        8.190, 6.738, 5.869, 4.755, 3.594, 2.416, 1.648, 0.975, 0.571, 0.314, 
        0.161, 0.068, 0.036, 0.019, 0.007, 0.004, 0.002, 0.001, 0.001, 0.000, 
        0.000
    ];

    /// Tristimulus weighting factors for Z (D50, 2-degree, 10nm, 380-780nm)
    pub const WZ_D50_2_10: [f32; 41] = [
        0.018, 0.067, 0.373, 1.306, 4.318, 8.921, 14.544, 17.196, 15.915, 11.320, 
        6.028, 2.014, 0.442, 0.056, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 
        0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 
        0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 
        0.000
    ];

    /// Sum of WY_D50_2_10 weights.
    pub const SUM_WY_D50_2_10: f32 = 100.000;
}

/// Bradford chromatic adaptation transform.
/// Converts XYZ from one illuminant to another using the Bradford cone response model.
///
/// Reference: Lindbloom (http://www.brucelindbloom.com/index.html?Eqn_ChromAdapt.html)
/// Note: The Bradford matrix used here is the "Sharp" variant commonly used in ICC profiles.
pub mod chromatic_adaptation {
    use super::XYZ;

    /// Apply Bradford transform to adapt XYZ from source to destination white point.
    ///
    /// # Example
    /// ```
    /// use spectro_rs::colorimetry::{XYZ, illuminant, chromatic_adaptation};
    /// let xyz_d50 = XYZ { x: 0.5, y: 0.5, z: 0.4 };
    /// let xyz_d65 = chromatic_adaptation::bradford_adapt(xyz_d50, illuminant::D50, illuminant::D65);
    /// ```
    #[expect(
        clippy::excessive_precision,
        reason = "The Bradford transform matrix requires high precision as defined in the ICC/Lindbloom standard"
    )]
    pub fn bradford_adapt(xyz: XYZ, src_wp: XYZ, dst_wp: XYZ) -> XYZ {
        // Bradford M matrix (XYZ to LMS cone response)
        // Source: Bruce Lindbloom, ICC Profile specification
        #[rustfmt::skip]
        let m = [
            [ 0.8951000,  0.2664000, -0.1614000],
            [-0.7502000,  1.7135000,  0.0367000],
            [ 0.0389000, -0.0685000,  1.0296000],
        ];
        // Inverse Bradford M matrix (computed to match M exactly)
        #[rustfmt::skip]
        let m_inv = [
            [ 0.9869929, -0.1470543,  0.1599627],
            [ 0.4323053,  0.5183603,  0.0492912],
            [-0.0085287,  0.0400428,  0.9684867],
        ];

        // Convert to LMS
        let src_lms = [
            m[0][0] * src_wp.x + m[0][1] * src_wp.y + m[0][2] * src_wp.z,
            m[1][0] * src_wp.x + m[1][1] * src_wp.y + m[1][2] * src_wp.z,
            m[2][0] * src_wp.x + m[2][1] * src_wp.y + m[2][2] * src_wp.z,
        ];
        let dst_lms = [
            m[0][0] * dst_wp.x + m[0][1] * dst_wp.y + m[0][2] * dst_wp.z,
            m[1][0] * dst_wp.x + m[1][1] * dst_wp.y + m[1][2] * dst_wp.z,
            m[2][0] * dst_wp.x + m[2][1] * dst_wp.y + m[2][2] * dst_wp.z,
        ];

        // Scaling factors
        let scale = [
            dst_lms[0] / src_lms[0],
            dst_lms[1] / src_lms[1],
            dst_lms[2] / src_lms[2],
        ];

        // Convert input XYZ to LMS
        let lms = [
            m[0][0] * xyz.x + m[0][1] * xyz.y + m[0][2] * xyz.z,
            m[1][0] * xyz.x + m[1][1] * xyz.y + m[1][2] * xyz.z,
            m[2][0] * xyz.x + m[2][1] * xyz.y + m[2][2] * xyz.z,
        ];

        // Scale LMS
        let lms_adapted = [lms[0] * scale[0], lms[1] * scale[1], lms[2] * scale[2]];

        // Convert back to XYZ
        XYZ {
            x: m_inv[0][0] * lms_adapted[0]
                + m_inv[0][1] * lms_adapted[1]
                + m_inv[0][2] * lms_adapted[2],
            y: m_inv[1][0] * lms_adapted[0]
                + m_inv[1][1] * lms_adapted[1]
                + m_inv[1][2] * lms_adapted[2],
            z: m_inv[2][0] * lms_adapted[0]
                + m_inv[2][1] * lms_adapted[1]
                + m_inv[2][2] * lms_adapted[2],
        }
    }
}

#[derive(Debug, Clone, Copy, PartialEq, serde::Serialize, serde::Deserialize)]
pub struct XYZ {
    pub x: f32,
    pub y: f32,
    pub z: f32,
}

#[derive(Debug, Clone, Copy, PartialEq, serde::Serialize, serde::Deserialize)]
pub struct Lab {
    pub l: f32,
    pub a: f32,
    pub b: f32,
}

/// LMS color space representing cone responses (Long, Medium, Short wavelengths).
/// Based on CIE 2015 physiologically-based sensitivity data.
#[derive(Debug, Clone, Copy, PartialEq, serde::Serialize, serde::Deserialize)]
pub struct LMS {
    pub l: f32,
    pub m: f32,
    pub s: f32,
}

/// Jzazbz: A modern perceptually uniform color space (Safdar et al., 2017).
/// Designed for HDR content with excellent uniformity across the entire
/// luminance range (0-10,000 nits). Euclidean distance in this space
/// directly represents perceptual difference.
#[derive(Debug, Clone, Copy, PartialEq, serde::Serialize, serde::Deserialize)]
pub struct Jzazbz {
    /// Lightness (0.0 = black, higher = brighter, no upper limit for HDR)
    pub jz: f32,
    /// Red-green opponent (negative = green, positive = red)
    pub az: f32,
    /// Blue-yellow opponent (negative = blue, positive = yellow)
    pub bz: f32,
}

pub fn calculate_cct(spd: &crate::spectrum::SpectralData) -> (f32, f32) {
    let xyz = spd.to_xyz();
    let sum = xyz.x + xyz.y + xyz.z;
    if sum == 0.0 {
        return (0.0, 0.0);
    }
    let x = xyz.x / sum;
    let y = xyz.y / sum;

    // McCamy's formula
    let n = (x - 0.3320) / (0.1858 - y);
    let cct = 449.0 * n.powi(3) + 3525.0 * n.powi(2) + 6823.3 * n + 5524.3;

    (cct, 0.0)
}

impl XYZ {
    /// Convert XYZ to CIE L*a*b* using the given white point.
    /// Uses precise CIE constants for continuity at the threshold.
    pub fn to_lab(&self, wp: XYZ) -> Lab {
        // CIE standard constants for continuity
        const EPSILON: f32 = 216.0 / 24389.0; // ≈ 0.008856
        const KAPPA: f32 = 24389.0 / 27.0; // ≈ 903.2963

        let f = |t: f32| -> f32 {
            if t > EPSILON {
                t.powf(1.0 / 3.0)
            } else {
                (KAPPA * t + 16.0) / 116.0
            }
        };

        let fx = f(self.x / wp.x);
        let fy = f(self.y / wp.y);
        let fz = f(self.z / wp.z);

        Lab {
            l: 116.0 * fy - 16.0,
            a: 500.0 * (fx - fy),
            b: 200.0 * (fy - fz),
        }
    }

    /// Convert XYZ (D65) to linear sRGB, then apply gamma correction.
    /// Returns clamped (r, g, b) values in [0, 255].
    ///
    /// **IMPORTANT**: This function assumes the input XYZ is referenced to D65.
    /// If your XYZ values are based on a different illuminant (e.g., D50 from
    /// an ICC profile), you must first use `chromatic_adaptation::bradford_adapt`
    /// to convert to D65 before calling this method.
    #[expect(
        clippy::excessive_precision,
        reason = "High precision is required for accurate sRGB gamut conversion"
    )]
    pub fn to_srgb(&self) -> (u8, u8, u8) {
        // XYZ to linear sRGB matrix (IEC 61966-2-1, D65 reference)
        let r_lin = 3.2404542 * self.x - 1.5371385 * self.y - 0.4985314 * self.z;
        let g_lin = -0.9692660 * self.x + 1.8760108 * self.y + 0.0415560 * self.z;
        let b_lin = 0.0556434 * self.x - 0.2040259 * self.y + 1.0572252 * self.z;

        // Gamma correction (sRGB companding)
        fn gamma(c: f32) -> f32 {
            if c <= 0.0031308 {
                12.92 * c
            } else {
                1.055 * c.powf(1.0 / 2.4) - 0.055
            }
        }

        let r = (gamma(r_lin).clamp(0.0, 1.0) * 255.0).round() as u8;
        let g = (gamma(g_lin).clamp(0.0, 1.0) * 255.0).round() as u8;
        let b = (gamma(b_lin).clamp(0.0, 1.0) * 255.0).round() as u8;

        (r, g, b)
    }

    /// A safer version of `to_srgb` that takes the current white point of the XYZ
    /// and automatically performs Bradford chromatic adaptation to D65 if needed.
    pub fn to_srgb_safe(&self, current_wp: XYZ) -> (u8, u8, u8) {
        if (self.x - current_wp.x).abs() < 1e-5
            && (self.y - current_wp.y).abs() < 1e-5
            && (self.z - current_wp.z).abs() < 1e-5
        {
            // Already matched or specialized logic could go here
        }

        if current_wp == illuminant::D65 {
            self.to_srgb()
        } else {
            let adapted = chromatic_adaptation::bradford_adapt(*self, current_wp, illuminant::D65);
            adapted.to_srgb()
        }
    }

    /// Convert XYZ (absolute, D65) to Jzazbz color space.
    /// Jzazbz (Safdar et al., 2017) is designed for HDR and provides
    /// excellent perceptual uniformity across the full luminance range.
    ///
    /// # Parameters
    /// * `luminance_scale`: Scaling factor for absolute luminance (cd/m²).
    ///   For SDR, typically 100.0 or 200.0.
    #[expect(
        clippy::excessive_precision,
        reason = "Jzazbz conversion matrices depend on exact coefficients from the original paper for perceptual uniformity"
    )]
    pub fn to_jzazbz(&self, luminance_scale: f32) -> Jzazbz {
        // Step 1: Absolute luminance scaling
        let x = self.x * luminance_scale;
        let y = self.y * luminance_scale;
        let z = self.z * luminance_scale;

        // Step 2: XYZ to LMS (M1 matrix from Safdar et al. 2017)
        let l = 0.41478972 * x + 0.57999905 * y + 0.01464805 * z;
        let m = -0.20151003 * x + 1.12064859 * y + 0.05310084 * z;
        let s = -0.01660078 * x + 0.26480015 * y + 0.66847986 * z;

        // Step 3: PQ transfer function (Perceptual Quantizer, ST 2084)
        fn pq(v: f32) -> f32 {
            let v_abs = (v.max(0.0) / 10000.0) as f64; // PQ normalized to 10,000 nits
            let n = 2610.0 / 16384.0;
            let p = 1.7 * (2523.0 / 32.0);
            let c1 = 3424.0 / 4096.0;
            let c2 = 2413.0 / 128.0;
            let c3 = 2392.0 / 128.0;

            let v_pow_n = v_abs.powf(n);
            (((c1 + c2 * v_pow_n) / (1.0 + c3 * v_pow_n)).powf(p)) as f32
        }

        let lp = pq(l);
        let mp = pq(m);
        let sp = pq(s);

        // Step 4: Izazbz
        let iz = 0.5 * lp + 0.5 * mp;
        let az = 3.524000 * lp - 4.066708 * mp + 0.542708 * sp;
        let bz = 0.199076 * lp + 1.096799 * mp - 1.295875 * sp;

        // Step 5: Jz (with white point offset for neutral alignment)
        let d = -0.56;
        let jz = ((1.0 + d) * iz) / (1.0 + d * iz) - 0.005605;

        Jzazbz {
            jz: jz.max(0.0),
            az,
            bz,
        }
    }

    /// Calculate Tristimulus XYZ from spectral reflectance using ASTM E308 weighting factors.
    /// Assumes D65/2° and 10nm sampling (380-780nm).
    ///
    /// # Notes
    /// - ASTM E308 weighting factors already include the D65 SPD and normalization.
    /// - Input reflectance should be in the range 0.0-1.0 (or 0-100%).
    /// - Output Y=100 corresponds to a perfect diffuser.
    pub fn from_reflectance_10nm(reflectance: &[f32; 41]) -> Self {
        let (x, y, z) = reflectance
            .iter()
            .zip(weighting::WX_D65_2_10.iter())
            .zip(weighting::WY_D65_2_10.iter())
            .zip(weighting::WZ_D65_2_10.iter())
            .fold(
                (0.0f32, 0.0f32, 0.0f32),
                |(x, y, z), (((r, wx), wy), wz)| (x + r * wx, y + r * wy, z + r * wz),
            );

        // ASTM E308 weights are pre-scaled so that sum(WY * R) = Y directly
        // when R is in 0-1 range. Multiply by 100 to get standard Y=100 scale.
        Self {
            x: x * 100.0 / weighting::SUM_WY_D65_2_10,
            y: y * 100.0 / weighting::SUM_WY_D65_2_10,
            z: z * 100.0 / weighting::SUM_WY_D65_2_10,
        }
    }

    /// Convert XYZ to CIE 1960 UCS (u, v) coordinates.
    /// This is used for CCT calculation and CRI reference illuminant alignment.
    pub fn to_uv_1960(&self) -> (f32, f32) {
        let denom = self.x + 15.0 * self.y + 3.0 * self.z;
        if denom.abs() < 1e-9 {
            return (0.0, 0.0);
        }
        let u = (4.0 * self.x) / denom;
        let v = (6.0 * self.y) / denom;
        (u, v)
    }

    /// Convert XYZ to CIE 1964 (U*, V*, W*) color space.
    /// `white_uv`: The (u, v) coordinates of the white point.
    pub fn to_uvw_1964(&self, white_uv: (f32, f32)) -> (f32, f32, f32) {
        let (u, v) = self.to_uv_1960();
        let w_star = 25.0 * self.y.powf(1.0 / 3.0) - 17.0;
        let u_star = 13.0 * w_star * (u - white_uv.0);
        let v_star = 13.0 * w_star * (v - white_uv.1);
        (u_star, v_star, w_star)
    }
}

pub mod generation {
    /// CIE Standard Illuminant A/D series generation.
    ///
    /// CIE S0, S1, S2 spectral power components for Daylight reconstruction (380nm-780nm, 10nm steps)
    pub const S0: [f32; 41] = [
        0.0, 0.0, 33.4, 37.4, 117.4, 117.8, 114.9, 115.9, 108.8, 109.3, 107.8, 104.8, 107.7, 104.4,
        104.0, 100.0, 96.0, 95.1, 89.1, 90.5, 90.3, 88.4, 84.0, 85.1, 81.9, 82.6, 84.9, 81.3, 71.9,
        74.3, 76.4, 63.3, 71.7, 77.0, 65.2, 47.7, 68.6, 65.0, 66.0, 61.0, 53.3,
    ];
    pub const S1: [f32; 41] = [
        0.0, 0.0, -1.1, -0.5, -0.7, -1.2, -2.6, -2.9, -2.8, -4.5, -6.1, -7.6, -9.7, -11.7, -12.2,
        -13.6, -12.0, -13.3, -12.9, -10.6, -11.6, -10.8, -8.1, -10.3, -11.0, -11.5, -10.8, -10.9,
        -8.8, -7.3, -12.9, -15.8, -15.1, -12.2, -10.2, -8.6, -12.0, -14.6, -15.1, -14.9, -13.7,
    ];
    pub const S2: [f32; 41] = [
        0.0, 0.0, -2.1, -1.9, -1.1, -2.2, -3.5, -3.5, -3.3, -2.0, -1.2, -1.1, -0.5, 0.2, 0.5, 2.1,
        3.2, 4.1, 4.7, 5.1, 6.7, 7.3, 8.6, 9.8, 10.2, 14.9, 18.1, 15.9, 16.8, 24.2, 31.7, 15.3,
        18.9, 21.2, 15.6, 8.3, 18.9, 14.6, 15.5, 15.4, 14.6,
    ];

    /// Generate Planckian radiator SPD (black-body) for a given CCT and wavelengths.
    pub fn generate_planckian(cct: f32, wavelengths: &[f32]) -> Vec<f32> {
        let c1 = 3.741771e-16;
        let c2 = 1.4388e-2;
        wavelengths
            .iter()
            .map(|&wl| {
                let wl_m = wl * 1e-9;
                if wl_m == 0.0 {
                    0.0
                } else {
                    c1 * wl_m.powi(-5) / ((c2 / (wl_m * cct)).exp() - 1.0)
                }
            })
            .collect()
    }

    /// Generate CIE Daylight SPD for a given CCT and wavelengths.
    pub fn generate_daylight(cct: f32, wavelengths: &[f32]) -> Vec<f32> {
        let x_d = if cct <= 7000.0 {
            -4.6070e9 / cct.powi(3) + 2.9678e6 / cct.powi(2) + 0.09911e3 / cct + 0.244063
        } else {
            -2.0064e9 / cct.powi(3) + 1.9018e6 / cct.powi(2) + 0.24748e3 / cct + 0.237040
        };

        let y_d = -3.000 * x_d * x_d + 2.870 * x_d - 0.275;

        let m1 = (-1.3515 - 1.7703 * x_d + 5.9114 * y_d) / (0.0241 + 0.2562 * x_d - 0.7341 * y_d);
        let m2 = (0.0300 - 31.4424 * x_d + 30.0717 * y_d) / (0.0241 + 0.2562 * x_d - 0.7341 * y_d);

        wavelengths
            .iter()
            .map(|&wl| {
                // Interpolate S0, S1, S2 from 10nm table (380-780, indices 0-40)
                let t = (wl - 380.0) / 10.0;
                let idx = t.floor() as i32;
                let x = t - idx as f32;

                let get_val = |table: &[f32; 41]| {
                    if idx < 0 {
                        table[0] // Clamp to 380nm value for < 380
                    } else if idx >= 40 {
                        table[40] // Clamp to 780nm value for > 780
                    } else {
                        let v0 = table[idx as usize];
                        let v1 = table[(idx + 1) as usize];
                        v0 + x * (v1 - v0)
                    }
                };

                let s0 = get_val(&S0);
                let s1 = get_val(&S1);
                let s2 = get_val(&S2);

                let val = s0 + m1 * s1 + m2 * s2;
                if val < 0.0 { 0.0 } else { val }
            })
            .collect()
    }
}

impl Lab {
    /// Convert Lab back to XYZ using the given white point.
    /// Uses precise CIE constants for continuity at the threshold.
    pub fn to_xyz(&self, wp: XYZ) -> XYZ {
        // CIE standard constants for continuity
        const EPSILON: f32 = 216.0 / 24389.0; // ≈ 0.008856
        const KAPPA: f32 = 24389.0 / 27.0; // ≈ 903.2963

        let fy = (self.l + 16.0) / 116.0;
        let fx = self.a / 500.0 + fy;
        let fz = fy - self.b / 200.0;

        let f_inv = |t: f32| -> f32 {
            let t3 = t.powi(3);
            if t3 > EPSILON {
                t3
            } else {
                (116.0 * t - 16.0) / KAPPA
            }
        };

        XYZ {
            x: wp.x * f_inv(fx),
            y: wp.y * f_inv(fy),
            z: wp.z * f_inv(fz),
        }
    }

    /// Convert Lab to sRGB via XYZ (using D65 white point).
    /// Returns clamped (r, g, b) values in [0, 255].
    pub fn to_srgb(&self) -> (u8, u8, u8) {
        self.to_xyz(illuminant::D65_2).to_srgb()
    }

    /// Calculates Delta E*ab (CIE 1976).
    pub fn delta_e_76(&self, other: &Lab) -> f32 {
        ((self.l - other.l).powi(2) + (self.a - other.a).powi(2) + (self.b - other.b).powi(2))
            .sqrt()
    }

    /// Calculates Delta E*00 (CIE 2000) using the Sharma (2005) reference implementation.
    /// This is the industry standard for perceptual color difference.
    pub fn delta_e_2000(&self, other: &Lab) -> f32 {
        // Weight factors (default = 1.0)
        let k_l = 1.0;
        let k_c = 1.0;
        let k_h = 1.0;

        let c1 = (self.a.powi(2) + self.b.powi(2)).sqrt();
        let c2 = (other.a.powi(2) + other.b.powi(2)).sqrt();
        let avg_c = (c1 + c2) / 2.0;

        // Calculate G and adjusted a'
        let g = 0.5 * (1.0 - (avg_c.powi(7) / (avg_c.powi(7) + 25.0f32.powi(7))).sqrt());
        let a1p = (1.0 + g) * self.a;
        let a2p = (1.0 + g) * other.a;

        let c1p = (a1p.powi(2) + self.b.powi(2)).sqrt();
        let c2p = (a2p.powi(2) + other.b.powi(2)).sqrt();

        // Key fix: atan2(b, a') - correct parameter order
        let get_hp = |b: f32, ap: f32| -> f32 {
            if b == 0.0 && ap == 0.0 {
                0.0
            } else {
                let h = b.atan2(ap).to_degrees();
                if h < 0.0 { h + 360.0 } else { h }
            }
        };
        let h1p = get_hp(self.b, a1p);
        let h2p = get_hp(other.b, a2p);

        // Calculate delta values
        let d_lp = other.l - self.l;
        let d_cp = c2p - c1p;

        let mut d_hp_deg = h2p - h1p;
        if c1p * c2p != 0.0 {
            if d_hp_deg.abs() > 180.0 {
                if h2p <= h1p {
                    d_hp_deg += 360.0;
                } else {
                    d_hp_deg -= 360.0;
                }
            }
        } else {
            d_hp_deg = 0.0;
        }
        let d_hp = 2.0 * (c1p * c2p).sqrt() * (d_hp_deg / 2.0).to_radians().sin();

        // Calculate averages
        let avg_lp = (self.l + other.l) / 2.0;
        let avg_cp = (c1p + c2p) / 2.0;

        let mut avg_hp = h1p + h2p;
        if c1p * c2p != 0.0 {
            if (h1p - h2p).abs() > 180.0 {
                if h1p + h2p < 360.0 {
                    avg_hp += 360.0;
                } else {
                    avg_hp -= 360.0;
                }
            }
            avg_hp /= 2.0;
        } else {
            avg_hp = h1p + h2p;
        }

        // T term
        let t = 1.0 - 0.17 * (avg_hp - 30.0).to_radians().cos()
            + 0.24 * (2.0 * avg_hp).to_radians().cos()
            + 0.32 * (3.0 * avg_hp + 6.0).to_radians().cos()
            - 0.20 * (4.0 * avg_hp - 63.0).to_radians().cos();

        // Key fix: SL denominator structure
        let s_l = 1.0 + (0.015 * (avg_lp - 50.0).powi(2)) / (20.0 + (avg_lp - 50.0).powi(2)).sqrt();
        let s_c = 1.0 + 0.045 * avg_cp;
        let s_h = 1.0 + 0.015 * avg_cp * t;

        // Rotation term RT
        let d_theta = 30.0 * (-((avg_hp - 275.0) / 25.0).powi(2)).exp();
        let rc = 2.0 * (avg_cp.powi(7) / (avg_cp.powi(7) + 25.0f32.powi(7))).sqrt();
        let rt = -rc * (2.0 * d_theta.to_radians()).sin();

        // Final Delta E 00
        ((d_lp / (k_l * s_l)).powi(2)
            + (d_cp / (k_c * s_c)).powi(2)
            + (d_hp / (k_h * s_h)).powi(2)
            + rt * (d_cp / (k_c * s_c)) * (d_hp / (k_h * s_h)))
            .sqrt()
    }

    /// Mix two Lab colors by a given ratio (0.0 = self, 1.0 = other).
    pub fn mix(&self, other: &Lab, ratio: f32) -> Lab {
        let ratio = ratio.clamp(0.0, 1.0);
        Lab {
            l: self.l * (1.0 - ratio) + other.l * ratio,
            a: self.a * (1.0 - ratio) + other.a * ratio,
            b: self.b * (1.0 - ratio) + other.b * ratio,
        }
    }

    /// Calculate chroma (C*) from a* and b*.
    pub fn chroma(&self) -> f32 {
        (self.a.powi(2) + self.b.powi(2)).sqrt()
    }

    /// Calculate hue angle (h°) in degrees [0, 360).
    pub fn hue(&self) -> f32 {
        let h = self.b.atan2(self.a).to_degrees();
        if h < 0.0 { h + 360.0 } else { h }
    }
}

impl Jzazbz {
    /// Calculate Delta Ez (Jzazbz Euclidean distance).
    /// Unlike CIEDE2000, this is a simple Euclidean distance that
    /// directly represents perceptual difference due to Jzazbz's
    /// excellent uniformity. Simpler and faster than Delta E 2000.
    pub fn delta_ez(&self, other: &Jzazbz) -> f32 {
        ((self.jz - other.jz).powi(2) + (self.az - other.az).powi(2) + (self.bz - other.bz).powi(2))
            .sqrt()
    }

    /// Calculate chroma (Cz) in Jzazbz space.
    pub fn chroma(&self) -> f32 {
        (self.az.powi(2) + self.bz.powi(2)).sqrt()
    }

    /// Calculate hue angle (hz) in degrees [0, 360).
    pub fn hue(&self) -> f32 {
        let h = self.bz.atan2(self.az).to_degrees();
        if h < 0.0 { h + 360.0 } else { h }
    }

    /// Mix two Jzazbz colors by a given ratio.
    pub fn mix(&self, other: &Jzazbz, ratio: f32) -> Jzazbz {
        let ratio = ratio.clamp(0.0, 1.0);
        Jzazbz {
            jz: self.jz * (1.0 - ratio) + other.jz * ratio,
            az: self.az * (1.0 - ratio) + other.az * ratio,
            bz: self.bz * (1.0 - ratio) + other.bz * ratio,
        }
    }
}

/// Color Rendering Index (CRI) and other light quality metrics.
pub mod metrics {
    use super::*;
    use crate::spectrum::SpectralData;

    /// Test Color Samples (TCS) for CRI calculation (380-780nm, 10nm).
    /// Only including TCS01-TCS08 (for Ra) and TCS09 (for R9).
    #[rustfmt::skip]
    pub const TCS: [[f32; 41]; 9] = [
        // TCS01: Light greyish red
        [0.22, 0.25, 0.26, 0.25, 0.24, 0.24, 0.23, 0.23, 0.22, 0.22, 0.21, 0.22, 0.22, 0.23, 0.23, 0.23, 0.24, 0.25, 0.27, 0.30, 0.34, 0.39, 0.42, 0.44, 0.45, 0.45, 0.45, 0.45, 0.45, 0.45, 0.46, 0.46, 0.46, 0.46, 0.46, 0.47, 0.47, 0.47, 0.47, 0.47, 0.47],
        // TCS02: Dark greyish yellow
        [0.07, 0.09, 0.11, 0.12, 0.12, 0.12, 0.12, 0.12, 0.13, 0.13, 0.14, 0.15, 0.17, 0.21, 0.24, 0.26, 0.27, 0.27, 0.28, 0.30, 0.32, 0.34, 0.34, 0.34, 0.34, 0.34, 0.34, 0.34, 0.34, 0.33, 0.33, 0.33, 0.33, 0.33, 0.32, 0.32, 0.32, 0.32, 0.32, 0.31, 0.31],
        // TCS03: Strong yellow green
        [0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.08, 0.09, 0.11, 0.15, 0.20, 0.24, 0.28, 0.34, 0.39, 0.40, 0.38, 0.35, 0.30, 0.26, 0.25, 0.24, 0.22, 0.22, 0.22, 0.23, 0.25, 0.29, 0.34, 0.39, 0.43, 0.46, 0.48, 0.49, 0.50, 0.51, 0.52, 0.52, 0.53, 0.54],
        // TCS04: Moderate yellowish green
        [0.07, 0.09, 0.12, 0.12, 0.13, 0.14, 0.14, 0.16, 0.19, 0.23, 0.28, 0.33, 0.37, 0.39, 0.39, 0.38, 0.34, 0.31, 0.28, 0.25, 0.21, 0.18, 0.16, 0.16, 0.15, 0.15, 0.15, 0.15, 0.16, 0.17, 0.17, 0.17, 0.17, 0.17, 0.18, 0.19, 0.19, 0.20, 0.21, 0.23, 0.25],
        // TCS05: Light bluish green
        [0.30, 0.31, 0.31, 0.32, 0.33, 0.34, 0.36, 0.38, 0.40, 0.42, 0.42, 0.41, 0.40, 0.39, 0.37, 0.35, 0.31, 0.28, 0.25, 0.22, 0.19, 0.19, 0.18, 0.18, 0.18, 0.18, 0.18, 0.19, 0.19, 0.20, 0.20, 0.20, 0.20, 0.20, 0.21, 0.21, 0.22, 0.22, 0.23, 0.24, 0.27],
        // TCS06: Light blue
        [0.15, 0.27, 0.41, 0.49, 0.52, 0.53, 0.54, 0.56, 0.55, 0.54, 0.52, 0.49, 0.45, 0.41, 0.36, 0.31, 0.25, 0.23, 0.23, 0.22, 0.22, 0.22, 0.22, 0.23, 0.24, 0.26, 0.27, 0.28, 0.28, 0.29, 0.30, 0.33, 0.35, 0.38, 0.40, 0.43, 0.45, 0.47, 0.49, 0.51, 0.53],
        // TCS07: Light violet
        [0.38, 0.52, 0.55, 0.56, 0.56, 0.56, 0.54, 0.51, 0.47, 0.43, 0.39, 0.34, 0.31, 0.30, 0.27, 0.26, 0.26, 0.26, 0.26, 0.25, 0.26, 0.28, 0.32, 0.36, 0.39, 0.41, 0.43, 0.44, 0.45, 0.47, 0.47, 0.48, 0.49, 0.50, 0.51, 0.53, 0.54, 0.55, 0.57, 0.58, 0.59],
        // TCS08: Light reddish purple
        [0.10, 0.17, 0.32, 0.46, 0.49, 0.48, 0.45, 0.43, 0.40, 0.37, 0.34, 0.31, 0.29, 0.28, 0.26, 0.25, 0.26, 0.27, 0.27, 0.28, 0.32, 0.38, 0.48, 0.57, 0.63, 0.66, 0.69, 0.70, 0.71, 0.71, 0.72, 0.72, 0.72, 0.72, 0.73, 0.73, 0.73, 0.73, 0.73, 0.73, 0.73],
        // TCS09: Strong red (R9)
        [0.066, 0.058, 0.052, 0.051, 0.050, 0.048, 0.046, 0.041, 0.035, 0.030, 0.028, 0.028, 0.030, 0.031, 0.032, 0.033, 0.041, 0.048, 0.060, 0.102, 0.190, 0.336, 0.505, 0.641, 0.717, 0.758, 0.781, 0.797, 0.809, 0.819, 0.828, 0.831, 0.835, 0.836, 0.838, 0.839, 0.839, 0.839, 0.839, 0.839, 0.838],
    ];

    /// Calculate CRI Ra and R9 for a given spectral power distribution.
    pub fn calculate_cri(spd: &SpectralData) -> (f32, f32) {
        let xyz = spd.to_xyz_emissive_2();
        let cct = xyz.to_cct();

        // 1. Generate reference SPD
        let ref_spd = if cct < 5000.0 {
            generate_planckian(cct)
        } else {
            generate_daylight(cct)
        };

        // 2. Calculate Ri for each TCS
        let mut ris = [0.0f32; 9];
        for i in 0..9 {
            ris[i] = calculate_ri(spd, &ref_spd, &TCS[i]);
        }

        // Ra is average of first 8
        let ra = ris[0..8].iter().sum::<f32>() / 8.0;
        let r9 = ris[8];

        (ra, r9)
    }

    fn calculate_ri(test_spd: &SpectralData, ref_spd: &[f32; 41], tcs: &[f32; 41]) -> f32 {
        let (xb, yb, zb) = Observer::CIE1931_2.get_cmfs();

        let calc_xyz = |spd_vals: &[f32], tcs_vals: &[f32]| -> XYZ {
            let mut x = 0.0;
            let mut y = 0.0;
            let mut z = 0.0;
            for i in 0..41 {
                let val = spd_vals[i] * tcs_vals[i];
                x += val * xb[i];
                y += val * yb[i];
                z += val * zb[i];
            }
            XYZ { x, y, z }
        };

        let test_xyz = calc_xyz(&test_spd.values, tcs);
        let ref_xyz = calc_xyz(ref_spd, tcs);

        // Simplified CRI calculation (skipping full Von Kries for brevity,
        // but using Lab Delta E as a proxy for Ri)
        // Ri = 100 - 4.6 * Delta E
        let test_lab = test_xyz.to_lab(test_spd.to_xyz_emissive_2());

        let mut ref_white_x = 0.0;
        let mut ref_white_y = 0.0;
        let mut ref_white_z = 0.0;
        for i in 0..41 {
            ref_white_x += ref_spd[i] * xb[i];
            ref_white_y += ref_spd[i] * yb[i];
            ref_white_z += ref_spd[i] * zb[i];
        }
        let ref_white = XYZ {
            x: ref_white_x,
            y: ref_white_y,
            z: ref_white_z,
        };
        let ref_lab = ref_xyz.to_lab(ref_white);

        let de = test_lab.delta_e_76(&ref_lab);
        100.0 - 4.6 * de
    }

    fn generate_planckian(cct: f32) -> [f32; 41] {
        let mut spd = [0.0f32; 41];
        let c1 = 3.741771e-16_f32;
        let c2 = 1.4388e-2_f32;
        for (i, val) in spd.iter_mut().enumerate() {
            let wl = (380 + i * 10) as f32 * 1e-9_f32;
            *val = c1 / (wl.powi(5) * ((c2 / (wl * cct)).exp() - 1.0));
        }
        spd
    }

    fn generate_daylight(cct: f32) -> [f32; 41] {
        // Simplified D-series generator
        let xd = if cct <= 7000.0 {
            -4.6070e9 / cct.powi(3) + 2.9678e6 / cct.powi(2) + 0.09911e3 / cct + 0.244063
        } else {
            -2.0064e9 / cct.powi(3) + 1.9018e6 / cct.powi(2) + 0.24748e3 / cct + 0.237040
        };

        let yd = -3.000 * xd * xd + 2.870 * xd - 0.275;

        let m1 = (-1.3515 - 1.7703 * xd + 5.9114 * yd) / (0.0241 + 0.2562 * xd - 0.7341 * yd);
        let m2 = (0.0300 - 31.4424 * xd + 30.0717 * yd) / (0.0241 + 0.2562 * xd - 0.7341 * yd);

        let mut spd = [0.0f32; 41];
        for (i, val) in spd.iter_mut().enumerate() {
            // Using D65 as a base and applying M1, M2 (simplified)
            *val = illuminant::spd::D65[i] * (1.0 + m1 * 0.01 + m2 * 0.01);
        }
        spd
    }
}

/// Color appearance and analysis utilities.
pub mod appearance {
    use super::{Lab, XYZ, illuminant};
    use crate::spectrum::SpectralData;

    /// Calculate Metamerism Index between two spectral samples.
    /// Compares how differently the samples appear under a test illuminant
    /// relative to a reference illuminant (typically D65).
    ///
    /// A high metamerism index means the samples look similar under one illuminant
    /// but different under another — a common issue in color matching.
    pub fn metamerism_index(
        sample1: &SpectralData,
        sample2: &SpectralData,
        ref_illuminant: XYZ,
        test_illuminant: XYZ,
    ) -> f32 {
        // Calculate Lab under reference illuminant
        let xyz1_ref = sample1.to_xyz();
        let xyz2_ref = sample2.to_xyz();
        let lab1_ref = XYZ {
            x: xyz1_ref.x / 100.0,
            y: xyz1_ref.y / 100.0,
            z: xyz1_ref.z / 100.0,
        }
        .to_lab(ref_illuminant);
        let lab2_ref = XYZ {
            x: xyz2_ref.x / 100.0,
            y: xyz2_ref.y / 100.0,
            z: xyz2_ref.z / 100.0,
        }
        .to_lab(ref_illuminant);

        // Adapt XYZ to test illuminant using Bradford
        let xyz1_test =
            super::chromatic_adaptation::bradford_adapt(xyz1_ref, illuminant::D65, test_illuminant);
        let xyz2_test =
            super::chromatic_adaptation::bradford_adapt(xyz2_ref, illuminant::D65, test_illuminant);

        let lab1_test = XYZ {
            x: xyz1_test.x / 100.0,
            y: xyz1_test.y / 100.0,
            z: xyz1_test.z / 100.0,
        }
        .to_lab(test_illuminant);
        let lab2_test = XYZ {
            x: xyz2_test.x / 100.0,
            y: xyz2_test.y / 100.0,
            z: xyz2_test.z / 100.0,
        }
        .to_lab(test_illuminant);

        // Delta E under reference
        let de_ref = lab1_ref.delta_e_2000(&lab2_ref);
        // Delta E under test
        let de_test = lab1_test.delta_e_2000(&lab2_test);

        // Metamerism index is the difference in color difference
        (de_test - de_ref).abs()
    }

    /// Simulate how a color appears under a different illuminant.
    /// Uses Bradford chromatic adaptation.
    pub fn simulate_illuminant(lab: &Lab, from: XYZ, to: XYZ) -> Lab {
        let xyz = lab.to_xyz(from);
        let adapted = super::chromatic_adaptation::bradford_adapt(xyz, from, to);
        adapted.to_lab(to)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_d65_white_point_from_weighting() {
        // Perfect diffuser (100% reflectance across all wavelengths)
        let reflectance = [1.0f32; 41];
        let xyz = XYZ::from_reflectance_10nm(&reflectance);

        // Expected D65 white point (2-degree)
        // X = 95.047
        // Y = 100.000
        // Z = 108.883
        assert!(
            (xyz.y - 100.0).abs() < 0.05,
            "Y should be 100, got {}",
            xyz.y
        );
        assert!(
            (xyz.x - 95.05).abs() < 0.05,
            "X should be ~95.05, got {}",
            xyz.x
        );
        assert!(
            (xyz.z - 108.88).abs() < 0.05,
            "Z should be ~108.88, got {}",
            xyz.z
        );
    }

    #[test]
    fn test_xyz_to_lab_d65() {
        let white = illuminant::D65;
        let lab = white.to_lab(white);

        // Lab of the white point itself should be (100, 0, 0)
        assert!((lab.l - 100.0).abs() < 1e-4);
        assert!(lab.a.abs() < 1e-4);
        assert!(lab.b.abs() < 1e-4);
    }
}