Atmo

def update_density_ratio(self) -> None:
    """Recompute density ratio for changed humidity."""
    self._density_ratio = Atmo.calculate_air_density(self._t0, self._p0, self.humidity) / cStandardDensityMetric

def temperature_at_altitude(self, altitude: float) -> float:
    """Interpolate temperature (°C) at altitude using lapse rate.

    Args:
        altitude: Altitude above mean sea level (ft).

    Returns:
        Temperature in degrees Celsius (bounded by model lower limit).
    """
    t = (altitude - self._a0) * cLapseRateKperFoot + self._t0
    if t < Atmo.cLowestTempC:
        t = Atmo.cLowestTempC
        warnings.warn(
            f"Temperature interpolated from altitude fell below minimum model limit. Bounded at {cLowestTempF}°F.",
            RuntimeWarning,
        )
    return t

def pressure_at_altitude(self, altitude: float) -> float:
    """Interpolate pressure (hPa) at altitude using barometric formula.

    Args:
        altitude: Altitude above mean sea level (ft).

    Returns:
        Pressure in hPa.
    """
    return self._p0 * math.pow(
        1 + cLapseRateKperFoot * (altitude - self._a0) / (self._t0 + cDegreesCtoK), cPressureExponent
    )

def get_density_and_mach_for_altitude(self, altitude: float) -> Tuple[float, float]:
    """Compute density ratio and Mach (fps) for the specified altitude.

    Uses lapse-rate interpolation unless altitude is within 30 ft of the base altitude,
        in which case the initial cached values are used for performance.

    Args:
        altitude: Altitude above mean sea level (ft).

    Returns:
        Tuple (density_ratio, mach_fps).
    """
    if math.fabs(self._a0 - altitude) < 30:  # fast path near base altitude
        return self._density_ratio, self._mach

    if altitude > 36089:  # troposphere limit ~36k ft
        warnings.warn(
            "Density request for altitude above modeled troposphere. Atmospheric model not valid here.",
            RuntimeWarning,
        )

    t_k = self.temperature_at_altitude(altitude) + cDegreesCtoK
    mach = Velocity.MPS(Atmo.machK(t_k)) >> Velocity.FPS
    p = self.pressure_at_altitude(altitude)
    density_delta = ((self._t0 + cDegreesCtoK) * p) / (self._p0 * t_k)
    density_ratio = self._density_ratio * density_delta
    # Alternative simplified exponential model (retained for reference):
    # density_ratio = self._density_ratio * math.exp(-(altitude - self._a0) / 34122)
    return density_ratio, mach

@staticmethod
def standard_temperature(altitude: Distance) -> Temperature:
    """ICAO standard temperature for altitude (valid to ~36,000 ft)."""
    return Temperature.Fahrenheit(cStandardTemperatureF + (altitude >> Distance.Foot) * cLapseRateImperial)

@staticmethod
def standard_pressure(altitude: Distance) -> Pressure:
    """ICAO standard pressure for altitude (valid to ~36,000 ft)."""
    return Pressure.hPa(
        cStandardPressureMetric
        * math.pow(
            1 + cLapseRateMetric * (altitude >> Distance.Meter) / (cStandardTemperatureC + cDegreesCtoK),
            cPressureExponent,
        )
    )

@staticmethod
def icao(
    altitude: Union[float, Distance] = 0,
    temperature: Optional[Temperature] = None,
    humidity: float = cStandardHumidity,
) -> Atmo:
    """Create a standard ICAO atmosphere at altitude.

    Args:
        altitude: Altitude (defaults to sea level).
        temperature: Optional override temperature (defaults to standard at altitude).
        humidity: Relative humidity (fraction or percent). Defaults to standard humidity.

    Returns:
        Atmo instance representing standard atmosphere at altitude.
    """
    altitude = PreferredUnits.distance(altitude)
    temperature = temperature or Atmo.standard_temperature(altitude)
    pressure = Atmo.standard_pressure(altitude)
    return Atmo(altitude, pressure, temperature, humidity)

@staticmethod
def machF(fahrenheit: float) -> float:
    """Mach 1 (fps) for given Fahrenheit temperature."""
    if fahrenheit < -cDegreesFtoR:
        bad_temp = fahrenheit
        fahrenheit = cLowestTempF
        warnings.warn(f"Invalid temperature: {bad_temp}°F. Adjusted to ({cLowestTempF}°F).", RuntimeWarning)
    return math.sqrt(fahrenheit + cDegreesFtoR) * cSpeedOfSoundImperial

@staticmethod
def machC(celsius: float) -> float:
    """Mach 1 (m/s) for given Celsius temperature."""
    if celsius < -cDegreesCtoK:
        bad_temp = celsius
        celsius = Atmo.cLowestTempC
        warnings.warn(f"Invalid temperature: {bad_temp}°C. Adjusted to ({celsius}°C).", RuntimeWarning)
    return Atmo.machK(celsius + cDegreesCtoK)

@staticmethod
def machK(kelvin: float) -> float:
    """Mach 1 (m/s) for given Kelvin temperature."""
    if kelvin < 0:
        bad_temp = kelvin
        kelvin = Atmo.cLowestTempC + cDegreesCtoK
        warnings.warn(f"Invalid temperature: {bad_temp}K. Adjusted to ({kelvin}K).", RuntimeWarning)
    return math.sqrt(kelvin) * cSpeedOfSoundMetric

@staticmethod
def calculate_air_density(t: float, p_hpa: float, humidity: float) -> float:
    """Air density from temperature (°C), pressure (hPa), and humidity.

    Args:
        t: Temperature in degrees Celsius.
        p_hpa: Pressure in hPa (hectopascals). Internally converted to Pa.
        humidity: Relative humidity (fraction or percent).

    Returns:
        Air density in kg/m³.

    Notes:
        - Divide result by `cDensityImperialToMetric` to get density in lb/ft³.
        - Source: CIPM-2007 (https://www.nist.gov/system/files/documents/calibrations/CIPM-2007.pdf)
    """
    R = 8.314472  # J/(mol·K), universal gas constant
    M_a = 28.96546e-3  # kg/mol, molar mass of dry air
    M_v = 18.01528e-3  # kg/mol, molar mass of water vapor

    def saturation_vapor_pressure(T):  # noqa: N802 (retain formula variable naming)
        A = [1.2378847e-5, -1.9121316e-2, 33.93711047, -6.3431645e3]
        return math.exp(A[0] * T**2 + A[1] * T + A[2] + A[3] / T)

    def enhancement_factor(p, T):  # noqa: N802
        alpha = 1.00062
        beta = 3.14e-8
        gamma = 5.6e-7
        return alpha + beta * p + gamma * T**2

    def compressibility_factor(p, T, x_v):  # noqa: N802
        a0 = 1.58123e-6
        a1 = -2.9331e-8
        a2 = 1.1043e-10
        b0 = 5.707e-6
        b1 = -2.051e-8
        c0 = 1.9898e-4
        c1 = -2.376e-6
        d = 1.83e-11
        e = -0.765e-8
        t_l = T - cDegreesCtoK
        Z = (
            1
            - (p / T) * (a0 + a1 * t_l + a2 * t_l**2 + (b0 + b1 * t_l) * x_v + (c0 + c1 * t_l) * x_v**2)
            + (p / T) ** 2 * (d + e * x_v**2)
        )
        return Z

    # Normalize humidity to fraction [0..1]
    rh = float(humidity)
    rh_frac = rh / 100.0 if rh > 1.0 else rh
    rh_frac = max(0.0, min(1.0, rh_frac))

    # Convert inputs for CIPM equations
    T_K = t + cDegreesCtoK  # Kelvin
    p = float(p_hpa) * 100.0  # hPa -> Pa

    # Calculation of saturated vapor pressure and enhancement factor
    p_sv = saturation_vapor_pressure(T_K)  # Pa (saturated vapor pressure)
    f = enhancement_factor(p, t)  # Enhancement factor (p in Pa, t in °C)

    # Partial pressure of water vapor and mole fraction
    p_v = rh_frac * f * p_sv  # Pa
    x_v = p_v / p  # Mole fraction of water vapor

    # Calculation of compressibility factor
    Z = compressibility_factor(p, T_K, x_v)
    return (p * M_a) / (Z * R * T_K) * (1.0 - x_v * (1.0 - M_v / M_a))