apsis-rs/crates/rf/apsis-atmosphere/src/p839.rs

95 lines
3.8 KiB
Rust

//! The rain height model of ITU-R Recommendation P.839.
//!
//! The height of the rain layer is taken as the height of the zero degree
//! isotherm above sea level. The recommendation tabulates that height on a
//! global grid, and the value at a given location is found by bilinear
//! interpolation between the four surrounding grid points.
//!
//! The rain height is the height at which the falling snow melts, below
//! which the precipitation is liquid. Above it, the attenuation is much
//! smaller, so it is the effective thickness of the attenuating layer and an
//! input to the path calculation of Recommendation P.618.
//!
//! Source: ITU-R Recommendation P.839-4, *Rain height model for prediction
//! methods*.
use crate::p839_table::ZERO_ISOTHERM_HEIGHT_METRES;
/// The grid resolution of the tabulated isotherm height, in degrees.
const STEP_DEGREES: f64 = 1.5;
/// The number of rows in the table, running from north to south.
const ROWS: usize = 121;
/// The number of columns in the table, running from east to east.
const COLUMNS: usize = 241;
/// The offset added to the isotherm height to obtain the rain height, in
/// kilometers.
///
/// The recommendation states that the rain height is the isotherm height
/// plus 0.36 km, which represents the depth of the melting layer below the
/// isotherm.
pub const RAIN_HEIGHT_OFFSET_KM: f64 = 0.36;
/// Returns the height of the zero degree isotherm at a location, in
/// kilometers above mean sea level.
///
/// The latitude is in degrees, positive north, and the longitude in degrees,
/// positive east. The value is interpolated bilinearly from the tabulated
/// grid, and longitudes are reduced into `[0, 360)` first so that any
/// convention is accepted.
///
/// Returns `None` when either coordinate is not finite, or when the latitude
/// is outside the range of the grid.
///
/// Source: ITU-R Recommendation P.839-4, equation 1 and the tabulated
/// isotherm height.
pub fn zero_isotherm_height(latitude: f64, longitude: f64) -> Option<f64> {
if !latitude.is_finite() || !longitude.is_finite() {
return None;
}
if !(-90.0..=90.0).contains(&latitude) {
return None;
}
// The rows run from north to south, so the row index grows as the
// latitude falls.
let row_position = (90.0 - latitude) / STEP_DEGREES;
let row_lower = (row_position.floor() as usize).min(ROWS - 1);
let row_upper = (row_lower + 1).min(ROWS - 1);
let row_fraction = row_position - row_lower as f64;
// The columns run from 0 to 360 degrees east, wrapping at the seam.
let longitude = longitude.rem_euclid(360.0);
let column_position = longitude / STEP_DEGREES;
let column_lower = (column_position.floor() as usize).min(COLUMNS - 1);
let column_upper = (column_lower + 1) % COLUMNS;
let column_fraction = column_position - column_lower as f64;
let sample =
|row: usize, column: usize| f64::from(ZERO_ISOTHERM_HEIGHT_METRES[row][column]) / 1000.0;
// Bilinear interpolation.
let top = sample(row_lower, column_lower)
+ (sample(row_lower, column_upper) - sample(row_lower, column_lower)) * column_fraction;
let bottom = sample(row_upper, column_lower)
+ (sample(row_upper, column_upper) - sample(row_upper, column_lower)) * column_fraction;
Some(top + (bottom - top) * row_fraction)
}
/// Returns the rain height at a location, in kilometers above mean sea level.
///
/// This is the zero degree isotherm height plus
/// [`RAIN_HEIGHT_OFFSET_KM`], as the recommendation defines it.
///
/// Returns `None` when the coordinates are not usable.
pub fn rain_height(latitude: f64, longitude: f64) -> Option<f64> {
Some(zero_isotherm_height(latitude, longitude)? + RAIN_HEIGHT_OFFSET_KM)
}
/// Returns the grid resolution of the tabulated isotherm height, in degrees.
pub const fn grid_step_degrees() -> f64 {
STEP_DEGREES
}