//! MovingBoundaryHX - Zone Discretization Heat Exchanger Component
//!
//! A heat exchanger component that discretizes the heat transfer area into
//! phase zones (superheated, two-phase, subcooled) for more accurate modeling
//! of refrigerant-side heat transfer.

use super::bphx_correlation::CorrelationSelector;
use super::bphx_geometry::BphxGeometry;
use super::eps_ntu::EpsNtuModel;
use super::exchanger::HeatExchanger;
use crate::state_machine::{CircuitId, OperationalState, StateManageable};
use crate::{
    Component, ComponentError, ConnectedPort, JacobianBuilder, ResidualVector, StateSlice,
};
use entropyk_core::{Enthalpy, MassFlow, Power};
use std::cell::{Cell, RefCell};
use std::sync::Arc;

/// Zone type for refrigerant-side phase classification
#[derive(Debug, Clone, Copy, PartialEq, Eq, Default)]
pub enum ZoneType {
    /// Superheated vapor zone (T > Tsat)
    Superheated,
    /// Two-phase zone (mixture of liquid and vapor)
    #[default]
    TwoPhase,
    /// Subcooled liquid zone (T < Tsat)
    Subcooled,
}

/// Zone boundary with relative position and zone type
#[derive(Debug, Clone)]
pub struct ZoneBoundary {
    /// Relative position along the heat exchanger (0.0 to 1.0)
    pub position: f64,
    /// Zone type at this boundary
    pub zone_type: ZoneType,
    /// UA value for this zone (W/K)
    pub ua: f64,
    /// Hot-side temperature at this boundary (K)
    pub t_hot: f64,
    /// Cold-side temperature at this boundary (K)
    pub t_cold: f64,
    /// Vapor quality at this boundary (0-1 for two-phase)
    pub quality: f64,
}

/// Zone discretization result containing all zones and summary data
#[derive(Debug, Clone, Default)]
pub struct ZoneDiscretization {
    /// List of zone boundaries (ordered by position)
    pub boundaries: Vec<ZoneBoundary>,
    /// Total UA (sum of all zone UAs) (W/K)
    pub total_ua: f64,
    /// Pinch temperature (minimum temperature difference) (K)
    pub pinch_temp: f64,
    /// Position of pinch point (relative, 0.0 to 1.0)
    pub pinch_position: f64,
}

/// Cache for MovingBoundaryHX calculations
#[derive(Debug, Clone, Default)]
pub struct MovingBoundaryCache {
    /// Whether the cache is valid and initialized
    pub valid: bool,
    /// Reference pressure for cache validity (Pa)
    pub p_ref: f64,
    /// Reference mass flow for cache validity (kg/s)
    pub m_ref: f64,
    /// Cached liquid saturation enthalpy (J/kg)
    pub h_sat_l: f64,
    /// Cached vapor saturation enthalpy (J/kg)
    pub h_sat_v: f64,
    /// Cached zone discretization result
    pub discretization: ZoneDiscretization,
}

impl MovingBoundaryCache {
    /// Checks if the cache remains valid given the current pressure and mass flow.
    /// Cache is valid if pressure deviates < 5% and mass flow deviates < 10%.
    pub fn is_valid_for(&self, p_current: f64, m_current: f64) -> bool {
        if !self.valid {
            return false;
        }

        let p_dev = (p_current - self.p_ref).abs() / self.p_ref.max(1e-10);
        let m_dev = (m_current - self.m_ref).abs() / self.m_ref.max(1e-10);

        p_dev < 0.05 && m_dev < 0.10
    }
}

/// MovingBoundaryHX - Zone discretization heat exchanger component
pub struct MovingBoundaryHX {
    inner: HeatExchanger<EpsNtuModel>,
    geometry: BphxGeometry,
    _correlation_selector: CorrelationSelector,
    refrigerant_id: String,
    secondary_fluid_id: String,
    fluid_backend: Option<Arc<dyn entropyk_fluids::FluidBackend>>,
    // Discretization parameters
    n_discretization: usize,
    cache: RefCell<MovingBoundaryCache>,

    // Internal state caching
    _last_htc: Cell<f64>,
    _last_validity_warning: Cell<bool>,
}

impl Default for MovingBoundaryHX {
    fn default() -> Self {
        Self::new()
    }
}

impl MovingBoundaryHX {
    /// Creates a new `MovingBoundaryHX` with default settings and 51 discretization points.
    pub fn new() -> Self {
        let geometry = BphxGeometry::from_dh_area(0.003, 0.5, 20);
        let model = EpsNtuModel::counter_flow(1000.0);

        Self {
            inner: HeatExchanger::new(model, "MovingBoundaryHX"),
            geometry,
            _correlation_selector: CorrelationSelector::default(),
            refrigerant_id: String::new(),
            secondary_fluid_id: String::new(),
            fluid_backend: None,
            n_discretization: 51,
            cache: RefCell::new(MovingBoundaryCache::default()),
            _last_htc: Cell::new(0.0),
            _last_validity_warning: Cell::new(false),
        }
    }

    /// Returns the number of discretization points.
    pub fn n_discretization(&self) -> usize {
        self.n_discretization
    }

    /// Sets the number of discretization points and returns self.
    pub fn with_discretization(mut self, n: usize) -> Self {
        self.n_discretization = n;
        self
    }

    /// Sets the geometry specification.
    pub fn with_geometry(mut self, geometry: BphxGeometry) -> Self {
        self.geometry = geometry;
        self
    }

    /// Sets the refrigerant fluid identifier.
    pub fn with_refrigerant(mut self, fluid: impl Into<String>) -> Self {
        self.refrigerant_id = fluid.into();
        self
    }

    /// Sets the secondary fluid identifier.
    pub fn with_secondary_fluid(mut self, fluid: impl Into<String>) -> Self {
        self.secondary_fluid_id = fluid.into();
        self
    }

    /// Attaches a fluid backend and returns self.
    pub fn with_fluid_backend(mut self, backend: Arc<dyn entropyk_fluids::FluidBackend>) -> Self {
        self.fluid_backend = Some(backend.clone());
        self.inner = self.inner.with_fluid_backend(backend);
        self
    }
}

impl Component for MovingBoundaryHX {
    fn n_equations(&self) -> usize {
        3
    }

    fn compute_residuals(
        &self,
        state: &StateSlice,
        residuals: &mut ResidualVector,
    ) -> Result<(), ComponentError> {
        let (p, m_refrig, t_sec_in, t_sec_out) = if let (Some(hot), Some(cold)) =
            (self.inner.hot_conditions(), self.inner.cold_conditions())
        {
            (
                hot.pressure_pa(),
                hot.mass_flow_kg_s(),
                cold.temperature_k(),
                cold.temperature_k() + 5.0,
            )
        } else {
            (500_000.0, 0.1, 300.0, 320.0)
        };

        // Extract enthalpies exactly as HeatExchanger does:
        let enthalpies = self.port_enthalpies(state)?;
        let h_in = enthalpies
            .get(0)
            .map(|h| h.to_joules_per_kg())
            .unwrap_or(400_000.0);
        let h_out = enthalpies
            .get(1)
            .map(|h| h.to_joules_per_kg())
            .unwrap_or(200_000.0);

        let mut cache = self.cache.borrow_mut();
        let use_cache = cache.is_valid_for(p, m_refrig);

        if !use_cache {
            let (disc, h_sat_l, h_sat_v) =
                self.identify_zones(h_in, h_out, p, t_sec_in, t_sec_out)?;
            cache.valid = true;
            cache.p_ref = p;
            cache.m_ref = m_refrig;
            cache.h_sat_l = h_sat_l;
            cache.h_sat_v = h_sat_v;
            cache.discretization = disc;
        }

        let total_ua = cache.discretization.total_ua;
        let base_ua = self.inner.ua_nominal();
        let custom_ua_scale = if base_ua > 0.0 {
            total_ua / base_ua
        } else {
            1.0
        };

        self.inner
            .compute_residuals_with_ua_scale(state, residuals, custom_ua_scale)
    }

    fn jacobian_entries(
        &self,
        state: &StateSlice,
        jacobian: &mut JacobianBuilder,
    ) -> Result<(), ComponentError> {
        self.inner.jacobian_entries(state, jacobian)
    }

    fn get_ports(&self) -> &[ConnectedPort] {
        self.inner.get_ports()
    }

    fn set_calib_indices(&mut self, indices: entropyk_core::CalibIndices) {
        self.inner.set_calib_indices(indices);
    }

    fn port_mass_flows(&self, state: &StateSlice) -> Result<Vec<MassFlow>, ComponentError> {
        self.inner.port_mass_flows(state)
    }

    fn port_enthalpies(&self, state: &StateSlice) -> Result<Vec<Enthalpy>, ComponentError> {
        self.inner.port_enthalpies(state)
    }

    fn energy_transfers(&self, state: &StateSlice) -> Option<(Power, Power)> {
        self.inner.energy_transfers(state)
    }

    fn set_fluid_backend_from_builder(&mut self, backend: std::sync::Arc<dyn entropyk_fluids::FluidBackend>) {
        if self.fluid_backend.is_none() {
            self.fluid_backend = Some(backend.clone());
            self.inner.set_fluid_backend_from_builder(backend);
        }
    }

    fn update_calib_factor(&mut self, factor: &str, value: f64) -> bool {
        self.inner.update_calib_factor(factor, value)
    }
}

impl StateManageable for MovingBoundaryHX {
    fn state(&self) -> OperationalState {
        self.inner.state()
    }

    fn set_state(&mut self, state: OperationalState) -> Result<(), ComponentError> {
        self.inner.set_state(state)
    }

    fn can_transition_to(&self, target: OperationalState) -> bool {
        self.inner.can_transition_to(target)
    }

    fn circuit_id(&self) -> &CircuitId {
        self.inner.circuit_id()
    }

    fn set_circuit_id(&mut self, circuit_id: CircuitId) {
        self.inner.set_circuit_id(circuit_id);
    }
}

impl MovingBoundaryHX {
    /// Identifies the phase zones along the heat exchanger and calculates boundaries.
    pub fn identify_zones(
        &self,
        h_refrig_in: f64,
        h_refrig_out: f64,
        p_refrig: f64,
        t_secondary_in: f64,
        t_secondary_out: f64,
    ) -> Result<(ZoneDiscretization, f64, f64), ComponentError> {
        let backend = self.fluid_backend.as_ref().ok_or_else(|| {
            ComponentError::CalculationFailed("No FluidBackend configured".to_string())
        })?;
        let fluid = entropyk_fluids::FluidId::new(&self.refrigerant_id);
        let p = entropyk_core::Pressure::from_pascals(p_refrig);

        let h_sat_l = backend
            .property(
                fluid.clone(),
                entropyk_fluids::Property::Enthalpy,
                entropyk_fluids::FluidState::from_px(p, entropyk_fluids::Quality::new(0.0)),
            )
            .map_err(|e| ComponentError::CalculationFailed(format!("h_sat_l failed: {}", e)))?;
        let h_sat_v = backend
            .property(
                fluid.clone(),
                entropyk_fluids::Property::Enthalpy,
                entropyk_fluids::FluidState::from_px(p, entropyk_fluids::Quality::new(1.0)),
            )
            .map_err(|e| ComponentError::CalculationFailed(format!("h_sat_v failed: {}", e)))?;

        let mut boundaries = Vec::new();

        // Calculate transition positions and types
        let is_condensing = h_refrig_in > h_refrig_out;

        // Add inlet boundary
        let inlet_type = if h_refrig_in > h_sat_v + 1e-3 {
            ZoneType::Superheated
        } else if h_refrig_in < h_sat_l - 1e-3 {
            ZoneType::Subcooled
        } else {
            ZoneType::TwoPhase
        };
        boundaries.push(self.create_boundary(
            0.0,
            h_refrig_in,
            p_refrig,
            inlet_type,
            t_secondary_in,
            h_sat_l,
            h_sat_v,
        )?);

        let (h_min, h_max) = if is_condensing {
            (h_refrig_out, h_refrig_in)
        } else {
            (h_refrig_in, h_refrig_out)
        };

        if h_min < h_sat_l && h_max > h_sat_l {
            let pos = (h_sat_l - h_refrig_in) / (h_refrig_out - h_refrig_in);
            let t_sec = t_secondary_in + pos * (t_secondary_out - t_secondary_in);
            // After sat_l, type is SC (if condensing) or TP (if evaporating)
            let post_type = if is_condensing {
                ZoneType::Subcooled
            } else {
                ZoneType::TwoPhase
            };
            boundaries.push(
                self.create_boundary(pos, h_sat_l, p_refrig, post_type, t_sec, h_sat_l, h_sat_v)?,
            );
        }

        if h_min < h_sat_v && h_max > h_sat_v {
            let pos = (h_sat_v - h_refrig_in) / (h_refrig_out - h_refrig_in);
            let t_sec = t_secondary_in + pos * (t_secondary_out - t_secondary_in);
            // After sat_v, type is TP (if condensing) or SH (if evaporating)
            let post_type = if is_condensing {
                ZoneType::TwoPhase
            } else {
                ZoneType::Superheated
            };
            boundaries.push(
                self.create_boundary(pos, h_sat_v, p_refrig, post_type, t_sec, h_sat_l, h_sat_v)?,
            );
        }

        // Add outlet boundary
        let outlet_type = if h_refrig_out > h_sat_v + 1e-3 {
            ZoneType::Superheated
        } else if h_refrig_out < h_sat_l - 1e-3 {
            ZoneType::Subcooled
        } else {
            ZoneType::TwoPhase
        };
        boundaries.push(self.create_boundary(
            1.0,
            h_refrig_out,
            p_refrig,
            outlet_type,
            t_secondary_out,
            h_sat_l,
            h_sat_v,
        )?);

        // Sort boundaries by position
        boundaries.sort_by(|a, b| a.position.partial_cmp(&b.position).unwrap());

        // Calculate UA for each zone
        let mut total_ua = 0.0;
        for i in 0..boundaries.len() - 1 {
            let ua_zone = self.compute_zone_ua(&boundaries[i], &boundaries[i + 1])?;
            boundaries[i].ua = ua_zone;
            total_ua += ua_zone;
        }

        let (pinch_temp, pinch_pos) = self.calculate_pinch(&boundaries);

        Ok((
            ZoneDiscretization {
                boundaries,
                total_ua,
                pinch_temp,
                pinch_position: pinch_pos,
            },
            h_sat_l,
            h_sat_v,
        ))
    }

    fn compute_zone_ua(&self, b1: &ZoneBoundary, b2: &ZoneBoundary) -> Result<f64, ComponentError> {
        let area_zone = self.geometry.area * (b2.position - b1.position);
        if area_zone <= 1e-10 {
            return Ok(0.0);
        }

        // Without access to fluid phase properties and geometry correlation,
        // we use a simplified approximation based on zone type.
        // A true implementation would query self.correlation_selector
        let h_refrig = match b1.zone_type {
            ZoneType::TwoPhase => 5000.0,   // Boiling or condensation
            ZoneType::Superheated => 500.0, // Vapor
            ZoneType::Subcooled => 1500.0,  // Liquid
        };
        let h_secondary = 5000.0; // Generally high for water/glycol

        let u_overall = 1.0 / (1.0 / h_refrig + 1.0 / h_secondary);

        Ok(u_overall * area_zone)
    }

    fn calculate_pinch(&self, boundaries: &[ZoneBoundary]) -> (f64, f64) {
        let mut min_dt = f64::MAX;
        let mut pinch_pos = 0.0;

        for b in boundaries {
            let dt = (b.t_hot - b.t_cold).abs();
            if dt < min_dt {
                min_dt = dt;
                pinch_pos = b.position;
            }
        }

        (min_dt, pinch_pos)
    }

    fn create_boundary(
        &self,
        pos: f64,
        h: f64,
        p: f64,
        zone_type: ZoneType,
        t_sec: f64,
        h_sat_l: f64,
        h_sat_v: f64,
    ) -> Result<ZoneBoundary, ComponentError> {
        let quality = if h_sat_v > h_sat_l {
            ((h - h_sat_l) / (h_sat_v - h_sat_l)).clamp(0.0, 1.0)
        } else {
            0.0
        };

        let t_refrig = if let Some(backend) = &self.fluid_backend {
            let fluid = entropyk_fluids::FluidId::new(&self.refrigerant_id);
            backend
                .property(
                    fluid,
                    entropyk_fluids::Property::Temperature,
                    entropyk_fluids::FluidState::from_ph(
                        entropyk_core::Pressure::from_pascals(p),
                        entropyk_core::Enthalpy::from_joules_per_kg(h),
                    ),
                )
                .map_err(|e| ComponentError::CalculationFailed(format!("T_refrig failed: {}", e)))?
        } else {
            300.0
        };

        Ok(ZoneBoundary {
            position: pos,
            zone_type,
            ua: 0.0,
            t_hot: if t_sec > t_refrig { t_sec } else { t_refrig },
            t_cold: if t_sec > t_refrig { t_refrig } else { t_sec },
            quality,
        })
    }
}

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

    #[test]
    fn test_zone_type_enum_exists() {
        let zone = ZoneType::Superheated;
        assert_eq!(zone, ZoneType::Superheated);

        let zone = ZoneType::TwoPhase;
        assert_eq!(zone, ZoneType::TwoPhase);

        let zone = ZoneType::Subcooled;
        assert_eq!(zone, ZoneType::Subcooled);
    }

    #[test]
    fn test_zone_boundary_struct_exists() {
        let boundary = ZoneBoundary {
            position: 0.5,
            zone_type: ZoneType::TwoPhase,
            ua: 1000.0,
            t_hot: 300.0,
            t_cold: 290.0,
            quality: 0.5,
        };

        assert!((boundary.position - 0.5).abs() < 1e-10);
        assert_eq!(boundary.zone_type, ZoneType::TwoPhase);
        assert!((boundary.ua - 1000.0).abs() < 1e-10);
    }

    #[test]
    fn test_moving_boundary_hx_with_fluids() {
        let hx = MovingBoundaryHX::new()
            .with_refrigerant("R410A")
            .with_secondary_fluid("Water");
        assert_eq!(hx.refrigerant_id, "R410A");
        assert_eq!(hx.secondary_fluid_id, "Water");
    }

    #[test]
    fn test_identify_zones_basic() {
        use entropyk_core::Pressure;
        use entropyk_fluids::{
            CriticalPoint, FluidError, FluidId, FluidResult, Phase, ThermoState,
        };

        struct MockBackend {
            h_sat_l: f64,
            h_sat_v: f64,
            t_sat: f64,
        }

        impl entropyk_fluids::FluidBackend for MockBackend {
            fn property(
                &self,
                _fluid: FluidId,
                property: entropyk_fluids::Property,
                state: entropyk_fluids::FluidState,
            ) -> FluidResult<f64> {
                match property {
                    entropyk_fluids::Property::Temperature => Ok(self.t_sat),
                    entropyk_fluids::Property::Enthalpy => {
                        let q = match state {
                            entropyk_fluids::FluidState::PressureQuality(_, q) => Some(q.value()),
                            _ => None,
                        };
                        match q {
                            Some(0.0) => Ok(self.h_sat_l),
                            Some(1.0) => Ok(self.h_sat_v),
                            _ => Ok(self.h_sat_v),
                        }
                    }
                    _ => Err(FluidError::UnsupportedProperty {
                        property: format!("{:?}", property),
                    }),
                }
            }

            fn critical_point(&self, fluid: FluidId) -> FluidResult<CriticalPoint> {
                Err(FluidError::NoCriticalPoint { fluid: fluid.0 })
            }

            fn is_fluid_available(&self, _fluid: &FluidId) -> bool {
                true
            }
            fn phase(
                &self,
                _fluid: FluidId,
                _state: entropyk_fluids::FluidState,
            ) -> FluidResult<Phase> {
                Ok(Phase::Unknown)
            }
            fn full_state(
                &self,
                _fluid: FluidId,
                _p: Pressure,
                _h: Enthalpy,
            ) -> FluidResult<ThermoState> {
                Err(FluidError::UnsupportedProperty {
                    property: "full_state".to_string(),
                })
            }
            fn list_fluids(&self) -> Vec<FluidId> {
                vec![]
            }
        }

        let backend = MockBackend {
            h_sat_l: 200_000.0,
            h_sat_v: 400_000.0,
            t_sat: 280.0,
        };

        let hx = MovingBoundaryHX::new()
            .with_refrigerant("R410A")
            .with_fluid_backend(Arc::new(backend));

        // Condensing: 450,000 (SH) -> 150,000 (SC)
        let result = hx.identify_zones(450_000.0, 150_000.0, 500_000.0, 300.0, 320.0);
        assert!(result.is_ok());
        let (disc, h_sat_l_res, h_sat_v_res) = result.unwrap();

        assert_eq!(h_sat_l_res, 200_000.0);
        assert_eq!(h_sat_v_res, 400_000.0);

        // Should have 4 boundaries: inlet(SH), sat_v(SH/TP), sat_l(TP/SC), outlet(SC)
        assert_eq!(disc.boundaries.len(), 4);
        assert_eq!(disc.boundaries[0].zone_type, ZoneType::Superheated);
        assert_eq!(disc.boundaries[1].zone_type, ZoneType::TwoPhase);
        assert_eq!(disc.boundaries[2].zone_type, ZoneType::Subcooled);
        assert_eq!(disc.boundaries[3].zone_type, ZoneType::Subcooled);

        // Total UA should be positive
        assert!(disc.total_ua > 0.0);
    }

    #[test]
    fn test_cache_is_valid_for() {
        let mut cache = MovingBoundaryCache {
            valid: true,
            p_ref: 100_000.0,
            m_ref: 1.0,
            h_sat_l: 100.0,
            h_sat_v: 200.0,
            discretization: ZoneDiscretization::default(),
        };

        // Identical
        assert!(cache.is_valid_for(100_000.0, 1.0));

        // P < 5% deviation (104,000 is 4%)
        assert!(cache.is_valid_for(104_000.0, 1.0));

        // P > 5% deviation (106,000 is 6%)
        assert!(!cache.is_valid_for(106_000.0, 1.0));

        // M < 10% deviation (1.09 is 9%)
        assert!(cache.is_valid_for(100_000.0, 1.09));

        // M > 10% deviation (1.11 is 11%)
        assert!(!cache.is_valid_for(100_000.0, 1.11));

        // Invalid if explicitly invalid
        cache.valid = false;
        assert!(!cache.is_valid_for(100_000.0, 1.0));
    }

    #[test]
    fn test_compute_residuals_uses_cache() {
        use crate::{Component, ResidualVector};
        use entropyk_core::Pressure;
        use entropyk_fluids::{
            CriticalPoint, FluidError, FluidId, FluidResult, Phase, ThermoState,
        };

        struct TrackingMockBackend {
            pub calls: std::sync::atomic::AtomicUsize,
        }

        impl entropyk_fluids::FluidBackend for TrackingMockBackend {
            fn property(
                &self,
                _fluid: FluidId,
                _property: entropyk_fluids::Property,
                _state: entropyk_fluids::FluidState,
            ) -> FluidResult<f64> {
                self.calls.fetch_add(1, std::sync::atomic::Ordering::SeqCst);
                Ok(100.0)
            }
            fn critical_point(&self, _fluid: FluidId) -> FluidResult<CriticalPoint> {
                Err(FluidError::NoCriticalPoint {
                    fluid: "".to_string(),
                })
            }
            fn is_fluid_available(&self, _fluid: &FluidId) -> bool {
                true
            }
            fn phase(
                &self,
                _fluid: FluidId,
                _state: entropyk_fluids::FluidState,
            ) -> FluidResult<Phase> {
                Ok(Phase::Unknown)
            }
            fn full_state(
                &self,
                _fluid: FluidId,
                _p: Pressure,
                _h: Enthalpy,
            ) -> FluidResult<ThermoState> {
                Err(FluidError::UnsupportedProperty {
                    property: "full_state".to_string(),
                })
            }
            fn list_fluids(&self) -> Vec<FluidId> {
                vec![]
            }
        }

        let backend = Arc::new(TrackingMockBackend {
            calls: std::sync::atomic::AtomicUsize::new(0),
        });

        let hx = MovingBoundaryHX::new()
            .with_refrigerant("R410A")
            .with_fluid_backend(backend.clone());

        let state = vec![500_000.0, 400_000.0];
        let mut residuals = vec![0.0; 3];

        // First call should calculate property (backend calls)
        let _ = hx.compute_residuals(&state, &mut residuals);
        let calls_first = backend.calls.load(std::sync::atomic::Ordering::SeqCst);
        assert!(calls_first > 0);

        // Second call with same state should use cache -> 0 new backend calls
        let _ = hx.compute_residuals(&state, &mut residuals);
        let calls_second = backend.calls.load(std::sync::atomic::Ordering::SeqCst);
        assert_eq!(calls_first, calls_second); // Calls remained the same because cache was used
    }

    #[test]
    fn test_performance_speedup() {
        use crate::{Component, ResidualVector};
        use entropyk_core::Pressure;
        use entropyk_fluids::{
            CriticalPoint, FluidError, FluidId, FluidResult, Phase, ThermoState,
        };
        use std::time::Instant;

        struct SlowMockBackend;

        impl entropyk_fluids::FluidBackend for SlowMockBackend {
            fn property(
                &self,
                _fluid: FluidId,
                _property: entropyk_fluids::Property,
                _state: entropyk_fluids::FluidState,
            ) -> FluidResult<f64> {
                // Simulate somewhat slow fluid property calculation
                std::thread::sleep(std::time::Duration::from_micros(10));
                Ok(100.0)
            }
            fn critical_point(&self, _fluid: FluidId) -> FluidResult<CriticalPoint> {
                Err(FluidError::NoCriticalPoint {
                    fluid: "".to_string(),
                })
            }
            fn is_fluid_available(&self, _fluid: &FluidId) -> bool {
                true
            }
            fn phase(
                &self,
                _fluid: FluidId,
                _state: entropyk_fluids::FluidState,
            ) -> FluidResult<Phase> {
                Ok(Phase::Unknown)
            }
            fn full_state(
                &self,
                _fluid: FluidId,
                _p: Pressure,
                _h: Enthalpy,
            ) -> FluidResult<ThermoState> {
                Err(FluidError::UnsupportedProperty {
                    property: "full_state".to_string(),
                })
            }
            fn list_fluids(&self) -> Vec<FluidId> {
                vec![]
            }
        }

        let backend = Arc::new(SlowMockBackend);

        let hx = MovingBoundaryHX::new()
            .with_refrigerant("R410A")
            .with_fluid_backend(backend.clone());

        let state = vec![500_000.0, 400_000.0];
        let mut residuals = vec![0.0; 3];

        // First run (no cache)
        let start = Instant::now();
        let _ = hx.compute_residuals(&state, &mut residuals);
        let duration_uncached = start.elapsed();

        // Second run (cached)
        let start = Instant::now();
        let _ = hx.compute_residuals(&state, &mut residuals);
        let duration_cached = start.elapsed();

        println!("Uncached duration: {:?}", duration_uncached);
        println!("Cached duration: {:?}", duration_cached);

        let speedup = duration_uncached.as_secs_f64() / duration_cached.as_secs_f64().max(1e-9);
        println!("Speedup multiplier: {:.1}x", speedup);

        assert!(duration_cached < duration_uncached);
    }
}