apsis-rs/crates/rf/apsis-atmosphere/tests/p835.rs

//! Tests for the ITU-R P.835 reference standard atmosphere.

use apsis_atmosphere::{
    geopotential_height, radio_refractive_index, standard_pressure, standard_temperature,
    standard_water_vapour_density,
};

const TOLERANCE: f64 = 1e-9;

#[test]
fn the_geopotential_height_is_below_the_geometric_height() {
    // The conversion compresses the height because gravity weakens with
    // altitude, so the geopotential height is always the smaller value.
    for geometric in [0.0, 1.0, 11.0, 50.0, 100.0] {
        let geopotential = geopotential_height(geometric);
        assert!(geopotential <= geometric, "height {geometric}");
    }
    // At the surface the two coincide.
    assert_eq!(geopotential_height(0.0), 0.0);
}

#[test]
fn the_surface_conditions_match_the_recommendation() {
    // The mean annual global reference atmosphere starts at 288.15 K and
    // 1013.25 hPa.
    assert!((standard_temperature(0.0) - 288.15).abs() < TOLERANCE);
    assert!((standard_pressure(0.0) - 1_013.25).abs() < TOLERANCE);
    assert!((standard_water_vapour_density(0.0, 7.5) - 7.5).abs() < TOLERANCE);
}

#[test]
fn the_temperature_profile_has_the_tabulated_segments() {
    // The troposphere cools at 6.5 K per kilometer of geopotential height,
    // so the rate is expressed against that height, not the geometric one.
    let expected = 288.15 - 6.5 * geopotential_height(1.0);
    assert!((standard_temperature(1.0) - expected).abs() < TOLERANCE);

    // The stratosphere is isothermal at 216.65 K between 11 and 20 km of
    // geopotential height. A geometric height of 15 km has a geopotential
    // height of about 14.97 km, which is inside that band.
    let geopotential = geopotential_height(15.0);
    assert!(
        (11.0..=20.0).contains(&geopotential),
        "height {geopotential}"
    );
    assert!((standard_temperature(15.0) - 216.65).abs() < TOLERANCE);

    // Above 20 km geopotential the temperature rises at 1 K per kilometer.
    // A geometric height of 30 km reaches the middle of that segment.
    let geopotential = geopotential_height(30.0);
    assert!(
        (20.0..32.0).contains(&geopotential),
        "height {geopotential}"
    );
    let expected = 216.65 + (geopotential - 20.0);
    assert!((standard_temperature(30.0) - expected).abs() < TOLERANCE);

    // The temperature is continuous across each segment boundary. The
    // boundaries are in geopotential height, so they are converted to the
    // geometric heights at which they occur before sampling either side.
    for geopotential in [11.0_f64, 20.0, 32.0, 47.0, 51.0, 71.0] {
        let radius = 6_356.766;
        let geometric = radius * geopotential / (radius - geopotential);
        let below = standard_temperature(geometric - 0.001);
        let above = standard_temperature(geometric + 0.001);
        assert!(
            (below - above).abs() < 0.01,
            "boundary {geopotential}: {below} against {above}"
        );
    }
}

#[test]
fn the_temperature_is_continuous_across_every_segment() {
    for boundary in [11.0_f64, 20.0, 32.0, 47.0, 51.0, 71.0, 84.852] {
        // The boundaries are in geopotential height, so a geometric height
        // just either side is used.
        let below = standard_temperature(boundary - 0.01);
        let above = standard_temperature(boundary + 0.01);
        assert!(
            (below - above).abs() < 0.2,
            "boundary {boundary}: {below} against {above}"
        );
    }
}

#[test]
fn the_pressure_falls_with_height() {
    let mut previous = f64::INFINITY;
    for height in [0.0, 1.0, 5.0, 11.0, 20.0, 32.0, 47.0, 51.0, 71.0, 84.852] {
        let pressure = standard_pressure(height);
        assert!(pressure < previous, "height {height} gave {pressure}");
        previous = pressure;
    }
}

#[test]
fn the_pressure_is_continuous_across_the_segments() {
    for boundary in [11.0_f64, 20.0, 32.0, 47.0, 51.0, 71.0, 84.852] {
        let below = standard_pressure(boundary - 0.01);
        let above = standard_pressure(boundary + 0.01);
        let relative = (below - above).abs() / below;
        assert!(
            relative < 0.01,
            "boundary {boundary}: {below} against {above}, relative {relative}"
        );
    }
}

#[test]
fn the_pressure_at_the_tropopause_matches_the_reference() {
    // The tropopause is at 11 km of geopotential height, which is a
    // geometric height of about 11.02 km. The pressure there is 226.32 hPa.
    let geometric = 11.019_068;
    assert!((geopotential_height(geometric) - 11.0).abs() < 1e-5);
    let pressure = standard_pressure(geometric);
    assert!((pressure - 226.32).abs() < 0.01, "pressure {pressure}");
}

#[test]
fn the_water_vapour_profile_is_exponential() {
    // The scale height is two kilometers, so the density halves about every
    // 1.386 km.
    let surface = 7.5;
    let at_two = standard_water_vapour_density(2.0, surface);
    assert!((at_two - surface / core::f64::consts::E).abs() < 1e-9);
    // The profile scales with the surface value.
    let doubled = standard_water_vapour_density(2.0, 15.0);
    assert!((doubled / at_two - 2.0).abs() < 1e-12);
}

#[test]
fn the_radio_refractive_index_at_the_surface_is_about_one() {
    // A refractivity of about three hundred N units gives an index about
    // 1.0003.
    let index = radio_refractive_index(1_013.25, 10.0, 288.15);
    assert!(index > 1.0 && index < 1.001, "index {index}");
    assert!((index - 1.000_315).abs() < 1e-5, "index {index}");
}

#[test]
fn the_radio_refractive_index_grows_with_humidity() {
    let dry = radio_refractive_index(1_013.25, 0.0, 288.15);
    let moist = radio_refractive_index(1_013.25, 20.0, 288.15);
    assert!(moist > dry);
}

#[test]
fn a_zero_temperature_is_handled_without_dividing_by_zero() {
    // The expression divides by the temperature, so a zero would give an
    // infinite value. The function returns the index of vacuum instead.
    assert_eq!(radio_refractive_index(1_013.25, 10.0, 0.0), 1.0);
}