Repository: pyfar/pyfar
Pull request #771 (2025-02-17 12:10:36)

Rebase: Speed of sound cramer
### Which issue(s) are closed by this pull request?

Contributes to #428 
replaces #738

### Changes proposed in this pull request:

- [x] same changes as in #738
- [x] rebase, to merge into ``speed_of_sound`` instead of ``air_attenuation``, to avoid more chaos
- [x] merge #732 first



Patch:
diff --git a/pyfar/constants/__init__.py b/pyfar/constants/__init__.py
index cca81ad0..b75e7f2b 100644
--- a/pyfar/constants/__init__.py
+++ b/pyfar/constants/__init__.py
@@ -8,6 +8,7 @@ from .constants import (
 
 from .speed_of_sound import (
     speed_of_sound_simple,
+    speed_of_sound_cramer,
     speed_of_sound_ideal_gas,
 )
 
@@ -18,8 +19,9 @@ from .constants import (
 __all__ = [
     'saturation_vapor_pressure_magnus',
     'speed_of_sound_simple',
-    'air_attenuation',
+    'speed_of_sound_cramer',
     'speed_of_sound_ideal_gas',
+    'air_attenuation',
 ]
 
 
diff --git a/pyfar/constants/speed_of_sound.py b/pyfar/constants/speed_of_sound.py
index 49cca5a4..1762b22a 100644
--- a/pyfar/constants/speed_of_sound.py
+++ b/pyfar/constants/speed_of_sound.py
@@ -144,3 +144,117 @@ def speed_of_sound_ideal_gas(
         gamma_a * R_a * temperature_kelvin * (
             1 + (alpha * (1 + delta - nu) - 1) * C),
     )
+
+
+def speed_of_sound_cramer(
+        temperature, relative_humidity, co2_ppm=425.19,
+        atmospheric_pressure=None,
+        ):
+    r"""
+    Calculate the speed of sound in air based on temperature, atmospheric
+    pressure, humidity and CO\ :sub:`2` concentration.
+
+    This implements Cramers method described in [#]_.
+
+    Parameters
+    ----------
+    temperature : float, array_like
+        Temperature in degree Celsius from 0°C to 30°C.
+    relative_humidity : float, array_like
+        Relative humidity in the range of 0 to 1.
+    co2_ppm : float, array_like, optional
+        CO\ :sub:`2` concentration in parts per million. The default is
+        425.19 ppm, based on [#]_. Value must be below
+        10 000 ppm (1%).
+    atmospheric_pressure : float, array_like, optional
+        Atmospheric pressure in Pascal, by default
+        :py:attr:`reference_atmospheric_pressure`.
+        Value must be between 75 000 Pa to 102000 Pa.
+
+    Returns
+    -------
+    speed_of_sound : float, array_like
+        Speed of sound in air in m/s.
+
+    References
+    ----------
+    .. [#] O. Cramer, “The variation of the specific heat ratio and the
+           speed of sound in air with temperature, pressure, humidity, and
+           CO\ :sub:`2` concentration,” The Journal of the Acoustical Society
+           of America, vol. 93, no. 5, pp. 2510-2516, May 1993,
+           doi: 10.1121/1.405827.
+    .. [#] Lan, X., Tans, P. and K.W. Thoning: Trends in globally-averaged
+           CO2 determined from NOAA Global Monitoring Laboratory measurements.
+           Version Friday, 14-Mar-2025 11:33:44 MDT
+           https://doi.org/10.15138/9N0H-ZH07
+
+    """
+    if atmospheric_pressure is None:
+        atmospheric_pressure = pf.constants.reference_atmospheric_pressure
+    # convert to array
+    temperature = np.array(temperature)
+    relative_humidity = np.array(relative_humidity)
+    atmospheric_pressure = np.array(atmospheric_pressure)
+    co2_ppm = np.array(co2_ppm)
+
+    # check inputs:
+    if np.any(temperature < 0) or np.any(temperature > 30):
+        raise ValueError("Temperature must be between 0°C and 30°C.")
+    if np.any(relative_humidity < 0) or np.any(relative_humidity > 1):
+        raise ValueError("Relative humidity must be between 0 and 1.")
+    if np.any(atmospheric_pressure < 75e3) or np.any(
+            atmospheric_pressure > 102e3):
+        raise ValueError(
+            "Atmospheric pressure must be between 75 000 Pa to 102 000 Pa.")
+    if np.any(co2_ppm < 0) or np.any(co2_ppm > .01e6):
+        raise ValueError(
+            "CO2 concentration (ppm) must be between 0 ppm to 10 000 ppm.")
+
+    # mole fraction of CO2 in air
+    x_c = co2_ppm * 1e-6
+    # pressure(atmospheric pressure) in Pa
+    p = atmospheric_pressure
+    # relative humidity as a fraction
+    h = relative_humidity
+    # temperature in Celsius
+    t = temperature
+    # thermodynamic temperature
+    T = t - pf.constants.absolute_zero_celsius
+
+    # enhancement factor f (Equation A2)
+    f = 1.00062+3.14e-8 * p + 5.6e-7*temperature**2
+    # saturation vapor pressure in air (Equation A3)
+    p_sv = np.exp(1.2811805e-5*T**2-1.9509874e-2*T+34.04926034-6.3536311e3/T)
+    # water vapor mole fraction (Equation A1)
+    x_w = h * f * p_sv / p
+
+    # check water mole fraction after Cramer
+    if np.any(x_w < 0) or np.any(x_w > 0.06):
+        raise ValueError("Water mole fraction must be between 0 and 0.06.")
+
+    # Coefficients according to Cramer (Table. III)
+    a0 = 331.5024
+    a1 = 0.603055
+    a2 = -0.000528
+    a3 = 51.471935
+    a4 = 0.1495874
+    a5 = -0.000782
+    a6 = -1.82e-7
+    a7 = 3.73e-8
+    a8 = -2.93e-10
+    a9 = -85.20931
+    a10 = -0.228525
+    a11 = 5.91e-5
+    a12 = -2.835149
+    a13 = -2.15e-13
+    a14 = 29.179762
+    a15 = 0.000486
+
+    # approximation for c according to Cramer (Eq. 15)
+    c1 = a0+a1*t+a2*t**2
+    c2 = (a3+a4*t+a5*t**2)*x_w
+    c3 = (a6+a7*t+a8*t**2)*p
+    c4 = (a9+a10*t+a11*t**2)*x_c
+    c5 = a12*x_w**2 + a13*p**2 + a14*x_c**2 + a15*x_w*p*x_c
+
+    return c1 + c2 + c3 + c4 + c5


Test patch:
diff --git a/tests/test_constants_speed_of_sound_cramer.py b/tests/test_constants_speed_of_sound_cramer.py
new file mode 100644
index 00000000..bdf9e772
--- /dev/null
+++ b/tests/test_constants_speed_of_sound_cramer.py
@@ -0,0 +1,147 @@
+import pytest
+import numpy as np
+import pyfar as pf
+import numpy.testing as npt
+import re
+
+
+@pytest.mark.parametrize(("temperature", "relative_humidity", "expected"), [
+    (0, 0.0047393, 331.5079),
+    (0, 0.30174, 331.627),
+    (0, 0.60032, 331.7063),
+    (0, 0.99526, 331.9048),
+    (10, 0.0063191, 337.4603),
+    (10, 0.30332, 337.6984),
+    (10, 0.6019, 337.8571),
+    (10, 0.9921, 338.1746),
+    (20, 0.0078989, 343.3333),
+    (20, 0.30174, 343.7302),
+    (20, 0.6019, 344.1667),
+    (20, 0.98736, 344.6825),
+    (30, 0.0047393, 349.1667),
+    (30, 0.30332, 349.881),
+    (30, 0.6019, 350.5556),
+    (30, 0.98894, 351.5079),
+    ])
+def test_speed_of_sound_cramer_figure_1_rel_hum(
+        temperature, relative_humidity, expected):
+    """Test speed_of_sound_cramer with figure 1 vs relative humidity."""
+
+    speed = pf.constants.speed_of_sound_cramer(
+        temperature, relative_humidity, 314)
+    npt.assert_almost_equal(speed, expected, 1)
+
+
+@pytest.mark.parametrize((
+        "temperature", "atmospheric_pressure", "expected"), [
+            (0, 75265.49, 331.4732),
+            (0, 85035.4, 331.3839),
+            (0, 95070.8, 331.3839),
+            (10, 75212.39, 337.4107),
+            (10, 85088.5, 337.4554),
+            (10, 95123.89, 337.5),
+            (20, 75212.39, 343.3036),
+            (20, 85088.5, 343.3482),
+            (20, 95017.7, 343.3929),
+            (30, 75053.1, 349.1000),
+            (30, 85088.5, 349.1518),
+            (30, 94858.41, 349.2411),
+    ])
+def test_speed_of_sound_cramer_figure_1_pressure(
+        temperature, atmospheric_pressure, expected):
+    """Test speed_of_sound_cramer with figure 1 vs atmospheric pressure."""
+    relative_humidity = 0
+
+    speed = pf.constants.speed_of_sound_cramer(
+        temperature, relative_humidity, 314, atmospheric_pressure)
+    npt.assert_almost_equal(speed, expected, 1)
+
+
+@pytest.mark.parametrize((
+        "temperature", "co2_percent", "expected"), [
+            (0, 0.006079, 331.4818),
+            (0, 0.29939, 331.2538),
+            (0, 0.6003, 331.0258),
+            (0, 0.9924, 330.6839),
+            (10, 0.0075988, 337.4848),
+            (10, 0.29939, 337.2188),
+            (10, 0.60182, 336.9529),
+            (10, 0.99088, 336.7249),
+            (20, 0.0091185, 343.4498),
+            (20, 0.30091, 343.1459),
+            (20, 0.6003, 342.9559),
+            (20, 0.98784, 342.576),
+            (30, 0.0091185, 349.035),
+            (30, 0.30091, 348.845),
+            (30, 0.60182, 348.617),
+            (30, 0.9848, 348.3511),
+    ])
+def test_speed_of_sound_cramer_figure_1_co2(
+        temperature, co2_percent, expected):
+    """Test speed_of_sound_cramer with figure 1 vs co2 concentration."""
+    co2_ppm = co2_percent / 100 * 1e6  # in parts per million
+    relative_humidity = 0
+
+    speed = pf.constants.speed_of_sound_cramer(
+        temperature, relative_humidity, co2_ppm=co2_ppm,
+        )
+    npt.assert_almost_equal(speed, expected, 1)
+
+
+@pytest.mark.parametrize((
+        "temperature", "relative_humidity", "atmospheric_pressure",
+        "co2_ppm", "expected"), [
+            (0, 0, 95070.8, 314, 331.3839),
+            (0, 0, 95070.8, [314], 331.3839),
+            (0, 0, [95070.8], [314], 331.3839),
+            (0, [0], 95070.8, [314], 331.3839),
+            ([0], 0, 95070.8, 314, 331.3839),
+            ([0], 0, (95070.8, ), 314, 331.3839),
+            ([0, 0], 0, 95070.8, 314, 331.3839),
+            (np.array([[0, 0], [0, 0]]), 0, 95070.8, 314, 331.3839),
+        ])
+def test_speed_of_sound_cramer_array_like(
+        temperature, relative_humidity, atmospheric_pressure,
+        co2_ppm, expected):
+    speed = pf.constants.speed_of_sound_cramer(
+        temperature, relative_humidity, co2_ppm, atmospheric_pressure,
+        )
+    npt.assert_almost_equal(speed, expected, 1)
+
+
+def test_speed_of_sound_cramer_invalid_temperature():
+    match = 'Temperature must be between 0°C and 30°C.'
+    with pytest.raises(ValueError, match=re.escape(match)):
+        pf.constants.speed_of_sound_cramer(
+            -10, .5)
+    with pytest.raises(ValueError, match=re.escape(match)):
+        pf.constants.speed_of_sound_cramer(
+            31, .5)
+
+
+def test_speed_of_sound_cramer_invalid_atmospheric_pressure():
+    match = 'Atmospheric pressure must be between 75 000 Pa to 102 000 Pa.'
+    with pytest.raises(ValueError, match=re.escape(match)):
+        pf.constants.speed_of_sound_cramer(0, .5, 341, 74999)
+    with pytest.raises(ValueError, match=re.escape(match)):
+        pf.constants.speed_of_sound_cramer(1, .5, 341, 102001)
+
+
+def test_speed_of_sound_cramer_invalid_relative_humidity():
+    match = 'Relative humidity must be between 0 and 1.'
+    with pytest.raises(ValueError, match=re.escape(match)):
+        pf.constants.speed_of_sound_cramer(
+            0, -0.1)
+    with pytest.raises(ValueError, match=re.escape(match)):
+        pf.constants.speed_of_sound_cramer(
+            1, 1.1)
+
+
+def test_speed_of_sound_cramer_invalid_co2_ppm():
+    match = 'CO2 concentration (ppm) must be between 0 ppm to 10 000 ppm.'
+    with pytest.raises(ValueError, match=re.escape(match)):
+        pf.constants.speed_of_sound_cramer(
+            0, 0, -1)
+    with pytest.raises(ValueError, match=re.escape(match)):
+        pf.constants.speed_of_sound_cramer(
+            1, .5, 10001)