Atmo

def update_density_ratio(self) -> None:
    """Update the density ratio based on current conditions."""
    self._density_ratio = Atmo.calculate_air_density(self._t0, self._p0, self.humidity) / cStandardDensityMetric

def temperature_at_altitude(self, altitude: float) -> float:
    """Temperature at altitude interpolated from zero conditions using lapse rate.

    Args:
        altitude: ASL in ft

    Returns:
        temperature in °C
    """
    t = (altitude - self._a0) * cLapseRateKperFoot + self._t0
    if t < Atmo.cLowestTempC:
        t = Atmo.cLowestTempC
        warnings.warn(f"Temperature interpolated from altitude fell below minimum temperature limit.  "
                      f"Model not accurate here.  Temperature bounded at cLowestTempF: {cLowestTempF}°F.",
                      RuntimeWarning)
    return t

def pressure_at_altitude(self, altitude: float) -> float:
    """Pressure at altitude interpolated from zero conditions using lapse rate.

    Ref: https://en.wikipedia.org/wiki/Barometric_formula#Pressure_equations

    Args:
        altitude: ASL in ft

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

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

    Ref: https://en.wikipedia.org/wiki/Barometric_formula#Density_equations

    Args:
        altitude: ASL in units of feet.
            Note: Altitude above 36,000 ft not modelled this way.

    Returns:
        density ratio and Mach 1 (fps) for the specified altitude
    """
    # Within 30 ft of initial altitude use initial values to save compute
    if math.fabs(self._a0 - altitude) < 30:
        mach = self._mach
        density_ratio = self._density_ratio
    else:
        if altitude > 36089:
            warnings.warn("Density request for altitude above troposphere."
                          " Atmospheric model not valid here.", RuntimeWarning)
        t = self.temperature_at_altitude(altitude) + cDegreesCtoK
        mach = Velocity.MPS(Atmo.machK(t)) >> Velocity.FPS
        p = self.pressure_at_altitude(altitude)
        density_delta = ((self._t0 + cDegreesCtoK) * p) / (self._p0 * t)
        density_ratio = self._density_ratio * density_delta
        # # Alternative simplified model:
        # # Ref https://en.wikipedia.org/wiki/Density_of_air#Exponential_approximation
        # # see doc/'Air Density Models.svg' for comparison
        # density_ratio = self._density_ratio * math.exp(-(altitude - self._a0) / 34122)
    return density_ratio, mach

@staticmethod
def standard_temperature(altitude: Distance) -> Temperature:
    """Calculate ICAO standard temperature for altitude.

    Note: This model is only valid up to the troposphere (~36,000 ft).

    Args:
        altitude: ASL in units of feet.

    Returns:
        ICAO standard temperature for altitude
    """
    return Temperature.Fahrenheit(cStandardTemperatureF
                                  + (altitude >> Distance.Foot) * cLapseRateImperial)

@staticmethod
def standard_pressure(altitude: Distance) -> Pressure:
    """Calculate ICAO standard pressure for altitude.

    Note: This model only valid up to troposphere (~36,000 ft).
    Ref: https://en.wikipedia.org/wiki/Barometric_formula#Pressure_equations

    Args:
        altitude: Distance above sea level (ASL)

    Returns:
        ICAO standard pressure for altitude
    """
    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 Atmo instance of standard ICAO atmosphere at given altitude.

    Note: This model is only valid up to the troposphere (~36,000 ft).

    Args:
        altitude: relative to sea level.  Default is sea level (0 ft).
        temperature: air temperature.  Default is standard temperature at altitude.

    Returns:
        Atmo instance of standard ICAO atmosphere at given altitude.
        If temperature not specified uses standard temperature.
    """
    altitude = PreferredUnits.distance(altitude)
    if temperature is None:
        temperature = Atmo.standard_temperature(altitude)
    pressure = Atmo.standard_pressure(altitude)

    return Atmo(altitude, pressure, temperature, humidity)

@staticmethod
def machF(fahrenheit: float) -> float:
    """Calculate Mach 1 in fps for given Fahrenheit temperature.

    Args:
        fahrenheit: Fahrenheit temperature

    Returns:
        Mach 1 in fps for given 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:
    """Calculate Mach 1 in mps for given Celsius temperature.

    Args:
        celsius: Celsius temperature

    Returns:
        Mach 1 in mps for given 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:
    """Calculate Mach 1 in mps for given Kelvin temperature.

    Args:
        kelvin: Kelvin temperature

    Returns:
        Mach 1 in mps for given temperature
    """
    return math.sqrt(kelvin) * cSpeedOfSoundMetric

@staticmethod
def calculate_air_density(t: float, p_hpa: float, humidity: float) -> float:
    """Calculate air density from temperature, pressure, and humidity.

    Args:
        t: Temperature in degrees Celsius.
        p_hpa: Pressure in hPa (hectopascals). Internally converted to Pa.
        humidity: Relative humidity. Accepts either fraction [0..1] or percent [0..100].

    Returns:
        Air density in kg/m³.

    Notes:
        - Divide result by cDensityImperialToMetric to get density in lb/ft³
        - Source: 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):
        # Calculation of saturated vapor pressure according to CIPM 2007
        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):
        # Calculation of enhancement factor according to CIPM 2007
        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):
        # Calculation of compressibility factor according to CIPM 2007
        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 = T - cDegreesCtoK
        Z = 1 - (p / T) * (a0 + a1 * t + a2 * t ** 2 + (b0 + b1 * t) * x_v + (c0 + c1 * t) * 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)

    # Density (kg/m^3) using moist air composition and compressibility factor
    density = (p * M_a) / (Z * R * T_K) * (1.0 - x_v * (1.0 - M_v / M_a))
    return density