//! Bilinear interpolation for 2D property tables.
//!
//! Provides C1-continuous interpolation suitable for solver Jacobian assembly.

use std::cmp::Ordering;

/// Performs bilinear interpolation on a 2D grid.
///
/// Given a rectangular grid with values at (p_idx, t_idx), interpolates
/// the value at (p, t) where p and t are in the grid's coordinate space.
/// Returns None if (p, t) is outside the grid bounds.
///
/// # Arguments
/// * `p_grid` - Pressure grid (must be sorted ascending)
/// * `t_grid` - Temperature grid (must be sorted ascending)
/// * `values` - 2D array [p_idx][t_idx], row-major
/// * `p` - Query pressure (Pa)
/// * `t` - Query temperature (K)
#[inline]
pub fn bilinear_interpolate(
    p_grid: &[f64],
    t_grid: &[f64],
    values: &[f64],
    p: f64,
    t: f64,
) -> Option<f64> {
    let n_p = p_grid.len();
    let n_t = t_grid.len();

    if n_p < 2 || n_t < 2 || values.len() != n_p * n_t {
        return None;
    }

    // Reject NaN to avoid panic in binary_search_by (Zero-Panic Policy)
    if !p.is_finite() || !t.is_finite() {
        return None;
    }

    // Find P indices (p_grid must be ascending)
    let p_idx = match p_grid.binary_search_by(|x| x.partial_cmp(&p).unwrap_or(Ordering::Equal)) {
        Ok(i) => {
            if i >= n_p - 1 {
                return None;
            }
            i
        }
        Err(i) => {
            if i == 0 || i >= n_p {
                return None;
            }
            i - 1
        }
    };

    // Find T indices
    let t_idx = match t_grid.binary_search_by(|x| x.partial_cmp(&t).unwrap_or(Ordering::Equal)) {
        Ok(i) => {
            if i >= n_t - 1 {
                return None;
            }
            i
        }
        Err(i) => {
            if i == 0 || i >= n_t {
                return None;
            }
            i - 1
        }
    };

    let p0 = p_grid[p_idx];
    let p1 = p_grid[p_idx + 1];
    let t0 = t_grid[t_idx];
    let t1 = t_grid[t_idx + 1];

    let dp = p1 - p0;
    let dt = t1 - t0;

    if dp <= 0.0 || dt <= 0.0 {
        return None;
    }

    let fp = (p - p0) / dp;
    let ft = (t - t0) / dt;

    // Clamp to [0,1] for edge cases
    let fp = fp.clamp(0.0, 1.0);
    let ft = ft.clamp(0.0, 1.0);

    let v00 = values[p_idx * n_t + t_idx];
    let v01 = values[p_idx * n_t + t_idx + 1];
    let v10 = values[(p_idx + 1) * n_t + t_idx];
    let v11 = values[(p_idx + 1) * n_t + t_idx + 1];

    let v0 = v00 * (1.0 - ft) + v01 * ft;
    let v1 = v10 * (1.0 - ft) + v11 * ft;

    Some(v0 * (1.0 - fp) + v1 * fp)
}

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

    #[test]
    fn test_bilinear_inside() {
        let p = [100000.0, 200000.0, 300000.0];
        let t = [250.0, 300.0, 350.0];
        let v = [
            1.0, 2.0, 3.0, // p=100k
            4.0, 5.0, 6.0, // p=200k
            7.0, 8.0, 9.0, // p=300k
        ];

        let result = bilinear_interpolate(&p, &t, &v, 200000.0, 300.0);
        assert!(result.is_some());
        assert!((result.unwrap() - 5.0).abs() < 1e-10);
    }

    #[test]
    fn test_bilinear_interpolated() {
        let p = [0.0, 1.0];
        let t = [0.0, 1.0];
        let v = [0.0, 1.0, 1.0, 2.0]; // v(0,0)=0, v(0,1)=1, v(1,0)=1, v(1,1)=2

        let result = bilinear_interpolate(&p, &t, &v, 0.5, 0.5);
        assert!(result.is_some());
        // At center: (0+1+1+2)/4 = 1.0
        assert!((result.unwrap() - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_bilinear_out_of_bounds() {
        let p = [100000.0, 200000.0];
        let t = [250.0, 300.0];
        let v = [1.0, 2.0, 3.0, 4.0];

        assert!(bilinear_interpolate(&p, &t, &v, 50000.0, 300.0).is_none());
        assert!(bilinear_interpolate(&p, &t, &v, 300000.0, 300.0).is_none());
    }

    #[test]
    fn test_bilinear_nan_rejected() {
        let p = [100000.0, 200000.0];
        let t = [250.0, 300.0];
        let v = [1.0, 2.0, 3.0, 4.0];

        assert!(bilinear_interpolate(&p, &t, &v, f64::NAN, 300.0).is_none());
        assert!(bilinear_interpolate(&p, &t, &v, 150000.0, f64::NAN).is_none());
        assert!(bilinear_interpolate(&p, &t, &v, f64::INFINITY, 300.0).is_none());
    }
}
