Repository: pyfar/pyfar (python)

Issue:
We need to implement calculations of commonly used constants
We need an implementation to calculate and store typically used constants, such as air density, speed of sound, etc.
Ideally, this could be implemented as a Dataclass or similar.

Pull request:
Reopen: Speed of sound ideal gas
### Which issue(s) are closed by this pull request?

Contributes to #428

### Changes proposed in this pull request:

- introduces magnus formula for saturation vapour pressure
- implement ideal gas speed of sound calculation using magnus Formula
- uses the new constants


Patch:
diff --git a/pyfar/constants/__init__.py b/pyfar/constants/__init__.py
index a429abd7..3bdca282 100644
--- a/pyfar/constants/__init__.py
+++ b/pyfar/constants/__init__.py
@@ -2,12 +2,19 @@
 
 from typing import Final
 
+from .constants import (
+    saturation_vapor_pressure_magnus,
+)
+
 from .speed_of_sound import (
     speed_of_sound_simple,
+    speed_of_sound_ideal_gas,
 )
 
 __all__ = [
+    'saturation_vapor_pressure_magnus',
     'speed_of_sound_simple',
+    'speed_of_sound_ideal_gas',
 ]
 
 
diff --git a/pyfar/constants/constants.py b/pyfar/constants/constants.py
new file mode 100644
index 00000000..8db6fb3d
--- /dev/null
+++ b/pyfar/constants/constants.py
@@ -0,0 +1,45 @@
+"""Constant calculation."""
+import numpy as np
+
+
+def saturation_vapor_pressure_magnus(temperature):
+    r"""
+    Calculate the saturation vapor pressure of water in Pascal using the
+    Magnus formula.
+
+    The Magnus formula is valid for temperatures between -45°C and 60°C [#]_.
+
+    .. math::
+
+        e_s = 610.94 \cdot \exp\left(\frac{17.625 \cdot T}{T + 243.04}\right)
+
+
+    Parameters
+    ----------
+    temperature : float, array_like
+        Temperature in degrees Celsius (°C).
+
+    Returns
+    -------
+    p_sat : float, array_like
+        Saturation vapor pressure in Pa.
+
+    References
+    ----------
+    .. [#] O. A. Alduchov and R. E. Eskridge, “Improved Magnus Form
+           Approximation of Saturation Vapor Pressure,” Journal of Applied
+           Meteorology and Climatology, vol. 35, no. 4, pp. 60-609, Apr. 1996
+    """
+
+    if not isinstance(temperature, (int, float, np.ndarray, list, tuple)):
+        raise TypeError(
+            'temperature must be a number or array of numbers')
+    if np.any(np.array(
+            temperature) < -45) or np.any(np.array(temperature) > 60):
+        raise ValueError("Temperature must be in the range of -45°C and 60°C.")
+    if isinstance(temperature, (np.ndarray, list, tuple)):
+        temperature = np.asarray(temperature, dtype=float)
+
+    # Eq. (21)
+    e_s = 6.1094 * np.exp((17.625 * temperature) / (temperature + 243.04))
+    return 100 * e_s
diff --git a/pyfar/constants/speed_of_sound.py b/pyfar/constants/speed_of_sound.py
index 75379772..49cca5a4 100644
--- a/pyfar/constants/speed_of_sound.py
+++ b/pyfar/constants/speed_of_sound.py
@@ -1,5 +1,6 @@
 """File containing all speed of sound calculation functions."""
 import numpy as np
+from . import constants
 import pyfar as pf
 
 
@@ -52,3 +53,94 @@ def speed_of_sound_simple(temperature):
     t_0 = pf.constants.absolute_zero_celsius
     c_ref = pf.constants.reference_speed_of_sound
     return c_ref*np.sqrt((temperature-t_0)/(t_ref-t_0))
+
+
+def speed_of_sound_ideal_gas(
+        temperature, relative_humidity, atmospheric_pressure=None,
+        saturation_vapor_pressure=None):
+    """Calculate speed of sound in air using the ideal gas law.
+
+    The speed of sound in air can be calculated based on chapter 6.3 in [#]_.
+    All input parameters must be broadcastable to the same shape.
+
+    Parameters
+    ----------
+    temperature : float, array_like
+        Temperature in degree Celsius.
+    relative_humidity : float, array_like
+        Relative humidity in the range of 0 to 1.
+    atmospheric_pressure : float, array_like, optional
+        Atmospheric pressure in Pascal, by default
+        :py:attr:`reference_atmospheric_pressure`.
+    saturation_vapor_pressure : float, array_like, optional
+        Saturation vapor pressure in Pascal.
+        If not given, the function
+        :py:func:`~pyfar.constants.saturation_vapor_pressure_magnus` is used.
+        Note that the valid temperature range therefore also depends on
+        :py:func:`~pyfar.constants.saturation_vapor_pressure_magnus`.
+
+    Returns
+    -------
+    speed_of_sound : float, array_like
+        Speed of sound in air in m/s.
+
+    References
+    ----------
+    .. [#] V. E. Ostashev and D. K. Wilson, Acoustics in Moving Inhomogeneous
+           Media, 2nd ed. London: CRC Press, 2015. doi: 10.1201/b18922.
+    """
+    if atmospheric_pressure is None:
+        atmospheric_pressure = pf.constants.reference_atmospheric_pressure
+    # check inputs
+    if not isinstance(temperature, (int, float, np.ndarray, list, tuple)):
+        raise TypeError(
+            'Temperature must be a number or array of numbers')
+    if not isinstance(
+            relative_humidity, (int, float, np.ndarray, list, tuple)):
+        raise TypeError(
+            'Relative humidity must be a number or array of numbers')
+    if not isinstance(
+            atmospheric_pressure, (int, float, np.ndarray, list, tuple)):
+        raise TypeError(
+            'Atmospheric pressure must be a number or array of numbers')
+    temperature = np.array(temperature, dtype=float)
+    temperature_kelvin = temperature - pf.constants.absolute_zero_celsius
+    if np.any(np.array(temperature_kelvin) < 0):
+        raise ValueError("Temperature must be above -273.15°C.")
+    if np.any(np.array(relative_humidity) < 0) or np.any(
+            np.array(relative_humidity) > 1):
+        raise ValueError("Relative humidity must be between 0 and 1.")
+    if np.any(np.array(atmospheric_pressure) <= 0):
+        raise ValueError("Atmospheric pressure must be greater than 0 Pa.")
+    P = np.array(atmospheric_pressure, dtype=float)  # Pa
+    relative_humidity = np.array(relative_humidity, dtype=float)
+
+    # Constants according to 6.3.2
+    R = 8.314  # J/(K mol)
+    gamma_a = 1.400
+    gamma_w = 1.330
+    mu_a = 28.97*1e-3  # kg/mol
+    mu_w = 18.02*1e-3  # kg/mol
+    R_a = R / mu_a
+
+    # partial pressure of water vapor in Pa
+    if saturation_vapor_pressure is None:
+        p = constants.saturation_vapor_pressure_magnus(temperature)  # Pa
+    else:
+        p = np.array(saturation_vapor_pressure, dtype=float)  # Pa
+    e = relative_humidity * p  # Pa
+
+    # next to Equation 6.69
+    alpha = mu_a / mu_w
+    # Equation 6.80
+    delta = (1 - (1/gamma_a)) / (1 - (1/gamma_w))
+    nu = (gamma_a - 1) / (gamma_w - 1)
+
+    # Equation 6.70
+    C = (e/P) / (alpha * (1 - e/P))
+
+    # Equation 6.84
+    return np.sqrt(
+        gamma_a * R_a * temperature_kelvin * (
+            1 + (alpha * (1 + delta - nu) - 1) * C),
+    )


Test patch:
diff --git a/tests/test_constants_speed_of_sound.py b/tests/test_constants_speed_of_sound.py
index 2a4d0e7c..1a253d0c 100644
--- a/tests/test_constants_speed_of_sound.py
+++ b/tests/test_constants_speed_of_sound.py
@@ -34,3 +34,109 @@ def test_speed_of_sound_simple_edge_cases():
     expected = 343.2 * np.sqrt((np.array(temperatures) + 273.15) / 293.15)
     result = pf.constants.speed_of_sound_simple(temperatures)
     npt.assert_array_equal(result, expected)
+
+
+def test_speed_of_sound_ideal_gas_typical_values():
+    temperature = 20.0  # Celsius
+    relative_humidity = 0  # 0%
+    result = pf.constants.speed_of_sound_ideal_gas(
+        temperature, relative_humidity)
+    expected = pf.constants.reference_speed_of_sound
+    npt.assert_almost_equal(result, expected, decimal=2)
+
+
+@pytest.mark.parametrize("temperature", [  # Celsius
+    [20, 20.0],
+    np.array([20, 20.0]),
+    (20, 20.0),
+    20,
+])
+@pytest.mark.parametrize("relative_humidity", [
+    [.5, .5],
+    np.array([.5, .5]),
+    (.5, .5),
+    .5,
+])
+@pytest.mark.parametrize("atmospheric_pressure", [  # Pa
+    [101325, 101325],
+    np.array([101325, 101325]),
+    (101325, 101325),
+    101325,
+])
+def test_speed_of_sound_ideal_gas_array_like(
+        temperature, relative_humidity, atmospheric_pressure):
+    result = pf.constants.speed_of_sound_ideal_gas(
+        temperature, relative_humidity, atmospheric_pressure)
+    npt.assert_almost_equal(result, 343.8278399)
+
+
+def test_speed_of_sound_ideal_gas_values_p_water():
+    temperature = 20.0  # Celsius
+    relative_humidity = 0.5  # 50%
+    atmospheric_pressure = pf.constants.reference_atmospheric_pressure  # Pa
+    p_water = pf.constants.saturation_vapor_pressure_magnus(temperature)
+    result = pf.constants.speed_of_sound_ideal_gas(
+        temperature, relative_humidity, atmospheric_pressure, p_water)
+    expected = pf.constants.speed_of_sound_ideal_gas(
+        temperature, relative_humidity, atmospheric_pressure)
+    npt.assert_allclose(result, expected)
+
+
+def test_speed_of_sound_ideal_gas_edge_temperature():
+    temperature = -20.0  # Celsius
+    relative_humidity = 0.0  # 0%
+    atmospheric_pressure = pf.constants.reference_atmospheric_pressure  # Pa
+    result = pf.constants.speed_of_sound_ideal_gas(
+        temperature, relative_humidity, atmospheric_pressure)
+    # calculated using the simplified ideal gas formula 'speed_of_sound_simple'
+    expected = 318.927  # m/s
+    npt.assert_allclose(result, expected, rtol=1e-4)
+
+
+def test_speed_of_sound_ideal_gas_humidity():
+    temperature = 20.0  # Celsius
+    relative_humidity = 1.0  # 100%
+    result = pf.constants.speed_of_sound_ideal_gas(
+        temperature, relative_humidity)
+    # calculated using the simplified ideal gas formula 'speed_of_sound_simple'
+    expected = pf.constants.reference_speed_of_sound
+    npt.assert_allclose(result, expected, rtol=1e-2)
+    assert result > expected
+
+
+def test_speed_of_sound_ideal_gas_edge_pressure():
+    temperature = 20.0  # Celsius
+    relative_humidity = 0.0  # 0%
+    atmospheric_pressure = 90000  # Pa
+    result = pf.constants.speed_of_sound_ideal_gas(
+        temperature, relative_humidity, atmospheric_pressure)
+    expected = pf.constants.reference_speed_of_sound
+    assert result < expected
+    npt.assert_allclose(result, expected, rtol=1e-2)
+
+
+def test_speed_of_sound_ideal_gas_invalid_temperature():
+    temperature = -300.0  # Celsius, out of valid range
+    relative_humidity = 0.0  # 0%
+    atmospheric_pressure = pf.constants.reference_atmospheric_pressure  # Pa
+    with pytest.raises(ValueError, match="Temperature must be"):
+        pf.constants.speed_of_sound_ideal_gas(
+            temperature, relative_humidity, atmospheric_pressure)
+
+
+def test_speed_of_sound_ideal_gas_invalid_humidity():
+    temperature = 20.0  # Celsius
+    relative_humidity = 1.5  # 150%, out of valid range
+    atmospheric_pressure = pf.constants.reference_atmospheric_pressure  # Pa
+    with pytest.raises(ValueError, match="Relative humidity must be"):
+        pf.constants.speed_of_sound_ideal_gas(
+            temperature, relative_humidity, atmospheric_pressure)
+
+
+def test_speed_of_sound_ideal_gas_invalid_pressure():
+    temperature = 20.0  # Celsius
+    relative_humidity = 0.5  # 50%
+    atmospheric_pressure = -50000  # Pa, out of valid range
+    with pytest.raises(ValueError, match="Atmospheric pressure must be"):
+        pf.constants.speed_of_sound_ideal_gas(
+            temperature, relative_humidity, atmospheric_pressure)
diff --git a/tests/test_constants_utils.py b/tests/test_constants_utils.py
new file mode 100644
index 00000000..f5c64992
--- /dev/null
+++ b/tests/test_constants_utils.py
@@ -0,0 +1,46 @@
+import pytest
+import numpy as np
+import pyfar as pf
+import numpy.testing as npt
+
+
+def test_saturation_vapor_pressure_scalar():
+    temperature = 25
+    expected = 3161.736
+    result = pf.constants.saturation_vapor_pressure_magnus(temperature)
+    npt.assert_allclose(result, expected, atol=0.001)
+
+
+def test_saturation_vapor_pressure_array():
+    temperature = np.array([0, 10, 20, 30])
+    expected = np.array([610.94, 1226.0206, 2333.4406, 4236.6503])
+    result = pf.constants.saturation_vapor_pressure_magnus(temperature)
+    npt.assert_allclose(result, expected, atol=0.001)
+
+
+def test_saturation_vapor_pressure_list():
+    temperature = [0, 10, 20, 30]
+    expected = [610.94, 1226.0206, 2333.4406, 4236.6503]
+    result = pf.constants.saturation_vapor_pressure_magnus(temperature)
+    npt.assert_allclose(result, expected, atol=0.001)
+
+
+def test_saturation_vapor_pressure_out_of_range_low():
+    with pytest.raises(
+            ValueError,
+            match="Temperature must be in the range of -45°C and 60°C."):
+        pf.constants.saturation_vapor_pressure_magnus(-50)
+
+
+def test_saturation_vapor_pressure_out_of_range_high():
+    with pytest.raises(
+            ValueError,
+            match="Temperature must be in the range of -45°C and 60°C."):
+        pf.constants.saturation_vapor_pressure_magnus(70)
+
+
+def test_saturation_vapor_pressure_invalid_type():
+    with pytest.raises(
+            TypeError,
+            match="temperature must be a number or array of numbers"):
+        pf.constants.saturation_vapor_pressure_magnus("invalid")