rf: add the ITU-R S.465 reference antenna pattern

crates/apsis-rs/Cargo.toml

@@ -169,3 +169,7 @@ required-features = ["rf"]
 [[example]]
 name = "dump_rain_path"
 required-features = ["rf"]
+
+[[example]]
+name = "dump_s465"
+required-features = ["rf"]

crates/apsis-rs/examples/dump_s465.rs

@@ -0,0 +1,36 @@
+//! Writes the ITU-R S.465 reference pattern computed by this workspace, for
+//! cross-checking against the recommendation.
+//!
+//! Usage:
+//!
+//! ```text
+//! cargo run -p apsis-rs --features rf --example dump_s465 \
+//!     > validation/data/expected/rust_s465.csv
+//! ```
+
+use apsis_antenna::{main_lobe_break_angle, reference_pattern_gain, reference_peak_gain};
+
+/// The antenna sizes sampled, as diameter to wavelength ratios.
+const RATIOS: &[f64] = &[50.0, 100.0, 200.0];
+
+/// The off-axis angles sampled, in degrees.
+const ANGLES_DEGREES: &[f64] = &[0.1, 0.5, 1.0, 5.0, 20.0, 36.0, 48.0, 90.0];
+
+fn main() {
+    println!("ratio,angle_deg,peak_dbi,break_deg,gain_dbi");
+    for &ratio in RATIOS {
+        let Some(peak) = reference_peak_gain(ratio) else {
+            continue;
+        };
+        let Some(break_angle) = main_lobe_break_angle(ratio) else {
+            continue;
+        };
+        for &angle_degrees in ANGLES_DEGREES {
+            let angle = angle_degrees.to_radians();
+            let Some(gain) = reference_pattern_gain(ratio, angle) else {
+                continue;
+            };
+            println!("{ratio},{angle_degrees},{peak:.9},{break_angle:.9},{gain:.9}");
+        }
+    }
+}

crates/rf/apsis-antenna/src/lib.rs

@@ -50,6 +50,11 @@
 
 mod parabolic;
 mod pattern;
+mod reference;
 
 pub use parabolic::{ParabolicAperture, beam_width, floor_angle};
 pub use pattern::{ElementPattern, Pattern, mask_gain, off_axis_attenuation};
+pub use reference::{
+    PEAK_GAIN_OFFSET_DB, diameter_over_wavelength_from_peak_gain, main_lobe_break_angle,
+    reference_beam_width, reference_pattern_decibel, reference_pattern_gain, reference_peak_gain,
+};

crates/rf/apsis-antenna/src/pattern.rs

@@ -7,11 +7,10 @@ use core::f64::consts::PI;
 /// Returns the gain of an antenna at an angle from boresight, given a peak
 /// gain and a half power beam width.
 ///
-/// This is the squared angle model used for interference analysis, in the
-/// form adopted by the Radio Regulations and by ITU-R Recommendation S.465
-/// for the reference pattern of a simple Earth station antenna. The gain
-/// falls as the square of the off-axis angle until it reaches the peak gain
-/// minus a fixed floor, and stays there beyond. In decibels,
+/// This is the squared angle model that several interference assessments
+/// use, expressed in terms of the half power beam width so that any antenna
+/// can be described by the two numbers a link budget has to hand. In
+/// decibels,
 ///
 /// ```text
 /// G(theta) = G_peak - min(12 (theta / theta_0)^2, G_floor)
@@ -23,6 +22,14 @@ use core::f64::consts::PI;
 /// illuminated aperture, which is a good approximation down to the first
 /// sidelobe.
 ///
+/// The coefficient twelve is chosen so that the attenuation reaches three
+/// decibels at `theta_0 / 2`, which is the definition of the half power beam
+/// width. A pattern taken from a recommendation should be used in place of
+/// this one where the recommendation applies: see
+/// [`reference_pattern_gain`](crate::reference_pattern_gain) for the ITU-R
+/// S.465 pattern, which differs from this model in the main lobe and also in
+/// the shape of its sidelobe region.
+///
 /// Returns `None` when the peak gain has no decibel form, or when the beam
 /// width or floor is not positive and finite.
 pub fn off_axis_attenuation(

crates/rf/apsis-antenna/src/reference.rs

@@ -0,0 +1,178 @@
+//! The ITU-R S.465 reference antenna pattern.
+//!
+//! Recommendation S.465 defines a reference pattern for an Earth station
+//! antenna, used to assess interference between satellite networks. The
+//! pattern has three regions: a main lobe that falls as the square of the
+//! off-axis angle, a sidelobe region that falls as the logarithm of the
+//! angle, and a floor.
+//!
+//! The pattern is selected only by the ratio of the antenna diameter to the
+//! wavelength. The recommendation gives one pattern for antennas with
+//! `D / lambda` of at least 100 and another for smaller antennas, whose main
+//! lobe is broader and whose sidelobe region therefore begins later.
+//!
+//! The main lobe expression is what distinguishes this pattern from the
+//! generic reference pattern of the Radio Regulations Appendix 8, which
+//! models the sidelobe region only and takes the main lobe width as an
+//! input. The two share the shape of the sidelobe region.
+//!
+//! Source: ITU-R Recommendation S.465, *Reference earth-station radiation
+//! pattern for use in coordination and interference assessment in the
+//! frequency range from 2 to about 30 GHz*.
+
+use apsis_rf::Decibel;
+use core::f64::consts::PI;
+
+/// The peak gain offset of the reference pattern, in decibels.
+///
+/// The recommendation gives the peak gain as `20 log10(D / lambda) + 7.7`,
+/// which corresponds to an aperture efficiency near 0.6. It is slightly
+/// below what a well adjusted dish reaches, because the reference pattern
+/// represents a typical rather than an optimal antenna.
+pub const PEAK_GAIN_OFFSET_DB: f64 = 7.7;
+
+/// The diameter to wavelength ratio at which the pattern changes form.
+///
+/// A larger antenna has a narrower main lobe, so its sidelobe region begins
+/// earlier and reaches the floor earlier.
+const LARGE_APERTURE_RATIO: f64 = 100.0;
+
+/// The constant and floor angle of the sidelobe region.
+///
+/// The constant is the value of `G + 25 log10(phi_deg)` in the sidelobe
+/// region, and the floor angle is the angle at which the region reaches its
+/// floor.
+fn sidelobe_branch(diameter_over_wavelength: f64) -> (f64, f64) {
+    if diameter_over_wavelength >= LARGE_APERTURE_RATIO {
+        (29.0, 36.0)
+    } else {
+        (32.0, 48.0)
+    }
+}
+
+/// Returns the peak gain of the reference pattern, in decibels relative to an
+/// isotropic radiator.
+///
+/// The expression is `20 log10(D / lambda) + 7.7`.
+///
+/// Returns `None` when the ratio is not positive and finite.
+pub fn reference_peak_gain(diameter_over_wavelength: f64) -> Option<f64> {
+    if !diameter_over_wavelength.is_finite() || diameter_over_wavelength <= 0.0 {
+        return None;
+    }
+    Some(20.0 * diameter_over_wavelength.log10() + PEAK_GAIN_OFFSET_DB)
+}
+
+/// Returns the angle at which the main lobe meets the sidelobe region, in
+/// degrees.
+///
+/// The two expressions intersect at this angle, and there is no closed form,
+/// so the value is found by bisection. The main lobe falls as the square of
+/// the angle while the sidelobe region falls as its logarithm, so the
+/// difference has a single zero in the interval searched, which makes
+/// bisection reliable.
+///
+/// The result lies between zero and the floor angle of the sidelobe region.
+pub fn main_lobe_break_angle(diameter_over_wavelength: f64) -> Option<f64> {
+    let peak = reference_peak_gain(diameter_over_wavelength)?;
+    let (head, floor_angle) = sidelobe_branch(diameter_over_wavelength);
+
+    let difference = |angle_degrees: f64| {
+        let main_lobe = peak - 2.5e-3 * (diameter_over_wavelength * angle_degrees).powi(2);
+        let sidelobe = head - 25.0 * angle_degrees.log10();
+        main_lobe - sidelobe
+    };
+
+    // The main lobe dominates at very small angles and the sidelobe
+    // expression dominates at the floor angle, so the zero lies between.
+    let mut low = 1.0e-9;
+    let mut high = floor_angle;
+    if difference(high) > 0.0 {
+        return Some(high);
+    }
+    for _ in 0..100 {
+        let middle = 0.5 * (low + high);
+        if difference(middle) > 0.0 {
+            low = middle;
+        } else {
+            high = middle;
+        }
+    }
+    Some(0.5 * (low + high))
+}
+
+/// Returns the gain of the S.465 reference pattern at an off-axis angle.
+///
+/// The angle is in radians. The three regions are
+///
+/// ```text
+/// G = G_max - 2.5e-3 (D/lambda phi_deg)^2   for 0 <= phi < phi_break
+/// G = head - 25 log10(phi_deg)              for phi_break <= phi < phi_floor
+/// G = -10                                   for phi_floor <= phi <= 180
+/// ```
+///
+/// where `phi_break` is the angle returned by [`main_lobe_break_angle`], the
+/// head is 29 for a large aperture and 32 for a small one, and the floor
+/// angle is 36 or 48 degrees accordingly.
+///
+/// Returns `None` when the ratio or the angle is not usable.
+pub fn reference_pattern_gain(diameter_over_wavelength: f64, off_axis_angle: f64) -> Option<f64> {
+    if !off_axis_angle.is_finite() || off_axis_angle < 0.0 {
+        return None;
+    }
+    let peak = reference_peak_gain(diameter_over_wavelength)?;
+    let break_angle = main_lobe_break_angle(diameter_over_wavelength)?;
+    let (head, floor_angle) = sidelobe_branch(diameter_over_wavelength);
+    let angle_degrees = off_axis_angle.to_degrees();
+
+    if angle_degrees < break_angle {
+        Some(peak - 2.5e-3 * (diameter_over_wavelength * angle_degrees).powi(2))
+    } else if angle_degrees <= floor_angle {
+        Some(head - 25.0 * angle_degrees.log10())
+    } else {
+        Some(-10.0)
+    }
+}
+
+/// Returns the gain of the S.465 reference pattern as a decibel value.
+///
+/// Returns `None` when the ratio or the angle is not usable.
+pub fn reference_pattern_decibel(
+    diameter_over_wavelength: f64,
+    off_axis_angle: f64,
+) -> Option<Decibel> {
+    let gain = reference_pattern_gain(diameter_over_wavelength, off_axis_angle)?;
+    Decibel::from_db(gain)
+        .is_finite()
+        .then_some(Decibel::from_db(gain))
+}
+
+/// Returns the diameter to wavelength ratio of an antenna with a given peak
+/// gain.
+///
+/// This inverts the peak gain expression, which is useful when a station is
+/// specified by its gain rather than its size.
+///
+/// Returns `None` when the gain is not finite.
+pub fn diameter_over_wavelength_from_peak_gain(peak_gain_dbi: f64) -> Option<f64> {
+    if !peak_gain_dbi.is_finite() {
+        return None;
+    }
+    Some(10f64.powf((peak_gain_dbi - PEAK_GAIN_OFFSET_DB) / 20.0))
+}
+
+/// Returns the half power beam width of the reference pattern main lobe, in
+/// radians.
+///
+/// The main lobe falls to three decibels at
+/// `phi = sqrt(3 / 2.5e-3) / (D / lambda)` degrees from boresight, so the
+/// full width between the two three decibel points is twice that value.
+///
+/// Returns `None` when the ratio is not positive and finite.
+pub fn reference_beam_width(diameter_over_wavelength: f64) -> Option<f64> {
+    if !diameter_over_wavelength.is_finite() || diameter_over_wavelength <= 0.0 {
+        return None;
+    }
+    let half_angle_degrees = (3.0_f64 / 2.5e-3).sqrt() / diameter_over_wavelength;
+    Some(2.0 * half_angle_degrees * PI / 180.0)
+}

crates/rf/apsis-antenna/tests/reference.rs

@@ -0,0 +1,198 @@
+//! Tests for the ITU-R S.465 reference antenna pattern.
+
+use apsis_antenna::{
+    PEAK_GAIN_OFFSET_DB, diameter_over_wavelength_from_peak_gain, main_lobe_break_angle,
+    reference_beam_width, reference_pattern_gain, reference_peak_gain,
+};
+
+const TOLERANCE: f64 = 1e-9;
+
+#[test]
+fn the_peak_gain_follows_the_recommendation() {
+    // The recommendation gives 20 log10(D / lambda) + 7.7.
+    for ratio in [10.0_f64, 50.0, 100.0, 200.0, 1_000.0] {
+        let expected = 20.0 * ratio.log10() + PEAK_GAIN_OFFSET_DB;
+        let actual = reference_peak_gain(ratio).expect("valid ratio");
+        assert!((actual - expected).abs() < TOLERANCE, "ratio {ratio}");
+    }
+    // A ratio of one hundred gives 47.7 dBi.
+    assert!((reference_peak_gain(100.0).expect("valid") - 47.7).abs() < TOLERANCE);
+}
+
+#[test]
+fn the_peak_gain_is_inverted_by_the_diameter_helper() {
+    for ratio in [10.0, 100.0, 500.0] {
+        let gain = reference_peak_gain(ratio).expect("valid");
+        let recovered = diameter_over_wavelength_from_peak_gain(gain).expect("finite gain");
+        assert!((recovered / ratio - 1.0).abs() < 1e-12, "ratio {ratio}");
+    }
+}
+
+#[test]
+fn the_pattern_at_boresight_is_the_peak_gain() {
+    for ratio in [50.0, 100.0, 200.0] {
+        let peak = reference_peak_gain(ratio).expect("valid");
+        let boresight = reference_pattern_gain(ratio, 0.0).expect("valid");
+        assert!((boresight - peak).abs() < TOLERANCE, "ratio {ratio}");
+    }
+}
+
+#[test]
+fn the_main_lobe_falls_as_the_square_of_the_angle() {
+    // Inside the main lobe the attenuation is 2.5e-3 (D/lambda phi)^2.
+    let ratio = 100.0;
+    let peak = reference_peak_gain(ratio).expect("valid");
+    for angle_degrees in [0.05_f64, 0.1, 0.2, 0.3] {
+        let gain = reference_pattern_gain(ratio, angle_degrees.to_radians()).expect("valid");
+        let expected = peak - 2.5e-3 * (ratio * angle_degrees).powi(2);
+        assert!(
+            (gain - expected).abs() < TOLERANCE,
+            "{angle_degrees} degrees: got {gain}, expected {expected}"
+        );
+    }
+}
+
+#[test]
+fn the_pattern_matches_recorded_values() {
+    // These are the values of the recommendation for a ratio of one hundred,
+    // computed independently from its equations.
+    let cases: [(f64, f64); 8] = [
+        (0.1, 47.45),
+        (0.5, 41.45),
+        (1.0, 29.0),
+        (5.0, 11.525_749_892),
+        (20.0, -3.525_749_892),
+        (36.0, -9.907_562_519),
+        (48.0, -10.0),
+        (90.0, -10.0),
+    ];
+    for (angle_degrees, expected) in cases {
+        let gain = reference_pattern_gain(100.0, angle_degrees.to_radians()).expect("valid");
+        assert!(
+            (gain - expected).abs() < 1e-6,
+            "{angle_degrees} degrees: got {gain}, expected {expected}"
+        );
+    }
+}
+
+#[test]
+fn the_pattern_matches_recorded_values_for_a_small_aperture() {
+    // A ratio of fifty uses the small aperture branch, with a head of 32 and
+    // a floor angle of 48 degrees.
+    // At exactly the floor angle the branch switches to the floor, so the
+    // last value of the falling branch is approached from below.
+    let cases: [(f64, f64); 6] = [
+        (0.1, 41.616_900_087),
+        (1.0, 35.429_400_087),
+        (5.0, 14.525_749_892),
+        (36.0, -6.907_562_519),
+        (47.99, -10.028_768_748_310_533),
+        (48.0, -10.0),
+    ];
+    for (angle_degrees, expected) in cases {
+        let gain = reference_pattern_gain(50.0, angle_degrees.to_radians()).expect("valid");
+        assert!(
+            (gain - expected).abs() < 1e-6,
+            "{angle_degrees} degrees: got {gain}, expected {expected}"
+        );
+    }
+}
+
+#[test]
+fn the_two_branches_meet_at_the_break_angle() {
+    // The main lobe and the sidelobe region must agree at the junction, or
+    // the pattern would jump.
+    for ratio in [50.0_f64, 100.0, 200.0, 500.0] {
+        let peak = reference_peak_gain(ratio).expect("valid");
+        let break_angle = main_lobe_break_angle(ratio).expect("valid");
+        let (head, _) = if ratio >= 100.0 {
+            (29.0, 36.0)
+        } else {
+            (32.0, 48.0)
+        };
+
+        let main_lobe = peak - 2.5e-3 * (ratio * break_angle).powi(2);
+        let sidelobe = head - 25.0 * break_angle.log10();
+        assert!(
+            (main_lobe - sidelobe).abs() < 1e-6,
+            "ratio {ratio}: {main_lobe} against {sidelobe}"
+        );
+    }
+}
+
+#[test]
+fn the_break_angle_shrinks_as_the_aperture_grows() {
+    // A larger antenna has a narrower main lobe, so the sidelobe region
+    // begins sooner.
+    let small = main_lobe_break_angle(50.0).expect("valid");
+    let medium = main_lobe_break_angle(100.0).expect("valid");
+    let large = main_lobe_break_angle(400.0).expect("valid");
+    assert!(small > medium, "{small} against {medium}");
+    assert!(medium > large, "{medium} against {large}");
+}
+
+#[test]
+fn the_pattern_floors_at_minus_ten_decibels() {
+    for ratio in [50.0, 100.0, 200.0] {
+        let gain = reference_pattern_gain(ratio, 90.0_f64.to_radians()).expect("valid");
+        assert!((gain + 10.0).abs() < TOLERANCE, "ratio {ratio}");
+        let beyond = reference_pattern_gain(ratio, core::f64::consts::PI).expect("valid");
+        assert!((beyond + 10.0).abs() < TOLERANCE, "ratio {ratio}");
+    }
+}
+
+#[test]
+fn the_pattern_never_rises_above_the_peak() {
+    let ratio = 100.0;
+    let peak = reference_peak_gain(ratio).expect("valid");
+    for tenth in 0..=100 {
+        let angle = tenth as f64 * 0.01;
+        let gain = reference_pattern_gain(ratio, angle).expect("valid");
+        assert!(gain <= peak + TOLERANCE, "angle {angle} gave {gain}");
+    }
+}
+
+#[test]
+fn the_pattern_falls_monotonically_across_the_main_lobe() {
+    let ratio = 100.0;
+    let mut previous = f64::INFINITY;
+    for step in 0..=80 {
+        let angle_degrees = step as f64 * 0.01;
+        let gain = reference_pattern_gain(ratio, angle_degrees.to_radians()).expect("valid");
+        assert!(gain <= previous + TOLERANCE, "angle {angle_degrees}");
+        previous = gain;
+    }
+}
+
+#[test]
+fn the_beam_width_follows_from_the_main_lobe_law() {
+    // Three decibels is reached at sqrt(3 / 2.5e-3) / (D/lambda) degrees
+    // from boresight, so the full width is twice that.
+    let ratio = 100.0;
+    let width = reference_beam_width(ratio).expect("valid");
+    let half_degrees = (3.0_f64 / 2.5e-3).sqrt() / ratio;
+    let expected = 2.0 * half_degrees.to_radians();
+    assert!((width - expected).abs() < 1e-15);
+
+    // The gain at half the beam width from boresight is three decibels down.
+    let gain = reference_pattern_gain(ratio, width / 2.0).expect("valid");
+    let peak = reference_peak_gain(ratio).expect("valid");
+    assert!(
+        (peak - gain - 3.0).abs() < 1e-9,
+        "peak - gain = {}",
+        peak - gain
+    );
+}
+
+#[test]
+fn invalid_arguments_are_rejected() {
+    assert_eq!(reference_peak_gain(0.0), None);
+    assert_eq!(reference_peak_gain(-1.0), None);
+    assert_eq!(reference_peak_gain(f64::NAN), None);
+    assert_eq!(reference_pattern_gain(100.0, -0.1), None);
+    assert_eq!(reference_pattern_gain(100.0, f64::NAN), None);
+    assert_eq!(reference_pattern_gain(0.0, 0.1), None);
+    assert_eq!(reference_beam_width(0.0), None);
+    assert_eq!(diameter_over_wavelength_from_peak_gain(f64::NAN), None);
+    assert_eq!(main_lobe_break_angle(0.0), None);
+}