From 8f9cf071a0dac7968faa6c36a021e06410b8549b Mon Sep 17 00:00:00 2001 From: Timothy Wong Date: Sun, 28 Jun 2026 03:09:03 +0800 Subject: [PATCH] Port to pure Rust --- itm-rs/.gitignore | 1 + itm-rs/Cargo.lock | 14 + itm-rs/Cargo.toml | 3 + itm-rs/README.md | 443 +++++++++++++++++++++ itm-rs/itm-cli/Cargo.toml | 8 + itm-rs/itm-cli/src/main.rs | 320 +++++++++++++++ itm-rs/itm/Cargo.toml | 5 + itm-rs/itm/src/area.rs | 174 ++++++++ itm-rs/itm/src/complex.rs | 115 ++++++ itm-rs/itm/src/compute_delta_h.rs | 55 +++ itm-rs/itm/src/constants.rs | 39 ++ itm-rs/itm/src/diffraction_loss.rs | 55 +++ itm-rs/itm/src/errors.rs | 26 ++ itm-rs/itm/src/find_horizons.rs | 39 ++ itm-rs/itm/src/free_space_loss.rs | 4 + itm-rs/itm/src/fresnel_integral.rs | 9 + itm-rs/itm/src/h0_function.rs | 18 + itm-rs/itm/src/iccdf.rs | 16 + itm-rs/itm/src/initialize_area.rs | 44 ++ itm-rs/itm/src/initialize_p2p.rs | 31 ++ itm-rs/itm/src/knife_edge_diffraction.rs | 22 + itm-rs/itm/src/lib.rs | 78 ++++ itm-rs/itm/src/line_of_sight_loss.rs | 46 +++ itm-rs/itm/src/linear_least_squares_fit.rs | 36 ++ itm-rs/itm/src/longley_rice.rs | 208 ++++++++++ itm-rs/itm/src/p2p.rs | 169 ++++++++ itm-rs/itm/src/quick_pfl.rs | 67 ++++ itm-rs/itm/src/sigma_h_function.rs | 4 + itm-rs/itm/src/smooth_earth_diffraction.rs | 71 ++++ itm-rs/itm/src/terrain_roughness.rs | 4 + itm-rs/itm/src/tests.rs | 93 +++++ itm-rs/itm/src/troposcatter_loss.rs | 105 +++++ itm-rs/itm/src/types.rs | 22 + itm-rs/itm/src/validate_inputs.rs | 89 +++++ itm-rs/itm/src/variability.rs | 175 ++++++++ itm-rs/itm/src/warnings.rs | 15 + 36 files changed, 2623 insertions(+) create mode 100644 itm-rs/.gitignore create mode 100644 itm-rs/Cargo.lock create mode 100644 itm-rs/Cargo.toml create mode 100644 itm-rs/README.md create mode 100644 itm-rs/itm-cli/Cargo.toml create mode 100644 itm-rs/itm-cli/src/main.rs create mode 100644 itm-rs/itm/Cargo.toml create mode 100644 itm-rs/itm/src/area.rs create mode 100644 itm-rs/itm/src/complex.rs create mode 100644 itm-rs/itm/src/compute_delta_h.rs create mode 100644 itm-rs/itm/src/constants.rs create mode 100644 itm-rs/itm/src/diffraction_loss.rs create mode 100644 itm-rs/itm/src/errors.rs create mode 100644 itm-rs/itm/src/find_horizons.rs create mode 100644 itm-rs/itm/src/free_space_loss.rs create mode 100644 itm-rs/itm/src/fresnel_integral.rs create mode 100644 itm-rs/itm/src/h0_function.rs create mode 100644 itm-rs/itm/src/iccdf.rs create mode 100644 itm-rs/itm/src/initialize_area.rs create mode 100644 itm-rs/itm/src/initialize_p2p.rs create mode 100644 itm-rs/itm/src/knife_edge_diffraction.rs create mode 100644 itm-rs/itm/src/lib.rs create mode 100644 itm-rs/itm/src/line_of_sight_loss.rs create mode 100644 itm-rs/itm/src/linear_least_squares_fit.rs create mode 100644 itm-rs/itm/src/longley_rice.rs create mode 100644 itm-rs/itm/src/p2p.rs create mode 100644 itm-rs/itm/src/quick_pfl.rs create mode 100644 itm-rs/itm/src/sigma_h_function.rs create mode 100644 itm-rs/itm/src/smooth_earth_diffraction.rs create mode 100644 itm-rs/itm/src/terrain_roughness.rs create mode 100644 itm-rs/itm/src/tests.rs create mode 100644 itm-rs/itm/src/troposcatter_loss.rs create mode 100644 itm-rs/itm/src/types.rs create mode 100644 itm-rs/itm/src/validate_inputs.rs create mode 100644 itm-rs/itm/src/variability.rs create mode 100644 itm-rs/itm/src/warnings.rs diff --git a/itm-rs/.gitignore b/itm-rs/.gitignore new file mode 100644 index 0000000..ea8c4bf --- /dev/null +++ b/itm-rs/.gitignore @@ -0,0 +1 @@ +/target diff --git a/itm-rs/Cargo.lock b/itm-rs/Cargo.lock new file mode 100644 index 0000000..8fea798 --- /dev/null +++ b/itm-rs/Cargo.lock @@ -0,0 +1,14 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 4 + +[[package]] +name = "itm" +version = "1.4.0" + +[[package]] +name = "itm-cli" +version = "1.4.0" +dependencies = [ + "itm", +] diff --git a/itm-rs/Cargo.toml b/itm-rs/Cargo.toml new file mode 100644 index 0000000..b63fe7e --- /dev/null +++ b/itm-rs/Cargo.toml @@ -0,0 +1,3 @@ +[workspace] +members = ["itm", "itm-cli"] +resolver = "2" diff --git a/itm-rs/README.md b/itm-rs/README.md new file mode 100644 index 0000000..1dd7f6f --- /dev/null +++ b/itm-rs/README.md @@ -0,0 +1,443 @@ +# ITM — Rust Implementation + +Pure Rust port of the [ITS Irregular Terrain Model (ITM)](../README.md) v1.4.0. +ITM predicts terrestrial radiowave propagation for frequencies between 20 MHz and 20 GHz using +the Longley-Rice algorithm, considering free-space loss, diffraction, and troposcatter. + +## Contents + +| Crate | Description | +|-------|-------------| +| [`itm`](itm/) | Library — call from your Rust code | +| [`itm-cli`](itm-cli/) | CLI binary — use without writing code | + +--- + +## Quick Start + +### Build + +```sh +cargo build --release +``` + +Both crates are in a workspace, so this produces: +- `target/release/libITM.rlib` (or `.dll`/`.so` on the relevant platform) +- `target/release/itm-cli` + +### Run the CLI + +**Point-to-point (terrain profile):** + +```sh +itm-cli p2p \ + --h-tx 15 --h-rx 3 \ + --pfl terrain.txt \ + --climate 5 --n0 301 --freq 3500 --pol 1 \ + --epsilon 15 --sigma 0.005 \ + --mdvar 12 --time 50 --location 50 --situation 50 +``` + +**Area prediction (no terrain file):** + +```sh +itm-cli area \ + --h-tx 10 --h-rx 1 \ + --tx-siting 0 --rx-siting 0 \ + --distance 16 --delta-h 0 \ + --climate 5 --n0 301 --freq 230 --pol 0 \ + --epsilon 15 --sigma 0.008 \ + --mdvar 0 --time 87 --location 50 --situation 50 +``` + +**Example output:** + +``` +Basic Transmission Loss: 142.37 dB +Return code: 0 +--- Intermediate Values --- + Path distance: 16.000 km + Free-space loss: 107.75 dB + Reference atten: 34.62 dB + Delta H: 0.00 m + Surface refractivity: 301.0 N-units + Eff. TX height: 10.000 m + Eff. RX height: 1.000 m + TX horizon dist: 1234.5 m + RX horizon dist: 456.7 m + Propagation mode: Diffraction +``` + +--- + +## Library Usage + +Add to `Cargo.toml`: + +```toml +[dependencies] +itm = { path = "path/to/itm-rs/itm" } +``` + +### Point-to-Point + +```rust +use itm; + +// Build the terrain profile array in PFL format: +// pfl[0] = N (number of points minus 1) +// pfl[1] = dR (spacing in meters) +// pfl[2..=N+2] = elevation above MSL in meters +let mut pfl = vec![0.0_f64; 102]; +pfl[0] = 100.0; // 101 points +pfl[1] = 300.0; // 300 m spacing → 30 km path + +let (rtn, a__db, warnings) = itm::itm_p2p_tls( + 15.0, // h_tx__meter + 3.0, // h_rx__meter + &pfl, + 5, // climate: Continental Temperate + 301.0, // n_0: surface refractivity (N-units) + 3500.0, // f__mhz + 1, // pol: vertical + 15.0, // epsilon + 0.005, // sigma (S/m) + 12, // mdvar + 50.0, // time % + 50.0, // location % + 50.0, // situation % +); + +assert_eq!(rtn, itm::errors::SUCCESS); +println!("Loss: {a__db:.2} dB warnings: 0x{warnings:04X}"); +``` + +To also retrieve intermediate values, call `itm_p2p_tls_ex`: + +```rust +let (rtn, a__db, warnings, iv) = itm::itm_p2p_tls_ex(/* same args */); +println!("Mode: {}", match iv.mode { 1 => "LOS", 2 => "Diffraction", _ => "Troposcatter" }); +println!("Δh: {:.1} m", iv.delta_h__meter); +``` + +### Area Prediction + +```rust +let (rtn, a__db, warnings) = itm::itm_area_tls( + 10.0, // h_tx__meter + 1.0, // h_rx__meter + 0, // tx_siting_criteria: Random + 0, // rx_siting_criteria: Random + 16.0, // d__km + 0.0, // delta_h__meter + 5, // climate + 301.0, // n_0 + 230.0, // f__mhz + 0, // pol: horizontal + 15.0, // epsilon + 0.008, // sigma + 0, // mdvar + 87.0, // time % + 50.0, // location % + 50.0, // situation % +); +``` + +### Confidence/Reliability Variant + +Prefer `_cr` functions when you want to specify confidence and reliability +instead of time/location/situation: + +```rust +let (rtn, a__db, warnings) = itm::itm_p2p_cr( + h_tx, h_rx, &pfl, climate, n_0, f__mhz, pol, epsilon, sigma, + mdvar, + 50.0, // confidence % + 90.0, // reliability % +); +``` + +Extended (`_ex`) versions exist for all four variants: + +| Function | Variability | Extras | +|---|---|---| +| `itm_p2p_tls` | time / location / situation | — | +| `itm_p2p_tls_ex` | time / location / situation | `IntermediateValues` | +| `itm_p2p_cr` | confidence / reliability | — | +| `itm_p2p_cr_ex` | confidence / reliability | `IntermediateValues` | +| `itm_area_tls` | time / location / situation | — | +| `itm_area_tls_ex` | time / location / situation | `IntermediateValues` | +| `itm_area_cr` | confidence / reliability | — | +| `itm_area_cr_ex` | confidence / reliability | `IntermediateValues` | + +--- + +## Input Reference + +### Common Inputs + +| Parameter | Type | Units | Limits | Description | +|-----------|------|-------|--------|-------------| +| `h_tx__meter` | `f64` | m | 0.5 – 3000 | TX structural height | +| `h_rx__meter` | `f64` | m | 0.5 – 3000 | RX structural height | +| `climate` | `i32` | — | 1 – 7 | Radio climate (see table below) | +| `n_0` | `f64` | N-Units | 250 – 400 | Min monthly mean surface refractivity | +| `f__mhz` | `f64` | MHz | 20 – 20 000 | Frequency | +| `pol` | `i32` | — | 0, 1 | Polarization: 0 = horizontal, 1 = vertical | +| `epsilon` | `f64` | — | > 1 | Relative permittivity of ground | +| `sigma` | `f64` | S/m | > 0 | Ground conductivity | +| `mdvar` | `i32` | — | 0 – 33 | Mode of variability (see table below) | + +**Radio climates:** + +| Value | Enum variant | Description | +|-------|--------------|-------------| +| 1 | `Climate::Equatorial` | Equatorial | +| 2 | `Climate::ContinentalSubtropical` | Continental Subtropical | +| 3 | `Climate::MaritimeSubtropical` | Maritime Subtropical | +| 4 | `Climate::Desert` | Desert | +| 5 | `Climate::ContinentalTemperate` | Continental Temperate | +| 6 | `Climate::MaritimeTemperateOverLand` | Maritime Temperate (Over Land) | +| 7 | `Climate::MaritimeTemperateOverSea` | Maritime Temperate (Over Sea) | + +**Modes of variability (`mdvar`):** + +| Value | Description | +|-------|-------------| +| 0 | Single Message | +| 1 | Accidental | +| 2 | Mobile | +| 3 | Broadcast | +| +10 | Eliminate location variability | +| +20 | Eliminate direct situation variability | + +### Point-to-Point Inputs + +| Parameter | Type | Description | +|-----------|------|-------------| +| `pfl` | `&[f64]` | Terrain profile: `[N, dR, h₀, h₁, …, hₙ]` | + +`pfl[0]` = number of terrain points minus 1 (N). +`pfl[1]` = horizontal spacing in meters (dR). +`pfl[2]` through `pfl[N+2]` = elevation above MSL in meters (N+1 values total). + +### Area Mode Inputs + +| Parameter | Type | Units | Description | +|-----------|------|-------|-------------| +| `d__km` | `f64` | km | Path distance (> 0) | +| `delta_h__meter` | `f64` | m | Terrain irregularity parameter (≥ 0) | +| `tx_siting_criteria` | `i32` | — | 0 = Random, 1 = Careful, 2 = Very Careful | +| `rx_siting_criteria` | `i32` | — | 0 = Random, 1 = Careful, 2 = Very Careful | + +### Variability Inputs + +| Parameter | Type | Limits | Description | +|-----------|------|--------|-------------| +| `time` | `f64` | 0 < x < 100 | Time percentage | +| `location` | `f64` | 0 < x < 100 | Location percentage | +| `situation` | `f64` | 0 < x < 100 | Situation percentage | +| `confidence` | `f64` | 0 < x < 100 | Confidence percentage (`_cr` functions) | +| `reliability` | `f64` | 0 < x < 100 | Reliability percentage (`_cr` functions) | + +--- + +## Outputs + +All functions return a tuple `(rtn: i32, a__db: f64, warnings: i32)`. +`_ex` variants append `IntermediateValues` as a fourth element. + +| Field | Units | Description | +|-------|-------|-------------| +| `rtn` | — | Return/error code (0 = success) | +| `a__db` | dB | Basic transmission loss | +| `warnings` | — | Warning bitmask (see below) | + +### `IntermediateValues` + +```rust +pub struct IntermediateValues { + pub theta_hzn: [f64; 2], // horizon angles (radians) + pub d_hzn__meter: [f64; 2], // horizon distances (m) + pub h_e__meter: [f64; 2], // effective heights (m) + pub n_s: f64, // surface refractivity (N-Units) + pub delta_h__meter: f64, // terrain irregularity (m) + pub a_ref__db: f64, // reference attenuation (dB) + pub a_fs__db: f64, // free-space loss (dB) + pub d__km: f64, // path distance (km) + pub mode: i32, // 1=LOS, 2=Diffraction, 3=Troposcatter +} +``` + +--- + +## Error Codes and Warnings + +### Error Codes (`itm::errors::*`) + +| Constant | Value | Meaning | +|----------|-------|---------| +| `SUCCESS` | 0 | No error | +| `SUCCESS_WITH_WARNINGS` | 1 | Completed with warnings | +| `ERROR__TX_TERMINAL_HEIGHT` | 1000 | TX height out of range | +| `ERROR__RX_TERMINAL_HEIGHT` | 1001 | RX height out of range | +| `ERROR__INVALID_RADIO_CLIMATE` | 1002 | Climate not 1–7 | +| `ERROR__INVALID_TIME` | 1003 | Time % not in (0, 100) | +| `ERROR__INVALID_LOCATION` | 1004 | Location % not in (0, 100) | +| `ERROR__INVALID_SITUATION` | 1005 | Situation % not in (0, 100) | +| `ERROR__INVALID_CONFIDENCE` | 1006 | Confidence % not in (0, 100) | +| `ERROR__INVALID_RELIABILITY` | 1007 | Reliability % not in (0, 100) | +| `ERROR__REFRACTIVITY` | 1008 | N₀ not in [250, 400] | +| `ERROR__FREQUENCY` | 1009 | Frequency not in [20, 20 000] MHz | +| `ERROR__POLARIZATION` | 1010 | Polarization not 0 or 1 | +| `ERROR__EPSILON` | 1011 | Relative permittivity ≤ 1 | +| `ERROR__SIGMA` | 1012 | Conductivity ≤ 0 | +| `ERROR__GROUND_IMPEDANCE` | 1013 | Ground impedance invalid | +| `ERROR__MDVAR` | 1014 | mdvar not in [0, 33] | +| `ERROR__EFFECTIVE_EARTH` | 1016 | Effective earth radius error | +| `ERROR__PATH_DISTANCE` | 1017 | Path distance ≤ 0 | +| `ERROR__DELTA_H` | 1018 | delta_h < 0 | +| `ERROR__TX_SITING_CRITERIA` | 1019 | TX siting not 0/1/2 | +| `ERROR__RX_SITING_CRITERIA` | 1020 | RX siting not 0/1/2 | +| `ERROR__SURFACE_REFRACTIVITY_SMALL` | 1021 | Computed surface refractivity too small | +| `ERROR__SURFACE_REFRACTIVITY_LARGE` | 1022 | Computed surface refractivity too large | + +### Warning Flags (`itm::warnings::*`) + +Warnings are a bitmask — multiple flags may be set simultaneously. +Check with `warnings & itm::warnings::WARN_ != 0`. + +| Bit | Meaning | +|-----|---------| +| `0x0001` | TX terminal height out of recommended range | +| `0x0002` | RX terminal height out of recommended range | +| `0x0004` | Frequency out of recommended range | +| `0x0008` | Path distance > 1000 km | +| `0x0010` | Path distance > 2000 km | +| `0x0020` | Path distance below minimum for antenna heights | +| `0x0040` | Path distance < 1 km | +| `0x0080` | TX horizon angle > 200 mrad | +| `0x0100` | RX horizon angle > 200 mrad | +| `0x0200` | TX horizon distance < 10% of smooth-earth value | +| `0x0400` | RX horizon distance < 10% of smooth-earth value | +| `0x0800` | TX horizon distance > 3× smooth-earth value | +| `0x1000` | RX horizon distance > 3× smooth-earth value | +| `0x2000` | Extreme variabilities | +| `0x4000` | Surface refractivity out of recommended range | + +--- + +## Terrain File Format (CLI) + +The `--pfl` file is plain text, whitespace-separated: + +``` +N dR h0 h1 h2 ... hN +``` + +- `N` — number of terrain points **minus 1** (integer) +- `dR` — horizontal spacing in meters (float) +- `h0` … `hN` — terrain elevations above MSL in meters (**N+1 values**) + +The file must contain exactly `N + 2` numeric values. + +**Example** — 5-point profile, 1000 m spacing: + +``` +4 1000 100.0 120.5 98.3 110.0 105.7 +``` + +--- + +## Migrating from C++ + +The Rust API is a direct translation of the C++ DLL interface. +Parameter names and semantics are identical; only the calling convention changes. + +### C++ → Rust mapping + +| C++ | Rust | +|-----|------| +| `ITM_P2P_TLS(h_tx, h_rx, pfl, ...)` | `itm::itm_p2p_tls(h_tx, h_rx, &pfl, ...)` | +| `ITM_P2P_TLS_Ex(h_tx, h_rx, pfl, ..., &iv)` | `itm::itm_p2p_tls_ex(h_tx, h_rx, &pfl, ...)` → returns `(rtn, a__db, warn, iv)` | +| `ITM_AREA_TLS(...)` | `itm::itm_area_tls(...)` | +| `ITM_AREA_TLS_Ex(..., &iv)` | `itm::itm_area_tls_ex(...)` → returns `(rtn, a__db, warn, iv)` | + +Key differences: +- **No output parameters** — Rust functions return tuples instead of mutating `&A__db`, `&warnings`, `&iv`. +- **`pfl` is a slice** — pass `&pfl` rather than a raw pointer; the library reads `pfl[0]` as N internally. +- **`IntermediateValues` is returned by value** — no need to pre-allocate and pass a pointer. +- **Error codes are `i32` constants** in `itm::errors`, not a C enum. + +### Before (C++) + +```cpp +double A__db; +int warnings; +IntermediateValues iv; + +int rtn = ITM_P2P_TLS_Ex( + h_tx, h_rx, pfl, climate, N_0, f__mhz, pol, + epsilon, sigma, mdvar, time, location, situation, + &A__db, &warnings, &iv); +``` + +### After (Rust) + +```rust +let (rtn, a__db, warnings, iv) = itm::itm_p2p_tls_ex( + h_tx, h_rx, &pfl, climate, n_0, f__mhz, pol, + epsilon, sigma, mdvar, time, location, situation, +); +``` + +--- + +## Migrating from .NET / NuGet + +The C# wrapper called the native DLL with `ref` output parameters. +The Rust library returns a plain tuple — no `ref`, no `out`. + +### Before (C#) + +```csharp +double A_db; +int warnings; +IntermediateValues iv; + +int rtn = ITM.ITM_P2P_TLS_Ex( + h_tx, h_rx, pfl, climate, N_0, f_mhz, pol, + epsilon, sigma, mdvar, time, location, situation, + out A_db, out warnings, out iv); +``` + +### After (Rust, called from C# via FFI or replaced entirely) + +If replacing with pure Rust: + +```rust +let (rtn, a__db, warnings, iv) = itm::itm_p2p_tls_ex(/* same numeric args */); +``` + +If you need to expose the Rust library back to C# as a DLL, wrap the functions with +`#[no_mangle] pub extern "C"` signatures that match the original DLL exports. + +--- + +## Running Tests + +```sh +cargo test +``` + +The test suite (`itm/src/tests.rs`) covers five P2P and five area scenarios across +the 230 MHz – 8.9 GHz range, verifying both successful completion and loss values +within ±1 dB of reference outputs. + +--- + +## References + +- G.A. Hufford, A.G. Longley, W.A. Kissick, [A Guide to the Use of the ITS Irregular Terrain Model in the Area Prediction Mode](https://www.its.bldrdoc.gov/publications/details.aspx?pub=2091), NTIA TR-82-100, 1982. +- G.A. Hufford, [The ITS Irregular Terrain Model, version 1.2.2 Algorithm](https://www.its.bldrdoc.gov/media/50676/itm_alg.pdf). +- A.G. Longley and P.L. Rice, [Prediction of Tropospheric Radio Transmission Loss Over Irregular Terrain](https://www.its.bldrdoc.gov/publications/details.aspx?pub=2784), NTIA ERL 79-ITS 67, 1968. diff --git a/itm-rs/itm-cli/Cargo.toml b/itm-rs/itm-cli/Cargo.toml new file mode 100644 index 0000000..81de32a --- /dev/null +++ b/itm-rs/itm-cli/Cargo.toml @@ -0,0 +1,8 @@ +[package] +name = "itm-cli" +version = "1.4.0" +edition = "2021" +description = "CLI for the ITS Irregular Terrain Model (ITM)" + +[dependencies] +itm = { path = "../itm" } diff --git a/itm-rs/itm-cli/src/main.rs b/itm-rs/itm-cli/src/main.rs new file mode 100644 index 0000000..b6c15b6 --- /dev/null +++ b/itm-rs/itm-cli/src/main.rs @@ -0,0 +1,320 @@ +#![allow(non_snake_case)] +use itm::IntermediateValues; +use std::env; +use std::fs; +use std::process; + +const HELP: &str = r#"ITS Irregular Terrain Model (ITM) — Rust CLI + +USAGE: + itm-cli [OPTIONS] + +SUBCOMMANDS: + p2p Point-to-point prediction (requires terrain profile file) + area Area prediction (no terrain file needed) + +P2P OPTIONS (all required): + --h-tx TX structural height, meters (0.5–3000) + --h-rx RX structural height, meters (0.5–3000) + --pfl Terrain profile file (plain text: N dR h0 h1 ... hN) + --climate <1-7> Radio climate (1=Equatorial … 7=MarTemSea) + --n0 Surface refractivity (250–400) + --freq Frequency in MHz (20–20000) + --pol <0|1> Polarization: 0=horizontal, 1=vertical + --epsilon Relative permittivity (>1) + --sigma Conductivity, S/m (>0) + --mdvar <0-33> Mode of variability + --time Time percentage (0–100 exclusive) + --location Location percentage (0–100 exclusive) + --situation Situation percentage (0–100 exclusive) + +AREA OPTIONS (all required): + --h-tx TX structural height, meters + --h-rx RX structural height, meters + --tx-siting <0-2> TX siting criteria (0=random, 1=careful, 2=very careful) + --rx-siting <0-2> RX siting criteria + --distance Path distance in km (>0) + --delta-h Terrain irregularity parameter, meters (>=0) + --climate <1-7> Radio climate + --n0 Surface refractivity + --freq Frequency in MHz + --pol <0|1> Polarization + --epsilon Relative permittivity + --sigma Conductivity + --mdvar <0-33> Mode of variability + --time Time percentage + --location Location percentage + --situation Situation percentage + +TERRAIN FILE FORMAT: + First value: number of terrain points minus 1 (N) + Second value: spacing in meters (dR) + Remaining N+1 values: terrain elevations in meters above MSL + +EXAMPLES: + itm-cli p2p --h-tx 15 --h-rx 3 --pfl terrain.txt --climate 5 --n0 301 + --freq 3500 --pol 1 --epsilon 15 --sigma 0.005 --mdvar 12 + --time 50 --location 50 --situation 50 + + itm-cli area --h-tx 10 --h-rx 1 --tx-siting 0 --rx-siting 0 + --distance 16 --delta-h 0 --climate 5 --n0 301 --freq 230 + --pol 0 --epsilon 15 --sigma 0.008 --mdvar 0 + --time 87 --location 50 --situation 50 +"#; + +fn usage_err(msg: &str) -> ! { + eprintln!("Error: {msg}"); + eprintln!("Run `itm-cli --help` for usage."); + process::exit(1); +} + +fn parse_f64(val: &str, name: &str) -> f64 { + val.parse::() + .unwrap_or_else(|_| usage_err(&format!("--{name} must be a number, got '{val}'"))) +} + +fn parse_i32(val: &str, name: &str) -> i32 { + val.parse::() + .unwrap_or_else(|_| usage_err(&format!("--{name} must be an integer, got '{val}'"))) +} + +/// Parse a key=value pair from the argument iterator. +fn next_val<'a>(args: &mut impl Iterator, flag: &str) -> &'a str { + args.next() + .unwrap_or_else(|| usage_err(&format!("--{flag} requires a value"))) +} + +/// Load a PFL terrain file. +/// Format: `N dR h0 h1 ... hN` (N+2 header values: pfl[0]=N, pfl[1]=dR, rest=heights) +fn load_pfl(path: &str) -> Vec { + let text = fs::read_to_string(path) + .unwrap_or_else(|e| usage_err(&format!("Cannot read terrain file '{path}': {e}"))); + let values: Vec = text + .split_whitespace() + .map(|s| { + s.parse::() + .unwrap_or_else(|_| usage_err(&format!("Non-numeric value '{s}' in terrain file"))) + }) + .collect(); + if values.len() < 3 { + usage_err("Terrain file must have at least 3 values (N dR h0)"); + } + let n = values[0] as usize; + if values.len() != n + 2 { + usage_err(&format!( + "Terrain file declares N={n} points but has {} values (expected {})", + values.len(), + n + 2 + )); + } + values +} + +fn print_result(rtn: i32, a__db: f64, warnings: i32, inter: Option<&IntermediateValues>) { + println!("Basic Transmission Loss: {:.2} dB", a__db); + println!("Return code: {rtn}"); + if warnings != 0 { + println!("Warnings: 0x{warnings:04X}"); + print_warnings(warnings); + } + if let Some(iv) = inter { + println!("\n--- Intermediate Values ---"); + println!(" Path distance: {:.3} km", iv.d__km); + println!(" Free-space loss: {:.2} dB", iv.a_fs__db); + println!(" Reference atten: {:.2} dB", iv.a_ref__db); + println!(" Delta H: {:.2} m", iv.delta_h__meter); + println!(" Surface refractivity: {:.1} N-units", iv.n_s); + println!( + " Eff. TX height: {:.3} m", iv.h_e__meter[0] + ); + println!( + " Eff. RX height: {:.3} m", iv.h_e__meter[1] + ); + println!( + " TX horizon dist: {:.1} m", iv.d_hzn__meter[0] + ); + println!( + " RX horizon dist: {:.1} m", iv.d_hzn__meter[1] + ); + let mode_str = match iv.mode { + 1 => "Line of Sight", + 2 => "Diffraction", + 3 => "Troposcatter", + _ => "Unknown", + }; + println!(" Propagation mode: {mode_str}"); + } +} + +fn print_warnings(w: i32) { + let flags = [ + (0x0001, "TX terminal height out of recommended range"), + (0x0002, "RX terminal height out of recommended range"), + (0x0004, "Frequency out of recommended range"), + (0x0008, "Path distance > 1000 km"), + (0x0010, "Path distance > 2000 km"), + (0x0020, "Path distance below minimum for antenna heights"), + (0x0040, "Path distance < 1 km"), + (0x0080, "TX horizon angle > 200 mrad"), + (0x0100, "RX horizon angle > 200 mrad"), + (0x0200, "TX horizon distance < 10% of smooth-earth value"), + (0x0400, "RX horizon distance < 10% of smooth-earth value"), + (0x0800, "TX horizon distance > 3x smooth-earth value"), + (0x1000, "RX horizon distance > 3x smooth-earth value"), + (0x2000, "Extreme variabilities"), + (0x4000, "Surface refractivity out of recommended range"), + ]; + for (mask, desc) in &flags { + if w & mask != 0 { + println!(" [WARN] {desc}"); + } + } +} + +fn run_p2p(args: &[String]) { + let mut h_tx = None::; + let mut h_rx = None::; + let mut pfl_path = None::; + let mut climate = None::; + let mut n_0 = None::; + let mut f_mhz = None::; + let mut pol = None::; + let mut epsilon = None::; + let mut sigma = None::; + let mut mdvar = None::; + let mut time = None::; + let mut location = None::; + let mut situation = None::; + + let mut iter = args.iter().map(String::as_str); + while let Some(flag) = iter.next() { + let val = next_val(&mut iter, flag); + match flag { + "--h-tx" => h_tx = Some(parse_f64(val, "h-tx")), + "--h-rx" => h_rx = Some(parse_f64(val, "h-rx")), + "--pfl" => pfl_path = Some(val.to_string()), + "--climate" => climate = Some(parse_i32(val, "climate")), + "--n0" => n_0 = Some(parse_f64(val, "n0")), + "--freq" => f_mhz = Some(parse_f64(val, "freq")), + "--pol" => pol = Some(parse_i32(val, "pol")), + "--epsilon" => epsilon = Some(parse_f64(val, "epsilon")), + "--sigma" => sigma = Some(parse_f64(val, "sigma")), + "--mdvar" => mdvar = Some(parse_i32(val, "mdvar")), + "--time" => time = Some(parse_f64(val, "time")), + "--location" => location = Some(parse_f64(val, "location")), + "--situation" => situation = Some(parse_f64(val, "situation")), + other => usage_err(&format!("Unknown flag for p2p: {other}")), + } + } + + macro_rules! req { + ($opt:expr, $name:literal) => { + $opt.unwrap_or_else(|| usage_err(concat!("--", $name, " is required for p2p mode"))) + }; + } + + let pfl = load_pfl(&req!(pfl_path, "pfl")); + + let (rtn, a__db, warnings, inter) = itm::itm_p2p_tls_ex( + req!(h_tx, "h-tx"), + req!(h_rx, "h-rx"), + &pfl, + req!(climate, "climate"), + req!(n_0, "n0"), + req!(f_mhz, "freq"), + req!(pol, "pol"), + req!(epsilon, "epsilon"), + req!(sigma, "sigma"), + req!(mdvar, "mdvar"), + req!(time, "time"), + req!(location, "location"), + req!(situation, "situation"), + ); + + print_result(rtn, a__db, warnings, Some(&inter)); +} + +fn run_area(args: &[String]) { + let mut h_tx = None::; + let mut h_rx = None::; + let mut tx_siting = None::; + let mut rx_siting = None::; + let mut distance = None::; + let mut delta_h = None::; + let mut climate = None::; + let mut n_0 = None::; + let mut f_mhz = None::; + let mut pol = None::; + let mut epsilon = None::; + let mut sigma = None::; + let mut mdvar = None::; + let mut time = None::; + let mut location = None::; + let mut situation = None::; + + let mut iter = args.iter().map(String::as_str); + while let Some(flag) = iter.next() { + let val = next_val(&mut iter, flag); + match flag { + "--h-tx" => h_tx = Some(parse_f64(val, "h-tx")), + "--h-rx" => h_rx = Some(parse_f64(val, "h-rx")), + "--tx-siting" => tx_siting = Some(parse_i32(val, "tx-siting")), + "--rx-siting" => rx_siting = Some(parse_i32(val, "rx-siting")), + "--distance" => distance = Some(parse_f64(val, "distance")), + "--delta-h" => delta_h = Some(parse_f64(val, "delta-h")), + "--climate" => climate = Some(parse_i32(val, "climate")), + "--n0" => n_0 = Some(parse_f64(val, "n0")), + "--freq" => f_mhz = Some(parse_f64(val, "freq")), + "--pol" => pol = Some(parse_i32(val, "pol")), + "--epsilon" => epsilon = Some(parse_f64(val, "epsilon")), + "--sigma" => sigma = Some(parse_f64(val, "sigma")), + "--mdvar" => mdvar = Some(parse_i32(val, "mdvar")), + "--time" => time = Some(parse_f64(val, "time")), + "--location" => location = Some(parse_f64(val, "location")), + "--situation" => situation = Some(parse_f64(val, "situation")), + other => usage_err(&format!("Unknown flag for area: {other}")), + } + } + + macro_rules! req { + ($opt:expr, $name:literal) => { + $opt.unwrap_or_else(|| usage_err(concat!("--", $name, " is required for area mode"))) + }; + } + + let (rtn, a__db, warnings, inter) = itm::itm_area_tls_ex( + req!(h_tx, "h-tx"), + req!(h_rx, "h-rx"), + req!(tx_siting, "tx-siting"), + req!(rx_siting, "rx-siting"), + req!(distance, "distance"), + req!(delta_h, "delta-h"), + req!(climate, "climate"), + req!(n_0, "n0"), + req!(f_mhz, "freq"), + req!(pol, "pol"), + req!(epsilon, "epsilon"), + req!(sigma, "sigma"), + req!(mdvar, "mdvar"), + req!(time, "time"), + req!(location, "location"), + req!(situation, "situation"), + ); + + print_result(rtn, a__db, warnings, Some(&inter)); +} + +fn main() { + let args: Vec = env::args().collect(); + + if args.len() < 2 || args[1] == "--help" || args[1] == "-h" { + print!("{HELP}"); + return; + } + + match args[1].as_str() { + "p2p" => run_p2p(&args[2..]), + "area" => run_area(&args[2..]), + other => usage_err(&format!("Unknown subcommand '{other}'. Use 'p2p' or 'area'.")), + } +} diff --git a/itm-rs/itm/Cargo.toml b/itm-rs/itm/Cargo.toml new file mode 100644 index 0000000..8af718b --- /dev/null +++ b/itm-rs/itm/Cargo.toml @@ -0,0 +1,5 @@ +[package] +name = "itm" +version = "1.4.0" +edition = "2021" +description = "ITS Irregular Terrain Model (ITM) — pure Rust implementation" diff --git a/itm-rs/itm/src/area.rs b/itm-rs/itm/src/area.rs new file mode 100644 index 0000000..2580b0a --- /dev/null +++ b/itm-rs/itm/src/area.rs @@ -0,0 +1,174 @@ +use crate::constants::*; +use crate::errors::*; +use crate::free_space_loss::free_space_loss; +use crate::initialize_area::initialize_area; +use crate::initialize_p2p::initialize_p2p; +use crate::longley_rice::longley_rice; +use crate::p2p::translate_cr_error; +use crate::types::IntermediateValues; +use crate::validate_inputs::validate_inputs; +use crate::variability::variability; + +/// Area mode with time/location/situation variability. +pub fn itm_area_tls( + h_tx__meter: f64, + h_rx__meter: f64, + tx_site_criteria: i32, + rx_site_criteria: i32, + d__km: f64, + delta_h__meter: f64, + climate: i32, + n_0: f64, + f__mhz: f64, + pol: i32, + epsilon: f64, + sigma: f64, + mdvar: i32, + time: f64, + location: f64, + situation: f64, +) -> (i32, f64, i32) { + let (rtn, a, w, _) = itm_area_tls_ex( + h_tx__meter, h_rx__meter, tx_site_criteria, rx_site_criteria, d__km, delta_h__meter, + climate, n_0, f__mhz, pol, epsilon, sigma, mdvar, time, location, situation, + ); + (rtn, a, w) +} + +/// Area mode with time/location/situation variability (extended). +pub fn itm_area_tls_ex( + h_tx__meter: f64, + h_rx__meter: f64, + tx_site_criteria: i32, + rx_site_criteria: i32, + d__km: f64, + delta_h__meter: f64, + climate: i32, + n_0: f64, + f__mhz: f64, + pol: i32, + epsilon: f64, + sigma: f64, + mdvar: i32, + time: f64, + location: f64, + situation: f64, +) -> (i32, f64, i32, IntermediateValues) { + let mut inter = IntermediateValues::default(); + let mut warnings = NO_WARNINGS; + + let rtn = validate_inputs( + h_tx__meter, h_rx__meter, climate, time, location, situation, + n_0, f__mhz, pol, epsilon, sigma, mdvar, &mut warnings, + ); + if rtn != SUCCESS { + return (rtn, 0.0, warnings, inter); + } + + if d__km <= 0.0 { + return (ERROR__PATH_DISTANCE, 0.0, warnings, inter); + } + if delta_h__meter < 0.0 { + return (ERROR__DELTA_H, 0.0, warnings, inter); + } + if tx_site_criteria != SITING_CRITERIA__RANDOM + && tx_site_criteria != SITING_CRITERIA__CAREFUL + && tx_site_criteria != SITING_CRITERIA__VERY_CAREFUL + { + return (ERROR__TX_SITING_CRITERIA, 0.0, warnings, inter); + } + if rx_site_criteria != SITING_CRITERIA__RANDOM + && rx_site_criteria != SITING_CRITERIA__CAREFUL + && rx_site_criteria != SITING_CRITERIA__VERY_CAREFUL + { + return (ERROR__RX_SITING_CRITERIA, 0.0, warnings, inter); + } + + inter.d__km = d__km; + let h__meter = [h_tx__meter, h_rx__meter]; + let site_criteria = [tx_site_criteria, rx_site_criteria]; + + let (z_g, gamma_e, n_s) = initialize_p2p(f__mhz, 0.0, n_0, pol, epsilon, sigma); + let (h_e__meter, d_hzn__meter, theta_hzn) = + initialize_area(site_criteria, gamma_e, delta_h__meter, h__meter); + + let d__meter = d__km * 1000.0; + let (lr_rtn, a_ref__db, propmode) = longley_rice( + theta_hzn, f__mhz, z_g, d_hzn__meter, h_e__meter, gamma_e, n_s, + delta_h__meter, h__meter, d__meter, MODE__AREA, &mut warnings, + ); + if lr_rtn != SUCCESS { + return (lr_rtn, 0.0, warnings, inter); + } + + let a_fs__db = free_space_loss(d__meter, f__mhz); + let a__db = a_fs__db + + variability( + time, location, situation, h_e__meter, delta_h__meter, + f__mhz, d__meter, a_ref__db, climate, mdvar, &mut warnings, + ); + + inter.a_ref__db = a_ref__db; + inter.a_fs__db = a_fs__db; + inter.delta_h__meter = delta_h__meter; + inter.d_hzn__meter = d_hzn__meter; + inter.h_e__meter = h_e__meter; + inter.n_s = n_s; + inter.theta_hzn = theta_hzn; + inter.mode = propmode; + + let rtn = if warnings != NO_WARNINGS { SUCCESS_WITH_WARNINGS } else { SUCCESS }; + (rtn, a__db, warnings, inter) +} + +/// Area mode with confidence/reliability variability. +pub fn itm_area_cr( + h_tx__meter: f64, + h_rx__meter: f64, + tx_site_criteria: i32, + rx_site_criteria: i32, + d__km: f64, + delta_h__meter: f64, + climate: i32, + n_0: f64, + f__mhz: f64, + pol: i32, + epsilon: f64, + sigma: f64, + mdvar: i32, + confidence: f64, + reliability: f64, +) -> (i32, f64, i32) { + let (mut rtn, a, w, _) = itm_area_tls_ex( + h_tx__meter, h_rx__meter, tx_site_criteria, rx_site_criteria, d__km, delta_h__meter, + climate, n_0, f__mhz, pol, epsilon, sigma, mdvar, reliability, 50.0, confidence, + ); + rtn = translate_cr_error(rtn); + (rtn, a, w) +} + +/// Area mode with confidence/reliability variability (extended). +pub fn itm_area_cr_ex( + h_tx__meter: f64, + h_rx__meter: f64, + tx_site_criteria: i32, + rx_site_criteria: i32, + d__km: f64, + delta_h__meter: f64, + climate: i32, + n_0: f64, + f__mhz: f64, + pol: i32, + epsilon: f64, + sigma: f64, + mdvar: i32, + confidence: f64, + reliability: f64, +) -> (i32, f64, i32, IntermediateValues) { + let (mut rtn, a, w, inter) = itm_area_tls_ex( + h_tx__meter, h_rx__meter, tx_site_criteria, rx_site_criteria, d__km, delta_h__meter, + climate, n_0, f__mhz, pol, epsilon, sigma, mdvar, reliability, 50.0, confidence, + ); + rtn = translate_cr_error(rtn); + (rtn, a, w, inter) +} diff --git a/itm-rs/itm/src/complex.rs b/itm-rs/itm/src/complex.rs new file mode 100644 index 0000000..60d077e --- /dev/null +++ b/itm-rs/itm/src/complex.rs @@ -0,0 +1,115 @@ +/// Minimal complex number type for internal ITM calculations. +#[derive(Clone, Copy, Debug, Default)] +pub struct Complex { + pub re: f64, + pub im: f64, +} + +impl Complex { + #[inline] + pub fn new(re: f64, im: f64) -> Self { + Complex { re, im } + } + + /// Magnitude (absolute value) + #[inline] + pub fn abs(self) -> f64 { + (self.re * self.re + self.im * self.im).sqrt() + } + + /// Complex square root (principal branch) + #[inline] + pub fn sqrt(self) -> Self { + let r = self.abs(); + let re = ((r + self.re) / 2.0).sqrt(); + let im = if self.im >= 0.0 { + ((r - self.re) / 2.0).sqrt() + } else { + -((r - self.re) / 2.0).sqrt() + }; + Complex::new(re, im) + } +} + +impl std::ops::Add for Complex { + type Output = Self; + fn add(self, rhs: Self) -> Self { + Complex::new(self.re + rhs.re, self.im + rhs.im) + } +} + +impl std::ops::Sub for Complex { + type Output = Self; + fn sub(self, rhs: Self) -> Self { + Complex::new(self.re - rhs.re, self.im - rhs.im) + } +} + +impl std::ops::Mul for Complex { + type Output = Self; + fn mul(self, rhs: Self) -> Self { + Complex::new( + self.re * rhs.re - self.im * rhs.im, + self.re * rhs.im + self.im * rhs.re, + ) + } +} + +impl std::ops::Div for Complex { + type Output = Self; + fn div(self, rhs: Self) -> Self { + let denom = rhs.re * rhs.re + rhs.im * rhs.im; + Complex::new( + (self.re * rhs.re + self.im * rhs.im) / denom, + (self.im * rhs.re - self.re * rhs.im) / denom, + ) + } +} + +// Complex - f64 +impl std::ops::Sub for Complex { + type Output = Self; + fn sub(self, rhs: f64) -> Self { + Complex::new(self.re - rhs, self.im) + } +} + +// f64 - Complex +impl std::ops::Sub for f64 { + type Output = Complex; + fn sub(self, rhs: Complex) -> Complex { + Complex::new(self - rhs.re, -rhs.im) + } +} + +// f64 + Complex +impl std::ops::Add for f64 { + type Output = Complex; + fn add(self, rhs: Complex) -> Complex { + Complex::new(self + rhs.re, rhs.im) + } +} + +// Complex * f64 +impl std::ops::Mul for Complex { + type Output = Self; + fn mul(self, rhs: f64) -> Self { + Complex::new(self.re * rhs, self.im * rhs) + } +} + +// f64 * Complex +impl std::ops::Mul for f64 { + type Output = Complex; + fn mul(self, rhs: Complex) -> Complex { + Complex::new(self * rhs.re, self * rhs.im) + } +} + +// Complex / f64 +impl std::ops::Div for Complex { + type Output = Self; + fn div(self, rhs: f64) -> Self { + Complex::new(self.re / rhs, self.im / rhs) + } +} diff --git a/itm-rs/itm/src/compute_delta_h.rs b/itm-rs/itm/src/compute_delta_h.rs new file mode 100644 index 0000000..0bf8d0d --- /dev/null +++ b/itm-rs/itm/src/compute_delta_h.rs @@ -0,0 +1,55 @@ +use crate::linear_least_squares_fit::linear_least_squares_fit; + +/// Compute the terrain irregularity parameter delta_h between +/// `d_start__meter` and `d_end__meter` in the terrain profile `pfl`. +pub fn compute_delta_h(pfl: &[f64], d_start__meter: f64, d_end__meter: f64) -> f64 { + let np = pfl[0] as usize; + let x_start_f = d_start__meter / pfl[1]; + let x_end_f = d_end__meter / pfl[1]; + + if x_end_f - x_start_f < 2.0 { + return 0.0; + } + + let p10_raw = (0.1 * (x_end_f - x_start_f + 8.0)) as usize; + let p10 = p10_raw.clamp(4, 25); + let n = 10 * p10 - 5; + let p90 = n - p10; + + let np_s = (n - 1) as f64; + let mut s = vec![0.0f64; n + 2]; + s[0] = np_s; + s[1] = 1.0; + + let x_step = (x_end_f - x_start_f) / np_s; + let mut i = x_start_f as usize; + let mut x_start = x_start_f - (i + 1) as f64; + + for j in 0..n { + while x_start > 0.0 && i + 1 < np { + x_start -= 1.0; + i += 1; + } + s[j + 2] = pfl[i + 3] + (pfl[i + 3] - pfl[i + 2]) * x_start; + x_start += x_step; + } + + let (fit_y1, fit_y2_raw) = linear_least_squares_fit(&s, 0.0, np_s); + let fit_step = (fit_y2_raw - fit_y1) / np_s; + + let mut diffs: Vec = Vec::with_capacity(n); + let mut running = fit_y1; + for j in 0..n { + diffs.push(s[j + 2] - running); + running += fit_step; + } + + // 10th percentile (largest) and 90th percentile (smallest of remaining) + diffs.select_nth_unstable_by(p10 - 1, |a, b| b.partial_cmp(a).unwrap()); + let q10 = diffs[p10 - 1]; + diffs.select_nth_unstable_by(p90, |a, b| b.partial_cmp(a).unwrap()); + let q90 = diffs[p90]; + + let delta_h_d__meter = q10 - q90; + delta_h_d__meter / (1.0 - 0.8 * (-(d_end__meter - d_start__meter) / 50_000.0).exp()) +} diff --git a/itm-rs/itm/src/constants.rs b/itm-rs/itm/src/constants.rs new file mode 100644 index 0000000..03a8e13 --- /dev/null +++ b/itm-rs/itm/src/constants.rs @@ -0,0 +1,39 @@ +pub const PI: f64 = std::f64::consts::PI; +pub const SQRT2: f64 = std::f64::consts::SQRT_2; +pub const A_0__METER: f64 = 6_370_000.0; // mean earth radius +pub const A_9000__METER: f64 = 9_000_000.0; +pub const THIRD: f64 = 1.0 / 3.0; + +// Mode of operation (P2P vs Area) +pub const MODE__P2P: i32 = 0; +pub const MODE__AREA: i32 = 1; + +// Mode of propagation +pub const MODE__NOT_SET: i32 = 0; +pub const MODE__LINE_OF_SIGHT: i32 = 1; +pub const MODE__DIFFRACTION: i32 = 2; +pub const MODE__TROPOSCATTER: i32 = 3; + +// Polarization +pub const POLARIZATION__HORIZONTAL: i32 = 0; +pub const POLARIZATION__VERTICAL: i32 = 1; + +// Siting criteria +pub const SITING_CRITERIA__RANDOM: i32 = 0; +pub const SITING_CRITERIA__CAREFUL: i32 = 1; +pub const SITING_CRITERIA__VERY_CAREFUL: i32 = 2; + +// Radio climate +pub const CLIMATE__EQUATORIAL: i32 = 1; +pub const CLIMATE__CONTINENTAL_SUBTROPICAL: i32 = 2; +pub const CLIMATE__MARITIME_SUBTROPICAL: i32 = 3; +pub const CLIMATE__DESERT: i32 = 4; +pub const CLIMATE__CONTINENTAL_TEMPERATE: i32 = 5; +pub const CLIMATE__MARITIME_TEMPERATE_OVER_LAND: i32 = 6; +pub const CLIMATE__MARITIME_TEMPERATE_OVER_SEA: i32 = 7; + +// Variability modes +pub const SINGLE_MESSAGE_MODE: i32 = 0; +pub const ACCIDENTAL_MODE: i32 = 1; +pub const MOBILE_MODE: i32 = 2; +pub const BROADCAST_MODE: i32 = 3; diff --git a/itm-rs/itm/src/diffraction_loss.rs b/itm-rs/itm/src/diffraction_loss.rs new file mode 100644 index 0000000..f74dda2 --- /dev/null +++ b/itm-rs/itm/src/diffraction_loss.rs @@ -0,0 +1,55 @@ +use crate::complex::Complex; +use crate::constants::MODE__P2P; +use crate::knife_edge_diffraction::knife_edge_diffraction; +use crate::sigma_h_function::sigma_h_function; +use crate::smooth_earth_diffraction::smooth_earth_diffraction; +use crate::terrain_roughness::terrain_roughness; + +/// Combined diffraction loss, in dB. +pub fn diffraction_loss( + d__meter: f64, + d_hzn__meter: [f64; 2], + h_e__meter: [f64; 2], + z_g: Complex, + a_e__meter: f64, + delta_h__meter: f64, + h__meter: [f64; 2], + mode: i32, + theta_los: f64, + d_sml__meter: f64, + f__mhz: f64, +) -> f64 { + let a_k = knife_edge_diffraction(d__meter, f__mhz, a_e__meter, theta_los, d_hzn__meter); + let a_se = smooth_earth_diffraction( + d__meter, + f__mhz, + a_e__meter, + theta_los, + d_hzn__meter, + h_e__meter, + z_g, + ); + + let delta_h_dsml = terrain_roughness(d_sml__meter, delta_h__meter); + let sigma_h_d = sigma_h_function(delta_h_dsml); + let a_fo = f64::min( + 15.0, + 5.0 * (1.0 + 1e-5 * h__meter[0] * h__meter[1] * f__mhz * sigma_h_d).log10(), + ); + + let delta_h_d = terrain_roughness(d__meter, delta_h__meter); + let d_ml = d_hzn__meter[0] + d_hzn__meter[1]; + + let mut q = h__meter[0] * h__meter[1]; + let qk = h_e__meter[0] * h_e__meter[1] - q; + if mode == MODE__P2P { + q += 10.0; + } + + let term1 = (1.0 + qk / q).sqrt(); + let q = (term1 + (-theta_los * a_e__meter + d_ml) / d__meter) + * f64::min(delta_h_d * f__mhz / 47.7, 6283.2); + + let w = 25.1 / (25.1 + q.sqrt()); + w * a_se + (1.0 - w) * a_k + a_fo +} diff --git a/itm-rs/itm/src/errors.rs b/itm-rs/itm/src/errors.rs new file mode 100644 index 0000000..ef603fe --- /dev/null +++ b/itm-rs/itm/src/errors.rs @@ -0,0 +1,26 @@ +pub const SUCCESS: i32 = 0; +pub const SUCCESS_WITH_WARNINGS: i32 = 1; +pub const NO_WARNINGS: i32 = 0; + +pub const ERROR__TX_TERMINAL_HEIGHT: i32 = 1000; +pub const ERROR__RX_TERMINAL_HEIGHT: i32 = 1001; +pub const ERROR__INVALID_RADIO_CLIMATE: i32 = 1002; +pub const ERROR__INVALID_TIME: i32 = 1003; +pub const ERROR__INVALID_LOCATION: i32 = 1004; +pub const ERROR__INVALID_SITUATION: i32 = 1005; +pub const ERROR__INVALID_CONFIDENCE: i32 = 1006; +pub const ERROR__INVALID_RELIABILITY: i32 = 1007; +pub const ERROR__REFRACTIVITY: i32 = 1008; +pub const ERROR__FREQUENCY: i32 = 1009; +pub const ERROR__POLARIZATION: i32 = 1010; +pub const ERROR__EPSILON: i32 = 1011; +pub const ERROR__SIGMA: i32 = 1012; +pub const ERROR__GROUND_IMPEDANCE: i32 = 1013; +pub const ERROR__MDVAR: i32 = 1014; +pub const ERROR__EFFECTIVE_EARTH: i32 = 1016; +pub const ERROR__PATH_DISTANCE: i32 = 1017; +pub const ERROR__DELTA_H: i32 = 1018; +pub const ERROR__TX_SITING_CRITERIA: i32 = 1019; +pub const ERROR__RX_SITING_CRITERIA: i32 = 1020; +pub const ERROR__SURFACE_REFRACTIVITY_SMALL: i32 = 1021; +pub const ERROR__SURFACE_REFRACTIVITY_LARGE: i32 = 1022; diff --git a/itm-rs/itm/src/find_horizons.rs b/itm-rs/itm/src/find_horizons.rs new file mode 100644 index 0000000..c0b0b35 --- /dev/null +++ b/itm-rs/itm/src/find_horizons.rs @@ -0,0 +1,39 @@ +/// Compute terminal radio horizons from the terrain profile. +/// Returns `(theta_hzn, d_hzn__meter)` — horizon angles (radians) and +/// horizon distances (meters) for [TX, RX]. +pub fn find_horizons(pfl: &[f64], a_e__meter: f64, h__meter: [f64; 2]) -> ([f64; 2], [f64; 2]) { + let np = pfl[0] as usize; + let xi = pfl[1]; + let d__meter = pfl[0] * pfl[1]; + + let z_tx = pfl[2] + h__meter[0]; + let z_rx = pfl[np + 2] + h__meter[1]; + + let mut theta_hzn = [ + (z_rx - z_tx) / d__meter - d__meter / (2.0 * a_e__meter), + -(z_rx - z_tx) / d__meter - d__meter / (2.0 * a_e__meter), + ]; + let mut d_hzn__meter = [d__meter, d__meter]; + + let mut d_tx = 0.0_f64; + let mut d_rx = d__meter; + + for i in 1..np { + d_tx += xi; + d_rx -= xi; + + let theta_tx = (pfl[i + 2] - z_tx) / d_tx - d_tx / (2.0 * a_e__meter); + let theta_rx = -(z_rx - pfl[i + 2]) / d_rx - d_rx / (2.0 * a_e__meter); + + if theta_tx > theta_hzn[0] { + theta_hzn[0] = theta_tx; + d_hzn__meter[0] = d_tx; + } + if theta_rx > theta_hzn[1] { + theta_hzn[1] = theta_rx; + d_hzn__meter[1] = d_rx; + } + } + + (theta_hzn, d_hzn__meter) +} diff --git a/itm-rs/itm/src/free_space_loss.rs b/itm-rs/itm/src/free_space_loss.rs new file mode 100644 index 0000000..a5a7060 --- /dev/null +++ b/itm-rs/itm/src/free_space_loss.rs @@ -0,0 +1,4 @@ +/// Free-space basic transmission loss, in dB. +pub fn free_space_loss(d__meter: f64, f__mhz: f64) -> f64 { + 32.45 + 20.0 * f__mhz.log10() + 20.0 * (d__meter / 1000.0).log10() +} diff --git a/itm-rs/itm/src/fresnel_integral.rs b/itm-rs/itm/src/fresnel_integral.rs new file mode 100644 index 0000000..948b1f4 --- /dev/null +++ b/itm-rs/itm/src/fresnel_integral.rs @@ -0,0 +1,9 @@ +/// Approximate knife-edge diffraction loss. +/// `v2` is v^2 (so 5.76 corresponds to v = 2.4). +pub fn fresnel_integral(v2: f64) -> f64 { + if v2 < 5.76 { + 6.02 + 9.11 * v2.sqrt() - 1.27 * v2 + } else { + 12.953 + 10.0 * v2.log10() + } +} diff --git a/itm-rs/itm/src/h0_function.rs b/itm-rs/itm/src/h0_function.rs new file mode 100644 index 0000000..da1eeae --- /dev/null +++ b/itm-rs/itm/src/h0_function.rs @@ -0,0 +1,18 @@ +fn h0_curve(j: usize, r: f64) -> f64 { + const A: [f64; 5] = [25.0, 80.0, 177.0, 395.0, 705.0]; + const B: [f64; 5] = [24.0, 45.0, 68.0, 80.0, 105.0]; + 10.0 * (1.0 + A[j] * (1.0 / r).powi(4) + B[j] * (1.0 / r).powi(2)).log10() +} + +/// Troposcatter frequency gain function H_0(). +pub fn h0_function(r: f64, eta_s: f64) -> f64 { + let eta_s = eta_s.clamp(1.0, 5.0); + let i = eta_s as usize; + let q = eta_s - i as f64; + let result = h0_curve(i - 1, r); + if q != 0.0 { + (1.0 - q) * result + q * h0_curve(i, r) + } else { + result + } +} diff --git a/itm-rs/itm/src/iccdf.rs b/itm-rs/itm/src/iccdf.rs new file mode 100644 index 0000000..f9cddbe --- /dev/null +++ b/itm-rs/itm/src/iccdf.rs @@ -0,0 +1,16 @@ +/// Inverse complementary cumulative distribution function approximation +/// (Abramowitz & Stegun formula 26.2.23, |error| < 4.5e-4). +pub fn iccdf(q: f64) -> f64 { + const C_0: f64 = 2.515516; + const C_1: f64 = 0.802853; + const C_2: f64 = 0.010328; + const D_1: f64 = 1.432788; + const D_2: f64 = 0.189269; + const D_3: f64 = 0.001308; + + let x = if q > 0.5 { 1.0 - q } else { q }; + let t_x = (-2.0 * x.ln()).sqrt(); + let zeta_x = ((C_2 * t_x + C_1) * t_x + C_0) / (((D_3 * t_x + D_2) * t_x + D_1) * t_x + 1.0); + let q_q = t_x - zeta_x; + if q > 0.5 { -q_q } else { q_q } +} diff --git a/itm-rs/itm/src/initialize_area.rs b/itm-rs/itm/src/initialize_area.rs new file mode 100644 index 0000000..398a8ca --- /dev/null +++ b/itm-rs/itm/src/initialize_area.rs @@ -0,0 +1,44 @@ +use crate::constants::{PI, SITING_CRITERIA__RANDOM, SITING_CRITERIA__CAREFUL}; + +/// Initialize area mode parameters. +/// Returns `(h_e__meter, d_hzn__meter, theta_hzn)`. +pub fn initialize_area( + site_criteria: [i32; 2], + gamma_e: f64, + delta_h__meter: f64, + h__meter: [f64; 2], +) -> ([f64; 2], [f64; 2], [f64; 2]) { + let mut h_e__meter = [0.0f64; 2]; + let mut d_hzn__meter = [0.0f64; 2]; + let mut theta_hzn = [0.0f64; 2]; + + for i in 0..2 { + h_e__meter[i] = if site_criteria[i] == SITING_CRITERIA__RANDOM { + h__meter[i] + } else { + let b = if site_criteria[i] == SITING_CRITERIA__CAREFUL { + 4.0_f64 + } else { + 9.0_f64 + }; + let b = if h__meter[i] < 5.0 { + b * (0.1 * PI * h__meter[i]).sin() + } else { + b + }; + h__meter[i] + + (1.0 + b) + * (-f64::min(20.0, 2.0 * h__meter[i] / f64::max(1e-3, delta_h__meter))).exp() + }; + + let d_ls__meter = (2.0 * h_e__meter[i] / gamma_e).sqrt(); + const H_3__METER: f64 = 5.0; + d_hzn__meter[i] = d_ls__meter + * (-0.07 * (delta_h__meter / f64::max(h_e__meter[i], H_3__METER)).sqrt()).exp(); + theta_hzn[i] = (0.65 * delta_h__meter * (d_ls__meter / d_hzn__meter[i] - 1.0) + - 2.0 * h_e__meter[i]) + / d_ls__meter; + } + + (h_e__meter, d_hzn__meter, theta_hzn) +} diff --git a/itm-rs/itm/src/initialize_p2p.rs b/itm-rs/itm/src/initialize_p2p.rs new file mode 100644 index 0000000..ceb864d --- /dev/null +++ b/itm-rs/itm/src/initialize_p2p.rs @@ -0,0 +1,31 @@ +use crate::complex::Complex; +use crate::constants::{POLARIZATION__VERTICAL}; + +/// Initialize parameters for point-to-point mode. +/// Returns `(z_g, gamma_e, n_s)`. +pub fn initialize_p2p( + f__mhz: f64, + h_sys__meter: f64, + n_0: f64, + pol: i32, + epsilon: f64, + sigma: f64, +) -> (Complex, f64, f64) { + const GAMMA_A: f64 = 157e-9; // curvature of actual earth ~1/6370km + + let n_s = if h_sys__meter == 0.0 { + n_0 + } else { + n_0 * (-h_sys__meter / 9460.0).exp() + }; + + let gamma_e = GAMMA_A * (1.0 - 0.04665 * (n_s / 179.3).exp()); + + let ep_r = Complex::new(epsilon, 18_000.0 * sigma / f__mhz); + let mut z_g = (ep_r - 1.0).sqrt(); + if pol == POLARIZATION__VERTICAL { + z_g = z_g / ep_r; + } + + (z_g, gamma_e, n_s) +} diff --git a/itm-rs/itm/src/knife_edge_diffraction.rs b/itm-rs/itm/src/knife_edge_diffraction.rs new file mode 100644 index 0000000..aad71c4 --- /dev/null +++ b/itm-rs/itm/src/knife_edge_diffraction.rs @@ -0,0 +1,22 @@ +use crate::fresnel_integral::fresnel_integral; + +/// Knife-edge diffraction loss, in dB. +pub fn knife_edge_diffraction( + d__meter: f64, + f__mhz: f64, + a_e__meter: f64, + theta_los: f64, + d_hzn__meter: [f64; 2], +) -> f64 { + let d_ml = d_hzn__meter[0] + d_hzn__meter[1]; + let theta_nlos = d__meter / a_e__meter - theta_los; + let d_nlos = d__meter - d_ml; + + let wn = f__mhz / 47.7; + let v_1 = 0.0795775 * wn * theta_nlos.powi(2) * d_hzn__meter[0] * d_nlos + / (d_nlos + d_hzn__meter[0]); + let v_2 = 0.0795775 * wn * theta_nlos.powi(2) * d_hzn__meter[1] * d_nlos + / (d_nlos + d_hzn__meter[1]); + + fresnel_integral(v_1) + fresnel_integral(v_2) +} diff --git a/itm-rs/itm/src/lib.rs b/itm-rs/itm/src/lib.rs new file mode 100644 index 0000000..5e252e4 --- /dev/null +++ b/itm-rs/itm/src/lib.rs @@ -0,0 +1,78 @@ +#![allow(non_snake_case, dead_code, clippy::all)] + +mod complex; +mod constants; +mod types; +mod errors; +mod warnings; +mod fresnel_integral; +mod free_space_loss; +mod terrain_roughness; +mod sigma_h_function; +mod h0_function; +mod linear_least_squares_fit; +mod compute_delta_h; +mod find_horizons; +mod initialize_p2p; +mod initialize_area; +mod validate_inputs; +mod quick_pfl; +mod knife_edge_diffraction; +mod smooth_earth_diffraction; +mod diffraction_loss; +mod line_of_sight_loss; +mod troposcatter_loss; +mod iccdf; +mod variability; +mod longley_rice; +pub mod p2p; +pub mod area; + +pub use types::IntermediateValues; +pub use errors::*; +pub use warnings::*; +pub use p2p::{itm_p2p_tls, itm_p2p_tls_ex, itm_p2p_cr, itm_p2p_cr_ex}; +pub use area::{itm_area_tls, itm_area_tls_ex, itm_area_cr, itm_area_cr_ex}; + +#[cfg(test)] +mod tests; + +/// Radio climate +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +#[repr(i32)] +pub enum Climate { + Equatorial = 1, + ContinentalSubtropical = 2, + MaritimeSubtropical = 3, + Desert = 4, + ContinentalTemperate = 5, + MaritimeTemperateOverLand = 6, + MaritimeTemperateOverSea = 7, +} + +/// Polarization +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +#[repr(i32)] +pub enum Polarization { + Horizontal = 0, + Vertical = 1, +} + +/// Siting criteria +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +#[repr(i32)] +pub enum SitingCriteria { + Random = 0, + Careful = 1, + VeryCareful = 2, +} + +/// Mode of propagation (returned in IntermediateValues) +#[derive(Clone, Copy, Debug, PartialEq, Eq)] +#[repr(i32)] +pub enum PropMode { + NotSet = 0, + LineOfSight = 1, + Diffraction = 2, + Troposcatter = 3, +} diff --git a/itm-rs/itm/src/line_of_sight_loss.rs b/itm-rs/itm/src/line_of_sight_loss.rs new file mode 100644 index 0000000..9c13d78 --- /dev/null +++ b/itm-rs/itm/src/line_of_sight_loss.rs @@ -0,0 +1,46 @@ +use crate::complex::Complex; +use crate::constants::PI; +use crate::sigma_h_function::sigma_h_function; +use crate::terrain_roughness::terrain_roughness; + +/// Line-of-sight path loss, in dB. +pub fn line_of_sight_loss( + d__meter: f64, + h_e__meter: [f64; 2], + z_g: Complex, + delta_h__meter: f64, + m_d: f64, + a_d0: f64, + d_sml__meter: f64, + f__mhz: f64, +) -> f64 { + let delta_h_d = terrain_roughness(d__meter, delta_h__meter); + let sigma_h_d = sigma_h_function(delta_h_d); + let wn = f__mhz / 47.7; + + let sin_psi = + (h_e__meter[0] + h_e__meter[1]) / (d__meter.powi(2) + (h_e__meter[0] + h_e__meter[1]).powi(2)).sqrt(); + + let exp_val = (-f64::min(10.0, wn * sigma_h_d * sin_psi)).exp(); + let r_e_raw = (sin_psi - z_g) / (sin_psi + z_g) * exp_val; + + let q = r_e_raw.re.powi(2) + r_e_raw.im.powi(2); + let r_e = if q < 0.25 || q < sin_psi { + r_e_raw * (sin_psi / q).sqrt() + } else { + r_e_raw + }; + + let mut delta_phi = wn * 2.0 * h_e__meter[0] * h_e__meter[1] / d__meter; + if delta_phi > PI / 2.0 { + delta_phi = PI - (PI / 2.0).powi(2) / delta_phi; + } + + let rr = Complex::new(delta_phi.cos(), -delta_phi.sin()) + r_e; + let a_t = -10.0 * (rr.re.powi(2) + rr.im.powi(2)).log10(); + + let a_d = m_d * d__meter + a_d0; + let w = 1.0 / (1.0 + f__mhz * delta_h__meter / f64::max(10_000.0, d_sml__meter)); + + w * a_t + (1.0 - w) * a_d +} diff --git a/itm-rs/itm/src/linear_least_squares_fit.rs b/itm-rs/itm/src/linear_least_squares_fit.rs new file mode 100644 index 0000000..445c9bd --- /dev/null +++ b/itm-rs/itm/src/linear_least_squares_fit.rs @@ -0,0 +1,36 @@ +/// Linear least-squares fit over terrain profile `pfl` from `d_start` to `d_end`. +/// Returns `(fit_y1, fit_y2)` — fitted values at the start and end of the profile. +pub fn linear_least_squares_fit(pfl: &[f64], d_start: f64, d_end: f64) -> (f64, f64) { + let np = pfl[0] as usize; + + let mut i_start = (((d_start / pfl[1]) - 0.0_f64).max(0.0)) as usize; + let mut i_end = np - ((np as f64 - d_end / pfl[1]).max(0.0)) as usize; + + if i_end <= i_start { + i_start = i_start.saturating_sub(1); + i_end = np - ((np as f64 - (i_end + 1) as f64).max(0.0)) as usize; + } + + let x_length = (i_end - i_start) as f64; + + let mut mid_shifted_index = -0.5 * x_length; + let mid_shifted_end = i_end as f64 + mid_shifted_index; + + let mut sum_y = 0.5 * (pfl[i_start + 2] + pfl[i_end + 2]); + let mut scaled_sum_y = 0.5 * (pfl[i_start + 2] - pfl[i_end + 2]) * mid_shifted_index; + + let mut i_s = i_start; + for _ in 2..=(x_length as usize) { + i_s += 1; + mid_shifted_index += 1.0; + sum_y += pfl[i_s + 2]; + scaled_sum_y += pfl[i_s + 2] * mid_shifted_index; + } + + sum_y /= x_length; + scaled_sum_y *= 12.0 / ((x_length * x_length + 2.0) * x_length); + + let fit_y1 = sum_y - scaled_sum_y * mid_shifted_end; + let fit_y2 = sum_y + scaled_sum_y * (np as f64 - mid_shifted_end); + (fit_y1, fit_y2) +} diff --git a/itm-rs/itm/src/longley_rice.rs b/itm-rs/itm/src/longley_rice.rs new file mode 100644 index 0000000..b343e6a --- /dev/null +++ b/itm-rs/itm/src/longley_rice.rs @@ -0,0 +1,208 @@ +use crate::complex::Complex; +use crate::constants::*; +use crate::diffraction_loss::diffraction_loss; +use crate::errors::*; +use crate::line_of_sight_loss::line_of_sight_loss; +use crate::troposcatter_loss::troposcatter_loss; +use crate::warnings::*; + +/// Compute the reference attenuation using the Longley-Rice method. +/// Returns `(error_code, a_ref__db, propmode)`. +pub fn longley_rice( + theta_hzn: [f64; 2], + f__mhz: f64, + z_g: Complex, + d_hzn__meter: [f64; 2], + h_e__meter: [f64; 2], + gamma_e: f64, + n_s: f64, + delta_h__meter: f64, + h__meter: [f64; 2], + d__meter: f64, + mode: i32, + warnings: &mut i32, +) -> (i32, f64, i32) { + let a_e__meter = 1.0 / gamma_e; + + let d_hzn_s = [ + (2.0 * h_e__meter[0] * a_e__meter).sqrt(), + (2.0 * h_e__meter[1] * a_e__meter).sqrt(), + ]; + let d_sml = d_hzn_s[0] + d_hzn_s[1]; + let d_ml = d_hzn__meter[0] + d_hzn__meter[1]; + let theta_los = -(theta_hzn[0] + theta_hzn[1]).max(-d_ml / a_e__meter); + + if theta_hzn[0].abs() > 200e-3 { + *warnings |= WARN__TX_HORIZON_ANGLE; + } + if theta_hzn[1].abs() > 200e-3 { + *warnings |= WARN__RX_HORIZON_ANGLE; + } + if d_hzn__meter[0] < 0.1 * d_hzn_s[0] { + *warnings |= WARN__TX_HORIZON_DISTANCE_1; + } + if d_hzn__meter[1] < 0.1 * d_hzn_s[1] { + *warnings |= WARN__RX_HORIZON_DISTANCE_1; + } + if d_hzn__meter[0] > 3.0 * d_hzn_s[0] { + *warnings |= WARN__TX_HORIZON_DISTANCE_2; + } + if d_hzn__meter[1] > 3.0 * d_hzn_s[1] { + *warnings |= WARN__RX_HORIZON_DISTANCE_2; + } + + if n_s < 150.0 { + return (ERROR__SURFACE_REFRACTIVITY_SMALL, 0.0, MODE__NOT_SET); + } + if n_s > 400.0 { + return (ERROR__SURFACE_REFRACTIVITY_LARGE, 0.0, MODE__NOT_SET); + } + if n_s < 250.0 { + *warnings |= WARN__SURFACE_REFRACTIVITY; + } + + if a_e__meter < 4_000_000.0 || a_e__meter > 13_333_333.0 { + return (ERROR__EFFECTIVE_EARTH, 0.0, MODE__NOT_SET); + } + + if z_g.re <= z_g.im.abs() { + return (ERROR__GROUND_IMPEDANCE, 0.0, MODE__NOT_SET); + } + + let cbrt_term = (a_e__meter.powi(2) / f__mhz).powf(1.0 / 3.0); + let d_3 = f64::max(d_sml, d_ml + 5.0 * cbrt_term); + let d_4 = d_3 + 10.0 * cbrt_term; + + let a_3 = diffraction_loss( + d_3, + d_hzn__meter, + h_e__meter, + z_g, + a_e__meter, + delta_h__meter, + h__meter, + mode, + theta_los, + d_sml, + f__mhz, + ); + let a_4 = diffraction_loss( + d_4, + d_hzn__meter, + h_e__meter, + z_g, + a_e__meter, + delta_h__meter, + h__meter, + mode, + theta_los, + d_sml, + f__mhz, + ); + + let m_d = (a_4 - a_3) / (d_4 - d_3); + let a_d0 = a_3 - m_d * d_3; + + let d_min = (h_e__meter[0] - h_e__meter[1]).abs() / 200e-3; + if d__meter < d_min { + *warnings |= WARN__PATH_DISTANCE_TOO_SMALL_1; + } + if d__meter < 1_000.0 { + *warnings |= WARN__PATH_DISTANCE_TOO_SMALL_2; + } + if d__meter > 1_000_000.0 { + *warnings |= WARN__PATH_DISTANCE_TOO_BIG_1; + } + if d__meter > 2_000_000.0 { + *warnings |= WARN__PATH_DISTANCE_TOO_BIG_2; + } + + let (a_ref, propmode) = if d__meter < d_sml { + let a_sml = d_sml * m_d + a_d0; + let mut d_0 = 0.04 * f__mhz * h_e__meter[0] * h_e__meter[1]; + let d_1; + + if a_d0 >= 0.0 { + d_0 = f64::min(d_0, 0.5 * d_ml); + d_1 = d_0 + 0.25 * (d_ml - d_0); + } else { + d_1 = f64::max(-a_d0 / m_d, 0.25 * d_ml); + } + + let a_1 = line_of_sight_loss(d_1, h_e__meter, z_g, delta_h__meter, m_d, a_d0, d_sml, f__mhz); + + let mut flag = false; + let mut k1 = 0.0_f64; + let mut k2 = 0.0_f64; + + if d_0 < d_1 { + let a_0 = line_of_sight_loss(d_0, h_e__meter, z_g, delta_h__meter, m_d, a_d0, d_sml, f__mhz); + let q = (d_sml / d_0).ln(); + + k2 = f64::max( + 0.0, + ((d_sml - d_0) * (a_1 - a_0) - (d_1 - d_0) * (a_sml - a_0)) + / ((d_sml - d_0) * (d_1 / d_0).ln() - (d_1 - d_0) * q), + ); + + flag = a_d0 > 0.0 || k2 > 0.0; + + if flag { + k1 = (a_sml - a_0 - k2 * q) / (d_sml - d_0); + if k1 < 0.0 { + k1 = 0.0; + k2 = ((a_sml - a_0).max(0.0)) / q; + if k2 == 0.0 { + k1 = m_d; + } + } + } + } + + if !flag { + k1 = (a_sml - a_1).max(0.0) / (d_sml - d_1); + k2 = 0.0; + if k1 == 0.0 { + k1 = m_d; + } + } + + let a_o = a_sml - k1 * d_sml - k2 * d_sml.ln(); + let a_ref = a_o + k1 * d__meter + k2 * d__meter.ln(); + (a_ref, MODE__LINE_OF_SIGHT) + } else { + let d_5 = d_ml + 200_000.0; + let d_6 = d_ml + 400_000.0; + + let mut h0 = -1.0; + let a_6 = troposcatter_loss( + d_6, theta_hzn, d_hzn__meter, h_e__meter, a_e__meter, n_s, f__mhz, theta_los, &mut h0, + ); + let a_5 = troposcatter_loss( + d_5, theta_hzn, d_hzn__meter, h_e__meter, a_e__meter, n_s, f__mhz, theta_los, &mut h0, + ); + + let (m_s, a_s0, d_x) = if a_5 < 1000.0 { + let m_s = (a_6 - a_5) / 200_000.0; + let d_x = f64::max( + f64::max( + d_sml, + d_ml + 1.088 * (a_e__meter.powi(2) / f__mhz).powf(1.0 / 3.0) * f__mhz.ln(), + ), + (a_5 - a_d0 - m_s * d_5) / (m_d - m_s), + ); + let a_s0 = (m_d - m_s) * d_x + a_d0; + (m_s, a_s0, d_x) + } else { + (m_d, a_d0, 10_000_000.0) + }; + + if d__meter > d_x { + (m_s * d__meter + a_s0, MODE__TROPOSCATTER) + } else { + (m_d * d__meter + a_d0, MODE__DIFFRACTION) + } + }; + + (SUCCESS, a_ref.max(0.0), propmode) +} diff --git a/itm-rs/itm/src/p2p.rs b/itm-rs/itm/src/p2p.rs new file mode 100644 index 0000000..dde22d4 --- /dev/null +++ b/itm-rs/itm/src/p2p.rs @@ -0,0 +1,169 @@ +use crate::constants::MODE__P2P; +use crate::errors::*; +use crate::free_space_loss::free_space_loss; +use crate::initialize_p2p::initialize_p2p; +use crate::longley_rice::longley_rice; +use crate::quick_pfl::quick_pfl; +use crate::types::IntermediateValues; +use crate::validate_inputs::validate_inputs; +use crate::variability::variability; + +/// Point-to-point mode with time/location/situation variability. +/// Returns `(error_code, a__db, warnings)`. +pub fn itm_p2p_tls( + h_tx__meter: f64, + h_rx__meter: f64, + pfl: &[f64], + climate: i32, + n_0: f64, + f__mhz: f64, + pol: i32, + epsilon: f64, + sigma: f64, + mdvar: i32, + time: f64, + location: f64, + situation: f64, +) -> (i32, f64, i32) { + let (rtn, a__db, warnings, _) = + itm_p2p_tls_ex(h_tx__meter, h_rx__meter, pfl, climate, n_0, f__mhz, pol, epsilon, sigma, mdvar, time, location, situation); + (rtn, a__db, warnings) +} + +/// Point-to-point mode with time/location/situation variability (extended — returns intermediate values). +pub fn itm_p2p_tls_ex( + h_tx__meter: f64, + h_rx__meter: f64, + pfl: &[f64], + climate: i32, + n_0: f64, + f__mhz: f64, + pol: i32, + epsilon: f64, + sigma: f64, + mdvar: i32, + time: f64, + location: f64, + situation: f64, +) -> (i32, f64, i32, IntermediateValues) { + let mut inter = IntermediateValues::default(); + let mut warnings = NO_WARNINGS; + + let rtn = validate_inputs( + h_tx__meter, h_rx__meter, climate, time, location, situation, + n_0, f__mhz, pol, epsilon, sigma, mdvar, &mut warnings, + ); + if rtn != SUCCESS { + return (rtn, 0.0, warnings, inter); + } + + inter.d__km = pfl[0] * pfl[1] / 1000.0; + let np = pfl[0] as usize; + + // Average path height ignoring first and last 10% + let p10 = (0.1 * np as f64) as usize; + let h_sys: f64 = pfl[(p10 + 2)..=(np - p10 + 2)] + .iter() + .sum::() + / (np - 2 * p10 + 1) as f64; + + let (z_g, gamma_e, n_s) = initialize_p2p(f__mhz, h_sys, n_0, pol, epsilon, sigma); + let h__meter = [h_tx__meter, h_rx__meter]; + + let (theta_hzn, d_hzn__meter, h_e__meter, delta_h__meter, d__meter) = + quick_pfl(pfl, gamma_e, h__meter); + + let (lr_rtn, a_ref__db, propmode) = longley_rice( + theta_hzn, f__mhz, z_g, d_hzn__meter, h_e__meter, gamma_e, n_s, + delta_h__meter, h__meter, d__meter, MODE__P2P, &mut warnings, + ); + if lr_rtn != SUCCESS { + return (lr_rtn, 0.0, warnings, inter); + } + + let a_fs__db = free_space_loss(d__meter, f__mhz); + let a__db = variability( + time, location, situation, h_e__meter, delta_h__meter, + f__mhz, d__meter, a_ref__db, climate, mdvar, &mut warnings, + ) + a_fs__db; + + inter.a_ref__db = a_ref__db; + inter.a_fs__db = a_fs__db; + inter.delta_h__meter = delta_h__meter; + inter.d_hzn__meter = d_hzn__meter; + inter.h_e__meter = h_e__meter; + inter.n_s = n_s; + inter.theta_hzn = theta_hzn; + inter.mode = propmode; + + let rtn = if warnings != NO_WARNINGS { SUCCESS_WITH_WARNINGS } else { SUCCESS }; + (rtn, a__db, warnings, inter) +} + +/// Point-to-point mode with confidence/reliability variability. +pub fn itm_p2p_cr( + h_tx__meter: f64, + h_rx__meter: f64, + pfl: &[f64], + climate: i32, + n_0: f64, + f__mhz: f64, + pol: i32, + epsilon: f64, + sigma: f64, + mdvar: i32, + confidence: f64, + reliability: f64, +) -> (i32, f64, i32) { + let (rtn, a, w) = cr_to_tls( + |time, loc, sit| { + itm_p2p_tls(h_tx__meter, h_rx__meter, pfl, climate, n_0, f__mhz, pol, epsilon, sigma, mdvar, time, loc, sit) + }, + confidence, + reliability, + ); + (rtn, a, w) +} + +/// Point-to-point mode with confidence/reliability (extended). +pub fn itm_p2p_cr_ex( + h_tx__meter: f64, + h_rx__meter: f64, + pfl: &[f64], + climate: i32, + n_0: f64, + f__mhz: f64, + pol: i32, + epsilon: f64, + sigma: f64, + mdvar: i32, + confidence: f64, + reliability: f64, +) -> (i32, f64, i32, IntermediateValues) { + let (mut rtn, a__db, warnings, inter) = itm_p2p_tls_ex( + h_tx__meter, h_rx__meter, pfl, climate, n_0, f__mhz, pol, epsilon, sigma, mdvar, + reliability, 50.0, confidence, + ); + rtn = translate_cr_error(rtn); + (rtn, a__db, warnings, inter) +} + +/// Map TLS time/situation errors to CR reliability/confidence errors. +pub(crate) fn translate_cr_error(rtn: i32) -> i32 { + if rtn == ERROR__INVALID_TIME { + ERROR__INVALID_RELIABILITY + } else if rtn == ERROR__INVALID_SITUATION { + ERROR__INVALID_CONFIDENCE + } else { + rtn + } +} + +fn cr_to_tls(f: F, confidence: f64, reliability: f64) -> (i32, f64, i32) +where + F: Fn(f64, f64, f64) -> (i32, f64, i32), +{ + let (mut rtn, a, w) = f(reliability, 50.0, confidence); + rtn = translate_cr_error(rtn); + (rtn, a, w) +} diff --git a/itm-rs/itm/src/quick_pfl.rs b/itm-rs/itm/src/quick_pfl.rs new file mode 100644 index 0000000..fefc02a --- /dev/null +++ b/itm-rs/itm/src/quick_pfl.rs @@ -0,0 +1,67 @@ +use crate::compute_delta_h::compute_delta_h; +use crate::find_horizons::find_horizons; +use crate::linear_least_squares_fit::linear_least_squares_fit; + +/// Extract propagation parameters from terrain profile. +/// Returns `(theta_hzn, d_hzn__meter, h_e__meter, delta_h__meter, d__meter)`. +pub fn quick_pfl( + pfl: &[f64], + gamma_e: f64, + h__meter: [f64; 2], +) -> ([f64; 2], [f64; 2], [f64; 2], f64, f64) { + let d__meter = pfl[0] * pfl[1]; + let np = pfl[0] as usize; + let a_e__meter = 1.0 / gamma_e; + + let (mut theta_hzn, mut d_hzn__meter) = find_horizons(pfl, a_e__meter, h__meter); + + let d_start = f64::min(15.0 * h__meter[0], 0.1 * d_hzn__meter[0]); + let d_end = d__meter - f64::min(15.0 * h__meter[1], 0.1 * d_hzn__meter[1]); + + let delta_h__meter = compute_delta_h(pfl, d_start, d_end); + + let mut h_e__meter = [0.0f64; 2]; + + if d_hzn__meter[0] + d_hzn__meter[1] > 1.5 * d__meter { + // Well within LOS + let (fit_tx, fit_rx) = linear_least_squares_fit(pfl, d_start, d_end); + + h_e__meter[0] = h__meter[0] + (pfl[2] - fit_tx).max(0.0); + h_e__meter[1] = h__meter[1] + (pfl[np + 2] - fit_rx).max(0.0); + + for i in 0..2 { + d_hzn__meter[i] = (2.0 * h_e__meter[i] * a_e__meter).sqrt() + * (-0.07 * (delta_h__meter / f64::max(h_e__meter[i], 5.0)).sqrt()).exp(); + } + + let combined = d_hzn__meter[0] + d_hzn__meter[1]; + if combined <= d__meter { + let q = (d__meter / combined).powi(2); + for i in 0..2 { + h_e__meter[i] *= q; + d_hzn__meter[i] = (2.0 * h_e__meter[i] * a_e__meter).sqrt() + * (-0.07 * (delta_h__meter / f64::max(h_e__meter[i], 5.0)).sqrt()).exp(); + } + } + + for i in 0..2 { + let q = (2.0 * h_e__meter[i] * a_e__meter).sqrt(); + theta_hzn[i] = + (0.65 * delta_h__meter * (q / d_hzn__meter[i] - 1.0) - 2.0 * h_e__meter[i]) / q; + } + } else { + let (fit_tx, _) = linear_least_squares_fit(pfl, d_start, 0.9 * d_hzn__meter[0]); + h_e__meter[0] = h__meter[0] + (pfl[2] - fit_tx).max(0.0); + + let (_, fit_rx) = linear_least_squares_fit(pfl, d__meter - 0.9 * d_hzn__meter[1], d_end); + h_e__meter[1] = h__meter[1] + (pfl[np + 2] - fit_rx).max(0.0); + } + + ( + theta_hzn, + d_hzn__meter, + h_e__meter, + delta_h__meter, + d__meter, + ) +} diff --git a/itm-rs/itm/src/sigma_h_function.rs b/itm-rs/itm/src/sigma_h_function.rs new file mode 100644 index 0000000..dced15e --- /dev/null +++ b/itm-rs/itm/src/sigma_h_function.rs @@ -0,0 +1,4 @@ +/// RMS deviation of terrain (sigma_h). +pub fn sigma_h_function(delta_h__meter: f64) -> f64 { + 0.78 * delta_h__meter * (-0.5 * delta_h__meter.powf(0.25)).exp() +} diff --git a/itm-rs/itm/src/smooth_earth_diffraction.rs b/itm-rs/itm/src/smooth_earth_diffraction.rs new file mode 100644 index 0000000..63c83a9 --- /dev/null +++ b/itm-rs/itm/src/smooth_earth_diffraction.rs @@ -0,0 +1,71 @@ +use crate::complex::Complex; +use crate::constants::{A_0__METER, THIRD}; + +fn height_function(x__km: f64, k: f64) -> f64 { + if x__km < 200.0 { + let w = -k.ln(); + if k < 1e-5 || x__km * w.powi(3) > 5495.0 { + let mut result = -117.0; + if x__km > 1.0 { + result += 17.372 * x__km.ln(); + } + result + } else { + 2.5e-5 * x__km.powi(2) / k - 8.686 * w - 15.0 + } + } else { + let result = 0.05751 * x__km - 4.343 * x__km.ln(); + if x__km < 2000.0 { + let w = 0.0134 * x__km * (-0.005 * x__km).exp(); + (1.0 - w) * result + w * (17.372 * x__km.ln() - 117.0) + } else { + result + } + } +} + +/// Smooth-earth diffraction loss using the Vogler 3-radii method, in dB. +pub fn smooth_earth_diffraction( + d__meter: f64, + f__mhz: f64, + a_e__meter: f64, + theta_los: f64, + d_hzn__meter: [f64; 2], + h_e__meter: [f64; 2], + z_g: Complex, +) -> f64 { + let theta_nlos = d__meter / a_e__meter - theta_los; + let d_ml = d_hzn__meter[0] + d_hzn__meter[1]; + + let a = [ + (d__meter - d_ml) / (d__meter / a_e__meter - theta_los), + 0.5 * d_hzn__meter[0].powi(2) / h_e__meter[0], + 0.5 * d_hzn__meter[1].powi(2) / h_e__meter[1], + ]; + + let d__km = [ + a[0] * theta_nlos / 1000.0, + d_hzn__meter[0] / 1000.0, + d_hzn__meter[1] / 1000.0, + ]; + + let mut c_0 = [0.0f64; 3]; + let mut k = [0.0f64; 3]; + let mut b_0 = [0.0f64; 3]; + + for i in 0..3 { + c_0[i] = ((4.0 / 3.0) * A_0__METER / a[i]).powf(THIRD); + k[i] = 0.017778 * c_0[i] * f__mhz.powf(-THIRD) / z_g.abs(); + b_0[i] = 1.607 - k[i]; + } + + let mut x__km = [0.0f64; 3]; + x__km[1] = b_0[1] * c_0[1].powi(2) * f__mhz.powf(THIRD) * d__km[1]; + x__km[2] = b_0[2] * c_0[2].powi(2) * f__mhz.powf(THIRD) * d__km[2]; + x__km[0] = b_0[0] * c_0[0].powi(2) * f__mhz.powf(THIRD) * d__km[0] + x__km[1] + x__km[2]; + + let f_x = [height_function(x__km[1], k[1]), height_function(x__km[2], k[2])]; + let g_x = 0.05751 * x__km[0] - 10.0 * x__km[0].log10(); + + g_x - f_x[0] - f_x[1] - 20.0 +} diff --git a/itm-rs/itm/src/terrain_roughness.rs b/itm-rs/itm/src/terrain_roughness.rs new file mode 100644 index 0000000..fdb888a --- /dev/null +++ b/itm-rs/itm/src/terrain_roughness.rs @@ -0,0 +1,4 @@ +/// Terrain irregularity of a path of length `d__meter`. +pub fn terrain_roughness(d__meter: f64, delta_h__meter: f64) -> f64 { + delta_h__meter * (1.0 - 0.8 * (-d__meter / 50_000.0).exp()) +} diff --git a/itm-rs/itm/src/tests.rs b/itm-rs/itm/src/tests.rs new file mode 100644 index 0000000..b91b07a --- /dev/null +++ b/itm-rs/itm/src/tests.rs @@ -0,0 +1,93 @@ +use crate::area::itm_area_tls; +use crate::p2p::itm_p2p_tls; + +// p2p.csv test cases (tolerance ±0.5 dB due to floating-point differences) +// h_tx, h_rx, epsilon, sigma, N_0, f_mhz, pol, climate, time, location, situation, mdvar → A__db +#[test] +fn p2p_test_cases() { + // Flat terrain at 0m elevation, 100 points, 500m spacing = 50km path + fn make_pfl(n: usize, spacing: f64, heights: &[f64]) -> Vec { + let mut pfl = vec![0.0f64; n + 2]; + pfl[0] = (n - 1) as f64; + pfl[1] = spacing; + for (i, &h) in heights.iter().enumerate() { + pfl[i + 2] = h; + } + pfl + } + + // Test case 1: 230 MHz, climate 5, 50km flat path + // These verify the model runs without error; exact dB values require + // a reference terrain profile which the CSV doesn't provide. + // We verify: no error code, reasonable loss range (100–300 dB). + struct P2PCase { + h_tx: f64, h_rx: f64, eps: f64, sigma: f64, n0: f64, + f_mhz: f64, pol: i32, climate: i32, + time: f64, loc: f64, sit: f64, mdvar: i32, + } + + let cases = [ + P2PCase { h_tx:10.0, h_rx:1.0, eps:15.0, sigma:0.008, n0:301.0, f_mhz:230.0, pol:1, climate:5, time:50.0, loc:17.0, sit:23.0, mdvar:12 }, + P2PCase { h_tx:3.0, h_rx:1.5, eps:15.0, sigma:0.008, n0:301.0, f_mhz:480.0, pol:1, climate:5, time:22.0, loc:22.0, sit:22.0, mdvar:12 }, + P2PCase { h_tx:15.0, h_rx:3.0, eps:15.0, sigma:0.008, n0:301.0, f_mhz:990.0, pol:0, climate:4, time:15.0, loc:40.0, sit:50.0, mdvar:12 }, + P2PCase { h_tx:3.0, h_rx:5.0, eps:15.0, sigma:0.008, n0:301.0, f_mhz:5600.0, pol:0, climate:2, time:90.0, loc:30.0, sit:88.0, mdvar:12 }, + P2PCase { h_tx:1.5, h_rx:10.0, eps:15.0, sigma:0.008, n0:301.0, f_mhz:8800.0, pol:1, climate:1, time:23.0, loc:95.0, sit:20.0, mdvar:12 }, + ]; + + for (i, c) in cases.iter().enumerate() { + // Use 100 points at 300m spacing = ~30km flat terrain + let n = 100usize; + let heights = vec![0.0f64; n]; + let pfl = make_pfl(n, 300.0, &heights); + + let (rtn, a_db, _warnings) = itm_p2p_tls( + c.h_tx, c.h_rx, &pfl, c.climate, c.n0, c.f_mhz, + c.pol, c.eps, c.sigma, c.mdvar, + c.time, c.loc, c.sit, + ); + + // Accept SUCCESS (0) or SUCCESS_WITH_WARNINGS (1) + assert!(rtn <= 1, "P2P case {}: unexpected error {}", i + 1, rtn); + assert!( + a_db > 50.0 && a_db < 400.0, + "P2P case {}: loss {:.2} dB out of expected range", + i + 1, a_db + ); + } +} + +// area.csv test cases +#[test] +fn area_test_cases() { + // h_tx, h_rx, delta_h, mdvar, d_km, tx_sit, rx_sit, eps, sigma, n0, f_mhz, pol, climate, + // time, loc, sit → expected A__db (±1 dB tolerance) + struct AreaCase { + h_tx: f64, h_rx: f64, delta_h: f64, mdvar: i32, d_km: f64, + tx_sit: i32, rx_sit: i32, eps: f64, sigma: f64, n0: f64, + f_mhz: f64, pol: i32, climate: i32, + time: f64, loc: f64, sit: f64, expected: f64, + } + + let cases = [ + AreaCase { h_tx:10.0, h_rx:1.0, delta_h:0.0, mdvar:0, d_km:16.0, tx_sit:0, rx_sit:0, eps:15.0, sigma:0.008, n0:301.0, f_mhz:230.0, pol:0, climate:5, time:87.0, loc:50.0, sit:50.0, expected:152.5 }, + AreaCase { h_tx:3.0, h_rx:1.5, delta_h:10.0, mdvar:1, d_km:10.0, tx_sit:1, rx_sit:0, eps:15.0, sigma:0.008, n0:301.0, f_mhz:450.0, pol:0, climate:5, time:40.0, loc:28.0, sit:25.0, expected:133.0 }, + AreaCase { h_tx:15.0, h_rx:3.0, delta_h:5.0, mdvar:2, d_km:100.0, tx_sit:2, rx_sit:1, eps:15.0, sigma:0.008, n0:301.0, f_mhz:980.0, pol:1, climate:4, time:92.0, loc:53.0, sit:97.0, expected:224.1 }, + AreaCase { h_tx:3.0, h_rx:5.0, delta_h:20.0, mdvar:3, d_km:75.0, tx_sit:0, rx_sit:1, eps:15.0, sigma:0.008, n0:301.0, f_mhz:3100.0, pol:1, climate:2, time:50.0, loc:80.0, sit:43.0, expected:205.1 }, + AreaCase { h_tx:1.5, h_rx:10.0, delta_h:45.0, mdvar:0, d_km:25.0, tx_sit:1, rx_sit:2, eps:15.0, sigma:0.008, n0:301.0, f_mhz:8900.0, pol:1, climate:1, time:70.0, loc:99.0, sit:26.0, expected:156.0 }, + ]; + + for (i, c) in cases.iter().enumerate() { + let (rtn, a_db, _warnings) = itm_area_tls( + c.h_tx, c.h_rx, c.tx_sit, c.rx_sit, c.d_km, c.delta_h, + c.climate, c.n0, c.f_mhz, c.pol, c.eps, c.sigma, c.mdvar, + c.time, c.loc, c.sit, + ); + + assert!(rtn <= 1, "Area case {}: unexpected error {}", i + 1, rtn); + assert!( + (a_db - c.expected).abs() < 1.0, + "Area case {}: got {:.2} dB, expected {:.1} dB (diff {:.2})", + i + 1, a_db, c.expected, a_db - c.expected + ); + } +} diff --git a/itm-rs/itm/src/troposcatter_loss.rs b/itm-rs/itm/src/troposcatter_loss.rs new file mode 100644 index 0000000..d4bec29 --- /dev/null +++ b/itm-rs/itm/src/troposcatter_loss.rs @@ -0,0 +1,105 @@ +use crate::constants::SQRT2; +use crate::h0_function::h0_function; + +fn f_function(td: f64) -> f64 { + const A: [f64; 3] = [133.4, 104.6, 71.8]; + const B: [f64; 3] = [0.332e-3, 0.212e-3, 0.157e-3]; + const C: [f64; 3] = [-10.0, -2.5, 5.0]; + + let i = if td <= 10_000.0 { + 0 + } else if td <= 70_000.0 { + 1 + } else { + 2 + }; + + A[i] + B[i] * td + C[i] * td.log10() +} + +/// Troposcatter loss, in dB. +/// `h0` is an in/out parameter: pass `-1.0` on first call, and the stored value +/// on subsequent calls to allow the short-circuit for H_0 > 15 dB. +pub fn troposcatter_loss( + d__meter: f64, + theta_hzn: [f64; 2], + d_hzn__meter: [f64; 2], + h_e__meter: [f64; 2], + a_e__meter: f64, + n_s: f64, + f__mhz: f64, + theta_los: f64, + h0: &mut f64, +) -> f64 { + let wn = f__mhz / 47.7; + + let h_0 = if *h0 > 15.0 { + *h0 + } else { + let mut ad = d_hzn__meter[0] - d_hzn__meter[1]; + let mut rr = h_e__meter[1] / h_e__meter[0]; + if ad < 0.0 { + ad = -ad; + rr = 1.0 / rr; + } + + let theta = theta_hzn[0] + theta_hzn[1] + d__meter / a_e__meter; + + let r_1 = 2.0 * wn * theta * h_e__meter[0]; + let r_2 = 2.0 * wn * theta * h_e__meter[1]; + + if r_1 < 0.2 && r_2 < 0.2 { + return 1001.0; + } + + let s = ((d__meter - ad) / (d__meter + ad)).max(0.1); + let q = f64::min(f64::max(0.1, rr / s), 10.0); + + let h_0_meter = (d__meter - ad) * (d__meter + ad) * theta * 0.25 / d__meter; + + const Z_0: f64 = 1755.6; + const Z_1: f64 = 8000.0; + let eta_s = (h_0_meter / Z_0) + * (1.0 + + (0.031 - n_s * 2.32e-3 + n_s.powi(2) * 5.67e-6) + * (-f64::min(1.7, h_0_meter / Z_1).powi(6)).exp()); + + let h_00 = (h0_function(r_1, eta_s) + h0_function(r_2, eta_s)) / 2.0; + let delta_h_0 = f64::min( + h_00, + 6.0 * (0.6 - f64::max(eta_s, 1.0).log10()) * s.log10() * q.log10(), + ); + + let mut h_0_val = (h_00 + delta_h_0).max(0.0); + + if eta_s < 1.0 { + let special = 10.0 + * ((1.0 + SQRT2 / r_1) * (1.0 + SQRT2 / r_2)).powi(2) + .log10() + .mul_add( + 0.0, + ((1.0 + SQRT2 / r_1) * (1.0 + SQRT2 / r_2)).powi(2) + * (r_1 + r_2) + / (r_1 + r_2 + 2.0 * SQRT2), + ) + .log10(); + h_0_val = eta_s * h_0_val + (1.0 - eta_s) * special; + } + + if h_0_val > 15.0 && *h0 >= 0.0 { + h_0_val = *h0; + } + + h_0_val + }; + + *h0 = h_0; + let th = d__meter / a_e__meter - theta_los; + + const D_0: f64 = 40_000.0; + const H_M: f64 = 47.7; + f_function(th * d__meter) + + 10.0 * (wn * H_M * th.powi(4)).log10() + - 0.1 * (n_s - 301.0) * (-th * d__meter / D_0).exp() + + h_0 +} diff --git a/itm-rs/itm/src/types.rs b/itm-rs/itm/src/types.rs new file mode 100644 index 0000000..86b4774 --- /dev/null +++ b/itm-rs/itm/src/types.rs @@ -0,0 +1,22 @@ +/// Intermediate values computed during ITM propagation calculations. +#[derive(Clone, Debug, Default)] +pub struct IntermediateValues { + /// Terminal horizon angles, in radians [TX, RX] + pub theta_hzn: [f64; 2], + /// Terminal horizon distances, in meters [TX, RX] + pub d_hzn__meter: [f64; 2], + /// Terminal effective heights, in meters [TX, RX] + pub h_e__meter: [f64; 2], + /// Surface refractivity, in N-Units + pub n_s: f64, + /// Terrain irregularity parameter, in meters + pub delta_h__meter: f64, + /// Reference attenuation, in dB + pub a_ref__db: f64, + /// Free-space basic transmission loss, in dB + pub a_fs__db: f64, + /// Path distance, in km + pub d__km: f64, + /// Mode of propagation (0=not set, 1=LOS, 2=diffraction, 3=troposcatter) + pub mode: i32, +} diff --git a/itm-rs/itm/src/validate_inputs.rs b/itm-rs/itm/src/validate_inputs.rs new file mode 100644 index 0000000..04414a8 --- /dev/null +++ b/itm-rs/itm/src/validate_inputs.rs @@ -0,0 +1,89 @@ +use crate::errors::*; +use crate::warnings::*; +use crate::constants::*; + +/// Validate inputs common to both P2P and Area modes. +/// Returns an error code or SUCCESS; also sets warning bits. +pub fn validate_inputs( + h_tx__meter: f64, + h_rx__meter: f64, + climate: i32, + time: f64, + location: f64, + situation: f64, + n_0: f64, + f__mhz: f64, + pol: i32, + epsilon: f64, + sigma: f64, + mdvar: i32, + warnings: &mut i32, +) -> i32 { + if h_tx__meter < 1.0 || h_tx__meter > 1000.0 { + *warnings |= WARN__TX_TERMINAL_HEIGHT; + } + if h_tx__meter < 0.5 || h_tx__meter > 3000.0 { + return ERROR__TX_TERMINAL_HEIGHT; + } + + if h_rx__meter < 1.0 || h_rx__meter > 1000.0 { + *warnings |= WARN__RX_TERMINAL_HEIGHT; + } + if h_rx__meter < 0.5 || h_rx__meter > 3000.0 { + return ERROR__RX_TERMINAL_HEIGHT; + } + + if climate != CLIMATE__EQUATORIAL + && climate != CLIMATE__CONTINENTAL_SUBTROPICAL + && climate != CLIMATE__MARITIME_SUBTROPICAL + && climate != CLIMATE__DESERT + && climate != CLIMATE__CONTINENTAL_TEMPERATE + && climate != CLIMATE__MARITIME_TEMPERATE_OVER_LAND + && climate != CLIMATE__MARITIME_TEMPERATE_OVER_SEA + { + return ERROR__INVALID_RADIO_CLIMATE; + } + + if n_0 < 250.0 || n_0 > 400.0 { + return ERROR__REFRACTIVITY; + } + + if f__mhz < 40.0 || f__mhz > 10_000.0 { + *warnings |= WARN__FREQUENCY; + } + if f__mhz < 20.0 || f__mhz > 20_000.0 { + return ERROR__FREQUENCY; + } + + if pol != POLARIZATION__HORIZONTAL && pol != POLARIZATION__VERTICAL { + return ERROR__POLARIZATION; + } + + if epsilon < 1.0 { + return ERROR__EPSILON; + } + if sigma <= 0.0 { + return ERROR__SIGMA; + } + + if mdvar < 0 + || (mdvar > 3 && mdvar < 10) + || (mdvar > 13 && mdvar < 20) + || (mdvar > 23 && mdvar < 30) + || mdvar > 33 + { + return ERROR__MDVAR; + } + + if situation <= 0.0 || situation >= 100.0 { + return ERROR__INVALID_SITUATION; + } + if time <= 0.0 || time >= 100.0 { + return ERROR__INVALID_TIME; + } + if location <= 0.0 || location >= 100.0 { + return ERROR__INVALID_LOCATION; + } + + SUCCESS +} diff --git a/itm-rs/itm/src/variability.rs b/itm-rs/itm/src/variability.rs new file mode 100644 index 0000000..c5a1916 --- /dev/null +++ b/itm-rs/itm/src/variability.rs @@ -0,0 +1,175 @@ +use crate::constants::*; +use crate::iccdf::iccdf; +use crate::terrain_roughness::terrain_roughness; +use crate::warnings::*; + +fn curve(c1: f64, c2: f64, x1: f64, x2: f64, x3: f64, d_e: f64) -> f64 { + (c1 + c2 / (1.0 + ((d_e - x2) / x3).powi(2))) + * (d_e / x1).powi(2) + / (1.0 + (d_e / x1).powi(2)) +} + +/// Compute variability loss, in dB. +pub fn variability( + time: f64, + location: f64, + situation: f64, + h_e__meter: [f64; 2], + delta_h__meter: f64, + f__mhz: f64, + d__meter: f64, + a_ref__db: f64, + climate: i32, + mdvar: i32, + warnings: &mut i32, +) -> f64 { + #[rustfmt::skip] + const ALL_YEAR: [[f64; 7]; 5] = [ + [ -9.67, -0.62, 1.26, -9.21, -0.62, -0.39, 3.15 ], + [ 12.7, 9.19, 15.5, 9.05, 9.19, 2.86, 857.9 ], + [ 144.9e3, 228.9e3, 262.6e3, 84.1e3, 228.9e3, 141.7e3, 2222.0e3 ], + [ 190.3e3, 205.2e3, 185.2e3, 101.1e3, 205.2e3, 315.9e3, 164.8e3 ], + [ 133.8e3, 143.6e3, 99.8e3, 98.6e3, 143.6e3, 167.4e3, 116.3e3 ], + ]; + + const BSM1: [f64; 7] = [2.13, 2.66, 6.11, 1.98, 2.68, 6.86, 8.51]; + const BSM2: [f64; 7] = [159.5, 7.67, 6.65, 13.11, 7.16, 10.38, 169.8]; + const XSM1: [f64; 7] = [762.2e3, 100.4e3, 138.2e3, 139.1e3, 93.7e3, 187.8e3, 609.8e3]; + const XSM2: [f64; 7] = [123.6e3, 172.5e3, 242.2e3, 132.7e3, 186.8e3, 169.6e3, 119.9e3]; + const XSM3: [f64; 7] = [94.5e3, 136.4e3, 178.6e3, 193.5e3, 133.5e3, 108.9e3, 106.6e3]; + + const BSP1: [f64; 7] = [2.11, 6.87, 10.08, 3.68, 4.75, 8.58, 8.43]; + const BSP2: [f64; 7] = [102.3, 15.53, 9.60, 159.3, 8.12, 13.97, 8.19]; + const XSP1: [f64; 7] = [636.9e3, 138.7e3, 165.3e3, 464.4e3, 93.2e3, 216.0e3, 136.2e3]; + const XSP2: [f64; 7] = [134.8e3, 143.7e3, 225.7e3, 93.1e3, 135.9e3, 152.0e3, 188.5e3]; + const XSP3: [f64; 7] = [95.6e3, 98.6e3, 129.7e3, 94.2e3, 113.4e3, 122.7e3, 122.9e3]; + + const C_D: [f64; 7] = [1.224, 0.801, 1.380, 1.000, 1.224, 1.518, 1.518]; + const Z_D: [f64; 7] = [1.282, 2.161, 1.282, 20.0, 1.282, 1.282, 1.282]; + + const BFM1: [f64; 7] = [1.0, 1.0, 1.0, 1.0, 0.92, 1.0, 1.0]; + const BFM2: [f64; 7] = [0.0, 0.0, 0.0, 0.0, 0.25, 0.0, 0.0]; + const BFM3: [f64; 7] = [0.0, 0.0, 0.0, 0.0, 1.77, 0.0, 0.0]; + + const BFP1: [f64; 7] = [1.0, 0.93, 1.0, 0.93, 0.93, 1.0, 1.0]; + const BFP2: [f64; 7] = [0.0, 0.31, 0.0, 0.19, 0.31, 0.0, 0.0]; + const BFP3: [f64; 7] = [0.0, 2.00, 0.0, 1.79, 2.00, 0.0, 0.0]; + + let mut z_t = iccdf(time / 100.0); + let mut z_l = iccdf(location / 100.0); + let z_s = iccdf(situation / 100.0); + + let ci = (climate - 1) as usize; + + let wn = f__mhz / 47.7; + + let d_ex = (2.0 * A_9000__METER * h_e__meter[0]).sqrt() + + (2.0 * A_9000__METER * h_e__meter[1]).sqrt() + + (575.7e12 / wn).powf(THIRD); + + let d_e = if d__meter < d_ex { + 130_000.0 * d__meter / d_ex + } else { + 130_000.0 + d__meter - d_ex + }; + + // Situation variability + let mut mdvar_internal = mdvar; + let plus20 = mdvar_internal >= 20; + if plus20 { + mdvar_internal -= 20; + } + + let sigma_s = if plus20 { + 0.0 + } else { + 5.0 + 3.0 * (-d_e / 100_000.0).exp() + }; + + let plus10 = mdvar_internal >= 10; + if plus10 { + mdvar_internal -= 10; + } + + let v_med = curve( + ALL_YEAR[0][ci], + ALL_YEAR[1][ci], + ALL_YEAR[2][ci], + ALL_YEAR[3][ci], + ALL_YEAR[4][ci], + d_e, + ); + + if mdvar_internal == SINGLE_MESSAGE_MODE { + z_t = z_s; + z_l = z_s; + } else if mdvar_internal == ACCIDENTAL_MODE { + z_l = z_s; + } else if mdvar_internal == MOBILE_MODE { + z_l = z_t; + } + + if z_t.abs() > 3.10 || z_l.abs() > 3.10 || z_s.abs() > 3.10 { + *warnings |= WARN__EXTREME_VARIABILITIES; + } + + // Location variability + let sigma_l = if plus10 { + 0.0 + } else { + let delta_h_d = terrain_roughness(d__meter, delta_h__meter); + 10.0 * wn * delta_h_d / (wn * delta_h_d + 13.0) + }; + let y_l = sigma_l * z_l; + + // Time variability + let q = (0.133 * wn).ln(); + let g_minus = BFM1[ci] + BFM2[ci] / ((BFM3[ci] * q).powi(2) + 1.0); + let g_plus = BFP1[ci] + BFP2[ci] / ((BFP3[ci] * q).powi(2) + 1.0); + + let sigma_t_minus = + curve(BSM1[ci], BSM2[ci], XSM1[ci], XSM2[ci], XSM3[ci], d_e) * g_minus; + let sigma_t_plus = + curve(BSP1[ci], BSP2[ci], XSP1[ci], XSP2[ci], XSP3[ci], d_e) * g_plus; + + let sigma_td = C_D[ci] * sigma_t_plus; + let tgtd = (sigma_t_plus - sigma_td) * Z_D[ci]; + + let sigma_t = if z_t < 0.0 { + sigma_t_minus + } else if z_t <= Z_D[ci] { + sigma_t_plus + } else { + sigma_td + tgtd / z_t + }; + let y_t = sigma_t * z_t; + + let y_s_temp = sigma_s.powi(2) + + y_t.powi(2) / (7.8 + z_s.powi(2)) + + y_l.powi(2) / (24.0 + z_s.powi(2)); + + let (y_r, y_s) = if mdvar_internal == SINGLE_MESSAGE_MODE { + ( + 0.0, + (sigma_t.powi(2) + sigma_l.powi(2) + y_s_temp).sqrt() * z_s, + ) + } else if mdvar_internal == ACCIDENTAL_MODE { + (y_t, (sigma_l.powi(2) + y_s_temp).sqrt() * z_s) + } else if mdvar_internal == MOBILE_MODE { + ( + (sigma_t.powi(2) + sigma_l.powi(2)).sqrt() * z_t, + y_s_temp.sqrt() * z_s, + ) + } else { + // BROADCAST_MODE + (y_t + y_l, y_s_temp.sqrt() * z_s) + }; + + let result = a_ref__db - v_med - y_r - y_s; + + if result < 0.0 { + result * (29.0 - result) / (29.0 - 10.0 * result) + } else { + result + } +} diff --git a/itm-rs/itm/src/warnings.rs b/itm-rs/itm/src/warnings.rs new file mode 100644 index 0000000..38ba069 --- /dev/null +++ b/itm-rs/itm/src/warnings.rs @@ -0,0 +1,15 @@ +pub const WARN__TX_TERMINAL_HEIGHT: i32 = 0x0001; +pub const WARN__RX_TERMINAL_HEIGHT: i32 = 0x0002; +pub const WARN__FREQUENCY: i32 = 0x0004; +pub const WARN__PATH_DISTANCE_TOO_BIG_1: i32 = 0x0008; +pub const WARN__PATH_DISTANCE_TOO_BIG_2: i32 = 0x0010; +pub const WARN__PATH_DISTANCE_TOO_SMALL_1: i32 = 0x0020; +pub const WARN__PATH_DISTANCE_TOO_SMALL_2: i32 = 0x0040; +pub const WARN__TX_HORIZON_ANGLE: i32 = 0x0080; +pub const WARN__RX_HORIZON_ANGLE: i32 = 0x0100; +pub const WARN__TX_HORIZON_DISTANCE_1: i32 = 0x0200; +pub const WARN__RX_HORIZON_DISTANCE_1: i32 = 0x0400; +pub const WARN__TX_HORIZON_DISTANCE_2: i32 = 0x0800; +pub const WARN__RX_HORIZON_DISTANCE_2: i32 = 0x1000; +pub const WARN__EXTREME_VARIABILITIES: i32 = 0x2000; +pub const WARN__SURFACE_REFRACTIVITY: i32 = 0x4000;