diff --git a/.github/workflows/basic.yml b/.github/workflows/basic.yml new file mode 100644 index 0000000..af032fd --- /dev/null +++ b/.github/workflows/basic.yml @@ -0,0 +1,65 @@ +name: Basic Validation + +on: + push: + pull_request: + workflow_dispatch: + schedule: + - cron: '26 3 20 * *' + +env: + CARGO_TERM_COLOR: always + RUSTFLAGS: "-Dwarnings" + +jobs: + build-and-test: + + name: Build and Test + + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v4 + - name: Prepare environment + run: | + rustup update stable + cargo clean + - name: Test with cargo + run: | + cargo test + cargo test --no-default-features + cargo test --features "double_precision" + cargo test --features "argdebug" + cargo test --features "array" + cargo test --features "parallel" + cargo test --features "double_precision, argdebug" + cargo test --features "double_precision, array" + cargo test --features "double_precision, parallel" + cargo test --features "argdebug, array" + cargo test --features "argdebug, parallel" + cargo test --features "array, parallel" + cargo test --features "double_precision, argdebug, array" + cargo test --features "double_precision, argdebug, parallel" + cargo test --features "double_precision, array, parallel" + cargo test --features "argdebug, array, parallel" + cargo test --all-features + cargo clean + - name: Check release builds + run: | + cargo build --release + cargo build --release --no-default-features + cargo build --release --features "double_precision" + cargo build --release --features "argdebug" + cargo build --release --features "array" + cargo build --release --features "parallel" + cargo build --release --features "double_precision, argdebug" + cargo build --release --features "double_precision, array" + cargo build --release --features "double_precision, parallel" + cargo build --release --features "argdebug, array" + cargo build --release --features "argdebug, parallel" + cargo build --release --features "array, parallel" + cargo build --release --features "double_precision, argdebug, array" + cargo build --release --features "double_precision, argdebug, parallel" + cargo build --release --features "double_precision, array, parallel" + cargo build --release --features "argdebug, array, parallel" + cargo build --release --all-features diff --git a/.github/workflows/clippy-check.yml b/.github/workflows/clippy-check.yml new file mode 100644 index 0000000..78dddbf --- /dev/null +++ b/.github/workflows/clippy-check.yml @@ -0,0 +1,19 @@ +on: pull_request + +name: Cargo Clippy + +env: + CARGO_TERM_COLOR: always + RUSTFLAGS: "-Dwarnings" + +jobs: + clippy_check: + name: Check + + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + - name: Clippy + run: | + rustup component add clippy + cargo clippy -- -D warnings diff --git a/.github/workflows/doc-coverage.yml b/.github/workflows/doc-coverage.yml new file mode 100644 index 0000000..6ed7ab3 --- /dev/null +++ b/.github/workflows/doc-coverage.yml @@ -0,0 +1,45 @@ +name: Cargo Doc + +on: + pull_request + +jobs: + pr-comment: + name: Check Coverage + runs-on: ubuntu-latest + steps: + - name: Checkout + uses: actions/checkout@v4 + - name: Install rust toolchain + uses: actions-rs/toolchain@v1 + with: + toolchain: nightly + - name: Fetch base + run: git fetch origin ${{ github.event.pull_request.base.sha }} + - name: Checkout base + run: git checkout ${{ github.event.pull_request.base.sha }} + - name: Calculate base doc coverage + uses: bewee/rustdoc-coverage-action@v1 + - name: Fetch head + run: git fetch origin ${{ github.event.pull_request.head.sha }} + - name: Checkout head + run: git checkout ${{ github.event.pull_request.head.sha }} + - name: Calculate doc coverage + id: coverage + uses: bewee/rustdoc-coverage-action@v1 + - name: Find Comment + uses: peter-evans/find-comment@v1 + id: fc + with: + issue-number: ${{ github.event.pull_request.number }} + comment-author: "github-actions[bot]" + body-includes: "## Documentation Coverage:" + - name: Create or update comment + uses: peter-evans/create-or-update-comment@v1 + with: + comment-id: ${{ steps.fc.outputs.comment-id }} + issue-number: ${{ github.event.pull_request.number }} + body: | + ## Documentation Coverage: + ${{ steps.coverage.outputs.table }} + edit-mode: replace \ No newline at end of file diff --git a/.github/workflows/rust.yml b/.github/workflows/rust.yml deleted file mode 100644 index 989d3fc..0000000 --- a/.github/workflows/rust.yml +++ /dev/null @@ -1,46 +0,0 @@ -name: cargo - -on: - push: - branches: [ main ] - pull_request: - branches: [ main ] - workflow_dispatch: - -env: - CARGO_TERM_COLOR: always - -jobs: - build: - - runs-on: ubuntu-latest - - steps: - - uses: actions/checkout@v2 - - name: Prepare environment - run: | - sudo apt-get update - sudo apt-get install clang - sudo apt-get install libclang1 - sudo apt-get install gnuplot - rustup update stable - cargo install cargo-criterion - cargo clean - - name: Build with cargo - run: | - cargo build --release - cargo build --release --features "double_precision" - cargo clean - - name: Check with clippy - run: | - cargo clippy -- -W clippy::pedantic - cargo clean - - name: Test with cargo - run: | - cargo test - cargo test --features "double_precision" - cargo clean - - name: Benchmark with criterion - run: | - cargo criterion - cargo criterion --features "double_precision" diff --git a/Cargo.toml b/Cargo.toml index 80b73e4..032334c 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,53 +1,69 @@ [package] name = "floccus" -version = "0.3.7" +version = "0.4.0-alpha.9" authors = ["Jakub Lewandowski "] -edition = "2021" +edition = "2024" description = "Formulae for air thermodynamic calculations" repository = "https://github.com/ScaleWeather/floccus" readme = "README.md" keywords = ["oceanography", "thermodynamics", "meteorology", "weather"] categories = ["mathematics", "science"] license = "Apache-2.0" -exclude = [ - ".github/*", -] +exclude = [".github/*"] [dependencies] -thiserror = "^1.0.30" -float-cmp = "^0.9.0" -floccus-proc = {version = "0.2.5", optional = true} -log = "^0.4.14" +thiserror = "2.0" +float-cmp = "0.10" +tracing = { version = "0.1", optional = true, default-features = false, features = [ + "attributes", +] } +ndarray = { version = "0.16", default-features = false, optional = true } +rayon = { version = "1.10", default-features = false, optional = true } +uom = { version = "0.37", default-features = false, features = [ + "si", + "std", + "autoconvert", + "f32", +] } + +[lib] +doctest = false [dev-dependencies] -criterion = "0.4.0" +criterion = "0.6" +itertools = "0.14" +log = "0.4" +testing_logger = "0.1" +floccus = { path = ".", features = ["double_precision"] } [features] -debug = ["floccus-proc"] -double_precision = [] +default = [] +argdebug = ["dep:tracing"] +double_precision = ["uom/f64"] +array = ["dep:ndarray"] +parallel = ["array", "ndarray/rayon", "dep:rayon"] [[bench]] -name = "virtual_temperature" +name = "equivalent_potential_temperature" harness = false [[bench]] -name = "vapour_pressure" +name = "mixing_ratio" harness = false [[bench]] -name = "vapour_pressure_deficit" -harness = false +name = "potential_temperature" [[bench]] -name = "mixing_ratio" +name = "relative_humidity" harness = false [[bench]] -name = "wet_bulb_temperature" +name = "saturation_mixing_ratio" harness = false [[bench]] -name = "relative_humidity" +name = "saturation_vapour_pressure" harness = false [[bench]] @@ -55,13 +71,21 @@ name = "specific_humidity" harness = false [[bench]] -name = "potential_temperature" +name = "vapour_pressure_deficit" harness = false [[bench]] -name = "equivalent_potential_temperature" +name = "vapour_pressure" +harness = false + +[[bench]] +name = "virtual_temperature" harness = false [[bench]] name = "wet_bulb_potential_temperature" harness = false + +[[bench]] +name = "wet_bulb_temperature" +harness = false diff --git a/README.md b/README.md index 7282a4e..a9c1017 100644 --- a/README.md +++ b/README.md @@ -1,9 +1,9 @@ # floccus -[![License](https://img.shields.io/github/license/ScaleWeather/floccus)](https://choosealicense.com/licenses/apache-2.0/) -[![Crates.io](https://img.shields.io/crates/v/floccus)](https://crates.io/crates/floccus) -[![dependency status](https://deps.rs/repo/github/ScaleWeather/floccus/status.svg)](https://deps.rs/repo/github/ScaleWeather/floccus) -[![GitHub Workflow Status](https://img.shields.io/github/actions/workflow/status/ScaleWeather/floccus/rust.yml?branch=main&label=cargo%20build)](https://github.com/ScaleWeather/floccus/actions) +[![Github Repository](https://img.shields.io/badge/Github-Repository-blue?style=flat-square&logo=github&color=blue)](https://github.com/ScaleWeather/floccus) +[![Crates.io](https://img.shields.io/crates/v/floccus?style=flat-square)](https://crates.io/crates/floccus) +[![License](https://img.shields.io/github/license/ScaleWeather/floccus?style=flat-square)](https://choosealicense.com/licenses/apache-2.0/) +[![dependency status](https://deps.rs/repo/github/ScaleWeather/floccus/status.svg?style=flat-square)](https://deps.rs/repo/github/ScaleWeather/floccus) Rust crate providing formulae for air thermodynamic calculations. @@ -58,6 +58,9 @@ To prevent any unexpected behavior, all functions check whether provided inputs Exact limits are specified in the documentation of each function. If the input is out of range the function will return an `InputError::OutOfRange` with erroneous input specified. +Each function also has `_unchecked` and `_validate` versions. The `_validate` version only checks the inputs with bounds defined for its "parent" function. +The `_unchecked` version performs only the calculation without any input checking. All "parent" functions simply call `_validate` and then `_unchecked`. + ## Debugging If additional information is needed about which function returns the error and why, `debug` feature can be enabled. @@ -66,6 +69,4 @@ information about the error. This feature potentially is not zero-cost so it is ## Benchmarks -Functions provided in this crate are intended for use in, i. a., numerical models. To provide the user information about performance overhead of each function all functions are benchmarked using [criterion.rs](https://bheisler.github.io/criterion.rs/book/index.html). Github Actions automatically runs all benchmarks. - -To check the latest benchmark results the newest workflow on [Github Actions page of floccus](https://github.com/ScaleWeather/floccus/actions). +Functions provided in this crate are intended for use in, i. a., numerical models. To provide the user information about performance overhead of each function all functions are can be benchmarked using [criterion.rs](https://bheisler.github.io/criterion.rs/book/index.html). diff --git a/benches/equivalent_potential_temperature.rs b/benches/equivalent_potential_temperature.rs index 419d3e4..97d740b 100644 --- a/benches/equivalent_potential_temperature.rs +++ b/benches/equivalent_potential_temperature.rs @@ -1,11 +1,59 @@ -use criterion::{Criterion, black_box, criterion_group, criterion_main}; -use floccus::equivalent_potential_temperature; +use criterion::{criterion_group, criterion_main, Criterion}; +use floccus::{formulas::equivalent_potential_temperature, formulas::Formula4}; -pub fn equivalent_potential_temperature_benchmark(c: &mut Criterion) { - c.bench_function("equivalent_potential_temperature::bryan1", |b| { - b.iter(|| equivalent_potential_temperature::bryan1(black_box(300.0), black_box(101325.0), black_box(3000.0))) +mod utils; +use utils::ReferenceValues; + +pub fn benchmark(c: &mut Criterion) { + let ref_norm = ReferenceValues::normal(); + + let mut group = c.benchmark_group("equivalent_potential_temperature"); + + group.bench_function("Bolton1", |b| { + b.iter(|| { + equivalent_potential_temperature::Bolton1::compute( + ref_norm.pres, + ref_norm.temp, + ref_norm.dwpt, + ref_norm.vapr, + ) + }) + }); + + group.bench_function("Bolton2", |b| { + b.iter(|| { + equivalent_potential_temperature::Bolton2::compute( + ref_norm.temp, + ref_norm.dwpt, + ref_norm.mxrt, + ref_norm.thet, + ) + }) + }); + + group.bench_function("Bryan1", |b| { + b.iter(|| { + equivalent_potential_temperature::Bryan1::compute( + ref_norm.temp, + ref_norm.mxrt, + ref_norm.rehu, + ref_norm.thet, + ) + }) + }); + + group.bench_function("Paluch1", |b| { + b.iter(|| { + equivalent_potential_temperature::Kerry1::compute( + ref_norm.temp, + ref_norm.pres, + ref_norm.mxrt, + ref_norm.rehu, + ) + }) }); + group.finish(); } -criterion_group!(benches, equivalent_potential_temperature_benchmark); -criterion_main!(benches); \ No newline at end of file +criterion_group!(benches, benchmark); +criterion_main!(benches); diff --git a/benches/mixing_ratio.rs b/benches/mixing_ratio.rs index f8f81a8..3a6cfbd 100644 --- a/benches/mixing_ratio.rs +++ b/benches/mixing_ratio.rs @@ -1,17 +1,19 @@ -use criterion::{Criterion, black_box, criterion_group, criterion_main}; -use floccus::mixing_ratio; +use criterion::{criterion_group, criterion_main, Criterion}; +use floccus::{formulas::mixing_ratio, formulas::Formula2}; -pub fn mixing_ratio_benchmark(c: &mut Criterion) { - c.bench_function("mixing_ratio::general1", |b| { - b.iter(|| mixing_ratio::general1(black_box(101325.0), black_box(3500.0))) - }); - c.bench_function("mixing_ratio::performance1", |b| { - b.iter(|| mixing_ratio::performance1(black_box(300.0), black_box(101325.0))) - }); - c.bench_function("mixing_ratio::accuracy1", |b| { - b.iter(|| mixing_ratio::accuracy1(black_box(300.0), black_box(101325.0))) +mod utils; +use utils::ReferenceValues; + +pub fn benchmark(c: &mut Criterion) { + let ref_norm = ReferenceValues::normal(); + + let mut group = c.benchmark_group("mixing_ratio"); + + group.bench_function("definition1", |b| { + b.iter(|| mixing_ratio::Definition1::compute(ref_norm.pres, ref_norm.vapr)) }); + group.finish(); } -criterion_group!(benches, mixing_ratio_benchmark); +criterion_group!(benches, benchmark); criterion_main!(benches); diff --git a/benches/potential_temperature.rs b/benches/potential_temperature.rs index 667ce1d..4b08b91 100644 --- a/benches/potential_temperature.rs +++ b/benches/potential_temperature.rs @@ -1,11 +1,21 @@ -use criterion::{Criterion, black_box, criterion_group, criterion_main}; -use floccus::potential_temperature; +use criterion::{criterion_group, criterion_main, Criterion}; +use floccus::{formulas::potential_temperature, formulas::Formula3}; -pub fn potential_temperature_benchmark(c: &mut Criterion) { - c.bench_function("potential_temperature::davies_jones1", |b| { - b.iter(|| potential_temperature::davies_jones1(black_box(300.0), black_box(101325.0), black_box(3000.0))) +mod utils; +use utils::ReferenceValues; + +pub fn benchmark(c: &mut Criterion) { + let ref_norm = ReferenceValues::normal(); + + let mut group = c.benchmark_group("potential_temperature"); + + group.bench_function("definition1", |b| { + b.iter(|| { + potential_temperature::Definition1::compute(ref_norm.temp, ref_norm.pres, ref_norm.vapr) + }) }); + group.finish(); } -criterion_group!(benches, potential_temperature_benchmark); -criterion_main!(benches); \ No newline at end of file +criterion_group!(benches, benchmark); +criterion_main!(benches); diff --git a/benches/relative_humidity.rs b/benches/relative_humidity.rs index 5dfcee4..5cee30a 100644 --- a/benches/relative_humidity.rs +++ b/benches/relative_humidity.rs @@ -1,27 +1,23 @@ -use criterion::{black_box, criterion_group, criterion_main, Criterion}; -use floccus::relative_humidity; +use criterion::{criterion_group, criterion_main, Criterion}; +use floccus::{formulas::relative_humidity, formulas::Formula2}; -pub fn relative_humidity_benchmark(c: &mut Criterion) { - c.bench_function("relative_humidity::general1", |b| { - b.iter(|| relative_humidity::general1(black_box(0.01064), black_box(0.01467))) - }); +mod utils; +use utils::ReferenceValues; - c.bench_function("relative_humidity::general2", |b| { - b.iter(|| relative_humidity::general2(black_box(1706.0), black_box(2339.0))) - }); +pub fn benchmark(c: &mut Criterion) { + let ref_norm = ReferenceValues::normal(); - c.bench_function("relative_humidity::general3", |b| { - b.iter(|| relative_humidity::general3(black_box(300.0), black_box(290.0))) - }); + let mut group = c.benchmark_group("relative_humidity"); - c.bench_function("relative_humidity::general4", |b| { - b.iter(|| relative_humidity::general4(black_box(300.0), black_box(290.0), black_box(101325.0))) + group.bench_function("definition1", |b| { + b.iter(|| relative_humidity::Definition1::compute(ref_norm.mxrt, ref_norm.smrt)) }); - c.bench_function("relative_humidity::general5", |b| { - b.iter(|| relative_humidity::general5(black_box(300.0), black_box(290.0), black_box(101325.0))) + group.bench_function("definition2", |b| { + b.iter(|| relative_humidity::Definition2::compute(ref_norm.vapr, ref_norm.savp)) }); + group.finish(); } -criterion_group!(benches, relative_humidity_benchmark); +criterion_group!(benches, benchmark); criterion_main!(benches); diff --git a/benches/saturation_mixing_ratio.rs b/benches/saturation_mixing_ratio.rs new file mode 100644 index 0000000..2eacc58 --- /dev/null +++ b/benches/saturation_mixing_ratio.rs @@ -0,0 +1,24 @@ +use criterion::{criterion_group, criterion_main, Criterion}; +use floccus::{formulas::saturation_mixing_ratio, formulas::Formula2}; + +mod utils; +use utils::ReferenceValues; + +pub fn benchmark(c: &mut Criterion) { + let ref_norm = ReferenceValues::normal(); + + let mut group = c.benchmark_group("saturation_mixing_ratio"); + + group.bench_function("definition1", |b| { + b.iter(|| saturation_mixing_ratio::Definition1::compute(ref_norm.pres, ref_norm.savp)) + }); + + group.bench_function("definition2", |b| { + b.iter(|| saturation_mixing_ratio::Definition2::compute(ref_norm.mxrt, ref_norm.rehu)) + }); + + group.finish(); +} + +criterion_group!(benches, benchmark); +criterion_main!(benches); diff --git a/benches/saturation_vapour_pressure.rs b/benches/saturation_vapour_pressure.rs new file mode 100644 index 0000000..f0c29a6 --- /dev/null +++ b/benches/saturation_vapour_pressure.rs @@ -0,0 +1,57 @@ +use criterion::{criterion_group, criterion_main, Criterion}; +use floccus::{formulas::saturation_vapour_pressure, formulas::{Formula1, Formula2}}; + +mod utils; +use utils::ReferenceValues; + +pub fn benchmark(c: &mut Criterion) { + let ref_norm = ReferenceValues::normal(); + let ref_freeze = ReferenceValues::freeze(); + + let mut group = c.benchmark_group("saturation_vapour_pressure"); + + group.bench_function("definition1", |b| { + b.iter(|| saturation_vapour_pressure::Definition1::compute(ref_norm.vapr, ref_norm.rehu)) + }); + + group.bench_function("tetens1", |b| { + b.iter(|| saturation_vapour_pressure::Tetens1::compute(ref_norm.temp)) + }); + + group.bench_function("buck1", |b| { + b.iter(|| saturation_vapour_pressure::Buck1::compute(ref_norm.temp, ref_norm.pres)) + }); + + group.bench_function("buck2", |b| { + b.iter(|| saturation_vapour_pressure::Buck2::compute(ref_freeze.temp, ref_freeze.pres)) + }); + + group.bench_function("buck3", |b| { + b.iter(|| saturation_vapour_pressure::Buck3::compute(ref_norm.temp, ref_norm.pres)) + }); + + group.bench_function("buck4", |b| { + b.iter(|| saturation_vapour_pressure::Buck4::compute(ref_freeze.temp, ref_freeze.pres)) + }); + + group.bench_function("buck3_simplified", |b| { + b.iter(|| saturation_vapour_pressure::Buck3Simplified::compute(ref_norm.temp)) + }); + + group.bench_function("buck4_simplified", |b| { + b.iter(|| saturation_vapour_pressure::Buck4Simplified::compute(ref_freeze.temp)) + }); + + group.bench_function("wexler1", |b| { + b.iter(|| saturation_vapour_pressure::Wexler1::compute(ref_norm.temp)) + }); + + group.bench_function("wexler2", |b| { + b.iter(|| saturation_vapour_pressure::Wexler2::compute(ref_freeze.temp)) + }); + + group.finish(); +} + +criterion_group!(benches, benchmark); +criterion_main!(benches); diff --git a/benches/specific_humidity.rs b/benches/specific_humidity.rs index fd90abc..78815e7 100644 --- a/benches/specific_humidity.rs +++ b/benches/specific_humidity.rs @@ -1,11 +1,20 @@ -use criterion::{Criterion, black_box, criterion_group, criterion_main}; -use floccus::specific_humidity; +use criterion::{criterion_group, criterion_main, Criterion}; +use floccus::{formulas::specific_humidity, formulas::Formula2}; -pub fn specific_humidity_benchmark(c: &mut Criterion) { - c.bench_function("specific_humidity::general1", |b| { - b.iter(|| specific_humidity::general1(black_box(3000.0), black_box(101325.0))) +mod utils; +use utils::ReferenceValues; + +pub fn benchmark(c: &mut Criterion) { + let ref_norm = ReferenceValues::normal(); + + let mut group = c.benchmark_group("specific_humidity"); + + group.bench_function("definition1", |b| { + b.iter(|| specific_humidity::Definition1::compute(ref_norm.vapr, ref_norm.pres)) }); + + group.finish(); } -criterion_group!(benches, specific_humidity_benchmark); -criterion_main!(benches); \ No newline at end of file +criterion_group!(benches, benchmark); +criterion_main!(benches); diff --git a/benches/utils.rs b/benches/utils.rs new file mode 100644 index 0000000..0cef4c6 --- /dev/null +++ b/benches/utils.rs @@ -0,0 +1,83 @@ +#![allow(unused)] + +use std::hint::black_box; +use floccus::quantities::{ + AtmosphericPressure, DewPointTemperature, DryBulbTemperature, EquivalentPotentialTemperature, + MixingRatio, PotentialTemperature, RelativeHumidity, SaturationMixingRatio, + SaturationVapourPressure, SpecificHumidity, VapourPressure, +}; +use uom::si::{pressure::pascal, ratio::ratio, thermodynamic_temperature::kelvin}; + +type Float = f64; + +pub(crate) const TEMP_NORM: Float = 300.0; +pub(crate) const DWPT_NORM: Float = 290.0; +pub(crate) const PRES_NORM: Float = 100000.0; +pub(crate) const VP_NORM: Float = 1919.4253257541593; +pub(crate) const SVP_NORM: Float = 3535.4235919263083; +pub(crate) const RH_NORM: Float = 0.5429124052171476; +pub(crate) const MR_NORM: Float = 0.012172079452423202; +pub(crate) const SMR_NROM: Float = 0.022419969290542845; +pub(crate) const SH_NORM: Float = 0.012025701656390478; +pub(crate) const THETAE_NORM: Float = 331.329289539998; +pub(crate) const THETA_NORM: Float = 301.66581400702955; + +pub(crate) const TEMP_FREEZ: Float = 260.0; +pub(crate) const DWPT_FREEZ: Float = 255.0; +pub(crate) const PRES_FREEZ: Float = 100000.0; +pub(crate) const VP_FREEZ: Float = 123.17937690212507; +pub(crate) const SVP_FREEZ: Float = 195.84980045970696; +pub(crate) const RH_FREEZ: Float = 0.6289481868911442; +pub(crate) const MR_FREEZ: Float = 0.0007670962389744638; +pub(crate) const SMR_FREEZ: Float = 0.0012196493367222787; +pub(crate) const SH_FREEZ: Float = 0.000766508253376156; +pub(crate) const THETAE_FREEZ: Float = 261.95287841149707; +pub(crate) const THETA_FREEZ: Float = 260.0915766593588; + +pub struct ReferenceValues { + pub temp: DryBulbTemperature, + pub pres: AtmosphericPressure, + pub dwpt: DewPointTemperature, + pub sphu: SpecificHumidity, + pub vapr: VapourPressure, + pub savp: SaturationVapourPressure, + pub rehu: RelativeHumidity, + pub mxrt: MixingRatio, + pub smrt: SaturationMixingRatio, + pub thte: EquivalentPotentialTemperature, + pub thet: PotentialTemperature, +} + +impl ReferenceValues { + pub fn normal() -> Self { + Self { + temp: black_box(DryBulbTemperature::new::(TEMP_NORM)), + pres: black_box(AtmosphericPressure::new::(PRES_NORM)), + dwpt: black_box(DewPointTemperature::new::(DWPT_NORM)), + sphu: black_box(SpecificHumidity::new::(SH_NORM)), + vapr: black_box(VapourPressure::new::(VP_NORM)), + savp: black_box(SaturationVapourPressure::new::(SVP_NORM)), + rehu: black_box(RelativeHumidity::new::(RH_NORM)), + mxrt: black_box(MixingRatio::new::(MR_NORM)), + smrt: black_box(SaturationMixingRatio::new::(SMR_NROM)), + thte: black_box(EquivalentPotentialTemperature::new::(THETAE_NORM)), + thet: black_box(PotentialTemperature::new::(THETA_NORM)), + } + } + + pub fn freeze() -> Self { + Self { + temp: black_box(DryBulbTemperature::new::(TEMP_FREEZ)), + pres: black_box(AtmosphericPressure::new::(PRES_FREEZ)), + dwpt: black_box(DewPointTemperature::new::(DWPT_FREEZ)), + sphu: black_box(SpecificHumidity::new::(SH_FREEZ)), + vapr: black_box(VapourPressure::new::(VP_FREEZ)), + savp: black_box(SaturationVapourPressure::new::(SVP_FREEZ)), + rehu: black_box(RelativeHumidity::new::(RH_FREEZ)), + mxrt: black_box(MixingRatio::new::(MR_FREEZ)), + smrt: black_box(SaturationMixingRatio::new::(SMR_FREEZ)), + thte: black_box(EquivalentPotentialTemperature::new::(THETAE_FREEZ)), + thet: black_box(PotentialTemperature::new::(THETA_FREEZ)), + } + } +} diff --git a/benches/vapour_pressure.rs b/benches/vapour_pressure.rs index 026530f..297913a 100644 --- a/benches/vapour_pressure.rs +++ b/benches/vapour_pressure.rs @@ -1,55 +1,60 @@ -use criterion::{black_box, criterion_group, criterion_main, Criterion}; -use floccus::vapour_pressure; +use criterion::{criterion_group, criterion_main, Criterion}; +use floccus::{formulas::vapour_pressure, formulas::{Formula1, Formula2}}; -pub fn vapour_pressure_benchmark(c: &mut Criterion) { - c.bench_function("vapour_pressure::general1", |b| { - b.iter(|| vapour_pressure::general1(black_box(0.022), black_box(101325.0))) - }); +mod utils; +use utils::ReferenceValues; + +pub fn benchmark(c: &mut Criterion) { + let ref_norm = ReferenceValues::normal(); + let ref_freeze = ReferenceValues::freeze(); + + let mut group = c.benchmark_group("vapour_pressure"); - c.bench_function("vapour_pressure::tetens1", |b| { - b.iter(|| vapour_pressure::tetens1(black_box(300.0))) + group.bench_function("definition1", |b| { + b.iter(|| vapour_pressure::Definition1::compute(ref_norm.sphu, ref_norm.pres)) }); - c.bench_function("vapour_pressure::buck1", |b| { - b.iter(|| vapour_pressure::buck1(black_box(300.0), black_box(101325.0))) + group.bench_function("definition2", |b| { + b.iter(|| vapour_pressure::Definition2::compute(ref_norm.savp, ref_norm.rehu)) }); - c.bench_function("vapour_pressure::buck2", |b| { - b.iter(|| vapour_pressure::buck2(black_box(250.0), black_box(101325.0))) + group.bench_function("tetens1", |b| { + b.iter(|| vapour_pressure::Tetens1::compute(ref_norm.dwpt)) }); - c.bench_function("vapour_pressure::buck3", |b| { - b.iter(|| vapour_pressure::buck3(black_box(300.0), black_box(101325.0))) + group.bench_function("buck1", |b| { + b.iter(|| vapour_pressure::Buck1::compute(ref_norm.dwpt, ref_norm.pres)) }); - c.bench_function("vapour_pressure::buck4", |b| { - b.iter(|| vapour_pressure::buck4(black_box(250.0), black_box(101325.0))) + group.bench_function("buck2", |b| { + b.iter(|| vapour_pressure::Buck2::compute(ref_freeze.dwpt, ref_freeze.pres)) }); - c.bench_function("vapour_pressure::buck3_simplified", |b| { - b.iter(|| vapour_pressure::buck3_simplified(black_box(300.0))) + group.bench_function("buck3", |b| { + b.iter(|| vapour_pressure::Buck3::compute(ref_norm.dwpt, ref_norm.pres)) }); - c.bench_function("vapour_pressure::buck4_simplified", |b| { - b.iter(|| vapour_pressure::buck4_simplified(black_box(250.0))) + group.bench_function("buck4", |b| { + b.iter(|| vapour_pressure::Buck4::compute(ref_freeze.dwpt, ref_freeze.pres)) }); - c.bench_function("vapour_pressure::saturation_specific1", |b| { - b.iter(|| vapour_pressure::saturation_specific1(black_box(3000.0), black_box(0.5))) + group.bench_function("buck3_simplified", |b| { + b.iter(|| vapour_pressure::Buck3Simplified::compute(ref_norm.dwpt)) }); - c.bench_function("vapour_pressure::saturation_specific2", |b| { - b.iter(|| vapour_pressure::saturation_specific2(black_box(3000.0), black_box(0.5))) + group.bench_function("buck4_simplified", |b| { + b.iter(|| vapour_pressure::Buck4Simplified::compute(ref_freeze.dwpt)) }); - c.bench_function("vapour_pressure::wexler1", |b| { - b.iter(|| vapour_pressure::wexler1(black_box(300.0))) + group.bench_function("wexler1", |b| { + b.iter(|| vapour_pressure::Wexler1::compute(ref_norm.dwpt)) }); - c.bench_function("vapour_pressure::wexler2", |b| { - b.iter(|| vapour_pressure::wexler2(black_box(250.0))) + group.bench_function("wexler2", |b| { + b.iter(|| vapour_pressure::Wexler2::compute(ref_freeze.dwpt)) }); + group.finish(); } -criterion_group!(benches, vapour_pressure_benchmark); +criterion_group!(benches, benchmark); criterion_main!(benches); diff --git a/benches/vapour_pressure_deficit.rs b/benches/vapour_pressure_deficit.rs index e376fc4..681fbf8 100644 --- a/benches/vapour_pressure_deficit.rs +++ b/benches/vapour_pressure_deficit.rs @@ -1,19 +1,19 @@ -use criterion::{Criterion, black_box, criterion_group, criterion_main}; -use floccus::vapour_pressure_deficit; +use criterion::{criterion_group, criterion_main, Criterion}; +use floccus::{formulas::vapour_pressure_deficit, formulas::Formula2}; -pub fn virtual_temperature_benchmark(c: &mut Criterion) { - c.bench_function("vapour_pressure_deficit::general1", |b| { - b.iter(|| vapour_pressure_deficit::general1(black_box(3000.0), black_box(3550.0))) - }); +mod utils; +use utils::ReferenceValues; - c.bench_function("vapour_pressure_deficit::general2", |b| { - b.iter(|| vapour_pressure_deficit::general2(black_box(300.0), black_box(290.0), black_box(101325.0))) - }); +pub fn benchmark(c: &mut Criterion) { + let ref_norm = ReferenceValues::normal(); + + let mut group = c.benchmark_group("vapour_pressure_deficit"); - c.bench_function("vapour_pressure_deficit::general3", |b| { - b.iter(|| vapour_pressure_deficit::general3(black_box(300.0), black_box(0.5), black_box(101325.0))) + group.bench_function("definition1", |b| { + b.iter(|| vapour_pressure_deficit::Definition1::compute(ref_norm.vapr, ref_norm.savp)) }); + group.finish(); } -criterion_group!(benches, virtual_temperature_benchmark); +criterion_group!(benches, benchmark); criterion_main!(benches); diff --git a/benches/virtual_temperature.rs b/benches/virtual_temperature.rs index d476a32..01ce84d 100644 --- a/benches/virtual_temperature.rs +++ b/benches/virtual_temperature.rs @@ -1,19 +1,29 @@ -use criterion::{Criterion, black_box, criterion_group, criterion_main}; -use floccus::virtual_temperature; +use criterion::{criterion_group, criterion_main, Criterion}; +use floccus::{formulas::virtual_temperature, formulas::{Formula2, Formula3}}; -pub fn virtual_temperature_benchmark(c: &mut Criterion) { - c.bench_function("virtual_temperature::general1", |b| { - b.iter(|| virtual_temperature::general1(black_box(300.0), black_box(0.022))) +mod utils; +use utils::ReferenceValues; + +pub fn benchmark(c: &mut Criterion) { + let ref_norm = ReferenceValues::normal(); + + let mut group = c.benchmark_group("virtual_temperature"); + + group.bench_function("definition1", |b| { + b.iter(|| virtual_temperature::Definition1::compute(ref_norm.temp, ref_norm.mxrt)) }); - c.bench_function("virtual_temperature::general2", |b| { - b.iter(|| virtual_temperature::general2(black_box(300.0), black_box(101325.0), black_box(3550.0))) + group.bench_function("definition2", |b| { + b.iter(|| { + virtual_temperature::Definition2::compute(ref_norm.temp, ref_norm.pres, ref_norm.vapr) + }) }); - c.bench_function("virtual_temperature::general3", |b| { - b.iter(|| virtual_temperature::general3(black_box(300.0), black_box(0.022))) + group.bench_function("definition3", |b| { + b.iter(|| virtual_temperature::Definition3::compute(ref_norm.temp, ref_norm.sphu)) }); + group.finish(); } -criterion_group!(benches, virtual_temperature_benchmark); +criterion_group!(benches, benchmark); criterion_main!(benches); diff --git a/benches/wet_bulb_potential_temperature.rs b/benches/wet_bulb_potential_temperature.rs index 104f841..38016e3 100644 --- a/benches/wet_bulb_potential_temperature.rs +++ b/benches/wet_bulb_potential_temperature.rs @@ -1,11 +1,19 @@ -use criterion::{Criterion, black_box, criterion_group, criterion_main}; -use floccus::wet_bulb_potential_temperature; +use criterion::{criterion_group, criterion_main, Criterion}; +use floccus::{formulas::wet_bulb_potential_temperature, formulas::Formula1}; -pub fn wet_bulb_potential_temperature_benchmark(c: &mut Criterion) { - c.bench_function("wet_bulb_potential_temperature::davies_jones1", |b| { - b.iter(|| wet_bulb_potential_temperature::davies_jones1(black_box(300.0))) +mod utils; +use utils::ReferenceValues; + +pub fn benchmark(c: &mut Criterion) { + let ref_norm = ReferenceValues::normal(); + + let mut group = c.benchmark_group("wet_bulb_potential_temperature"); + + group.bench_function("DaviesJones1", |b| { + b.iter(|| wet_bulb_potential_temperature::DaviesJones1::compute(ref_norm.thte)) }); + group.finish(); } -criterion_group!(benches, wet_bulb_potential_temperature_benchmark); -criterion_main!(benches); \ No newline at end of file +criterion_group!(benches, benchmark); +criterion_main!(benches); diff --git a/benches/wet_bulb_temperature.rs b/benches/wet_bulb_temperature.rs index 928e673..3df907a 100644 --- a/benches/wet_bulb_temperature.rs +++ b/benches/wet_bulb_temperature.rs @@ -1,11 +1,19 @@ -use criterion::{Criterion, black_box, criterion_group, criterion_main}; -use floccus::wet_bulb_temperature; +use criterion::{criterion_group, criterion_main, Criterion}; +use floccus::{formulas::wet_bulb_temperature, formulas::Formula2}; -pub fn wet_bulb_temperature_benchmark(c: &mut Criterion) { - c.bench_function("wet_bulb_temperature::stull1", |b| { - b.iter(|| wet_bulb_temperature::stull1(black_box(300.0), black_box(0.5))) +mod utils; +use utils::ReferenceValues; + +pub fn benchmark(c: &mut Criterion) { + let ref_norm = ReferenceValues::normal(); + + let mut group = c.benchmark_group("wet_bulb_temperature"); + + group.bench_function("stull1", |b| { + b.iter(|| wet_bulb_temperature::Stull1::compute(ref_norm.temp, ref_norm.rehu)) }); + group.finish(); } -criterion_group!(benches, wet_bulb_temperature_benchmark); +criterion_group!(benches, benchmark); criterion_main!(benches); diff --git a/src/constants.rs b/src/constants.rs index 1b9c7cb..8245ec9 100644 --- a/src/constants.rs +++ b/src/constants.rs @@ -1,48 +1,140 @@ -//!Module containing physical constants - -use crate::Float; - -///Temperature of 0 Celsius in `K` -pub const ZERO_CELSIUS: Float = 273.15; - -///Gravitational acceleration in `m s^-2` -pub const G: Float = 9.80665; - -///Universal gas constant in `J K^-1 mol^-1` -pub const R: Float = 8.314_462_618_153_24; - -///Molar mass of dry air in `kg mol^-1` (ECMWF, 2020) -pub const M_D: Float = 0.028_964_4; - -///Molar mass of water vapour in `kg mol^-1` -pub const M_V: Float = 0.018_015_283_3; - -///Specific heat capacity of dry air at constant pressure in `J kg^-1 K^-1` (ECMWF, 2020) -pub const C_P: Float = 1004.709; - -///Specific heat capacity of dry air at constant volume in `J kg^-1 K^-1` (ECMWF, 2020) -pub const C_V: Float = 717.6493; - -///Specific heat capacity of water vapour at constant pressure in `J kg^-1 K^-1` (ECMWF, 2020) -pub const C_PV: Float = 1846.1; - -///Specific heat capacity of water vapour at constant volume in `J kg^-1 K^-1` (ECMWF, 2020) -pub const C_VV: Float = 1384.575; - -///Specific heat capacity of liquid water in `J kg^-1 K^-1` (ECMWF, 2020) -pub const C_L: Float = 4218.0; - -///Specific heat capacity of solid water in `J kg^-1 K^-1` (ECMWF, 2020) -pub const C_S: Float = 2106.0; - -///Mass latent heat of vapourization of water in `J kg^1` (ECMWF, 2020) -pub const L_V: Float = 2_500_800.0; - -///Ratio of molar masses of dry air and water vapour in `no unit` -pub const EPSILON: Float = M_V / M_D; - -///Specific gas constant for dry air in `J kg^-1 K^-1` -pub const R_D: Float = R / M_D; - -///Specific gas constant for water vapour in `J kg^-1 K^-1` -pub const R_V: Float = R / M_V; +//! Module containing physical constants + +use std::marker::PhantomData; + +use crate::Storage; + +/// Gravitational acceleration of Earth +pub const G: Storage::Acceleration = Storage::Acceleration { + dimension: PhantomData, + units: PhantomData, + value: 9.806_65, +}; + +/// Universal gas constant +pub const R: Storage::MolarHeatCapacity = Storage::MolarHeatCapacity { + dimension: PhantomData, + units: PhantomData, + value: 8.314_462_618_153_24, +}; + +/// Molar mass of dry air (ECMWF, 2020) +pub const M_D: Storage::MolarMass = Storage::MolarMass { + dimension: PhantomData, + units: PhantomData, + value: 0.028_964_4, +}; + +/// Molar mass of water vapour +pub const M_V: Storage::MolarMass = Storage::MolarMass { + dimension: PhantomData, + units: PhantomData, + value: 0.018_015_283_3, +}; + +/// Specific heat capacity of dry air at constant pressure (ECMWF, 2020) +pub const C_P: Storage::SpecificHeatCapacity = Storage::SpecificHeatCapacity { + dimension: PhantomData, + units: PhantomData, + value: 1004.709, +}; + +/// Specific heat capacity of dry air at constant volume (ECMWF, 2020) +pub const C_V: Storage::SpecificHeatCapacity = Storage::SpecificHeatCapacity { + dimension: PhantomData, + units: PhantomData, + value: 717.6493, +}; + +/// Specific heat capacity of water vapour at constant pressure (ECMWF, 2020) +pub const C_PV: Storage::SpecificHeatCapacity = Storage::SpecificHeatCapacity { + dimension: PhantomData, + units: PhantomData, + value: 1846.1, +}; + +/// Specific heat capacity of water vapour at constant volume (ECMWF, 2020) +pub const C_VV: Storage::SpecificHeatCapacity = Storage::SpecificHeatCapacity { + dimension: PhantomData, + units: PhantomData, + value: 1384.575, +}; + +/// Specific heat capacity of liquid water (ECMWF, 2020) +pub const C_L: Storage::SpecificHeatCapacity = Storage::SpecificHeatCapacity { + dimension: PhantomData, + units: PhantomData, + value: 4218.0, +}; + +/// Specific heat capacity of solid water (ECMWF, 2020) +pub const C_S: Storage::SpecificHeatCapacity = Storage::SpecificHeatCapacity { + dimension: PhantomData, + units: PhantomData, + value: 2106.0, +}; + +/// Specific latent heat of vapourization of water (ECMWF, 2020) +pub const L_V: Storage::AvailableEnergy = Storage::AvailableEnergy { + dimension: PhantomData, + units: PhantomData, + value: 2_500_800.0, +}; + +/// Specific gas constant for dry air +pub const R_D: Storage::SpecificHeatCapacity = Storage::SpecificHeatCapacity { + dimension: PhantomData, + units: PhantomData, + value: R.value / M_D.value, +}; + +/// Specific gas constant for water vapour +pub const R_V: Storage::SpecificHeatCapacity = Storage::SpecificHeatCapacity { + dimension: PhantomData, + units: PhantomData, + value: R.value / M_V.value, +}; + +// Internal Constants (commonly appearing in formulas to use them with uom units) + +/// Reference pressure level of 1000hPa. Used in adiabatic processes +pub(crate) const REF_PRES: Storage::Pressure = Storage::Pressure { + dimension: PhantomData, + units: PhantomData, + value: 100_000., +}; + +/// Ratio of molar masses of dry air and water vapour +pub(crate) const EPSILON: Storage::Ratio = Storage::Ratio { + dimension: PhantomData, + units: PhantomData, + value: M_V.value / M_D.value, +}; + +/// Ratio of specific gas constant and specific heat capacity for dry air +pub(crate) const KAPPA: Storage::Ratio = Storage::Ratio { + dimension: PhantomData, + units: PhantomData, + value: R_D.value / C_P.value, +}; + +/// Inverse of KAPPA +pub(crate) const LAMBDA: Storage::Ratio = Storage::Ratio { + dimension: PhantomData, + units: PhantomData, + value: C_P.value / R_D.value, +}; + +pub(crate) const DIMLESS_ONE: Storage::Ratio = Storage::Ratio { + dimension: PhantomData, + units: PhantomData, + value: 1.0, +}; + +/// Useful to convert `TemperatureInterval` into `ThermodynamicTemperature` +pub(crate) const ZERO_KELVIN: Storage::ThermodynamicTemperature = + Storage::ThermodynamicTemperature { + dimension: PhantomData, + units: PhantomData, + value: 0.0, + }; diff --git a/src/equivalent_potential_temperature.rs b/src/equivalent_potential_temperature.rs deleted file mode 100644 index cd737bb..0000000 --- a/src/equivalent_potential_temperature.rs +++ /dev/null @@ -1,226 +0,0 @@ -//!Functions to calculate equivalent potential temperature of air in K. -use crate::constants::{C_L, R_V}; -use crate::Float; -use crate::{ - constants::{C_P, EPSILON, L_V, R_D}, - errors::InputError, - mixing_ratio, potential_temperature, relative_humidity, vapour_pressure, -}; - -#[cfg(feature = "debug")] -use floccus_proc::logerr; - -///Most accuarte formula for computing equivalent potential temperature of unsaturated air from -///temperature, pressure and vapour pressure. -/// -///Implementation of this formula assumes no liquid or solid water in the air parcel. -/// -///Provided in Emmanuel, Kerry (1994). Atmospheric Convection. Oxford University Press. -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `temperature` range: 253K - 324K\ -///Valid `pressure` range: 100Pa - 150000Pa\ -///Valid `vapour_pressure` range: 0Pa - 10000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn general1( - temperature: Float, - pressure: Float, - vapour_pressure: Float, -) -> Result { - if !(253.0..=324.0).contains(&temperature) { - return Err(InputError::OutOfRange(String::from("temperature"))); - } - - if !(20000.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - if !(0.0..=10_000.0).contains(&vapour_pressure) { - return Err(InputError::OutOfRange(String::from("vapour_pressure"))); - } - - let p0 = 100_000.0; - - let mixing_ratio = mixing_ratio::general1(pressure, vapour_pressure)?; - let saturation_vapour_pressure = vapour_pressure::buck1(temperature, pressure)?; - - let relative_humidity = - relative_humidity::general2(vapour_pressure, saturation_vapour_pressure)?; - - let result = temperature - * (p0 / pressure).powf(R_D / (C_P + mixing_ratio * C_L)) - * relative_humidity.powf((-mixing_ratio * R_V) / (C_P + mixing_ratio * C_L)) - * ((L_V * mixing_ratio) / (temperature * (C_P + mixing_ratio * C_L))).exp(); - - Ok(result) -} - -///Formula for computing equivalent potential temperature of unsaturated air from -///temperature, pressure and vapour pressure. -/// -///Derived by G. H. Bryan (2008) [(doi:10.1175/2008MWR2593.1)](https://doi.org/10.1175/2008MWR2593.1) -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `temperature` range: 253K - 324K\ -///Valid `pressure` range: 100Pa - 150000Pa\ -///Valid `vapour_pressure` range: 0Pa - 10000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn bryan1( - temperature: Float, - pressure: Float, - vapour_pressure: Float, -) -> Result { - if !(253.0..=324.0).contains(&temperature) { - return Err(InputError::OutOfRange(String::from("temperature"))); - } - - if !(20000.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - if !(0.0..=10_000.0).contains(&vapour_pressure) { - return Err(InputError::OutOfRange(String::from("vapour_pressure"))); - } - - let kappa = R_D / C_P; - - let potential_temperature = - potential_temperature::davies_jones1(temperature, pressure, vapour_pressure)?; - - let saturation_vapour_pressure = vapour_pressure::buck3(temperature, pressure)?; - let relative_humidity = - relative_humidity::general2(vapour_pressure, saturation_vapour_pressure)?; - - let mixing_ratio = mixing_ratio::general1(pressure, vapour_pressure)?; - - let result = potential_temperature - * relative_humidity.powf((-kappa) * (mixing_ratio / EPSILON)) - * ((L_V * mixing_ratio) / (C_P * temperature)).exp(); - - Ok(result) -} - -///Approximate formula for computing equivalent potential temperature of unsaturated air from -///temperature, pressure and dewpoint. -/// -///Derived by D. Bolton (1980) -///[(doi:10.1175/1520-0493(1980)108<1046:TCOEPT>2.0.CO;2)](https://doi.org/10.1175/1520-0493(1980)108%3C1046:TCOEPT%3E2.0.CO;2) -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `pressure` range: 100Pa - 150000Pa\ -///Valid `temperature` range: 253K - 324K\ -///Valid `dewpoint` range: 253K - 324K -#[cfg_attr(feature = "debug", logerr)] -pub fn bolton1(pressure: Float, temperature: Float, dewpoint: Float) -> Result { - if !(20000.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - if !(253.0..=324.0).contains(&temperature) { - return Err(InputError::OutOfRange(String::from("temperature"))); - } - - if !(253.0..=324.0).contains(&dewpoint) { - return Err(InputError::OutOfRange(String::from("dewpoint"))); - } - - let kappa = R_D / C_P; - - let vapour_pressure = vapour_pressure::buck3(dewpoint, pressure)?; - let mixing_ratio = mixing_ratio::general1(pressure, vapour_pressure)?; - - let lcl_temp = - (1.0 / ((1.0 / (dewpoint - 56.0)) + ((temperature / dewpoint).ln() / 800.0))) + 56.0; - - let theta_dl = temperature - * (100000.0 / (pressure - vapour_pressure)).powf(kappa) - * (temperature / lcl_temp).powf(0.28 * mixing_ratio); - - let result = theta_dl - * (((3036.0 / lcl_temp) - 1.78) * mixing_ratio * (1.0 + 0.448 * mixing_ratio)).exp(); - - Ok(result) -} - -#[cfg(test)] -mod tests { - use crate::{ - equivalent_potential_temperature, - tests_framework::{self, Argument}, - }; - - #[test] - fn general1() { - assert!(tests_framework::test_with_3args( - &equivalent_potential_temperature::general1, - Argument { - name: "temperature", - def_val: 300.0, - range: [253.0, 324.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [20000.0, 150_000.0] - }, - Argument { - name: "vapour_pressure", - def_val: 991.189131, - range: [0.0, 10_000.0] - }, - 315.23724970376776 - )); - } - - #[test] - fn bryan1() { - assert!(tests_framework::test_with_3args( - &equivalent_potential_temperature::bryan1, - Argument { - name: "temperature", - def_val: 300.0, - range: [253.0, 324.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [20000.0, 150_000.0] - }, - Argument { - name: "vapour_pressure", - def_val: 991.189131, - range: [0.0, 10_000.0] - }, - 316.52762026634014 - )); - } - - #[test] - fn bolton1() { - assert!(tests_framework::test_with_3args( - &equivalent_potential_temperature::bolton1, - Argument { - name: "pressure", - def_val: 101325.0, - range: [20000.0, 150_000.0] - }, - Argument { - name: "temperature", - def_val: 300.0, - range: [253.0, 324.0] - }, - Argument { - name: "dewpoint", - def_val: 280.0, - range: [253.0, 324.0] - }, - 317.3855211897774 - )); - } -} diff --git a/src/errors.rs b/src/errors.rs index 2f60b7a..301e46c 100644 --- a/src/errors.rs +++ b/src/errors.rs @@ -2,10 +2,9 @@ use thiserror::Error; -#[derive(Error, Debug, PartialEq, Eq)] +#[derive(Error, Debug, PartialEq, Eq, Clone)] ///Error enum returned when provided input will cause function to return erronous result ///eg. `Inf` or `NaN` - pub enum InputError { #[error("Value of {0} out of a reasonable range.")] ///Error returned when provided input is out of reasonable range. @@ -28,7 +27,7 @@ pub enum InputError { /// ///If you find that in your use case input ranges are too narrow you should first look for a more relevant formula. ///If such formula does not exist do not hesitate to create an issue in Github repository. - OutOfRange(String), + OutOfRange(&'static str), ///Error returned when provided set of arguments will result in invalid output. ///Contains detailed information about the error. @@ -46,8 +45,10 @@ pub enum InputError { /// ///This error should be handled on case-to-case basis, as it can be returned by functions ///for different reasons. Check the documentation of function that you use to learn more - ///about when this error can appear. - #[error("Provided arguments result in erronous output. - Check documentation of the function and change one of arguments. Details: {0}")] - IncorrectArgumentSet(String), + ///about when this error can appear. + #[error( + "Provided arguments result in erronous output. + Check documentation of the function and change one of arguments. Details: {0}" + )] + IncorrectArgumentSet(&'static str), } diff --git a/src/formulas.rs b/src/formulas.rs new file mode 100644 index 0000000..d76aa30 --- /dev/null +++ b/src/formulas.rs @@ -0,0 +1,18 @@ +pub mod equivalent_potential_temperature; +pub mod mixing_ratio; +pub mod potential_temperature; +pub mod relative_humidity; +pub mod saturation_mixing_ratio; +pub mod saturation_vapour_pressure; +pub mod specific_humidity; +pub mod vapour_pressure; +pub mod vapour_pressure_deficit; +pub mod virtual_temperature; +pub mod wet_bulb_potential_temperature; +pub mod wet_bulb_temperature; + +mod traits; +pub mod isobaric_equivalent_temperature; +pub mod adiabatic_equivalent_temperature; + +pub use traits::{Formula1, Formula2, Formula3, Formula4}; diff --git a/src/formulas/adiabatic_equivalent_temperature.rs b/src/formulas/adiabatic_equivalent_temperature.rs new file mode 100644 index 0000000..9322bc1 --- /dev/null +++ b/src/formulas/adiabatic_equivalent_temperature.rs @@ -0,0 +1,151 @@ +//! Functions to calculate adiabatic potential temperature of unsaturated air +//! +//! Adiabatic equivalent temperature (also known as pseudoequivalent temperature) is a temperature that +//! an air parcel would have after undergoing the following process: dry-adiabatic expansion until saturated; +//! pseudoadiabatic expansion until all moisture is precipitated out; dry- adiabatic compression to the initial +//! pressure. This is the equivalent temperature as read from a thermodynamic chart and is always greater than +//! the isobaric equivalent temperature. +//! ([AMETSOC Glossary](https://glossary.ametsoc.org/wiki/Equivalent_temperature)). + +use uom::si::pressure::pascal; +use uom::si::ratio::ratio; +use uom::si::thermodynamic_temperature::kelvin; + +use crate::constants::{C_P, KAPPA, L_V, REF_PRES, ZERO_KELVIN}; +use crate::errors::InputError; +use crate::formulas::Formula2; +use crate::quantities::{ + AdiabaticEquivalentTemperature, AtmosphericPressure, DryBulbTemperature, + EquivalentPotentialTemperature, MixingRatio, QuantityHelpers, +}; + +type FormulaQuantity = AdiabaticEquivalentTemperature; + +/// Formula for computing adiabatic equivalent temperature from dry-bulb temperature and mixing ratio +/// +/// Valid `temperature` range: 253K - 324K +/// +/// Valid `mixing_ratio` range: 0.000_000_1 - 2.0 +pub struct Definition1; + +impl Formula2 for Definition1 { + #[inline] + fn validate_inputs_internal( + temperature: DryBulbTemperature, + mixing_ratio: MixingRatio, + ) -> Result<(), InputError> { + temperature.check_range_si(253.0, 324.0)?; + mixing_ratio.check_range_si(0.000_000_1, 2.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + temperature: DryBulbTemperature, + mixing_ratio: MixingRatio, + ) -> FormulaQuantity { + let temp_ae = temperature.0 * ((L_V * mixing_ratio.0) / (C_P * temperature.0)).exp(); + + AdiabaticEquivalentTemperature(temp_ae + ZERO_KELVIN) + } +} + +/// Formula for computing adiabatic equivalent temperature from +/// atmospheric pressure and equivalent potential temperature. +/// +/// Cited in [Davies-Jones (2008)](https://doi.org/10.1175/2007MWR2224.1) +/// +/// Valid `equivalent_potential_temperature` range: 173K - 373K +/// +/// Valid `pressure` range: 100Pa - 150000Pa +pub struct Definition2; + +impl Formula2 + for Definition2 +{ + #[inline] + fn validate_inputs_internal( + equivalent_potential_temperature: EquivalentPotentialTemperature, + pressure: AtmosphericPressure, + ) -> Result<(), InputError> { + equivalent_potential_temperature.check_range_si(173., 373.)?; + pressure.check_range_si(100., 150_000.)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + equivalent_potential_temperature: EquivalentPotentialTemperature, + pressure: AtmosphericPressure, + ) -> FormulaQuantity { + let p = pressure.get::(); + let p0 = REF_PRES.get::(); + let kappa = KAPPA.get::(); + let theta_e = equivalent_potential_temperature.get::(); + + let pi = (p / p0).powf(kappa); + let temp_ae = theta_e * pi; + + AdiabaticEquivalentTemperature::new::(temp_ae) + } +} + +#[cfg(test)] +mod tests { + + use crate::tests::{test_with_2args, testing_traits::ReferenceAtmosphere, Argument}; + + use super::*; + + #[test] + fn definition1_norm() { + test_with_2args::( + Argument::new([253., 324.]), + Argument::new([0.000_000_1, 2.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } + + #[test] + fn definition2_norm() { + test_with_2args::< + FormulaQuantity, + EquivalentPotentialTemperature, + AtmosphericPressure, + Definition2, + >( + Argument::new([173., 373.]), + Argument::new([100., 150_000.]), + ReferenceAtmosphere::Normal, + 1., + ); + } + + #[test] + fn definition1_freez() { + test_with_2args::( + Argument::new([253., 324.]), + Argument::new([0.000_000_1, 2.0]), + ReferenceAtmosphere::Freezing, + 1e-12, + ); + } + + #[test] + fn definition2_freez() { + test_with_2args::< + FormulaQuantity, + EquivalentPotentialTemperature, + AtmosphericPressure, + Definition2, + >( + Argument::new([173., 373.]), + Argument::new([100., 150_000.]), + ReferenceAtmosphere::Freezing, + 1e-1, + ); + } +} diff --git a/src/formulas/equivalent_potential_temperature.rs b/src/formulas/equivalent_potential_temperature.rs new file mode 100644 index 0000000..0af55c5 --- /dev/null +++ b/src/formulas/equivalent_potential_temperature.rs @@ -0,0 +1,407 @@ +//! Functions to calculate equivalent potential temperature of air +//! +//! Equivalent potential eemperature is a thermodynamic quantity, with its natural logarithm proportional +//! to the entropy of moist air, that is conserved in a reversible moist +//! adiabatic process ([AMETSOC Glossary](https://glossary.ametsoc.org/wiki/Equivalent_potential_temperature)). + +use float_cmp::approx_eq; +use uom::si::available_energy::joule_per_kilogram; +use uom::si::pressure::pascal; +use uom::si::ratio::ratio; +use uom::si::specific_heat_capacity::joule_per_kilogram_kelvin; +use uom::si::thermodynamic_temperature::kelvin; + +use crate::constants::{C_L, C_P, EPSILON, KAPPA, L_V, REF_PRES, R_D, R_V}; +use crate::errors::InputError; +use crate::quantities::{ + AtmosphericPressure, DewPointTemperature, DryBulbTemperature, EquivalentPotentialTemperature, + MixingRatio, PotentialTemperature, QuantityHelpers, RelativeHumidity, VapourPressure, +}; +use crate::{ + formulas::{mixing_ratio, Formula2, Formula4}, + Float, +}; + +type FormulaQuantity = EquivalentPotentialTemperature; + +/// Most accurate formula for computing equivalent potential temperature of unsaturated air from +/// temperature, pressure and mixing ratio and relative humidity. +/// +/// Implementation of this formula assumes no liquid or solid water in the air parcel. +/// +/// Appears in Emmanuel, Kerry (1994). Atmospheric Convection. Oxford University Press. +/// as equation 4.5.11 +/// +/// Valid `temperature` range: 253K - 324K +/// +/// Valid `pressure` range: 100Pa - 150000Pa +/// +/// Valid `mixing_ratio` range: 0.000_000_1 - 2.0 +/// +/// Valid `relative_humidity` range: 0.000_000_1 - 2.0 +pub struct Kerry1; + +impl + Formula4< + FormulaQuantity, + DryBulbTemperature, + AtmosphericPressure, + MixingRatio, + RelativeHumidity, + > for Kerry1 +{ + #[inline] + fn validate_inputs_internal( + temperature: DryBulbTemperature, + pressure: AtmosphericPressure, + mixing_ratio: MixingRatio, + relative_humidity: RelativeHumidity, + ) -> Result<(), InputError> { + temperature.check_range_si(253.0, 324.0)?; + pressure.check_range_si(100.0, 150_000.0)?; + mixing_ratio.check_range_si(0.000_000_1, 2.0)?; + relative_humidity.check_range_si(0.000_000_1, 2.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + temperature: DryBulbTemperature, + pressure: AtmosphericPressure, + mixing_ratio: MixingRatio, + relative_humidity: RelativeHumidity, + ) -> EquivalentPotentialTemperature { + let temperature = temperature.0.get::(); + let pressure = pressure.0.get::(); + let mixing_ratio = mixing_ratio.0.get::(); + let relative_humidity = relative_humidity.0.get::(); + + let p0 = REF_PRES.get::(); + let r_d = R_D.get::(); + let r_v = R_V.get::(); + let l_v = L_V.get::(); + let c_p = C_P.get::(); + let c_l = C_L.get::(); + + let result = temperature + * (p0 / pressure).powf(r_d / (c_p + mixing_ratio * c_l)) + * relative_humidity.powf((-mixing_ratio * r_v) / (c_p + mixing_ratio * c_l)) + * ((l_v * mixing_ratio) / (temperature * (c_p + mixing_ratio * c_l))).exp(); + + EquivalentPotentialTemperature::new::(result) + } +} + +/// Formula for computing equivalent potential temperature of unsaturated air from +/// temperature, mixing ratio, relative humidity and potential temperature. +/// +/// Derived by G. H. Bryan (2008) [(doi:10.1175/2008MWR2593.1)](https://doi.org/10.1175/2008MWR2593.1) +/// +/// Valid `temperature` range: 253K - 324K +/// +/// Valid `mixing_ratio` range: 0.000_000_1 - 2.0 +/// +/// Valid `relative_humidity` range: 0.000_000_1 - 2.0 +/// +/// Valid `potential_temperature` range: 253K - 324K +pub struct Bryan1; + +impl + Formula4< + FormulaQuantity, + DryBulbTemperature, + MixingRatio, + RelativeHumidity, + PotentialTemperature, + > for Bryan1 +{ + #[inline] + fn validate_inputs_internal( + temperature: DryBulbTemperature, + mixing_ratio: MixingRatio, + relative_humidity: RelativeHumidity, + potential_temperature: PotentialTemperature, + ) -> Result<(), InputError> { + temperature.check_range_si(253.0, 324.0)?; + mixing_ratio.check_range_si(0.000_000_1, 2.0)?; + relative_humidity.check_range_si(0.000_000_1, 2.0)?; + potential_temperature.check_range_si(253.0, 324.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + temperature: DryBulbTemperature, + mixing_ratio: MixingRatio, + relative_humidity: RelativeHumidity, + potential_temperature: PotentialTemperature, + ) -> EquivalentPotentialTemperature { + let temperature = temperature.0.get::(); + let mixing_ratio = mixing_ratio.0.get::(); + let relative_humidity = relative_humidity.0.get::(); + let potential_temperature = potential_temperature.0.get::(); + + let kappa = KAPPA.get::(); + let l_v = L_V.get::(); + let c_p = C_P.get::(); + let epsilon = EPSILON.get::(); + + let result = potential_temperature + * relative_humidity.powf((-kappa) * (mixing_ratio / epsilon)) + * ((l_v * mixing_ratio) / (c_p * temperature)).exp(); + + EquivalentPotentialTemperature::new::(result) + } +} + +/// Approximate formula for computing equivalent potential temperature of unsaturated air from +/// temperature, pressure, dewpoint and vapour pressure. +/// +/// Derived by D. Bolton (1980) +/// [(doi:10.1175/1520-0493(1980)108<1046:TCOEPT>2.0.CO;2)](https://doi.org/10.1175/1520-0493(1980)108%3C1046:TCOEPT%3E2.0.CO;2) +/// +/// Valid `pressure` range: 100Pa - 150000Pa +/// +/// Valid `temperature` range: 253K - 324K +/// +/// Valid `dewpoint` range: 253K - 324K +/// +/// Valid `vapour_pressure` range: 0Pa - 50000Pa +pub struct Bolton1; + +impl + Formula4< + FormulaQuantity, + AtmosphericPressure, + DryBulbTemperature, + DewPointTemperature, + VapourPressure, + > for Bolton1 +{ + #[inline] + fn validate_inputs_internal( + pressure: AtmosphericPressure, + temperature: DryBulbTemperature, + dewpoint: DewPointTemperature, + vapour_pressure: VapourPressure, + ) -> Result<(), InputError> { + pressure.check_range_si(100.0, 150_000.0)?; + temperature.check_range_si(253.0, 324.0)?; + dewpoint.check_range_si(253.0, 324.0)?; + vapour_pressure.check_range_si(0.0, 50_000.0)?; + + if approx_eq!( + Float, + pressure.get_si_value(), + vapour_pressure.get_si_value(), + ulps = 2 + ) { + return Err(InputError::IncorrectArgumentSet( + "pressure must be greater than vapour pressure", + )); + } + + if vapour_pressure.0 > pressure.0 { + return Err(InputError::IncorrectArgumentSet( + "pressure must be greater than vapour pressure", + )); + } + + if dewpoint.0 > temperature.0 { + return Err(InputError::IncorrectArgumentSet( + "dewpoint must be less than temperature" + )); + } + + let mixing_ratio = mixing_ratio::Definition1::compute_unchecked(pressure, vapour_pressure); + + mixing_ratio.check_range_si(0.000_000_1, 2.0).map_err(|_| InputError::IncorrectArgumentSet( + "pressure and vapour_pressure must give mixing_ratio less than 2 so cannot be close to each other", + ))?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + pressure: AtmosphericPressure, + temperature: DryBulbTemperature, + dewpoint: DewPointTemperature, + vapour_pressure: VapourPressure, + ) -> EquivalentPotentialTemperature { + let mixing_ratio = mixing_ratio::Definition1::compute_unchecked(pressure, vapour_pressure); + + let pressure = pressure.0.get::(); + let temperature = temperature.0.get::(); + let dewpoint = dewpoint.0.get::(); + let mixing_ratio = mixing_ratio.0.get::(); + let vapour_pressure = vapour_pressure.0.get::(); + let p0 = REF_PRES.get::(); + + let kappa = KAPPA.get::(); + + // technically LCL and Theta-DL could be extracted to separate functions + + let lcl_temp = + (1.0 / ((1.0 / (dewpoint - 56.0)) + ((temperature / dewpoint).ln() / 800.0))) + 56.0; + + let theta_dl = temperature + * (p0 / (pressure - vapour_pressure)).powf(kappa) + * (temperature / lcl_temp).powf(0.28 * mixing_ratio); + + let result = theta_dl + * (((3036.0 / lcl_temp) - 1.78) * mixing_ratio * (1.0 + 0.448 * mixing_ratio)).exp(); + + EquivalentPotentialTemperature::new::(result) + } +} + +/// Approximate formula for computing equivalent potential temperature of unsaturated air from +/// temperature, dewpoint, mixing ratio and potential temperature. +/// +/// Derived by D. Bolton (1980) +/// [(doi:10.1175/1520-0493(1980)108<1046:TCOEPT>2.0.CO;2)](https://doi.org/10.1175/1520-0493(1980)108%3C1046:TCOEPT%3E2.0.CO;2) +/// +/// Valid `temperature` range: 253K - 324K +/// +/// Valid `dewpoint` range: 253K - 324K +/// +/// Valid `mixing_ratio` range: 0.000_000_1 - 2.0 +/// +/// Valid `potential_temperature` range: 253K - 324K +pub struct Bolton2; + +impl + Formula4< + EquivalentPotentialTemperature, + DryBulbTemperature, + DewPointTemperature, + MixingRatio, + PotentialTemperature, + > for Bolton2 +{ + fn compute_unchecked( + temperature: DryBulbTemperature, + dewpoint: DewPointTemperature, + mixing_ratio: MixingRatio, + potential_temperature: PotentialTemperature, + ) -> EquivalentPotentialTemperature { + let temperature = temperature.0.get::(); + let dewpoint = dewpoint.0.get::(); + let mixing_ratio = mixing_ratio.0.get::(); + let theta = potential_temperature.0.get::(); + + let lcl_temp = + (1.0 / ((1.0 / (dewpoint - 56.0)) + ((temperature / dewpoint).ln() / 800.0))) + 56.0; + + let result = theta + * (((3.376 / lcl_temp) - 0.00254) + * (1000.0 * mixing_ratio) + * (1.0 + (0.81 * mixing_ratio))) + .exp(); + + EquivalentPotentialTemperature::new::(result) + } + + fn validate_inputs_internal( + temperature: DryBulbTemperature, + dewpoint: DewPointTemperature, + mixing_ratio: MixingRatio, + potential_temperature: PotentialTemperature, + ) -> Result<(), InputError> { + temperature.check_range_si(253.0, 324.0)?; + dewpoint.check_range_si(253.0, 324.0)?; + mixing_ratio.check_range_si(0.000_000_1, 2.0)?; + potential_temperature.check_range_si(253.0, 324.0)?; + + Ok(()) + } +} + +#[cfg(test)] +mod tests { + use crate::{ + quantities::{AtmosphericPressure, DryBulbTemperature, VapourPressure}, + tests::{test_with_4args, testing_traits::ReferenceAtmosphere, Argument}, + }; + + use super::*; + + #[test] + fn paluch1() { + test_with_4args::< + FormulaQuantity, + DryBulbTemperature, + AtmosphericPressure, + MixingRatio, + RelativeHumidity, + Kerry1, + >( + Argument::new([253.0, 324.0]), + Argument::new([100.0, 150_000.0]), + Argument::new([0.000_000_1, 2.0]), + Argument::new([0.000_000_1, 2.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } + + #[test] + fn bryan1() { + test_with_4args::< + FormulaQuantity, + DryBulbTemperature, + MixingRatio, + RelativeHumidity, + PotentialTemperature, + Bryan1, + >( + Argument::new([253.0, 324.0]), + Argument::new([0.000_000_1, 2.0]), + Argument::new([0.000_000_1, 2.0]), + Argument::new([253.0, 324.0]), + ReferenceAtmosphere::Normal, + 1e1, + ); + } + + #[test] + fn bolton1() { + test_with_4args::< + FormulaQuantity, + AtmosphericPressure, + DryBulbTemperature, + DewPointTemperature, + VapourPressure, + Bolton1, + >( + Argument::new([100.0, 150_000.0]), + Argument::new([253.0, 324.0]), + Argument::new([253.0, 324.0]), + Argument::new([0.0, 50_000.0]), + ReferenceAtmosphere::Normal, + 1e1, + ); + } + + #[test] + fn bolton2() { + test_with_4args::< + FormulaQuantity, + DryBulbTemperature, + DewPointTemperature, + MixingRatio, + PotentialTemperature, + Bolton2, + >( + Argument::new([253.0, 324.0]), + Argument::new([253.0, 324.0]), + Argument::new([0.000_000_1, 2.0]), + Argument::new([253.0, 324.0]), + ReferenceAtmosphere::Normal, + 1e1, + ); + } +} diff --git a/src/formulas/isobaric_equivalent_temperature.rs b/src/formulas/isobaric_equivalent_temperature.rs new file mode 100644 index 0000000..cfce124 --- /dev/null +++ b/src/formulas/isobaric_equivalent_temperature.rs @@ -0,0 +1,66 @@ +//! Functions to calculate isobaric potential temperature of unsaturated air +//! +//! +//! Isobaric equivalent temperature is a temperature that an air parcel would have if all water +//! vapor were condensed at constant pressure and the enthalpy released from the vapor used to heat the air. +//! This process is physically impossible in the atmosphere +//! ([AMETSOC Glossary](https://glossary.ametsoc.org/wiki/Equivalent_temperature)). + +use crate::constants::{C_P, DIMLESS_ONE, L_V, ZERO_KELVIN}; +use crate::errors::InputError; +use crate::formulas::Formula2; +use crate::quantities::{ + DryBulbTemperature, IsobaricEquivalentTemperature, MixingRatio, QuantityHelpers, +}; + +type FormulaQuantity = IsobaricEquivalentTemperature; + +/// Formula for computing isobaric equivalent temperature from dry-bulb temperature and mixing ratio +/// +/// Valid `temperature` range: 253K - 324K +/// +/// Valid `mixing_ratio` range: 0.000_000_1 - 2.0 +/// +pub struct Definition1; + +impl Formula2 for Definition1 { + #[inline] + fn validate_inputs_internal( + temperature: DryBulbTemperature, + mixing_ratio: MixingRatio, + ) -> Result<(), InputError> { + temperature.check_range_si(253.0, 324.0)?; + mixing_ratio.check_range_si(0.000_000_1, 2.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + temperature: DryBulbTemperature, + mixing_ratio: MixingRatio, + ) -> FormulaQuantity { + let temp_ie = + temperature.0 * (DIMLESS_ONE + ((L_V * mixing_ratio.0) / (C_P * temperature.0))); + + IsobaricEquivalentTemperature(ZERO_KELVIN + temp_ie) + } +} + +#[cfg(test)] +mod tests { + + use crate::tests::{test_with_2args, testing_traits::ReferenceAtmosphere, Argument}; + + use super::*; + + #[test] + fn definition1() { + test_with_2args::( + Argument::new([253., 324.]), + Argument::new([0.000_000_1, 2.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } +} diff --git a/src/formulas/mixing_ratio.rs b/src/formulas/mixing_ratio.rs new file mode 100644 index 0000000..cfda1ef --- /dev/null +++ b/src/formulas/mixing_ratio.rs @@ -0,0 +1,79 @@ +//! Functions to calculate mixing ratio of water vapour in unsaturated air +//! +//! Mixing ratio is the ratio of the mass of a variable atmospheric constituent to the mass +//! of dry air ([AMETSOC Glossary](https://glossary.ametsoc.org/wiki/Mixing_ratio)). + +use crate::Float; +use crate::formulas::Formula2; +use crate::quantities::{AtmosphericPressure, MixingRatio, QuantityHelpers, VapourPressure}; +use crate::{constants::EPSILON, errors::InputError}; +use float_cmp::approx_eq; + +type FormulaQuantity = MixingRatio; + +/// Formula for computing mixing ratio of unsaturated air from air pressure and vapour pressure +/// +/// Valid `pressure` range: 100Pa - 150000Pa +/// +/// Valid `vapour_pressure` range: 0Pa - 10000Pa +/// +/// Returns [`InputError::IncorrectArgumentSet`] when: +/// +/// - inputs are equal and division by 0 would occur +/// - `vapour_pressure` is greater than `pressure` +pub struct Definition1; + +impl Formula2 for Definition1 { + #[inline] + fn validate_inputs_internal( + pressure: AtmosphericPressure, + vapour_pressure: VapourPressure, + ) -> Result<(), InputError> { + pressure.check_range_si(100.0, 150_000.0)?; + vapour_pressure.check_range_si(0.0, 50_000.0)?; + + if vapour_pressure.0 > pressure.0 { + return Err(InputError::IncorrectArgumentSet( + "vapour_pressure cannot be greater than pressure", + )); + } + + if approx_eq!( + Float, + pressure.get_si_value(), + vapour_pressure.get_si_value(), + ulps = 2 + ) { + return Err(InputError::IncorrectArgumentSet( + "pressure and vapour_pressure cannot be equal", + )); + } + Ok(()) + } + + #[inline] + fn compute_unchecked( + pressure: AtmosphericPressure, + vapour_pressure: VapourPressure, + ) -> MixingRatio { + MixingRatio(EPSILON * (vapour_pressure.0 / (pressure.0 - vapour_pressure.0))) + } +} + +#[cfg(test)] +mod tests { + + use crate::tests::{Argument, test_with_2args, testing_traits::ReferenceAtmosphere}; + + use super::*; + + #[test] + fn definition1() { + test_with_2args::( + Argument::new([100.0, 150_000.0]), + Argument::new([0.0, 50_000.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } +} diff --git a/src/formulas/potential_temperature.rs b/src/formulas/potential_temperature.rs new file mode 100644 index 0000000..9a3e861 --- /dev/null +++ b/src/formulas/potential_temperature.rs @@ -0,0 +1,110 @@ +//! Functions to calculate potential temperature of dry air +//! +//! The temperature that an unsaturated parcel of dry air would have if brought +//! adiabatically and reversibly from its initial state to a +//! standard pressure, p0 = 100 kPa ([AMETSOC Glossary](https://glossary.ametsoc.org/wiki/Potential_temperature)). + +use crate::constants::{KAPPA, REF_PRES}; +use crate::errors::InputError; +use crate::formulas::Formula3; +use crate::quantities::{ + AtmosphericPressure, DryBulbTemperature, PotentialTemperature, QuantityHelpers, VapourPressure, +}; +use crate::Float; +use float_cmp::approx_eq; +use uom::si::pressure::pascal; +use uom::si::ratio::ratio; +use uom::si::thermodynamic_temperature::kelvin; + +type FormulaQuantity = PotentialTemperature; + +/// Formula for computing potential temperature of dry air from temperature, pressure and vapour pressure. +/// +/// Provided in by R. Davies-Jones (2009) [(doi:10.1175/2009MWR2774.1)](https://doi.org/10.1175/2009MWR2774.1) +/// +/// Valid `temperature` range: 253K - 324K +/// +/// Valid `pressure` range: 100Pa - 150000Pa +/// +/// Valid `vapour_pressure` range: 0Pa - 10000Pa +/// +/// Returns [`InputError::IncorrectArgumentSet`] when `pressure` and `vapour_pressure` are equal, +/// in which case division by 0 occurs. +/// +/// Returns [`InputError::IncorrectArgumentSet`] when `pressure` is lower than `vapour_pressure`, +/// in which case floating-point exponentation of negative number occurs. +pub struct Definition1; + +impl Formula3 + for Definition1 +{ + #[inline] + fn validate_inputs_internal( + temperature: DryBulbTemperature, + pressure: AtmosphericPressure, + vapour_pressure: VapourPressure, + ) -> Result<(), InputError> { + temperature.check_range_si(253.0, 324.0)?; + pressure.check_range_si(100.0, 150_000.0)?; + vapour_pressure.check_range_si(0.0, 10_000.0)?; + + if approx_eq!( + Float, + pressure.get_si_value(), + vapour_pressure.get_si_value(), + ulps = 2 + ) { + return Err(InputError::IncorrectArgumentSet( + "pressure and vapour_pressure cannot be equal", + )); + } + + if vapour_pressure.0 > pressure.0 { + return Err(InputError::IncorrectArgumentSet( + "vapour_pressure cannot be greater or equal to pressure", + )); + } + + Ok(()) + } + + #[inline] + fn compute_unchecked( + temperature: DryBulbTemperature, + pressure: AtmosphericPressure, + vapour_pressure: VapourPressure, + ) -> PotentialTemperature { + let temperature = temperature.0.get::(); + let pressure = pressure.0.get::(); + let vapour_pressure = vapour_pressure.0.get::(); + let p0 = REF_PRES.get::(); + + let kappa = KAPPA.get::(); + let result = temperature * (p0 / (pressure - vapour_pressure)).powf(kappa); + + PotentialTemperature::new::(result) + } +} + +#[cfg(test)] +mod tests { + use crate::tests::{test_with_3args, testing_traits::ReferenceAtmosphere, Argument}; + + use super::*; + #[test] + fn definition1() { + test_with_3args::< + FormulaQuantity, + DryBulbTemperature, + AtmosphericPressure, + VapourPressure, + Definition1, + >( + Argument::new([253.0, 324.0]), + Argument::new([100.0, 150_000.0]), + Argument::new([0.0, 10_000.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } +} diff --git a/src/formulas/relative_humidity.rs b/src/formulas/relative_humidity.rs new file mode 100644 index 0000000..cfc21fa --- /dev/null +++ b/src/formulas/relative_humidity.rs @@ -0,0 +1,98 @@ +//! Functions to calculate relative humidity + +use crate::errors::InputError; +use crate::formulas::Formula2; +use crate::quantities::{ + MixingRatio, QuantityHelpers, RelativeHumidity, SaturationMixingRatio, + SaturationVapourPressure, VapourPressure, +}; + +type FormulaQuantity = RelativeHumidity; + +/// Formula for computing relative humidity from mixing ratio and saturation mixing ratio. +/// Can be used interchangeably with [`general2`]. +/// +/// By the definition of mixing ratio, this formula is mathematically equivalent of +/// formula used in [`general2`]. +/// +/// Valid `mixing_ratio` range: 0.00001 - 10.0 +/// +/// Valid `saturation_mixing_ratio` range: 0.00001 - 10.0 +pub struct Definition1; + +impl Formula2 for Definition1 { + #[inline] + fn validate_inputs_internal( + mixing_ratio: MixingRatio, + saturation_mixing_ratio: SaturationMixingRatio, + ) -> Result<(), InputError> { + mixing_ratio.check_range_si(0.00001, 10.0)?; + saturation_mixing_ratio.check_range_si(0.00001, 10.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + mixing_ratio: MixingRatio, + saturation_mixing_ratio: SaturationMixingRatio, + ) -> RelativeHumidity { + RelativeHumidity(mixing_ratio.0 / saturation_mixing_ratio.0) + } +} + +/// Formula for computing relative humidity from vapour pressure and saturation vapour pressure. +/// Can be used interchangeably with [`general1`]. +/// +/// Valid `vapour_pressure` range: 0Pa - 50000Pa +/// +/// Valid `saturation_vapour_pressure` range: 0Pa - 50000Pa +pub struct Definition2; + +impl Formula2 for Definition2 { + #[inline] + fn validate_inputs_internal( + vapour_pressure: VapourPressure, + saturation_vapour_pressure: SaturationVapourPressure, + ) -> Result<(), InputError> { + vapour_pressure.check_range_si(0.0, 50_000.0)?; + saturation_vapour_pressure.check_range_si(0.1, 50_000.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + vapour_pressure: VapourPressure, + saturation_vapour_pressure: SaturationVapourPressure, + ) -> RelativeHumidity { + RelativeHumidity(vapour_pressure.0 / saturation_vapour_pressure.0) + } +} + +#[cfg(test)] +mod tests { + use crate::tests::{test_with_2args, testing_traits::ReferenceAtmosphere, Argument}; + + use super::*; + + #[test] + fn definition1() { + test_with_2args::( + Argument::new([0.00001, 10.0]), + Argument::new([0.00001, 10.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } + + #[test] + fn definition2() { + test_with_2args::( + Argument::new([0.0, 50_000.0]), + Argument::new([0.1, 50_000.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } +} diff --git a/src/formulas/saturation_mixing_ratio.rs b/src/formulas/saturation_mixing_ratio.rs new file mode 100644 index 0000000..c76eed2 --- /dev/null +++ b/src/formulas/saturation_mixing_ratio.rs @@ -0,0 +1,124 @@ +//! Functions to calculate saturation mixing ratio of unsaturated air +//! +//! Saturation mixing ration is the value of the mixing ratio of saturated air at the +//! given temperature and pressure ([AMETSOC Glossary](https://glossary.ametsoc.org/wiki/Saturation_mixing_ratio)). + +use crate::Float; +use crate::formulas::Formula2; +use crate::quantities::{ + AtmosphericPressure, MixingRatio, QuantityHelpers, RelativeHumidity, SaturationMixingRatio, + SaturationVapourPressure, +}; +use crate::{constants::EPSILON, errors::InputError}; +use float_cmp::approx_eq; + +type FormulaQuantity = SaturationMixingRatio; + +/// Formula for computing saturation mixing ratio of unsaturated air from air pressure and vapour pressure +/// +/// Valid `pressure` range: 100Pa - 150000Pa +/// +/// Valid `saturation_vapour_pressure` range: 0Pa - 10000Pa +/// +/// Returns [`InputError::IncorrectArgumentSet`] when inputs are equal and division by 0 would occur. +pub struct Definition1; + +impl Formula2 for Definition1 { + #[inline] + fn validate_inputs_internal( + pressure: AtmosphericPressure, + saturation_vapour_pressure: SaturationVapourPressure, + ) -> Result<(), InputError> { + pressure.check_range_si(100.0, 150_000.0)?; + saturation_vapour_pressure.check_range_si(0.0, 50_000.0)?; + + if saturation_vapour_pressure.0 > pressure.0 { + return Err(InputError::IncorrectArgumentSet( + "saturation_vapour_pressure cannot be greater than pressure", + )); + } + + if approx_eq!( + Float, + pressure.get_si_value(), + saturation_vapour_pressure.get_si_value(), + ulps = 2 + ) { + return Err(InputError::IncorrectArgumentSet( + "pressure and saturation_vapour_pressure cannot be equal", + )); + } + Ok(()) + } + + #[inline] + fn compute_unchecked( + pressure: AtmosphericPressure, + saturation_vapour_pressure: SaturationVapourPressure, + ) -> SaturationMixingRatio { + SaturationMixingRatio( + EPSILON * (saturation_vapour_pressure.0 / (pressure.0 - saturation_vapour_pressure.0)), + ) + } +} + +/// Formula for computing saturation mixing ratio of unsaturated air from +/// mixing ratio and relative humidity. +/// +/// Valid `mixing_ratio` range: 0.000_000_000_1 - 1.0 +/// +/// Valid `relative_humditity` range: 0.000_000_000_1 - 2.0 +pub struct Definition2; + +impl Formula2 for Definition2 { + #[inline] + fn validate_inputs_internal( + mixing_ratio: MixingRatio, + relative_humidity: RelativeHumidity, + ) -> Result<(), InputError> { + mixing_ratio.check_range_si(0.000_000_000_1, 1.0)?; + relative_humidity.check_range_si(0.000_000_000_1, 2.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + mixing_ratio: MixingRatio, + relative_humidity: RelativeHumidity, + ) -> SaturationMixingRatio { + SaturationMixingRatio(mixing_ratio.0 / relative_humidity.0) + } +} + +#[cfg(test)] +mod tests { + use crate::tests::{Argument, test_with_2args, testing_traits::ReferenceAtmosphere}; + + use super::*; + + #[test] + fn definition1() { + test_with_2args::< + FormulaQuantity, + AtmosphericPressure, + SaturationVapourPressure, + Definition1, + >( + Argument::new([100.0, 150_000.0]), + Argument::new([0.0, 50_000.0]), + ReferenceAtmosphere::Normal, + 1e-2, + ); + } + + #[test] + fn definition2() { + test_with_2args::( + Argument::new([0.000_000_000_1, 1.0]), + Argument::new([0.000_000_000_1, 2.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } +} diff --git a/src/formulas/saturation_vapour_pressure.rs b/src/formulas/saturation_vapour_pressure.rs new file mode 100644 index 0000000..1ba141f --- /dev/null +++ b/src/formulas/saturation_vapour_pressure.rs @@ -0,0 +1,531 @@ +//! Formulae to calculate saturation vapour pressure of the unsaturated air +//! +//! The vapor pressure of a system, at a given temperature, for which the vapor +//! of a substance is in equilibrium with a plane surface of that substance's pure +//! liquid or solid phase; that is, the vapor pressure of a system that has attained +//! saturation but not supersaturation ([AMETSOC Glossary](https://glossary.ametsoc.org/wiki/Saturation_vapor_pressure)). + +use crate::Float; +use crate::Storage::Pressure; +use crate::errors::InputError; +use crate::formulas::{Formula1, Formula2}; +use crate::quantities::{ + AtmosphericPressure, DryBulbTemperature, QuantityHelpers, RelativeHumidity, + SaturationVapourPressure, VapourPressure, +}; + +use uom::si::pressure::{hectopascal, kilopascal, pascal}; +use uom::si::thermodynamic_temperature::{degree_celsius, kelvin}; + +type FormulaQuantity = SaturationVapourPressure; + +/// Formula for computing saturation vapour pressure from vapour pressure and relative humidity. +/// +/// Valid `vapour_pressure` range: 0Pa - 50000Pa +/// +/// Valid `relative_humidity` range: 0.00001 - 2.0 +pub struct Definition1; + +impl Formula2 for Definition1 { + #[inline] + fn validate_inputs_internal( + vapour_pressure: VapourPressure, + relative_humidity: RelativeHumidity, + ) -> Result<(), InputError> { + vapour_pressure.check_range_si(0.0, 50_000.0)?; + relative_humidity.check_range_si(0.00001, 2.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + vapour_pressure: VapourPressure, + relative_humidity: RelativeHumidity, + ) -> SaturationVapourPressure { + SaturationVapourPressure(vapour_pressure.0 / relative_humidity.0) + } +} + +/// Formula for saturation computing saturation vapour pressure from dewpoint temperature and pressure. +/// Should be used for air over water when accuracy is desired. +/// +/// Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). +/// +/// Valid `dewpoint temperature` range: 232K - 324K +/// +/// Valid `atmospheric pressure` range: 100Pa - 150000Pa +pub struct Buck1; + +impl Formula2 for Buck1 { + #[inline] + fn validate_inputs_internal( + temperature: DryBulbTemperature, + pressure: AtmosphericPressure, + ) -> Result<(), InputError> { + temperature.check_range_si(232.0, 324.0)?; + pressure.check_range_si(100.0, 150_000.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + temperature: DryBulbTemperature, + pressure: AtmosphericPressure, + ) -> SaturationVapourPressure { + let dewpoint = temperature.0.get::(); + let pressure = pressure.0.get::(); + + let lower_a = 6.1121; + let lower_b = 18.729; + let lower_c = 257.87; + let lower_d = 227.3; + + let upper_a = 0.000_72; + let upper_b = 0.000_003_2; + let upper_c = 0.000_000_000_59; + + let lower_e = + lower_a * (((lower_b - (dewpoint / lower_d)) * dewpoint) / (dewpoint + lower_c)).exp(); + let lower_f = 1.0 + upper_a + (pressure * (upper_b + (upper_c * dewpoint * dewpoint))); + + let result = Pressure::new::(lower_e * lower_f); + + SaturationVapourPressure(result) + } +} + +/// Formula for computing saturation vapour pressure from dewpoint temperature and pressure. +/// Should be used for air over ice when accuracy is desired. +/// +/// Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). +/// +/// Valid `dewpoint` range: 193K - 274K +/// +/// Valid `pressure` range: 100Pa - 150000Pa +pub struct Buck2; + +impl Formula2 for Buck2 { + #[inline] + fn validate_inputs_internal( + temperature: DryBulbTemperature, + pressure: AtmosphericPressure, + ) -> Result<(), InputError> { + temperature.check_range_si(193.0, 274.0)?; + pressure.check_range_si(100.0, 150_000.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + temperature: DryBulbTemperature, + pressure: AtmosphericPressure, + ) -> SaturationVapourPressure { + let dewpoint = temperature.0.get::(); + let pressure = pressure.0.get::(); + + let lower_a = 6.1115; + let lower_b = 23.036; + let lower_c = 279.82; + let lower_d = 333.7; + + let upper_a = 0.000_22; + let upper_b = 0.000_003_83; + let upper_c = 0.000_000_000_64; + + let lower_e = + lower_a * (((lower_b - (dewpoint / lower_d)) * dewpoint) / (dewpoint + lower_c)).exp(); + let lower_f = 1.0 + upper_a + (pressure * (upper_b + (upper_c * dewpoint * dewpoint))); + + let result = Pressure::new::(lower_e * lower_f); + + SaturationVapourPressure(result) + } +} + +/// Formula for computing saturation vapour pressure from dewpoint temperature and pressure. +/// Should be used for air over water for general use. +/// +/// Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). +/// +/// Valid `dewpoint` range: 253K - 324K +/// +/// Valid `pressure` range: 100Pa - 150000Pa +pub struct Buck3; + +impl Formula2 for Buck3 { + #[inline] + fn validate_inputs_internal( + temperature: DryBulbTemperature, + pressure: AtmosphericPressure, + ) -> Result<(), InputError> { + temperature.check_range_si(253.0, 324.0)?; + pressure.check_range_si(100.0, 150_000.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + temperature: DryBulbTemperature, + pressure: AtmosphericPressure, + ) -> SaturationVapourPressure { + let dewpoint = temperature.0.get::(); + let pressure = pressure.0.get::(); + + let lower_a = 6.1121; + let lower_b = 17.502; + let lower_c = 240.97; + + let upper_a = 0.000_7; + let upper_b = 0.000_003_46; + + let lower_e = lower_a * ((lower_b * dewpoint) / (dewpoint + lower_c)).exp(); + let lower_f = 1.0 + upper_a + (pressure * upper_b); + + let result = Pressure::new::(lower_e * lower_f); + + SaturationVapourPressure(result) + } +} + +/// Formula for computing saturation vapour pressure from dewpoint temperature. +/// Simplified version of [`buck3`]. Very popular in meteorological sources. +/// +/// Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). +/// +/// Valid `dewpoint` range: 253K - 324K +pub struct Buck3Simplified; + +impl Formula1 for Buck3Simplified { + #[inline] + fn validate_inputs_internal(temperature: DryBulbTemperature) -> Result<(), InputError> { + temperature.check_range_si(253.0, 324.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked(temperature: DryBulbTemperature) -> SaturationVapourPressure { + let dewpoint = temperature.0.get::(); + + let lower_a = 6.1121; + let lower_b = 17.502; + let lower_c = 240.97; + + let lower_e = lower_a * ((lower_b * dewpoint) / (dewpoint + lower_c)).exp(); + + let result = Pressure::new::(lower_e); + + SaturationVapourPressure(result) + } +} + +/// Formula for computing vapour saturation pressure from dewpoint temperature and pressure. +/// Should be used for air over ice for general use. +/// +/// Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). +/// +/// Valid `dewpoint` range: 223K - 274K +/// +/// Valid `pressure` range: 100Pa - 150000Pa +pub struct Buck4; + +impl Formula2 for Buck4 { + #[inline] + fn validate_inputs_internal( + temperature: DryBulbTemperature, + pressure: AtmosphericPressure, + ) -> Result<(), InputError> { + temperature.check_range_si(223.0, 274.0)?; + pressure.check_range_si(100.0, 150_000.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + temperature: DryBulbTemperature, + pressure: AtmosphericPressure, + ) -> SaturationVapourPressure { + let dewpoint = temperature.0.get::(); + let pressure = pressure.0.get::(); + + let lower_a = 6.1115; + let lower_b = 22.452; + let lower_c = 272.55; + + let upper_a = 0.000_3; + let upper_b = 0.000_004_18; + + let lower_e = lower_a * ((lower_b * dewpoint) / (dewpoint + lower_c)).exp(); + let lower_f = 1.0 + upper_a + (pressure * upper_b); + + let result = Pressure::new::(lower_e * lower_f); + + SaturationVapourPressure(result) + } +} + +/// Formula for computing saturation vapour pressure from dewpoint temperature. +/// Simplified version of [`buck4`], analogical to [`buck3_simplified`]. +/// +/// Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). +/// +/// Valid `dewpoint` range: 223K - 274K +pub struct Buck4Simplified; + +impl Formula1 for Buck4Simplified { + #[inline] + fn validate_inputs_internal(temperature: DryBulbTemperature) -> Result<(), InputError> { + temperature.check_range_si(223.0, 274.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked(temperature: DryBulbTemperature) -> SaturationVapourPressure { + let dewpoint = temperature.0.get::(); + + let lower_a = 6.1115; + let lower_b = 22.452; + let lower_c = 272.55; + + let lower_e = lower_a * ((lower_b * dewpoint) / (dewpoint + lower_c)).exp(); + + let result = Pressure::new::(lower_e); + + SaturationVapourPressure(result) + } +} + +/// Formula for computing saturation vapour pressure over water from dewpoint temperature. +/// Should be used for temperatures above 273K. +/// +/// Derived by O. Tetens (1930). +/// +/// Valid `dewpoint` range: 273K - 353K +pub struct Tetens1; + +impl Formula1 for Tetens1 { + #[inline] + fn validate_inputs_internal(temperature: DryBulbTemperature) -> Result<(), InputError> { + temperature.check_range_si(273.0, 353.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked(temperature: DryBulbTemperature) -> SaturationVapourPressure { + let dewpoint = temperature.0.get::(); + + let lower_a = 0.61078; + let lower_b = 17.27; + let lower_c = 237.3; + + let result = lower_a * ((lower_b * dewpoint) / (dewpoint + lower_c)).exp(); + + let result = Pressure::new::(result); + + SaturationVapourPressure(result) + } +} + +/// Formula for computing saturation vapour pressure over water from dewpoint temperature. +/// Should be used when accuracy is required as it is +/// computationally expensive. +/// +/// Derived by A. Wexler (1976) [(doi: 10.6028/jres.080A.071)](https://dx.doi.org/10.6028%2Fjres.080A.071). +/// +/// Valid `dewpoint` range: 273K - 374K +pub struct Wexler1; + +impl Formula1 for Wexler1 { + #[inline] + fn validate_inputs_internal(temperature: DryBulbTemperature) -> Result<(), InputError> { + temperature.check_range_si(273.0, 374.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked(temperature: DryBulbTemperature) -> SaturationVapourPressure { + let dewpoint = temperature.get_si_value(); + + // constants from the paper + let g: [Float; 8] = [ + -2991.2729, + -6017.0128, + 18.876_438_54, + -0.028_354_721, + 0.000_017_838_3, + -0.000_000_000_841_504_17, + 0.000_000_000_000_444_125_43, + 2.858_487, + ]; + + let ln_p = g + .iter() + .enumerate() + .take(7) + .fold(g[7] * dewpoint.ln(), |acc, (i, &x)| { + acc + x * dewpoint.powi(i as i32 - 2) + }); + + let result = Pressure::new::(ln_p.exp()); + + SaturationVapourPressure(result) + } +} + +/// Formula for computing saturation vapour over ice pressure from dewpoint temperature. +/// Should be used when accuracy is required as it is +/// computationally expensive. +/// +/// Derived by A. Wexler (1977) [(doi: 10.6028/jres.081A.003)](https://dx.doi.org/10.6028%2Fjres.081A.003). +/// +/// Valid `dewpoint` range: 173K - 274K +pub struct Wexler2; + +impl Formula1 for Wexler2 { + #[inline] + fn validate_inputs_internal(temperature: DryBulbTemperature) -> Result<(), InputError> { + temperature.check_range_si(173.0, 274.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked(temperature: DryBulbTemperature) -> SaturationVapourPressure { + let dewpoint = temperature.0.get::(); + + // constants from the paper + let big_k: [Float; 6] = [ + -5865.3696, + 22.241_033, + 0.013_749_042, + -0.000_034_031_77, + 0.000_000_026_967_687, + 0.691_865_1, + ]; + + let ln_p = big_k + .iter() + .enumerate() + .take(5) + .fold(big_k[5] * dewpoint.ln(), |acc, (i, &x)| { + acc + x * dewpoint.powi(i as i32 - 1) + }); + + let result = Pressure::new::(ln_p.exp()); + + SaturationVapourPressure(result) + } +} + +#[cfg(test)] +mod tests { + use crate::{ + quantities::{AtmosphericPressure, DryBulbTemperature, RelativeHumidity, VapourPressure}, + tests::{Argument, test_with_1arg, test_with_2args, testing_traits::ReferenceAtmosphere}, + }; + + use super::*; + + #[test] + fn definition1() { + test_with_2args::( + Argument::new([0.0, 50_000.0]), + Argument::new([0.00001, 2.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } + + #[test] + fn buck1() { + test_with_2args::( + Argument::new([232.0, 324.0]), + Argument::new([100.0, 150_000.0]), + ReferenceAtmosphere::Normal, + 1e2, + ); + } + + #[test] + fn buck2() { + test_with_2args::( + Argument::new([193.0, 274.0]), + Argument::new([100.0, 150_000.0]), + ReferenceAtmosphere::Freezing, + 1e0, + ); + } + + #[test] + fn buck3() { + test_with_2args::( + Argument::new([253.0, 324.0]), + Argument::new([100.0, 150_000.0]), + ReferenceAtmosphere::Normal, + 1e2, + ); + } + + #[test] + fn buck4() { + test_with_2args::( + Argument::new([223.0, 274.0]), + Argument::new([100.0, 150_000.0]), + ReferenceAtmosphere::Freezing, + 1e0, + ); + } + + #[test] + fn buck3_simplified() { + test_with_1arg::( + Argument::new([253.0, 324.0]), + ReferenceAtmosphere::Normal, + 1e1, + ); + } + + #[test] + fn buck4_simplified() { + test_with_1arg::( + Argument::new([223.0, 274.0]), + ReferenceAtmosphere::Freezing, + 1e-1, + ); + } + + #[test] + fn tetens1() { + test_with_1arg::( + Argument::new([273.0, 353.0]), + ReferenceAtmosphere::Normal, + 1e1, + ); + } + + #[test] + fn wexler1() { + test_with_1arg::( + Argument::new([273.0, 374.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } + + #[test] + fn wexler2() { + test_with_1arg::( + Argument::new([173.0, 274.0]), + ReferenceAtmosphere::Freezing, + 1e-12, + ); + } +} diff --git a/src/formulas/specific_humidity.rs b/src/formulas/specific_humidity.rs new file mode 100644 index 0000000..f2bd48b --- /dev/null +++ b/src/formulas/specific_humidity.rs @@ -0,0 +1,65 @@ +//! Functions to calculate specific humidity of air +//! +//! Specific humidity (or moisture content) is the ratio of the mass +//! of water vapor to the total mass of the air parcel [Wikipedia](https://en.wikipedia.org/wiki/Humidity#Specific_humidity). +//! +//! Specific humidity is approximately equal to mixing ratio. + +use crate::constants::DIMLESS_ONE; +use crate::formulas::Formula2; +use crate::quantities::{AtmosphericPressure, QuantityHelpers, SpecificHumidity, VapourPressure}; +use crate::{constants::EPSILON, errors::InputError}; + +type FormulaQuantity = SpecificHumidity; + +/// Formula for computing specific humidity from vapour pressure and pressure. +/// Reverse function of [`vapour_pressure::general1`](crate::vapour_pressure::general1). +/// This function is theoretical not empirical. +/// +/// Provided by [Rogers & Yau (1989)](https://www.elsevier.com/books/a-short-course-in-cloud-physics/yau/978-0-08-057094-5). +/// +/// Valid `vapour_pressure` range: 0Pa - 50000OPa +/// +/// Valid `pressure` range: 100Pa - 150000Pa +pub struct Definition1; + +impl Formula2 for Definition1 { + #[inline] + fn validate_inputs_internal( + vapour_pressure: VapourPressure, + pressure: AtmosphericPressure, + ) -> Result<(), InputError> { + vapour_pressure.check_range_si(0.0, 50_000.0)?; + pressure.check_range_si(100.0, 150_000.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + vapour_pressure: VapourPressure, + pressure: AtmosphericPressure, + ) -> SpecificHumidity { + let result = EPSILON + * (vapour_pressure.0 / (pressure.0 - (vapour_pressure.0 * (DIMLESS_ONE - EPSILON)))); + + SpecificHumidity(result) + } +} + +#[cfg(test)] +mod tests { + use crate::tests::{test_with_2args, testing_traits::ReferenceAtmosphere, Argument}; + + use super::*; + + #[test] + fn definition1() { + test_with_2args::( + Argument::new([0.0, 50_000.0]), + Argument::new([100.0, 150_000.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } +} diff --git a/src/formulas/traits.rs b/src/formulas/traits.rs new file mode 100644 index 0000000..3c3a9dd --- /dev/null +++ b/src/formulas/traits.rs @@ -0,0 +1,411 @@ +#![allow(missing_docs)] +#![allow(clippy::missing_errors_doc)] + +use crate::{errors::InputError, quantities::ThermodynamicQuantity}; +#[cfg(feature = "array")] +use ndarray::{Array, ArrayView, Dimension, FoldWhile, Zip}; +#[cfg(feature = "parallel")] +use rayon::iter::{IndexedParallelIterator, IntoParallelIterator, ParallelIterator}; +#[cfg(feature = "argdebug")] +use tracing::instrument; + +pub trait Formula1 { + fn compute_unchecked(i1: I1) -> O; + + fn validate_inputs_internal(i1: I1) -> Result<(), InputError>; + + /// This function exist only to have tracing available + /// Hopefully compiler optimises it properly + #[inline] + #[cfg_attr(feature = "argdebug", instrument(level = "trace"))] + fn validate_inputs(i1: I1) -> Result<(), InputError> { + Self::validate_inputs_internal(i1) + } + + #[inline] + fn compute(i1: I1) -> Result { + Self::validate_inputs(i1)?; + Ok(Self::compute_unchecked(i1)) + } + + fn compute_vec(i1: &[I1]) -> Result, InputError> { + i1.iter().map(|&i1| Self::compute(i1)).collect() + } + + #[cfg(feature = "array")] + fn compute_ndarray<'a, D: Dimension + Copy, A: Into>>( + i1: A, + ) -> Result, InputError> + where + I1: 'a, + { + let i1: ArrayView = i1.into(); + + Zip::from(i1) + .fold_while(Ok(()), |_, &i1| match Self::validate_inputs(i1) { + Ok(_) => FoldWhile::Continue(Ok(())), + Err(e) => FoldWhile::Done(Err(e)), + }) + .into_inner()?; + + Ok(Zip::from(i1).map_collect(|&i1| Self::compute_unchecked(i1))) + } + + #[cfg(feature = "parallel")] + fn compute_vec_parallel(i1: &[I1]) -> Result, InputError> { + i1.into_par_iter().map(|&i1| Self::compute(i1)).collect() + } + + #[cfg(feature = "parallel")] + fn compute_ndarray_parallel<'a, D: Dimension + Copy, A: Into>>( + i1: A, + ) -> Result, InputError> + where + I1: 'a, + { + let i1: ArrayView = i1.into(); + + Zip::from(i1).par_fold( + || Ok(()), + |_, &i1| Self::validate_inputs(i1), + |a, b| a.and(b), + )?; + + Ok(Zip::from(i1).par_map_collect(|&a| Self::compute_unchecked(a))) + } +} + +pub trait Formula2 { + fn compute_unchecked(i1: I1, i2: I2) -> O; + + fn validate_inputs_internal(i1: I1, i2: I2) -> Result<(), InputError>; + + #[inline] + #[cfg_attr(feature = "argdebug", instrument(level = "trace"))] + fn validate_inputs(i1: I1, i2: I2) -> Result<(), InputError> { + Self::validate_inputs_internal(i1, i2) + } + + #[inline] + fn compute(i1: I1, i2: I2) -> Result { + Self::validate_inputs(i1, i2)?; + Ok(Self::compute_unchecked(i1, i2)) + } + + fn compute_vec(i1: &[I1], i2: &[I2]) -> Result, InputError> { + i1.iter() + .zip(i2.iter()) + .map(|(&i1, &i2)| Self::compute(i1, i2)) + .collect() + } + + #[cfg(feature = "array")] + fn compute_ndarray< + 'a, + D: Dimension + Copy, + A1: Into>, + A2: Into>, + >( + i1: A1, + i2: A2, + ) -> Result, InputError> + where + I1: 'a, + I2: 'a, + { + let i1: ArrayView = i1.into(); + let i2: ArrayView = i2.into(); + + Zip::from(i1) + .and(i2) + .fold_while(Ok(()), |_, &i1, &i2| match Self::validate_inputs(i1, i2) { + Ok(_) => FoldWhile::Continue(Ok(())), + Err(e) => FoldWhile::Done(Err(e)), + }) + .into_inner()?; + + Ok(Zip::from(i1) + .and(i2) + .map_collect(|&i1, &i2| Self::compute_unchecked(i1, i2))) + } + + #[cfg(feature = "parallel")] + fn compute_vec_parallel(i1: &[I1], i2: &[I2]) -> Result, InputError> { + i1.into_par_iter() + .zip(i2) + .map(|(&i1, &i2)| Self::compute(i1, i2)) + .collect() + } + + #[cfg(feature = "parallel")] + fn compute_ndarray_parallel< + 'a, + D: Dimension + Copy, + A1: Into>, + A2: Into>, + >( + i1: A1, + i2: A2, + ) -> Result, InputError> + where + I1: 'a, + I2: 'a, + { + let i1: ArrayView = i1.into(); + let i2: ArrayView = i2.into(); + + Zip::from(i1).and(i2).par_fold( + || Ok(()), + |_, &i1, &i2| Self::validate_inputs(i1, i2), + |a, b| a.and(b), + )?; + + Ok(Zip::from(i1) + .and(i2) + .par_map_collect(|&i1, &i2| Self::compute_unchecked(i1, i2))) + } +} + +pub trait Formula3< + O: ThermodynamicQuantity, + I1: ThermodynamicQuantity, + I2: ThermodynamicQuantity, + I3: ThermodynamicQuantity, +> +{ + fn compute_unchecked(i1: I1, i2: I2, i3: I3) -> O; + + fn validate_inputs_internal(i1: I1, i2: I2, i3: I3) -> Result<(), InputError>; + + #[inline] + #[cfg_attr(feature = "argdebug", instrument(level = "trace"))] + fn validate_inputs(i1: I1, i2: I2, i3: I3) -> Result<(), InputError> { + Self::validate_inputs_internal(i1, i2, i3) + } + + #[inline] + fn compute(i1: I1, i2: I2, i3: I3) -> Result { + Self::validate_inputs(i1, i2, i3)?; + Ok(Self::compute_unchecked(i1, i2, i3)) + } + + fn compute_vec(i1: &[I1], i2: &[I2], i3: &[I3]) -> Result, InputError> { + i1.iter() + .zip(i2.iter()) + .zip(i3.iter()) + .map(|((&i1, &i2), &i3)| Self::compute(i1, i2, i3)) + .collect() + } + + #[cfg(feature = "array")] + fn compute_ndarray< + 'a, + D: Dimension + Copy, + A1: Into>, + A2: Into>, + A3: Into>, + >( + i1: A1, + i2: A2, + i3: A3, + ) -> Result, InputError> + where + I1: 'a, + I2: 'a, + I3: 'a, + { + let i1: ArrayView = i1.into(); + let i2: ArrayView = i2.into(); + let i3: ArrayView = i3.into(); + + Zip::from(i1) + .and(i2) + .and(i3) + .fold_while(Ok(()), |_, &i1, &i2, &i3| { + match Self::validate_inputs(i1, i2, i3) { + Ok(_) => FoldWhile::Continue(Ok(())), + Err(e) => FoldWhile::Done(Err(e)), + } + }) + .into_inner()?; + + Ok(Zip::from(i1) + .and(i2) + .and(i3) + .map_collect(|&i1, &i2, &i3| Self::compute_unchecked(i1, i2, i3))) + } + + #[cfg(feature = "parallel")] + fn compute_vec_parallel(i1: &[I1], i2: &[I2], i3: &[I3]) -> Result, InputError> { + i1.into_par_iter() + .zip(i2) + .zip(i3) + .map(|((&i1, &i2), &i3)| Self::compute(i1, i2, i3)) + .collect() + } + + #[cfg(feature = "parallel")] + fn compute_ndarray_parallel< + 'a, + D: Dimension + Copy, + A1: Into>, + A2: Into>, + A3: Into>, + >( + i1: A1, + i2: A2, + i3: A3, + ) -> Result, InputError> + where + I1: 'a, + I2: 'a, + I3: 'a, + { + let i1: ArrayView = i1.into(); + let i2: ArrayView = i2.into(); + let i3: ArrayView = i3.into(); + + Zip::from(i1).and(i2).and(i3).par_fold( + || Ok(()), + |_, &i1, &i2, &i3| Self::validate_inputs(i1, i2, i3), + |a, b| a.and(b), + )?; + + Ok(Zip::from(i1) + .and(i2) + .and(i3) + .par_map_collect(|&i1, &i2, &i3| Self::compute_unchecked(i1, i2, i3))) + } +} + +pub trait Formula4< + O: ThermodynamicQuantity, + I1: ThermodynamicQuantity, + I2: ThermodynamicQuantity, + I3: ThermodynamicQuantity, + I4: ThermodynamicQuantity, +> +{ + fn compute_unchecked(i1: I1, i2: I2, i3: I3, i4: I4) -> O; + + fn validate_inputs_internal(i1: I1, i2: I2, i3: I3, i4: I4) -> Result<(), InputError>; + + #[inline] + #[cfg_attr(feature = "argdebug", instrument(level = "trace"))] + fn validate_inputs(i1: I1, i2: I2, i3: I3, i4: I4) -> Result<(), InputError> { + Self::validate_inputs_internal(i1, i2, i3, i4) + } + + #[inline] + fn compute(i1: I1, i2: I2, i3: I3, i4: I4) -> Result { + Self::validate_inputs(i1, i2, i3, i4)?; + Ok(Self::compute_unchecked(i1, i2, i3, i4)) + } + + fn compute_vec(i1: &[I1], i2: &[I2], i3: &[I3], i4: &[I4]) -> Result, InputError> { + i1.iter() + .zip(i2.iter()) + .zip(i3.iter()) + .zip(i4.iter()) + .map(|(((&i1, &i2), &i3), &i4)| Self::compute(i1, i2, i3, i4)) + .collect() + } + + #[cfg(feature = "array")] + fn compute_ndarray< + 'a, + D: Dimension + Copy, + A1: Into>, + A2: Into>, + A3: Into>, + A4: Into>, + >( + i1: A1, + i2: A2, + i3: A3, + i4: A4, + ) -> Result, InputError> + where + I1: 'a, + I2: 'a, + I3: 'a, + I4: 'a, + { + let i1: ArrayView = i1.into(); + let i2: ArrayView = i2.into(); + let i3: ArrayView = i3.into(); + let i4: ArrayView = i4.into(); + + Zip::from(i1) + .and(i2) + .and(i3) + .and(i4) + .fold_while( + Ok(()), + |_, &i1, &i2, &i3, &i4| match Self::validate_inputs(i1, i2, i3, i4) { + Ok(_) => FoldWhile::Continue(Ok(())), + Err(e) => FoldWhile::Done(Err(e)), + }, + ) + .into_inner()?; + + Ok(Zip::from(i1) + .and(i2) + .and(i3) + .and(i4) + .map_collect(|&i1, &i2, &i3, &i4| Self::compute_unchecked(i1, i2, i3, i4))) + } + + #[cfg(feature = "parallel")] + fn compute_vec_parallel( + i1: &[I1], + i2: &[I2], + i3: &[I3], + i4: &[I4], + ) -> Result, InputError> { + i1.into_par_iter() + .zip(i2) + .zip(i3) + .zip(i4) + .map(|(((&i1, &i2), &i3), &i4)| Self::compute(i1, i2, i3, i4)) + .collect() + } + + #[cfg(feature = "parallel")] + fn compute_ndarray_parallel< + 'a, + D: Dimension + Copy, + A1: Into>, + A2: Into>, + A3: Into>, + A4: Into>, + >( + i1: A1, + i2: A2, + i3: A3, + i4: A4, + ) -> Result, InputError> + where + I1: 'a, + I2: 'a, + I3: 'a, + I4: 'a, + { + let i1: ArrayView = i1.into(); + let i2: ArrayView = i2.into(); + let i3: ArrayView = i3.into(); + let i4: ArrayView = i4.into(); + + Zip::from(i1).and(i2).and(i3).and(i4).par_fold( + || Ok(()), + |_, &i1, &i2, &i3, &i4| Self::validate_inputs(i1, i2, i3, i4), + |a, b| a.and(b), + )?; + + Ok(Zip::from(i1) + .and(i2) + .and(i3) + .and(i4) + .par_map_collect(|&i1, &i2, &i3, &i4| Self::compute_unchecked(i1, i2, i3, i4))) + } +} diff --git a/src/formulas/vapour_pressure.rs b/src/formulas/vapour_pressure.rs new file mode 100644 index 0000000..af30b4d --- /dev/null +++ b/src/formulas/vapour_pressure.rs @@ -0,0 +1,578 @@ +#![allow(clippy::needless_range_loop)] +#![allow(clippy::cast_possible_truncation)] +#![allow(clippy::cast_possible_wrap)] + +//! Formulae to calculate partial vapour pressure of the unsaturated air. + +use crate::constants::DIMLESS_ONE; +use crate::formulas::{Formula1, Formula2}; +use crate::quantities::{ + AtmosphericPressure, DewPointTemperature, QuantityHelpers, RelativeHumidity, + SaturationVapourPressure, SpecificHumidity, VapourPressure, +}; +use crate::Float; +use crate::Storage::Pressure; +use crate::{constants::EPSILON, errors::InputError}; + +use uom::si::pressure::{hectopascal, kilopascal, pascal}; +use uom::si::thermodynamic_temperature::{degree_celsius, kelvin}; + +type FormulaQuantity = VapourPressure; + +/// Formula for computing vapour pressure from specific humidity and pressure. +/// This function is theoretical not empirical. +/// +/// Provided by [Rogers & Yau (1989)](https://www.elsevier.com/books/a-short-course-in-cloud-physics/yau/978-0-08-057094-5). +/// +/// Valid `specific humidity` range: 0.00001 - 2.0 +/// +/// Valid `atmospheric pressure` range: 100Pa - 150000Pa +pub struct Definition1; + +impl Formula2 for Definition1 { + #[inline] + fn validate_inputs_internal( + specific_humidity: SpecificHumidity, + pressure: AtmosphericPressure, + ) -> Result<(), InputError> { + specific_humidity.check_range_si(0.00001, 2.0)?; + pressure.check_range_si(100.0, 150_000.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + specific_humidity: SpecificHumidity, + pressure: AtmosphericPressure, + ) -> VapourPressure { + let specific_humidity = specific_humidity.0; + let pressure = pressure.0; + + let result = -((pressure * specific_humidity) + / ((specific_humidity * (EPSILON - DIMLESS_ONE)) - EPSILON)); + + VapourPressure(result) + } +} + +/// Formula for computing vapour pressure from saturation vapour pressure and relative humidity. +/// +/// Valid `saturation_vapour_pressure` range: 0Pa - 50000Pa +/// +/// Valid `relative_humidity` range: 0.0 - 2.0 +pub struct Definition2; + +impl Formula2 for Definition2 { + #[inline] + fn validate_inputs_internal( + saturation_vapour_pressure: SaturationVapourPressure, + relative_humidity: RelativeHumidity, + ) -> Result<(), InputError> { + relative_humidity.check_range_si(0.0, 2.0)?; + saturation_vapour_pressure.check_range_si(0.0, 50_000.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + saturation_vapour_pressure: SaturationVapourPressure, + relative_humidity: RelativeHumidity, + ) -> VapourPressure { + let result = saturation_vapour_pressure.0 * relative_humidity.0; + + VapourPressure(result) + } +} + +/// Formula for computing vapour pressure from dewpoint temperature and pressure. +/// Should be used for air over water when accuracy is desired. +/// +/// Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). +/// +/// Valid `dewpoint temperature` range: 232K - 324K +/// +/// Valid `atmospheric pressure` range: 100Pa - 150000Pa +pub struct Buck1; + +impl Formula2 for Buck1 { + #[inline] + fn validate_inputs_internal( + dewpoint: DewPointTemperature, + pressure: AtmosphericPressure, + ) -> Result<(), InputError> { + dewpoint.check_range_si(232.0, 324.0)?; + pressure.check_range_si(100.0, 150_000.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + dewpoint: DewPointTemperature, + pressure: AtmosphericPressure, + ) -> VapourPressure { + let dewpoint = dewpoint.0.get::(); + let pressure = pressure.0.get::(); + + let lower_a = 6.1121; + let lower_b = 18.729; + let lower_c = 257.87; + let lower_d = 227.3; + + let upper_a = 0.000_72; + let upper_b = 0.000_003_2; + let upper_c = 0.000_000_000_59; + + let lower_e = + lower_a * (((lower_b - (dewpoint / lower_d)) * dewpoint) / (dewpoint + lower_c)).exp(); + let lower_f = 1.0 + upper_a + (pressure * (upper_b + (upper_c * dewpoint * dewpoint))); + + let result = Pressure::new::(lower_e * lower_f); + + VapourPressure(result) + } +} + +/// Formula for computing vapour pressure from dewpoint temperature and pressure. +/// Should be used for air over ice when accuracy is desired. +/// +/// Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). +/// +/// Valid `dewpoint` range: 193K - 274K +/// +/// Valid `pressure` range: 100Pa - 150000Pa +pub struct Buck2; + +impl Formula2 for Buck2 { + #[inline] + fn validate_inputs_internal( + dewpoint: DewPointTemperature, + pressure: AtmosphericPressure, + ) -> Result<(), InputError> { + dewpoint.check_range_si(193.0, 274.0)?; + pressure.check_range_si(100.0, 150_000.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + dewpoint: DewPointTemperature, + pressure: AtmosphericPressure, + ) -> VapourPressure { + let dewpoint = dewpoint.0.get::(); + let pressure = pressure.0.get::(); + + let lower_a = 6.1115; + let lower_b = 23.036; + let lower_c = 279.82; + let lower_d = 333.7; + + let upper_a = 0.000_22; + let upper_b = 0.000_003_83; + let upper_c = 0.000_000_000_64; + + let lower_e = + lower_a * (((lower_b - (dewpoint / lower_d)) * dewpoint) / (dewpoint + lower_c)).exp(); + let lower_f = 1.0 + upper_a + (pressure * (upper_b + (upper_c * dewpoint * dewpoint))); + + let result = Pressure::new::(lower_e * lower_f); + + VapourPressure(result) + } +} + +/// Formula for computing vapour pressure from dewpoint temperature and pressure. +/// Should be used for air over water for general use. +/// +/// Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). +/// +/// Valid `dewpoint` range: 253K - 324K +/// +/// Valid `pressure` range: 100Pa - 150000Pa +pub struct Buck3; + +impl Formula2 for Buck3 { + #[inline] + fn validate_inputs_internal( + dewpoint: DewPointTemperature, + pressure: AtmosphericPressure, + ) -> Result<(), InputError> { + dewpoint.check_range_si(253.0, 324.0)?; + pressure.check_range_si(100.0, 150_000.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + dewpoint: DewPointTemperature, + pressure: AtmosphericPressure, + ) -> VapourPressure { + let dewpoint = dewpoint.0.get::(); + let pressure = pressure.0.get::(); + + let lower_a = 6.1121; + let lower_b = 17.502; + let lower_c = 240.97; + + let upper_a = 0.000_7; + let upper_b = 0.000_003_46; + + let lower_e = lower_a * ((lower_b * dewpoint) / (dewpoint + lower_c)).exp(); + let lower_f = 1.0 + upper_a + (pressure * upper_b); + + let result = Pressure::new::(lower_e * lower_f); + + VapourPressure(result) + } +} + +/// Formula for computing vapour pressure from dewpoint temperature. +/// Simplified version of [`buck3`]. Very popular in meteorological sources. +/// +/// Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). +/// +/// Valid `dewpoint` range: 253K - 324K +pub struct Buck3Simplified; + +impl Formula1 for Buck3Simplified { + #[inline] + fn validate_inputs_internal(dewpoint: DewPointTemperature) -> Result<(), InputError> { + dewpoint.check_range_si(253.0, 324.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked(dewpoint: DewPointTemperature) -> VapourPressure { + let dewpoint = dewpoint.0.get::(); + + let lower_a = 6.1121; + let lower_b = 17.502; + let lower_c = 240.97; + + let lower_e = lower_a * ((lower_b * dewpoint) / (dewpoint + lower_c)).exp(); + + let result = Pressure::new::(lower_e); + + VapourPressure(result) + } +} + +/// Formula for computing vapour pressure from dewpoint temperature and pressure. +/// Should be used for air over ice for general use. +/// +/// Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). +/// +/// Valid `dewpoint` range: 223K - 274K +/// +/// Valid `pressure` range: 100Pa - 150000Pa +pub struct Buck4; + +impl Formula2 for Buck4 { + #[inline] + fn validate_inputs_internal( + dewpoint: DewPointTemperature, + pressure: AtmosphericPressure, + ) -> Result<(), InputError> { + dewpoint.check_range_si(223.0, 274.0)?; + pressure.check_range_si(100.0, 150_000.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + dewpoint: DewPointTemperature, + pressure: AtmosphericPressure, + ) -> VapourPressure { + let dewpoint = dewpoint.0.get::(); + let pressure = pressure.0.get::(); + + let lower_a = 6.1115; + let lower_b = 22.452; + let lower_c = 272.55; + + let upper_a = 0.000_3; + let upper_b = 0.000_004_18; + + let lower_e = lower_a * ((lower_b * dewpoint) / (dewpoint + lower_c)).exp(); + let lower_f = 1.0 + upper_a + (pressure * upper_b); + + let result = Pressure::new::(lower_e * lower_f); + + VapourPressure(result) + } +} + +/// Formula for computing vapour pressure from dewpoint temperature. +/// Simplified version of [`buck4`], analogical to [`buck3_simplified`]. +/// +/// Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). +/// +/// Valid `dewpoint` range: 223K - 274K +pub struct Buck4Simplified; + +impl Formula1 for Buck4Simplified { + #[inline] + fn validate_inputs_internal(dewpoint: DewPointTemperature) -> Result<(), InputError> { + dewpoint.check_range_si(223.0, 274.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked(dewpoint: DewPointTemperature) -> VapourPressure { + let dewpoint = dewpoint.0.get::(); + + let lower_a = 6.1115; + let lower_b = 22.452; + let lower_c = 272.55; + + let lower_e = lower_a * ((lower_b * dewpoint) / (dewpoint + lower_c)).exp(); + + let result = Pressure::new::(lower_e); + + VapourPressure(result) + } +} + +/// Formula for computing vapour pressure over water from dewpoint temperature. +/// Should be used for temperatures above 273K. +/// +/// Derived by O. Tetens (1930). +/// +/// Valid `dewpoint` range: 273K - 353K +pub struct Tetens1; + +impl Formula1 for Tetens1 { + #[inline] + fn validate_inputs_internal(dewpoint: DewPointTemperature) -> Result<(), InputError> { + dewpoint.check_range_si(273.0, 353.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked(dewpoint: DewPointTemperature) -> VapourPressure { + let dewpoint = dewpoint.0.get::(); + + let lower_a = 0.61078; + let lower_b = 17.27; + let lower_c = 237.3; + + let result = lower_a * ((lower_b * dewpoint) / (dewpoint + lower_c)).exp(); + + let result = Pressure::new::(result); + + VapourPressure(result) + } +} + +/// Formula for computing vapour pressure over water from dewpoint temperature. +/// Should be used when accuracy is required as it is +/// computationally expensive. +/// +/// Derived by A. Wexler (1976) [(doi: 10.6028/jres.080A.071)](https://dx.doi.org/10.6028%2Fjres.080A.071). +/// +/// Valid `dewpoint` range: 273K - 374K +pub struct Wexler1; + +impl Formula1 for Wexler1 { + #[inline] + fn validate_inputs_internal(dewpoint: DewPointTemperature) -> Result<(), InputError> { + dewpoint.check_range_si(273.0, 374.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked(dewpoint: DewPointTemperature) -> VapourPressure { + let dewpoint = dewpoint.get_si_value(); + + // constants from the paper + let g: [Float; 8] = [ + -2991.2729, + -6017.0128, + 18.876_438_54, + -0.028_354_721, + 0.000_017_838_3, + -0.000_000_000_841_504_17, + 0.000_000_000_000_444_125_43, + 2.858_487, + ]; + + let mut ln_p = g[7] * dewpoint.ln(); + + for i in 0..=6 { + ln_p += g[i] * dewpoint.powi(i as i32 - 2); + } + + let result = Pressure::new::(ln_p.exp()); + + VapourPressure(result) + } +} + +/// Formula for computing vapour over ice pressure from dewpoint temperature. +/// Should be used when accuracy is required as it is +/// computationally expensive. +/// +/// Derived by A. Wexler (1977) [(doi: 10.6028/jres.081A.003)](https://dx.doi.org/10.6028%2Fjres.081A.003). +/// +/// Valid `dewpoint` range: 173K - 274K +pub struct Wexler2; + +impl Formula1 for Wexler2 { + #[inline] + fn validate_inputs_internal(dewpoint: DewPointTemperature) -> Result<(), InputError> { + dewpoint.check_range_si(173.0, 274.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked(dewpoint: DewPointTemperature) -> VapourPressure { + let dewpoint = dewpoint.0.get::(); + + // constants from the paper + let big_k: [Float; 6] = [ + -5865.3696, + 22.241_033, + 0.013_749_042, + -0.000_034_031_77, + 0.000_000_026_967_687, + 0.691_865_1, + ]; + + let mut ln_p = big_k[5] * dewpoint.ln(); + + for j in 0..=4 { + ln_p += big_k[j] * dewpoint.powi(j as i32 - 1); + } + + let result = Pressure::new::(ln_p.exp()); + + VapourPressure(result) + } +} + +#[cfg(test)] +mod tests { + use crate::{ + quantities::{ + AtmosphericPressure, RelativeHumidity, SaturationVapourPressure, SpecificHumidity, + }, + tests::{test_with_1arg, test_with_2args, testing_traits::ReferenceAtmosphere, Argument}, + }; + + use super::*; + + #[test] + fn definition1() { + test_with_2args::( + Argument::new([0.00001, 2.0]), + Argument::new([100.0, 150_000.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } + + #[test] + fn definition2() { + test_with_2args::( + Argument::new([0.0, 50_000.0]), + Argument::new([0.0, 2.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } + + #[test] + fn buck1() { + test_with_2args::( + Argument::new([232.0, 324.0]), + Argument::new([100.0, 150_000.0]), + ReferenceAtmosphere::Normal, + 1e1, + ); + } + + #[test] + fn buck2() { + test_with_2args::( + Argument::new([193.0, 274.0]), + Argument::new([100.0, 150_000.0]), + ReferenceAtmosphere::Freezing, + 1e0, + ); + } + + #[test] + fn buck3() { + test_with_2args::( + Argument::new([253.0, 324.0]), + Argument::new([100.0, 150_000.0]), + ReferenceAtmosphere::Normal, + 1e1, + ); + } + + #[test] + fn buck4() { + test_with_2args::( + Argument::new([223.0, 274.0]), + Argument::new([100.0, 150_000.0]), + ReferenceAtmosphere::Freezing, + 1e0, + ); + } + + #[test] + fn buck3_simplified() { + test_with_1arg::( + Argument::new([253.0, 324.0]), + ReferenceAtmosphere::Normal, + 1e0, + ); + } + + #[test] + fn buck4_simplified() { + test_with_1arg::( + Argument::new([223.0, 274.0]), + ReferenceAtmosphere::Freezing, + 1e-1, + ); + } + + #[test] + fn tetens1() { + test_with_1arg::( + Argument::new([273.0, 353.0]), + ReferenceAtmosphere::Normal, + 1e0, + ); + } + + #[test] + fn wexler1() { + test_with_1arg::( + Argument::new([273.0, 374.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } + + #[test] + fn wexler2() { + test_with_1arg::( + Argument::new([173.0, 274.0]), + ReferenceAtmosphere::Freezing, + 1e-12, + ); + } +} diff --git a/src/formulas/vapour_pressure_deficit.rs b/src/formulas/vapour_pressure_deficit.rs new file mode 100644 index 0000000..a9bf4f2 --- /dev/null +++ b/src/formulas/vapour_pressure_deficit.rs @@ -0,0 +1,58 @@ +//! Functions to calculate vapour pressure deficit +//! +//! Vapour-pressure deficit, is the difference (deficit) between +//! the amount of moisture in the air and how much moisture the air can hold +//! when it is saturated ([Wikipedia](https://en.wikipedia.org/wiki/Vapour-pressure_deficit)). + +use crate::errors::InputError; +use crate::formulas::Formula2; +use crate::quantities::{ + QuantityHelpers, SaturationVapourPressure, VapourPressure, VapourPressureDeficit, +}; + +type FormulaQuantity = VapourPressureDeficit; + +/// Formula for computing vapour pressure deficit from vapour pressure and saturation vapour pressure +/// +/// Valid `vapour_pressure` range: 0Pa - 50000Pa +/// +/// Valid `saturation_vapour_pressure` range: 0Pa - 50000Pa +pub struct Definition1; + +impl Formula2 for Definition1 { + #[inline] + fn validate_inputs_internal( + vapour_pressure: VapourPressure, + saturation_vapour_pressure: SaturationVapourPressure, + ) -> Result<(), InputError> { + vapour_pressure.check_range_si(0.0, 50_000.0)?; + saturation_vapour_pressure.check_range_si(0.0, 50_000.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + vapour_pressure: VapourPressure, + saturation_vapour_pressure: SaturationVapourPressure, + ) -> VapourPressureDeficit { + VapourPressureDeficit(saturation_vapour_pressure.0 - vapour_pressure.0) + } +} + +#[cfg(test)] +mod tests { + use crate::tests::{test_with_2args, testing_traits::ReferenceAtmosphere, Argument}; + + use super::*; + + #[test] + fn definition1() { + test_with_2args::( + Argument::new([0.0, 50_000.0]), + Argument::new([0.0, 50_000.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } +} diff --git a/src/formulas/virtual_temperature.rs b/src/formulas/virtual_temperature.rs new file mode 100644 index 0000000..d492e04 --- /dev/null +++ b/src/formulas/virtual_temperature.rs @@ -0,0 +1,166 @@ +//! Functions to calculate virtual temperature of air +//! +//! In atmospheric thermodynamics, the virtual temperature of a moist air parcel is the temperature +//! at which a theoretical dry air parcel would have a total pressure and density equal +//! to the moist parcel of air ([Wikipedia](https://en.wikipedia.org/wiki/Virtual_temperature)). + +use crate::constants::{DIMLESS_ONE, EPSILON, ZERO_KELVIN}; +use crate::errors::InputError; +use crate::formulas::{Formula2, Formula3}; +use crate::quantities::{ + AtmosphericPressure, DryBulbTemperature, MixingRatio, QuantityHelpers, SpecificHumidity, + VapourPressure, VirtualTemperature, +}; + +type FormulaQuantity = VirtualTemperature; + +/// Formula for computing virtual temperature from temperature and mixing ratio. +/// +/// Valid `temperature` range: 173K - 373K +/// +/// Valid `mixing_ratio` range: 0.0000000001 - 0.5 +pub struct Definition1; + +impl Formula2 for Definition1 { + #[inline] + fn validate_inputs_internal( + temperature: DryBulbTemperature, + mixing_ratio: MixingRatio, + ) -> Result<(), InputError> { + temperature.check_range_si(173.0, 354.0)?; + mixing_ratio.check_range_si(0.000_000_000_1, 0.5)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + temperature: DryBulbTemperature, + mixing_ratio: MixingRatio, + ) -> VirtualTemperature { + let result = temperature.0 + * ((mixing_ratio.0 + EPSILON) / (EPSILON * (DIMLESS_ONE + mixing_ratio.0))); + + // this is necessary because result is TemperatureInterval + let result = ZERO_KELVIN + result; + + VirtualTemperature(result) + } +} + +/// Formula for computing virtual temperature from air temperature, pressure and vapour pressure. +/// +/// Valid `temperature` range: 173K - 373K +/// +/// Valid `pressure` range: 100Pa - 150000Pa +/// +/// Valid `vapour_pressure` range: 0Pa - 10000Pa +pub struct Definition2; + +impl Formula3 + for Definition2 +{ + #[inline] + fn validate_inputs_internal( + temperature: DryBulbTemperature, + pressure: AtmosphericPressure, + vapour_pressure: VapourPressure, + ) -> Result<(), InputError> { + temperature.check_range_si(173.0, 354.0)?; + pressure.check_range_si(100.0, 150_000.0)?; + vapour_pressure.check_range_si(0.0, 10_000.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + temperature: DryBulbTemperature, + pressure: AtmosphericPressure, + vapour_pressure: VapourPressure, + ) -> VirtualTemperature { + let result = temperature.0 + / (DIMLESS_ONE - ((vapour_pressure.0 / pressure.0) * (DIMLESS_ONE - EPSILON))); + let result = ZERO_KELVIN + result; + + VirtualTemperature(result) + } +} + +///Formula for computing virtual temperature from air temperature and specific humidity. +/// +///Valid `temperature` range: 173K - 373K +/// +///Valid `specific_humidity` range: 100Pa - 150000Pa +pub struct Definition3; + +impl Formula2 for Definition3 { + #[inline] + fn validate_inputs_internal( + temperature: DryBulbTemperature, + specific_humidity: SpecificHumidity, + ) -> Result<(), InputError> { + temperature.check_range_si(173.0, 354.0)?; + specific_humidity.check_range_si(0.000_000_001, 2.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + temperature: DryBulbTemperature, + specific_humidity: SpecificHumidity, + ) -> VirtualTemperature { + let result = temperature.0 + * (DIMLESS_ONE + (specific_humidity.0 * ((DIMLESS_ONE / EPSILON) - DIMLESS_ONE))); + let result = ZERO_KELVIN + result; + + VirtualTemperature(result) + } +} + +#[cfg(test)] +mod tests { + use crate::tests::{ + test_with_2args, test_with_3args, testing_traits::ReferenceAtmosphere, Argument, + }; + + use super::*; + + #[test] + fn definition1() { + test_with_2args::( + Argument::new([173.0, 354.0]), + Argument::new([0.000_000_000_1, 0.5]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } + + #[test] + fn definition2() { + test_with_3args::< + FormulaQuantity, + DryBulbTemperature, + AtmosphericPressure, + VapourPressure, + Definition2, + >( + Argument::new([173.0, 354.0]), + Argument::new([100.0, 150_000.0]), + Argument::new([0.0, 10_000.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } + + #[test] + fn definition3() { + test_with_2args::( + Argument::new([173.0, 354.0]), + Argument::new([0.000000001, 2.0]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } +} diff --git a/src/formulas/wet_bulb_potential_temperature.rs b/src/formulas/wet_bulb_potential_temperature.rs new file mode 100644 index 0000000..cad1ed4 --- /dev/null +++ b/src/formulas/wet_bulb_potential_temperature.rs @@ -0,0 +1,119 @@ +//!Functions to calculate wet bulb potential temperature of unsaturated air in K. + +use uom::si::ratio::ratio; +use uom::si::thermodynamic_temperature::{degree_celsius, kelvin}; + +use crate::constants::LAMBDA; +use crate::errors::InputError; +use crate::formulas::Formula1; +use crate::quantities::{ + EquivalentPotentialTemperature, QuantityHelpers, WetBulbPotentialTemperature, +}; +use crate::Storage; + +type FormulaQuantity = WetBulbPotentialTemperature; + +/// Formula for computing wet bulb potential temperature from equivalent potential temperature. +/// +/// Derived by R. Davies-Jones (2008) [(doi:10.1175/2007MWR2224.1)](https://doi.org/10.1175/2007MWR2224.1). +/// This formula is a part of three formuale covering a wide range of theta_e inputs. +/// The other two formulas have not been (yet) implemented. +/// +/// Valid `equivalent_potential_temperature` range: 257K - 377K +pub struct DaviesJones1; + +impl Formula1 for DaviesJones1 { + #[inline] + fn validate_inputs_internal( + equivalent_potential_temperature: EquivalentPotentialTemperature, + ) -> Result<(), InputError> { + equivalent_potential_temperature.check_range_si(257.0, 377.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + equivalent_potential_temperature: EquivalentPotentialTemperature, + ) -> WetBulbPotentialTemperature { + let lambda = LAMBDA.get::(); + let equivalent_potential_temperature = equivalent_potential_temperature.0.get::(); + let result = 45.114 - 51.489 * (273.15 / equivalent_potential_temperature).powf(lambda); + + let result = Storage::ThermodynamicTemperature::new::(result); + + WetBulbPotentialTemperature(result) + } +} + +/// Formula for computing wet bulb potential temperature from equivalent potential temperature. +/// +/// Derived by R. Davies-Jones (2008) [(doi:10.1175/2007MWR2224.1)](https://doi.org/10.1175/2007MWR2224.1). +/// This is a very accurate rational-function approximation of [`DaviesJones1`] and two other (unimplemented) +/// formulas. +/// +/// Valid `equivalent_potential_temperature` range: 174K - 377K +pub struct DaviesJones2; + +impl Formula1 for DaviesJones2 { + #[inline] + fn validate_inputs_internal( + equivalent_potential_temperature: EquivalentPotentialTemperature, + ) -> Result<(), InputError> { + equivalent_potential_temperature.check_range_si(174.0, 377.0)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + equivalent_potential_temperature: EquivalentPotentialTemperature, + ) -> WetBulbPotentialTemperature { + let theta_e = equivalent_potential_temperature.get::(); + + let x = theta_e / 273.15; + let a0 = 7.101_574; + let a1 = -20.682_08; + let a2 = 16.111_82; + let a3 = 2.574_631; + let a4 = -5.205_688; + let b1 = -3.552_497; + let b2 = 3.781_782; + let b3 = -0.689_965_5; + let b4 = -0.592_934_0; + + let exponent = (a0 + a1 * x + a2 * x.powi(2) + a3 * x.powi(3) + a4 * x.powi(4)) + / (1. + b1 * x + b2 * x.powi(2) + b3 * x.powi(3) + b4 * x.powi(4)); + let theta_w = theta_e - 273.15 - exponent.exp(); + + WetBulbPotentialTemperature::new::(theta_w) + } +} + +#[cfg(test)] +mod tests { + use crate::{ + quantities::EquivalentPotentialTemperature, + tests::{test_with_1arg, testing_traits::ReferenceAtmosphere, Argument}, + }; + + use super::*; + + #[test] + fn davies_jones1() { + test_with_1arg::( + Argument::new([257.0, 377.0]), + ReferenceAtmosphere::Normal, + 1e-2, + ); + } + + #[test] + fn davies_jones2() { + test_with_1arg::( + Argument::new([174.0, 377.0]), + ReferenceAtmosphere::Normal, + 1e-1, + ); + } +} diff --git a/src/formulas/wet_bulb_temperature.rs b/src/formulas/wet_bulb_temperature.rs new file mode 100644 index 0000000..ed3f2b6 --- /dev/null +++ b/src/formulas/wet_bulb_temperature.rs @@ -0,0 +1,237 @@ +//! Functions to calculate wet bulb temperature of unsaturated air + +use uom::si::pressure::pascal; +use uom::si::ratio::{percent, ratio}; +use uom::si::thermodynamic_temperature::{degree_celsius, kelvin}; + +use crate::constants::{EPSILON, KAPPA, LAMBDA, REF_PRES}; +use crate::errors::InputError; +use crate::formulas::Formula2; +use crate::quantities::{ + AdiabaticEquivalentTemperature, AtmosphericPressure, DryBulbTemperature, QuantityHelpers, + RelativeHumidity, WetBulbTemperature, +}; + +type FormulaQuantity = WetBulbTemperature; + +/// Formula for computing wet bulb temperature from dry bulb temperature and relative humidity. +/// +/// Derived by R. Stull (2011) [(doi:10.1175/JAMC-D-11-0143.1)](https://doi.org/10.1175/JAMC-D-11-0143.1) +/// Created with use of gene-expression programming. +/// +/// Result error is within -1K to +0.65K, with mean absolute error of 0.28K +/// +/// Valid `temperature` range: 253K - 324K\ +/// +/// Valid `relative_humidity` range: 0.05 - 0.99 +pub struct Stull1; + +impl Formula2 for Stull1 { + #[inline] + fn validate_inputs_internal( + temperature: DryBulbTemperature, + relative_humidity: RelativeHumidity, + ) -> Result<(), InputError> { + temperature.check_range_si(253.0, 324.0)?; + relative_humidity.check_range_si(0.05, 0.99)?; + + Ok(()) + } + + #[inline] + fn compute_unchecked( + temperature: DryBulbTemperature, + relative_humidity: RelativeHumidity, + ) -> WetBulbTemperature { + let temperature = temperature.0.get::(); + let relative_humidity = relative_humidity.0.get::(); + + let result = (temperature * (0.151_977 * (relative_humidity + 8.313_659).sqrt()).atan()) + + (temperature + relative_humidity).atan() + - (relative_humidity - 1.676_331).atan() + + (0.003_918_38 * relative_humidity.powf(1.5) * (0.023_101 * relative_humidity).atan()) + - 4.686_035; + + WetBulbTemperature::new::(result) + } +} + +/// Formula for computing wet bulb temperature pressure from adiabatic equivalent temperature, +/// saturation mixing ratio and pressure. +/// +/// Derived by R. Davies-Jones (2008) [(doi:10.1175/2007MWR2224.1)](https://doi.org/10.1175/2007MWR2224.1). +/// +/// Note that this formula conditionally selects equation used for computation. +/// Thus it has discontinuities along the input range and has some performance overhead. +/// +/// Valid `equivalent_temperature` range: 174K - 377K +/// +/// Valid `pressure` range: 100Pa - 150000Pa +pub struct DaviesJones1; + +impl Formula2 + for DaviesJones1 +{ + #[inline] + fn validate_inputs_internal( + equivalent_temperature: AdiabaticEquivalentTemperature, + pressure: AtmosphericPressure, + ) -> Result<(), InputError> { + equivalent_temperature.check_range_si(174., 377.)?; + pressure.check_range_si(100., 150_000.)?; + Ok(()) + } + + #[inline] + fn compute_unchecked( + equivalent_temperature: AdiabaticEquivalentTemperature, + pressure: AtmosphericPressure, + ) -> WetBulbTemperature { + let lambda = LAMBDA.get::(); + let kappa = KAPPA.get::(); + let t_e = equivalent_temperature.get::(); + let p = pressure.get::(); + let p0 = REF_PRES.get::(); + let pi = (p / p0).powf(kappa); + + let indicator = (273.15 / t_e).powf(lambda); + let d_pi = (0.1859 * (pi / p0) + 0.6512).powi(-1); + + let k1_pi = || -38.5 * pi.powi(2) + 137.81 * pi - 53.737; + let k2_pi = || -4.392 * pi.powi(2) + 56.831 * pi - 0.384; + + let eq_4_9 = || k1_pi() - k2_pi() * indicator; + let eq_4_10 = || (k1_pi() - 1.21) - (k2_pi() - 1.21) * indicator; + let eq_4_11 = || { + (k1_pi() - 2.66) - (k2_pi() - 1.21) * indicator + 0.58 * (273.15 / t_e).powf(-lambda) + }; + + let eq_4_8 = || { + let eps = EPSILON.get::(); + let big_a = 2675.0; + let a = 17.67; + let b = 243.5; + let e_s = 6.112 * (a * (t_e - 273.15) / (t_e - 273.15 + b)).exp(); //effectively Tetens + let r_s = (eps * e_s) / ((p / 100.0) - e_s); + let d_ln_es_d = a * b / (t_e - 273.15 + b).powi(2); + + t_e - 273.15 - ((big_a * r_s) / (1.0 + (big_a * r_s * d_ln_es_d))) + }; + + let t_w; + + if indicator > d_pi { + t_w = eq_4_8(); + #[cfg(test)] + log::info!("eq 4.8"); + } else if indicator >= 1.0 && indicator <= d_pi { + t_w = eq_4_9(); + #[cfg(test)] + log::info!("eq 4.9"); + } else if (0.4..1.0).contains(&indicator) { + t_w = eq_4_10(); + #[cfg(test)] + log::info!("eq 4.10"); + } else { + t_w = eq_4_11(); + #[cfg(test)] + log::info!("eq 4.11"); + } + + WetBulbTemperature::new::(t_w) + } +} + +#[cfg(test)] +mod tests { + use float_cmp::assert_approx_eq; + use uom::si::pressure::hectopascal; + + use crate::tests::{test_with_2args, testing_traits::ReferenceAtmosphere, Argument}; + + use super::*; + + #[test] + fn stull1() { + test_with_2args::( + Argument::new([253.0, 324.0]), + Argument::new([0.05, 0.99]), + ReferenceAtmosphere::Normal, + 1e1, + ); + } + + #[test] + fn daviesjones1_norm() { + test_with_2args::< + FormulaQuantity, + AdiabaticEquivalentTemperature, + AtmosphericPressure, + DaviesJones1, + >( + Argument::new([174.0, 377.0]), + Argument::new([100., 150_000.]), + ReferenceAtmosphere::Normal, + 1e-12, + ); + } + + #[test] + fn daviesjones1_freez() { + test_with_2args::< + FormulaQuantity, + AdiabaticEquivalentTemperature, + AtmosphericPressure, + DaviesJones1, + >( + Argument::new([174.0, 377.0]), + Argument::new([100., 150_000.]), + ReferenceAtmosphere::Freezing, + 1e-12, + ); + } + + #[test] + fn daviesjones1_manual() { + testing_logger::setup(); + // eq 4.8 + let pressure = AtmosphericPressure::new::(500.0); + let t_e = AdiabaticEquivalentTemperature::new::(240.); + let wbt = DaviesJones1::compute_unchecked(t_e, pressure); + testing_logger::validate(|captured_logs| { + assert_eq!(captured_logs.len(), 1); + assert_eq!(captured_logs[0].body, "eq 4.8"); + }); + assert_approx_eq!(f64, 238.88007446399956, wbt.get_si_value()); + + // eq 4.9 + let pressure = AtmosphericPressure::new::(850.0); + let t_e = AdiabaticEquivalentTemperature::new::(260.); + let wbt = DaviesJones1::compute_unchecked(t_e, pressure); + testing_logger::validate(|captured_logs| { + assert_eq!(captured_logs.len(), 1); + assert_eq!(captured_logs[0].body, "eq 4.9"); + }); + assert_approx_eq!(f64, 256.61913438141625, wbt.get_si_value()); + + // eq 4.10 + let pressure = AtmosphericPressure::new::(1000.0); + let t_e = AdiabaticEquivalentTemperature::new::(330.); + let wbt = DaviesJones1::compute_unchecked(t_e, pressure); + testing_logger::validate(|captured_logs| { + assert_eq!(captured_logs.len(), 1); + assert_eq!(captured_logs[0].body, "eq 4.10"); + }); + assert_approx_eq!(f64, 291.2797584152462, wbt.get_si_value()); + + // eq 4.11 + let pressure = AtmosphericPressure::new::(1000.0); + let t_e = AdiabaticEquivalentTemperature::new::(370.); + let wbt = DaviesJones1::compute_unchecked(t_e, pressure); + testing_logger::validate(|captured_logs| { + assert_eq!(captured_logs.len(), 1); + assert_eq!(captured_logs[0].body, "eq 4.11"); + }); + assert_approx_eq!(f64, 300.1638142668842, wbt.get_si_value()); + } +} diff --git a/src/lib.rs b/src/lib.rs index 76ebfae..6ccfc52 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,3 +1,10 @@ +#![warn(clippy::pedantic)] +#![warn(missing_docs)] +#![warn(clippy::cargo)] +#![allow(clippy::excessive_precision)] +#![allow(clippy::must_use_candidate)] +#![allow(missing_docs)] // TODO: Remove it + //! Crate providing formulae for air thermodynamic calculations. //! //! # How to use @@ -48,6 +55,9 @@ //! Exact limits are specified in the documentation of each function. //! If the input is out of range the function will return an [`InputError::OutOfRange`](errors::InputError::OutOfRange) with erronous input specified. //! +//! Each function also has `_unchecked` and `_validate` versions. The `_validate` version only checks the inputs with bounds defined for its "parent" function. +//! The `_unchecked` version performs only the calculation without any input checking. All "parent" functions simply call `_validate` and then `_unchecked`. +//! //! # Units //! //! This crate uses basic SI units in the interface. @@ -72,21 +82,23 @@ //! information about the error. This feature potentially is not zero-cost so it is optional. pub mod constants; -pub mod equivalent_potential_temperature; -pub mod errors; -pub mod mixing_ratio; -pub mod potential_temperature; -pub mod relative_humidity; -pub mod specific_humidity; -mod tests_framework; -pub mod vapour_pressure; -pub mod vapour_pressure_deficit; -pub mod virtual_temperature; -pub mod wet_bulb_potential_temperature; -pub mod wet_bulb_temperature; +mod errors; +pub mod formulas; +pub mod quantities; + +pub use errors::InputError; + +#[cfg(test)] +mod tests; #[cfg(not(feature = "double_precision"))] type Float = f32; #[cfg(feature = "double_precision")] type Float = f64; + +#[cfg(not(feature = "double_precision"))] +pub(crate) use uom::si::f32 as Storage; + +#[cfg(feature = "double_precision")] +pub(crate) use uom::si::f64 as Storage; diff --git a/src/mixing_ratio.rs b/src/mixing_ratio.rs deleted file mode 100644 index 1d8f045..0000000 --- a/src/mixing_ratio.rs +++ /dev/null @@ -1,152 +0,0 @@ -//!Functions to calculate mixing ratio of air in kg*kg^-1. -//! -//!To calculate saturation mixing ratio input dry-bulb temperature in place of dewpoint -//!or saturation vapour pressure in place of vapour pressure. - -use crate::{constants::EPSILON, errors::InputError, vapour_pressure}; -use float_cmp::approx_eq; -use crate::Float; - -#[cfg(feature="debug")] -use floccus_proc::logerr; - -///Formula for computing mixing ratio of unsaturated air from air pressure and vapour pressure -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `pressure` range: 100Pa - 150000Pa\ -///Valid `vapour_pressure` range: 0Pa - 10000Pa -/// -///Returns [`InputError::IncorrectArgumentSet`] when inputs are equal, in which -///case division by 0 occurs. -#[cfg_attr(feature = "debug", logerr)] -pub fn general1(pressure: Float, vapour_pressure: Float) -> Result { - //validate inputs - if !(100.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - if !(0.0..=50_000.0).contains(&vapour_pressure) { - return Err(InputError::OutOfRange(String::from("vapour_pressure"))); - } - - if approx_eq!(Float, pressure, vapour_pressure, ulps = 2) { - return Err(InputError::IncorrectArgumentSet(String::from( - "pressure and vapour_pressure cannot be equal", - ))); - } - - let result = EPSILON * (vapour_pressure / (pressure - vapour_pressure)); - Ok(result) -} - -///Formula for computing mixing ratio of unsaturated air from dewpoint temperature and pressure. -///Optimised by performance. -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `dewpoint` range: 273K - 353K\ -///Valid `pressure` range: 100Pa - 150000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn performance1(dewpoint: Float, pressure: Float) -> Result { - //validate inputs - if !(273.0..=353.0).contains(&dewpoint) { - return Err(InputError::OutOfRange(String::from("dewpoint"))); - } - - if !(100.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - let vapour_pressure = vapour_pressure::tetens1(dewpoint)?; - let result = general1(pressure, vapour_pressure)?; - Ok(result) -} - -///Formula for computing mixing ratio of unsaturated air from dewpoint temperature and pressure. -///Optimised by accuracy. -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `dewpoint` range: 232K - 324K\ -///Valid `pressure` range: 100Pa - 150000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn accuracy1(dewpoint: Float, pressure: Float) -> Result { - //validate inputs - if !(232.0..=324.0).contains(&dewpoint) { - return Err(InputError::OutOfRange(String::from("dewpoint"))); - } - - if !(100.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - let vapour_pressure = vapour_pressure::buck1(dewpoint, pressure)?; - let result = general1(pressure, vapour_pressure)?; - Ok(result) -} - -#[cfg(test)] -mod tests { - use crate::{ - mixing_ratio, - tests_framework::{self, Argument}, - }; - - #[test] - fn general1() { - assert!(tests_framework::test_with_2args( - &mixing_ratio::general1, - Argument { - name: "pressure", - def_val: 101325.0, - range: [100.0, 150_000.0] - }, - Argument { - name: "vapour_pressure", - def_val: 3500.0, - range: [0.0, 50_000.0] - }, - 0.022253316630823517 - )); - } - - #[test] - fn performance1() { - assert!(tests_framework::test_with_2args( - &mixing_ratio::performance1, - Argument { - name: "dewpoint", - def_val: 300.0, - range: [273.0, 353.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [100.0, 150_000.0] - }, - 0.022477100514593465 - )); - } - - #[test] - fn accuracy1() { - assert!(tests_framework::test_with_2args( - &mixing_ratio::accuracy1, - Argument { - name: "dewpoint", - def_val: 300.0, - range: [232.0, 324.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [100.0, 150_000.0] - }, - 0.022587116896465847 - )); - } -} diff --git a/src/potential_temperature.rs b/src/potential_temperature.rs deleted file mode 100644 index 66610f8..0000000 --- a/src/potential_temperature.rs +++ /dev/null @@ -1,95 +0,0 @@ -//!Functions to calculate potential temperature of dry air in K. - -use float_cmp::approx_eq; -use crate::Float; -use crate::{ - constants::{C_P, R_D}, - errors::InputError, -}; - -#[cfg(feature="debug")] -use floccus_proc::logerr; - -///Formula for computing potential temperature of dry air from temperature, pressure and vapour pressure. -/// -///Provided by R. Davies-Jones (2009) [(doi:10.1175/2009MWR2774.1)](https://doi.org/10.1175/2009MWR2774.1) -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `temperature` range: 253K - 324K\ -///Valid `pressure` range: 100Pa - 150000Pa\ -///Valid `vapour_pressure` range: 0Pa - 10000Pa -/// -///Returns [`InputError::IncorrectArgumentSet`] when `pressure` and `vapour_pressure` are equal, -///in which case division by 0 occurs. -/// -///Returns [`InputError::IncorrectArgumentSet`] when `pressure` is lower than `vapour_pressure`, -///in which case floating-point exponentation of negative number occurs. -#[cfg_attr(feature = "debug", logerr)] -pub fn davies_jones1( - temperature: Float, - pressure: Float, - vapour_pressure: Float, -) -> Result { - if !(253.0..=324.0).contains(&temperature) { - return Err(InputError::OutOfRange(String::from("temperature"))); - } - - if !(100.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - if !(0.0..=10_000.0).contains(&vapour_pressure) { - return Err(InputError::OutOfRange(String::from("vapour_pressure"))); - } - - if approx_eq!(Float, pressure, vapour_pressure, ulps = 2) { - return Err(InputError::IncorrectArgumentSet(String::from( - "pressure and vapour_pressure cannot be equal", - ))); - } - - if vapour_pressure > pressure { - return Err(InputError::IncorrectArgumentSet(String::from( - "vapour_pressure cannot be higher than pressure", - ))); - } - - let kappa = R_D / C_P; - - let result = temperature * (100_000.0 / (pressure - vapour_pressure)).powf(kappa); - - Ok(result) -} - -#[cfg(test)] -mod tests { - use crate::{ - tests_framework::{self, Argument}, - potential_temperature, - }; - - #[test] - fn davies_jones1() { - assert!(tests_framework::test_with_3args( - &potential_temperature::davies_jones1, - Argument { - name: "temperature", - def_val: 300.0, - range: [253.0, 324.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [100.0, 150_000.0] - }, - Argument { - name: "vapour_pressure", - def_val: 3000.0, - range: [0.0, 10_000.0] - }, - 301.45136519081666 - )); - } -} diff --git a/src/quantities.rs b/src/quantities.rs new file mode 100644 index 0000000..6334efd --- /dev/null +++ b/src/quantities.rs @@ -0,0 +1,101 @@ +#![allow(missing_docs)] + +use crate::{Float, errors::InputError}; +use std::any::type_name; +use std::fmt::Debug; + +pub trait ThermodynamicQuantity: + Debug + Clone + Copy + PartialEq + PartialOrd + Default + Send + Sync +{ +} + +pub(crate) trait QuantityHelpers: ThermodynamicQuantity { + fn get_si_value(&self) -> Float; + + fn name() -> &'static str { + type_name::() + } + + #[inline] + fn check_range_si(&self, lower_bound: Float, upper_bound: Float) -> Result<(), InputError> { + if !(lower_bound..=upper_bound).contains(&self.get_si_value()) { + return Err(InputError::OutOfRange(Self::name())); + } + + Ok(()) + } +} + +macro_rules! define_quantity { + ($quantity:ident, $storage:ident, $uom_module:ident, $si_unit:ident) => { + #[derive(Debug, Clone, Copy, PartialEq, PartialOrd, Default)] + pub struct $quantity(pub crate::Storage::$storage); + + impl $quantity { + pub fn get(&self) -> Float + where + T: uom::si::$uom_module::Unit + uom::si::$uom_module::Conversion, + { + self.0.get::() + } + + pub fn new(value: Float) -> Self + where + T: uom::si::$uom_module::Unit + uom::si::$uom_module::Conversion, + { + Self(crate::Storage::$storage::new::(value)) + } + } + + impl ThermodynamicQuantity for $quantity {} + + impl QuantityHelpers for $quantity { + fn get_si_value(&self) -> Float { + self.get::() + } + } + }; +} + +macro_rules! define_temperature { + ($quantity:ident) => { + define_quantity!( + $quantity, + ThermodynamicTemperature, + thermodynamic_temperature, + kelvin + ); + }; +} + +macro_rules! define_pressure { + ($quantity:ident) => { + define_quantity!($quantity, Pressure, pressure, pascal); + }; +} + +macro_rules! define_ratio { + ($quantity:ident) => { + define_quantity!($quantity, Ratio, ratio, ratio); + }; +} + +define_temperature!(DryBulbTemperature); +define_temperature!(WetBulbTemperature); +define_temperature!(DewPointTemperature); +define_temperature!(VirtualTemperature); +define_temperature!(IsobaricEquivalentTemperature); +define_temperature!(AdiabaticEquivalentTemperature); +define_temperature!(PotentialTemperature); +define_temperature!(EquivalentPotentialTemperature); +define_temperature!(WetBulbPotentialTemperature); + +define_pressure!(AtmosphericPressure); +define_pressure!(VapourPressure); +define_pressure!(SaturationVapourPressure); +define_pressure!(VapourPressureDeficit); + +define_ratio!(MixingRatio); +define_ratio!(SaturationMixingRatio); +define_ratio!(SpecificHumidity); +define_ratio!(RelativeHumidity); diff --git a/src/relative_humidity.rs b/src/relative_humidity.rs deleted file mode 100644 index c65895c..0000000 --- a/src/relative_humidity.rs +++ /dev/null @@ -1,250 +0,0 @@ -//!Functions to calculate relative humidity in %/100 - -use crate::{errors::InputError, mixing_ratio, vapour_pressure}; -use crate::Float; - -#[cfg(feature="debug")] -use floccus_proc::logerr; - -///Formula for computing relative humidity from mixing ratio and saturation mixing ratio. -///Can be used interchangeably with [`general2`]. -/// -///By the definition of mixing ratio, this formula is mathematically equivalent of -///formula used in [`general2`]. -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `mixing_ratio` range: 0.00001 - 0.5\ -///Valid `saturation_mixing_ratio` range: 0.00001 - 0.5 -#[cfg_attr(feature = "debug", logerr)] -pub fn general1(mixing_ratio: Float, saturation_mixing_ratio: Float) -> Result { - if !(0.00001..=10.0).contains(&mixing_ratio) { - return Err(InputError::OutOfRange(String::from("mixing_ratio"))); - } - - if !(0.00001..=10.0).contains(&saturation_mixing_ratio) { - return Err(InputError::OutOfRange(String::from( - "saturation_mixing_ratio", - ))); - } - - Ok(mixing_ratio / saturation_mixing_ratio) -} - -///Formula for computing relative humidity from vapour pressure and saturation vapour pressure. -///Can be used interchangeably with [`general1`]. -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `vapour_pressure` range: 0Pa - 10000Pa -///Valid `saturation_vapour_pressure` range: 0Pa - 10000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn general2(vapour_pressure: Float, saturation_vapour_pressure: Float) -> Result { - if !(0.0..=50_000.0).contains(&vapour_pressure) { - return Err(InputError::OutOfRange(String::from("vapour_pressure"))); - } - - if !(0.1..=50_000.0).contains(&saturation_vapour_pressure) { - return Err(InputError::OutOfRange(String::from( - "saturation_vapour_pressure", - ))); - } - - Ok(vapour_pressure / saturation_vapour_pressure) -} - -///Formula for computing relative humidity from temperature and dewpoint using [`tetens1`](vapour_pressure::tetens1) -///function for vapour pressure calculation -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `temperature` range: 273K - 353K -///Valid `dewpoint` range: 273K - 353K -#[cfg_attr(feature = "debug", logerr)] -pub fn general3(temperature: Float, dewpoint: Float) -> Result { - if !(273.0..=353.0).contains(&temperature) { - return Err(InputError::OutOfRange(String::from("temperature"))); - } - - if !(273.0..=353.0).contains(&temperature) { - return Err(InputError::OutOfRange(String::from("temperature"))); - } - - let vapour_pressure = vapour_pressure::tetens1(dewpoint)?; - let saturation_vapour_pressure = vapour_pressure::tetens1(temperature)?; - let result = general2(vapour_pressure, saturation_vapour_pressure)?; - - Ok(result) -} - -///Formula for computing relative humidity from temperature, dewpoint and pressure using [`buck3`](vapour_pressure::buck3) -///function for vapour pressure calculation -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `temperature` range: 253K - 324K\ -///Valid `dewpoint` range: 253K - 324K\ -///Valid `pressure` range: 100Pa - 150000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn general4(temperature: Float, dewpoint: Float, pressure: Float) -> Result { - if !(253.0..=324.0).contains(&temperature) { - return Err(InputError::OutOfRange(String::from("temperature"))); - } - - if !(253.0..=324.0).contains(&temperature) { - return Err(InputError::OutOfRange(String::from("temperature"))); - } - - if !(100.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - let vapour_pressure = vapour_pressure::buck3(dewpoint, pressure)?; - let saturation_vapour_pressure = vapour_pressure::buck3(temperature, pressure)?; - let result = general2(vapour_pressure, saturation_vapour_pressure)?; - - Ok(result) -} - -///Formula for computing relative humidity from temperature, dewpoint and pressure using [`accuracy1`](mixing_ratio::accuracy1) -///function for mixing ratio calculation -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `temperature` range: 232K - 324K\ -///Valid `dewpoint` range: 232K - 324K\ -///Valid `pressure` range: 100Pa - 150000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn general5(temperature: Float, dewpoint: Float, pressure: Float) -> Result { - if !(232.0..=314.0).contains(&temperature) { - return Err(InputError::OutOfRange(String::from("temperature"))); - } - - if !(232.0..=314.0).contains(&dewpoint) { - return Err(InputError::OutOfRange(String::from("dewpoint"))); - } - - if !(10000.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - let mixing_ratio = mixing_ratio::accuracy1(dewpoint, pressure)?; - let saturation_mixing_ratio = mixing_ratio::accuracy1(temperature, pressure)?; - //println!("{} {}", mixing_ratio, saturation_mixing_ratio); - let result = general1(mixing_ratio, saturation_mixing_ratio)?; - - Ok(result) -} - -#[cfg(test)] -mod tests { - use crate::{ - relative_humidity, - tests_framework::{self, Argument}, - }; - - #[test] - fn general1() { - assert!(tests_framework::test_with_2args( - &relative_humidity::general1, - Argument { - name: "mixing_ratio", - def_val: 0.01064, - range: [0.00001, 10.0] - }, - Argument { - name: "saturation_mixing_ratio", - def_val: 0.01467, - range: [0.00001, 10.0] - }, - 0.7252897068847989 - )); - } - - #[test] - fn general2() { - assert!(tests_framework::test_with_2args( - &relative_humidity::general2, - Argument { - name: "vapour_pressure", - def_val: 1706.0, - range: [0.0, 50_000.0] - }, - Argument { - name: "saturation_vapour_pressure", - def_val: 2339.0, - range: [0.1, 50_000.0] - }, - 0.7293715262932877 - )); - } - - #[test] - fn general3() { - assert!(tests_framework::test_with_2args( - &relative_humidity::general3, - Argument { - name: "temperature", - def_val: 300.0, - range: [273.0, 353.0] - }, - Argument { - name: "dewpoint", - def_val: 290.0, - range: [273.0, 353.0] - }, - 0.5431069897660531 - )); - } - - #[test] - fn general4() { - assert!(tests_framework::test_with_3args( - &relative_humidity::general4, - Argument { - name: "temperature", - def_val: 300.0, - range: [253.0, 324.0] - }, - Argument { - name: "dewpoint", - def_val: 290.0, - range: [253.0, 324.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [100.0, 150_000.0] - }, - 0.5429224562155812 - )); - } - - #[test] - fn general5() { - assert!(tests_framework::test_with_3args( - &relative_humidity::general5, - Argument { - name: "temperature", - def_val: 300.0, - range: [232.0, 314.0] - }, - Argument { - name: "dewpoint", - def_val: 290.0, - range: [232.0, 314.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [10000.0, 150_000.0] - }, - 0.5338747953552858 - )); - } -} diff --git a/src/specific_humidity.rs b/src/specific_humidity.rs deleted file mode 100644 index 4f4af45..0000000 --- a/src/specific_humidity.rs +++ /dev/null @@ -1,63 +0,0 @@ -//!Functions to calculate specific humidity of air in kg*kg^-1. -//! -//!Specific humidity (or moisture content) is the ratio of the mass -//!of water vapor to the total mass of the air parcel [Wikipedia](https://en.wikipedia.org/wiki/Humidity#Specific_humidity). -//! -//!Specific humidity is approximately equal to mixing ratio. - -use crate::{constants::EPSILON, errors::InputError}; -use crate::Float; - -#[cfg(feature="debug")] -use floccus_proc::logerr; - -///Formula for computing specific humidity from vapour pressure and pressure. -///Reverse function of [`vapour_pressure::general1`](crate::vapour_pressure::general1). -///This function is theoretical not empirical. -/// -///Provided by [Rogers & Yau (1989)](https://www.elsevier.com/books/a-short-course-in-cloud-physics/yau/978-0-08-057094-5). -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `vapour_pressure` range: 0Pa - 50000OPa\, -///Valid `pressure` range: 100Pa - 150000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn general1(vapour_pressure: Float, pressure: Float) -> Result { - if !(0.0..=50_000.0).contains(&vapour_pressure) { - return Err(InputError::OutOfRange(String::from("vapour_pressure"))); - } - - if !(100.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - let result = EPSILON * (vapour_pressure / (pressure - (vapour_pressure * (1.0 - EPSILON)))); - Ok(result) -} - -#[cfg(test)] -mod tests { - use crate::{ - specific_humidity, - tests_framework::{self, Argument}, - }; - - #[test] - fn general1() { - assert!(tests_framework::test_with_2args( - &specific_humidity::general1, - Argument { - name: "vapour_pressure", - def_val: 3000.0, - range: [0.0, 50_000.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [100.0, 150_000.0] - }, - 0.018623845512674677 - )); - } -} diff --git a/src/tests.rs b/src/tests.rs new file mode 100644 index 0000000..78835ff --- /dev/null +++ b/src/tests.rs @@ -0,0 +1,55 @@ +mod four_arg; +mod one_arg; +pub(crate) mod reference_values; +pub(crate) mod testing_traits; +mod three_arg; +mod two_arg; + +use crate::{Float, InputError}; +use float_cmp::assert_approx_eq; +use std::marker::PhantomData; + +pub use self::four_arg::test_with_4args; +pub use self::one_arg::test_with_1arg; +pub use self::three_arg::test_with_3args; +pub use self::two_arg::test_with_2args; + +use self::testing_traits::{ReferenceAtmosphere, TestingQuantity}; + +#[derive(Copy, Clone, Debug)] +pub struct Argument { + pub range: [Float; 2], + _quantity: PhantomData, +} + +impl Argument { + pub fn new(range: [Float; 2]) -> Self { + Self { + range, + _quantity: PhantomData, + } + } + + pub fn quantity_name(&self) -> &str { + I::name() + } + + pub fn ref_val(&self, atm: ReferenceAtmosphere) -> I { + I::ref_val_si(atm) + } +} + +fn check_result(result: T, atm: ReferenceAtmosphere, eps: Float) { + let expected = T::ref_val_si(atm).get_si_value(); + let result = result.get_si_value(); + + assert_approx_eq!(Float, result, expected, epsilon = eps); +} + +pub fn check_range_error(result: InputError, expected_name: &str) { + if let InputError::OutOfRange(name) = result { + assert_eq!(name, expected_name); + } else { + panic!("wrong error type") + } +} diff --git a/src/tests/four_arg.rs b/src/tests/four_arg.rs new file mode 100644 index 0000000..4563d94 --- /dev/null +++ b/src/tests/four_arg.rs @@ -0,0 +1,289 @@ +use float_cmp::assert_approx_eq; +use itertools::multiunzip; +#[cfg(feature = "array")] +use ndarray::Array1; + +use super::Argument; +use super::check_result; +use super::testing_traits::{ReferenceAtmosphere, TestingQuantity}; +use crate::Float; +use crate::errors::InputError; +use crate::formulas::Formula4; +use crate::tests::check_range_error; +use std::mem::discriminant; + +pub fn test_with_4args< + O: TestingQuantity, + I1: TestingQuantity, + I2: TestingQuantity, + I3: TestingQuantity, + I4: TestingQuantity, + F: Formula4, +>( + arg1: Argument, + arg2: Argument, + arg3: Argument, + arg4: Argument, + atm: ReferenceAtmosphere, + eps: Float, +) { + let ref_result = F::compute( + arg1.ref_val(atm), + arg2.ref_val(atm), + arg3.ref_val(atm), + arg4.ref_val(atm), + ) + .unwrap(); + check_result(ref_result, atm, eps); + + let results = vec![ + F::compute( + arg1.ref_val(atm), + arg2.ref_val(atm), + arg3.ref_val(atm), + arg4.ref_val(atm), + ), + F::compute( + I1::new_si(-9999.0), + arg2.ref_val(atm), + arg3.ref_val(atm), + arg4.ref_val(atm), + ), + F::compute( + arg1.ref_val(atm), + I2::new_si(-9999.0), + arg3.ref_val(atm), + arg4.ref_val(atm), + ), + F::compute( + arg1.ref_val(atm), + arg2.ref_val(atm), + I3::new_si(-9999.0), + arg4.ref_val(atm), + ), + F::compute( + arg1.ref_val(atm), + arg2.ref_val(atm), + arg3.ref_val(atm), + I4::new_si(-9999.0), + ), + F::compute( + arg1.ref_val(atm), + I2::new_si(-9999.0), + I3::new_si(-9999.0), + arg4.ref_val(atm), + ), + F::compute( + I1::new_si(-9999.0), + I2::new_si(-9999.0), + arg3.ref_val(atm), + arg4.ref_val(atm), + ), + F::compute( + I1::new_si(-9999.0), + arg2.ref_val(atm), + I3::new_si(-9999.0), + arg4.ref_val(atm), + ), + F::compute( + I1::new_si(-9999.0), + arg2.ref_val(atm), + arg3.ref_val(atm), + I4::new_si(-9999.0), + ), + F::compute( + I1::new_si(-9999.0), + I2::new_si(-9999.0), + I3::new_si(-9999.0), + I4::new_si(-9999.0), + ), + ]; + + for result in results { + if let Ok(result) = result { + assert!(result.get_si_value().is_finite()); + } + } + + //the third promise of the crate is to always return finite f64 + //if all inputs are within the range + //the only allowed error is InccorectArgumentsSet as it can occur + //for values within valid range + for arg1_itr in 0..=20 { + for arg2_itr in 0..=20 { + for arg3_itr in 0..=20 { + for arg4_itr in 0..=20 { + let arg1_tmp = (((arg1.range[1] - arg1.range[0]) / 20.0) * arg1_itr as Float) + + arg1.range[0]; + let arg2_tmp = (((arg2.range[1] - arg2.range[0]) / 20.0) * arg2_itr as Float) + + arg2.range[0]; + let arg3_tmp = (((arg3.range[1] - arg3.range[0]) / 20.0) * arg3_itr as Float) + + arg3.range[0]; + let arg4_tmp = (((arg4.range[1] - arg4.range[0]) / 20.0) * arg4_itr as Float) + + arg4.range[0]; + + let arg1_tmp = I1::new_si(arg1_tmp); + let arg2_tmp = I2::new_si(arg2_tmp); + let arg3_tmp = I3::new_si(arg3_tmp); + let arg4_tmp = I4::new_si(arg4_tmp); + + let result = F::compute(arg1_tmp, arg2_tmp, arg3_tmp, arg4_tmp); + + match result.clone() { + Ok(r) => assert!(r.get_si_value().is_finite()), + + Err(e) => assert_eq!( + discriminant(&InputError::IncorrectArgumentSet("")), + discriminant(&e) + ), + } + } + } + } + } + + //the fourth promise of the crate is to return an error with + //erronous variable name when input is out of range + let result = F::compute( + I1::new_si(arg1.range[0] - 0.1), + arg2.ref_val(atm), + arg3.ref_val(atm), + arg4.ref_val(atm), + ) + .unwrap_err(); + check_range_error(result, arg1.quantity_name()); + let result = F::compute( + I1::new_si(arg1.range[1] + 0.1), + arg2.ref_val(atm), + arg3.ref_val(atm), + arg4.ref_val(atm), + ) + .unwrap_err(); + check_range_error(result, arg1.quantity_name()); + + let result = F::compute( + arg1.ref_val(atm), + I2::new_si(arg2.range[0] - 0.1), + arg3.ref_val(atm), + arg4.ref_val(atm), + ) + .unwrap_err(); + check_range_error(result, arg2.quantity_name()); + let result = F::compute( + arg1.ref_val(atm), + I2::new_si(arg2.range[1] + 0.1), + arg3.ref_val(atm), + arg4.ref_val(atm), + ) + .unwrap_err(); + check_range_error(result, arg2.quantity_name()); + + let result = F::compute( + arg1.ref_val(atm), + arg2.ref_val(atm), + I3::new_si(arg3.range[0] - 0.1), + arg4.ref_val(atm), + ) + .unwrap_err(); + check_range_error(result, arg3.quantity_name()); + let result = F::compute( + arg1.ref_val(atm), + arg2.ref_val(atm), + I3::new_si(arg3.range[1] + 0.1), + arg4.ref_val(atm), + ) + .unwrap_err(); + check_range_error(result, arg3.quantity_name()); + + let result = F::compute( + arg1.ref_val(atm), + arg2.ref_val(atm), + arg3.ref_val(atm), + I4::new_si(arg4.range[0] - 0.1), + ) + .unwrap_err(); + check_range_error(result, arg4.quantity_name()); + let result = F::compute( + arg1.ref_val(atm), + arg2.ref_val(atm), + arg3.ref_val(atm), + I4::new_si(arg4.range[1] + 0.1), + ) + .unwrap_err(); + check_range_error(result, arg4.quantity_name()); + + let arg_vecs = (-10..=10).map(|i| i as Float / 1000.0).map(|i| { + ( + I1::new_si((1.0 + i) * arg1.ref_val(atm).get_si_value()), + I2::new_si((1.0 + i) * arg2.ref_val(atm).get_si_value()), + I3::new_si((1.0 + i) * arg3.ref_val(atm).get_si_value()), + I4::new_si((1.0 + i) * arg4.ref_val(atm).get_si_value()), + ) + }); + + let arg_vecs: (Vec<_>, Vec<_>, Vec<_>, Vec<_>) = multiunzip(arg_vecs); + + #[cfg(feature = "array")] + let arg_arrs = ( + Array1::from(arg_vecs.0.clone()), + Array1::from(arg_vecs.1.clone()), + Array1::from(arg_vecs.2.clone()), + Array1::from(arg_vecs.3.clone()), + ); + + let result_vec = F::compute_vec(&arg_vecs.0, &arg_vecs.1, &arg_vecs.2, &arg_vecs.3).unwrap(); + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_vec[10].get_si_value(), + ulps = 4 + ); + + #[cfg(feature = "array")] + let result_arr = + F::compute_ndarray(&arg_arrs.0, &arg_arrs.1, &arg_arrs.2, &arg_arrs.3).unwrap(); + #[cfg(feature = "array")] + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_arr[10].get_si_value(), + ulps = 4 + ); + + #[cfg(feature = "parallel")] + let result_vec = + F::compute_vec_parallel(&arg_vecs.0, &arg_vecs.1, &arg_vecs.2, &arg_vecs.3).unwrap(); + #[cfg(feature = "parallel")] + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_vec[10].get_si_value(), + ulps = 4 + ); + + #[cfg(feature = "parallel")] + let result_arr = + F::compute_ndarray_parallel(&arg_arrs.0, &arg_arrs.1, &arg_arrs.2, &arg_arrs.3).unwrap(); + #[cfg(feature = "parallel")] + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_arr[10].get_si_value(), + ulps = 4 + ); + + let result_imperial = F::compute( + arg1.ref_val(atm).imperial(), + arg2.ref_val(atm).imperial(), + arg3.ref_val(atm).imperial(), + arg4.ref_val(atm).imperial(), + ) + .unwrap(); + + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_imperial.get_si_value(), + epsilon = 1e-12 + ); +} diff --git a/src/tests/one_arg.rs b/src/tests/one_arg.rs new file mode 100644 index 0000000..ef3e7f4 --- /dev/null +++ b/src/tests/one_arg.rs @@ -0,0 +1,123 @@ +use float_cmp::assert_approx_eq; +#[cfg(feature = "array")] +use ndarray::Array1; + +use super::Argument; +use super::check_result; +use super::testing_traits::{ReferenceAtmosphere, TestingQuantity}; +use crate::Float; +use crate::errors::InputError; +use crate::formulas::Formula1; +use crate::tests::check_range_error; +use std::mem::discriminant; + +pub fn test_with_1arg>( + arg1: Argument, + atm: ReferenceAtmosphere, + eps: Float, +) { + //the first promise of the crate is that returned value + //is calculated correctly + let ref_result = F::compute(arg1.ref_val(atm)).unwrap(); + check_result(ref_result, atm, eps); + + // the second promise of the crate is to never return NaN or Inf + // here we check several edge cases for that + // the function can only return a finite number or InputError + // check for error implicit as Result<> ensures that if value + // is not Ok() then it is Err() + // here we don't care about error being correct + let results = vec![ + F::compute(arg1.ref_val(atm)), + F::compute(I1::new_si(-9999.0)), + ]; + + for result in results { + if let Ok(result) = result { + assert!(result.get_si_value().is_finite()); + } + } + + //the third promise of the crate is to always return finite f64 + //if all inputs are within the range + //the only allowed error is InccorectArgumentsSet as it can occur + //for values within valid range + for arg1_itr in 0..=100 { + let arg1_tmp = + (((arg1.range[1] - arg1.range[0]) / 100.0) * arg1_itr as Float) + arg1.range[0]; + + let arg1_tmp = I1::new_si(arg1_tmp); + + let result = F::compute(arg1_tmp); + + match result { + Ok(r) => assert!(r.get_si_value().is_finite()), + Err(e) => assert_eq!( + discriminant(&InputError::IncorrectArgumentSet("")), + discriminant(&e) + ), + } + } + + //the fourth promise of the crate is to return an error with + //erronous variable name when input is out of range + let result = F::compute(I1::new_si(arg1.range[0] - 0.1)).unwrap_err(); + check_range_error(result, arg1.quantity_name()); + let result = F::compute(I1::new_si(arg1.range[1] + 0.1)).unwrap_err(); + check_range_error(result, arg1.quantity_name()); + + let arg_vecs: Vec<_> = (-10..=10) + .map(|i| i as Float / 1000.0) + .map(|i| I1::new_si((1.0 + i) * arg1.ref_val(atm).get_si_value())) + .collect(); + + #[cfg(feature = "array")] + let arg_arrs = Array1::from(arg_vecs.clone()); + + let result_vec = F::compute_vec(&arg_vecs).unwrap(); + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_vec[10].get_si_value(), + ulps = 4 + ); + + #[cfg(feature = "array")] + let result_arr = F::compute_ndarray(&arg_arrs).unwrap(); + #[cfg(feature = "array")] + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_arr[10].get_si_value(), + ulps = 4 + ); + + #[cfg(feature = "parallel")] + let result_vec = F::compute_vec_parallel(&arg_vecs).unwrap(); + #[cfg(feature = "parallel")] + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_vec[10].get_si_value(), + ulps = 4 + ); + + #[cfg(feature = "parallel")] + let result_arr = F::compute_ndarray_parallel(&arg_arrs).unwrap(); + #[cfg(feature = "parallel")] + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_arr[10].get_si_value(), + ulps = 4 + ); + + let result_imperial = F::compute(arg1.ref_val(atm).imperial()).unwrap(); + + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_imperial.get_si_value(), + epsilon = 1e-12 + ); +} diff --git a/src/tests/reference_values.rs b/src/tests/reference_values.rs new file mode 100644 index 0000000..b8284dc --- /dev/null +++ b/src/tests/reference_values.rs @@ -0,0 +1,39 @@ +use std::f64::NAN; + +use crate::Float; + +pub(crate) const TEMP_NORM: Float = 300.0; +pub(crate) const DWPT_NORM: Float = 290.0; +pub(crate) const PRES_NORM: Float = 100000.0; +pub(crate) const VP_NORM: Float = 1919.4253257541593; +pub(crate) const SVP_NORM: Float = 3535.4235919263083; +pub(crate) const RH_NORM: Float = 0.5429124052171476; +pub(crate) const MR_NORM: Float = 0.012172079452423202; +pub(crate) const SMR_NROM: Float = 0.022419969290542845; +pub(crate) const VPD_NORM: Float = 1615.998266172149; +pub(crate) const SH_NORM: Float = 0.012025701656390478; +pub(crate) const THETAE_NORM: Float = 331.329289539998; +pub(crate) const THETA_NORM: Float = 301.66581400702955; +pub(crate) const THETAW_NORM: Float = 292.0717306393948; +pub(crate) const WBT_NORM: Float = 291.79619238702986; +pub(crate) const VRT_NORM: Float = 302.1926517941886; +pub(crate) const TIE_NORM: Float = 330.29726646682764; +pub(crate) const TAE_NORM: Float = 331.8799684989062; + +pub(crate) const TEMP_FREEZ: Float = 260.0; +pub(crate) const DWPT_FREEZ: Float = 255.0; +pub(crate) const PRES_FREEZ: Float = 100000.0; +pub(crate) const VP_FREEZ: Float = 123.17937690212507; +pub(crate) const SVP_FREEZ: Float = 195.84980045970696; +pub(crate) const RH_FREEZ: Float = 0.6289481868911442; +pub(crate) const MR_FREEZ: Float = 0.0007670962389744638; +pub(crate) const SMR_FREEZ: Float = 0.0012196493367222787; +pub(crate) const VPD_FREEZ: Float = 72.67042355758188; +pub(crate) const SH_FREEZ: Float = 0.000766508253376156; +pub(crate) const THETAE_FREEZ: Float = 261.95287841149707; +pub(crate) const THETA_FREEZ: Float = 260.0915766593588; +pub(crate) const THETAW_FREEZ: Float = 258.6611332391296; +pub(crate) const WBT_FREEZ: Float = 258.4257262913491; +pub(crate) const VRT_FREEZ: Float = 260.12112343315795; +pub(crate) const TIE_FREEZ: Float = NAN; +pub(crate) const TAE_FREEZ: Float = 261.9163911760276; diff --git a/src/tests/testing_traits.rs b/src/tests/testing_traits.rs new file mode 100644 index 0000000..225748d --- /dev/null +++ b/src/tests/testing_traits.rs @@ -0,0 +1,343 @@ +use uom::si::{ + pressure::{pascal, pound_force_per_square_inch}, + ratio::{per_mille, percent, ratio}, + thermodynamic_temperature::{degree_fahrenheit, kelvin}, +}; + +use super::reference_values::*; +use crate::{quantities::*, Float}; + +#[derive(Copy, Clone, Debug)] +pub(crate) enum ReferenceAtmosphere { + Normal, + Freezing, +} + +pub(crate) trait TestingQuantity: QuantityHelpers { + fn new_si(value: Float) -> Self; + fn imperial(&self) -> Self; + fn ref_val_si(atm: ReferenceAtmosphere) -> Self; +} + +impl TestingQuantity for DryBulbTemperature { + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } + + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(TEMP_NORM), + ReferenceAtmosphere::Freezing => Self::new::(TEMP_FREEZ), + } + } +} + +impl TestingQuantity for DewPointTemperature { + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } + + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(DWPT_NORM), + ReferenceAtmosphere::Freezing => Self::new::(DWPT_FREEZ), + } + } +} + +impl TestingQuantity for EquivalentPotentialTemperature { + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } + + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(THETAE_NORM), + ReferenceAtmosphere::Freezing => Self::new::(THETAE_FREEZ), + } + } +} + +impl TestingQuantity for AtmosphericPressure { + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } + + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(PRES_NORM), + ReferenceAtmosphere::Freezing => Self::new::(PRES_FREEZ), + } + } +} + +impl TestingQuantity for VapourPressure { + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } + + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(VP_NORM), + ReferenceAtmosphere::Freezing => Self::new::(VP_FREEZ), + } + } +} + +impl TestingQuantity for SaturationVapourPressure { + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(SVP_NORM), + ReferenceAtmosphere::Freezing => Self::new::(SVP_FREEZ), + } + } + + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } +} + +impl TestingQuantity for MixingRatio { + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(MR_NORM), + ReferenceAtmosphere::Freezing => Self::new::(MR_FREEZ), + } + } + + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } +} + +impl TestingQuantity for SaturationMixingRatio { + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } + + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(SMR_NROM), + ReferenceAtmosphere::Freezing => Self::new::(SMR_FREEZ), + } + } +} + +impl TestingQuantity for RelativeHumidity { + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } + + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(RH_NORM), + ReferenceAtmosphere::Freezing => Self::new::(RH_FREEZ), + } + } +} + +impl TestingQuantity for SpecificHumidity { + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } + + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(SH_NORM), + ReferenceAtmosphere::Freezing => Self::new::(SH_FREEZ), + } + } +} + +impl TestingQuantity for VapourPressureDeficit { + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } + + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(VPD_NORM), + ReferenceAtmosphere::Freezing => Self::new::(VPD_FREEZ), + } + } +} + +impl TestingQuantity for VirtualTemperature { + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } + + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(VRT_NORM), + ReferenceAtmosphere::Freezing => Self::new::(VRT_FREEZ), + } + } +} + +impl TestingQuantity for WetBulbTemperature { + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } + + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(WBT_NORM), + ReferenceAtmosphere::Freezing => Self::new::(WBT_FREEZ), + } + } +} + +impl TestingQuantity for WetBulbPotentialTemperature { + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } + + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(THETAW_NORM), + ReferenceAtmosphere::Freezing => Self::new::(THETAW_FREEZ), + } + } +} + +impl TestingQuantity for PotentialTemperature { + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } + + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(THETA_NORM), + ReferenceAtmosphere::Freezing => Self::new::(THETA_FREEZ), + } + } +} + +impl TestingQuantity for IsobaricEquivalentTemperature { + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } + + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(TIE_NORM), + ReferenceAtmosphere::Freezing => Self::new::(TIE_FREEZ), + } + } +} + +impl TestingQuantity for AdiabaticEquivalentTemperature { + fn new_si(value: Float) -> Self { + Self::new::(value) + } + + fn imperial(&self) -> Self { + let value = self.0.get::(); + + Self::new::(value) + } + + fn ref_val_si(atm: ReferenceAtmosphere) -> Self { + match atm { + ReferenceAtmosphere::Normal => Self::new::(TAE_NORM), + ReferenceAtmosphere::Freezing => Self::new::(TAE_FREEZ), + } + } +} diff --git a/src/tests/three_arg.rs b/src/tests/three_arg.rs new file mode 100644 index 0000000..b91f8eb --- /dev/null +++ b/src/tests/three_arg.rs @@ -0,0 +1,206 @@ +use float_cmp::assert_approx_eq; +use itertools::multiunzip; +#[cfg(feature = "array")] +use ndarray::Array1; + +use super::Argument; +use super::check_result; +use super::testing_traits::{ReferenceAtmosphere, TestingQuantity}; +use crate::Float; +use crate::errors::InputError; +use crate::formulas::Formula3; +use crate::tests::check_range_error; +use std::mem::discriminant; + +pub fn test_with_3args< + O: TestingQuantity, + I1: TestingQuantity, + I2: TestingQuantity, + I3: TestingQuantity, + F: Formula3, +>( + arg1: Argument, + arg2: Argument, + arg3: Argument, + atm: ReferenceAtmosphere, + eps: Float, +) { + //the first promise of the crate is that returned value + //is calculated correctly + let ref_result = F::compute(arg1.ref_val(atm), arg2.ref_val(atm), arg3.ref_val(atm)).unwrap(); + check_result(ref_result, atm, eps); + + // the second promise of the crate is to never return NaN or Inf + // here we check several edge cases for that + // the function can only return a finite number or InputError + // check for error implicit as Result<> ensures that if value + // is not Ok() then it is Err() + // here we don't care about error being correct + let results = vec![ + F::compute(arg1.ref_val(atm), arg2.ref_val(atm), arg3.ref_val(atm)), + F::compute(I1::new_si(-9999.0), arg2.ref_val(atm), arg3.ref_val(atm)), + F::compute(arg1.ref_val(atm), I2::new_si(-9999.0), arg3.ref_val(atm)), + F::compute(arg1.ref_val(atm), arg2.ref_val(atm), I3::new_si(-9999.0)), + F::compute(arg1.ref_val(atm), I2::new_si(-9999.0), arg3.ref_val(atm)), + F::compute(arg1.ref_val(atm), I2::new_si(-9999.0), I3::new_si(-9999.0)), + F::compute(I1::new_si(-9999.0), arg2.ref_val(atm), I3::new_si(-9999.0)), + F::compute( + I1::new_si(-9999.0), + I2::new_si(-9999.0), + I3::new_si(-9999.0), + ), + ]; + + for result in results { + if let Ok(result) = result { + assert!(result.get_si_value().is_finite()); + } + } + + //the third promise of the crate is to always return finite f64 + //if all inputs are within the range + //the only allowed error is InccorectArgumentsSet as it can occur + //for values within valid range + for arg1_itr in 0..=100 { + for arg2_itr in 0..=100 { + for arg3_itr in 0..=100 { + let arg1_tmp = + (((arg1.range[1] - arg1.range[0]) / 100.0) * arg1_itr as Float) + arg1.range[0]; + let arg2_tmp = + (((arg2.range[1] - arg2.range[0]) / 100.0) * arg2_itr as Float) + arg2.range[0]; + let arg3_tmp = + (((arg3.range[1] - arg3.range[0]) / 100.0) * arg3_itr as Float) + arg3.range[0]; + + let arg1_tmp = I1::new_si(arg1_tmp); + let arg2_tmp = I2::new_si(arg2_tmp); + let arg3_tmp = I3::new_si(arg3_tmp); + + let result = F::compute(arg1_tmp, arg2_tmp, arg3_tmp); + + match result { + Ok(r) => assert!(r.get_si_value().is_finite()), + Err(e) => assert_eq!( + discriminant(&InputError::IncorrectArgumentSet("")), + discriminant(&e) + ), + } + } + } + } + + //the fourth promise of the crate is to return an error with + //erronous variable name when input is out of range + let result = F::compute( + I1::new_si(arg1.range[0] - 0.1), + arg2.ref_val(atm), + arg3.ref_val(atm), + ) + .unwrap_err(); + check_range_error(result, arg1.quantity_name()); + let result = F::compute( + I1::new_si(arg1.range[1] + 0.1), + arg2.ref_val(atm), + arg3.ref_val(atm), + ) + .unwrap_err(); + check_range_error(result, arg1.quantity_name()); + + let result = F::compute( + arg1.ref_val(atm), + I2::new_si(arg2.range[0] - 0.1), + arg3.ref_val(atm), + ) + .unwrap_err(); + check_range_error(result, arg2.quantity_name()); + let result = F::compute( + arg1.ref_val(atm), + I2::new_si(arg2.range[1] + 0.1), + arg3.ref_val(atm), + ) + .unwrap_err(); + check_range_error(result, arg2.quantity_name()); + + let result = F::compute( + arg1.ref_val(atm), + arg2.ref_val(atm), + I3::new_si(arg3.range[0] - 0.1), + ) + .unwrap_err(); + check_range_error(result, arg3.quantity_name()); + let result = F::compute( + arg1.ref_val(atm), + arg2.ref_val(atm), + I3::new_si(arg3.range[1] + 0.1), + ) + .unwrap_err(); + check_range_error(result, arg3.quantity_name()); + + let arg_vecs = (-10..=10).map(|i| i as Float / 1000.0).map(|i| { + ( + I1::new_si((1.0 + i) * arg1.ref_val(atm).get_si_value()), + I2::new_si((1.0 + i) * arg2.ref_val(atm).get_si_value()), + I3::new_si((1.0 + i) * arg3.ref_val(atm).get_si_value()), + ) + }); + + let arg_vecs: (Vec<_>, Vec<_>, Vec<_>) = multiunzip(arg_vecs); + + #[cfg(feature = "array")] + let arg_arrs = ( + Array1::from(arg_vecs.0.clone()), + Array1::from(arg_vecs.1.clone()), + Array1::from(arg_vecs.2.clone()), + ); + + let result_vec = F::compute_vec(&arg_vecs.0, &arg_vecs.1, &arg_vecs.2).unwrap(); + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_vec[10].get_si_value(), + ulps = 4 + ); + + #[cfg(feature = "array")] + let result_arr = F::compute_ndarray(&arg_arrs.0, &arg_arrs.1, &arg_arrs.2).unwrap(); + #[cfg(feature = "array")] + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_arr[10].get_si_value(), + ulps = 4 + ); + + #[cfg(feature = "parallel")] + let result_vec = F::compute_vec_parallel(&arg_vecs.0, &arg_vecs.1, &arg_vecs.2).unwrap(); + #[cfg(feature = "parallel")] + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_vec[10].get_si_value(), + ulps = 4 + ); + + #[cfg(feature = "parallel")] + let result_arr = F::compute_ndarray_parallel(&arg_arrs.0, &arg_arrs.1, &arg_arrs.2).unwrap(); + #[cfg(feature = "parallel")] + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_arr[10].get_si_value(), + ulps = 4 + ); + + let result_imperial = F::compute( + arg1.ref_val(atm).imperial(), + arg2.ref_val(atm).imperial(), + arg3.ref_val(atm).imperial(), + ) + .unwrap(); + + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_imperial.get_si_value(), + epsilon = 1e-12 + ); +} diff --git a/src/tests/two_arg.rs b/src/tests/two_arg.rs new file mode 100644 index 0000000..2bad9d2 --- /dev/null +++ b/src/tests/two_arg.rs @@ -0,0 +1,148 @@ +use float_cmp::assert_approx_eq; +#[cfg(feature = "array")] +use ndarray::Array1; + +use super::Argument; +use super::check_result; +use super::testing_traits::{ReferenceAtmosphere, TestingQuantity}; +use crate::Float; +use crate::errors::InputError; +use crate::formulas::Formula2; +use crate::tests::check_range_error; +use std::mem::discriminant; + +pub fn test_with_2args< + O: TestingQuantity, + I1: TestingQuantity, + I2: TestingQuantity, + F: Formula2, +>( + arg1: Argument, + arg2: Argument, + atm: ReferenceAtmosphere, + eps: Float, +) { + //the first promise of the crate is that returned value + //is calculated correctly + let ref_result = F::compute(arg1.ref_val(atm), arg2.ref_val(atm)).unwrap(); + check_result(ref_result, atm, eps); + + // the second promise of the crate is to never return NaN or Inf + // here we check several edge cases for that + // the function can only return a finite number or InputError + // check for error implicit as Result<> ensures that if value + // is not Ok() then it is Err() + // here we don't care about error being correct + let results = vec![ + F::compute(arg1.ref_val(atm), arg2.ref_val(atm)), + F::compute(I1::new_si(-9999.0), arg2.ref_val(atm)), + F::compute(arg1.ref_val(atm), I2::new_si(-9999.0)), + F::compute(I1::new_si(-9999.0), I2::new_si(-9999.0)), + ]; + + for result in results { + if let Ok(result) = result { + assert!(result.get_si_value().is_finite()); + } + } + + for arg1_itr in 0..=100 { + for arg2_itr in 0..=100 { + let arg1_tmp = + (((arg1.range[1] - arg1.range[0]) / 100.0) * arg1_itr as Float) + arg1.range[0]; + let arg2_tmp = + (((arg2.range[1] - arg2.range[0]) / 100.0) * arg2_itr as Float) + arg2.range[0]; + + let arg1_tmp = I1::new_si(arg1_tmp); + let arg2_tmp = I2::new_si(arg2_tmp); + + let result = F::compute(arg1_tmp, arg2_tmp); + + match result { + Ok(r) => assert!(r.get_si_value().is_finite()), + Err(e) => assert_eq!( + discriminant(&InputError::IncorrectArgumentSet("")), + discriminant(&e) + ), + } + } + } + + //the fourth promise of the crate is to return an error with + //erronous variable name when input is out of range + let result = F::compute(I1::new_si(arg1.range[0] - 0.1), arg2.ref_val(atm)).unwrap_err(); + check_range_error(result, arg1.quantity_name()); + let result = F::compute(I1::new_si(arg1.range[1] + 0.1), arg2.ref_val(atm)).unwrap_err(); + check_range_error(result, arg1.quantity_name()); + + let result = F::compute(arg1.ref_val(atm), I2::new_si(arg2.range[0] - 0.1)).unwrap_err(); + check_range_error(result, arg2.quantity_name()); + let result = F::compute(arg1.ref_val(atm), I2::new_si(arg2.range[1] + 0.1)).unwrap_err(); + check_range_error(result, arg2.quantity_name()); + + let arg_vecs: (Vec<_>, Vec<_>) = (-10..=10) + .map(|i| i as Float / 1000.0) + .map(|i| { + ( + I1::new_si((1.0 + i) * arg1.ref_val(atm).get_si_value()), + I2::new_si((1.0 + i) * arg2.ref_val(atm).get_si_value()), + ) + }) + .unzip(); + + dbg!(&arg_vecs); + + #[cfg(feature = "array")] + let arg_arrs = ( + Array1::from(arg_vecs.0.clone()), + Array1::from(arg_vecs.1.clone()), + ); + + let result_vec = F::compute_vec(&arg_vecs.0, &arg_vecs.1).unwrap(); + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_vec[10].get_si_value(), + ulps = 4 + ); + + #[cfg(feature = "array")] + let result_arr = F::compute_ndarray(&arg_arrs.0, &arg_arrs.1).unwrap(); + #[cfg(feature = "array")] + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_arr[10].get_si_value(), + ulps = 4 + ); + + #[cfg(feature = "parallel")] + let result_vec = F::compute_vec_parallel(&arg_vecs.0, &arg_vecs.1).unwrap(); + #[cfg(feature = "parallel")] + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_vec[10].get_si_value(), + ulps = 4 + ); + + #[cfg(feature = "parallel")] + let result_arr = F::compute_ndarray_parallel(&arg_arrs.0, &arg_arrs.1).unwrap(); + #[cfg(feature = "parallel")] + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_arr[10].get_si_value(), + ulps = 4 + ); + + let result_imperial = + F::compute(arg1.ref_val(atm).imperial(), arg2.ref_val(atm).imperial()).unwrap(); + + assert_approx_eq!( + Float, + ref_result.get_si_value(), + result_imperial.get_si_value(), + epsilon = 1e-12 + ); +} diff --git a/src/tests_framework.rs b/src/tests_framework.rs deleted file mode 100644 index ed35883..0000000 --- a/src/tests_framework.rs +++ /dev/null @@ -1,198 +0,0 @@ -use crate::errors::InputError; -use crate::Float; -use float_cmp::assert_approx_eq; -use std::mem::discriminant; - -#[derive(Copy, Clone, Debug, Default)] -pub struct Argument { - pub name: &'static str, - pub def_val: Float, - pub range: [Float; 2], -} - -//due to a bug [https://github.com/rust-lang/rust/issues/46379] -//cargo flags those functions as a dead code even though -//they are used in tests -#[allow(dead_code)] -//this function should work as a reference for other test functions below -pub fn test_with_2args( - tested_function: &dyn Fn(Float, Float) -> Result, - arg1: Argument, - arg2: Argument, - expected_result: Float, -) -> bool { - //the first promise of the crate is that returned value - //is calculated correctly - let result = tested_function(arg1.def_val, arg2.def_val).unwrap(); - assert_approx_eq!(Float, result, expected_result, epsilon = 0.01); - - //the second promise of the crate is to never return NaN or Inf - //here we check several edge cases for that - //the function can only return a finite number or InputError - //check for error implicit as Result<> ensures that if value - //is not Ok() then it is Err() - //here we don't care about error being correct - let results = vec![ - tested_function(0.0, arg2.def_val), - tested_function(arg1.def_val, 0.0), - tested_function(0.0, 0.0), - ]; - - for result in results { - if result.is_ok() { - assert!(result.unwrap().is_finite()); - } - } - - //the third promise of the crate is to always return finite f64 - //if all inputs are within the range - //the only allowed error is InccorectArgumentsSet as it can occur - //for values within valid range - for arg1_itr in 0..=100 { - for arg2_itr in 0..=100 { - let arg1_tmp = - (((arg1.range[1] - arg1.range[0]) / 100.0) * arg1_itr as Float) + arg1.range[0]; - let arg2_tmp = - (((arg2.range[1] - arg2.range[0]) / 100.0) * arg2_itr as Float) + arg2.range[0]; - - let result = tested_function(arg1_tmp, arg2_tmp); - - if result.is_err() { - assert!( - discriminant(&InputError::IncorrectArgumentSet(String::from(""))) - == discriminant(&result.unwrap_err()) - ); - } else { - assert!(result.unwrap().is_finite()); - } - } - } - - //the fourth promise of the crate is to return an error with - //erronous variable name when input is out of range - let expected = InputError::OutOfRange(String::from(arg1.name)); - let result = tested_function(arg1.range[0] - 0.1, arg2.def_val).unwrap_err(); - assert_eq!(result, expected); - let result = tested_function(arg1.range[1] + 0.1, arg2.def_val).unwrap_err(); - assert_eq!(result, expected); - - let expected = InputError::OutOfRange(String::from(arg2.name)); - let result = tested_function(arg1.def_val, arg2.range[0] - 0.1).unwrap_err(); - assert_eq!(result, expected); - let result = tested_function(arg1.def_val, arg2.range[1] + 0.1).unwrap_err(); - assert_eq!(result, expected); - - true -} - -#[allow(dead_code)] -pub fn test_with_1arg( - tested_function: &dyn Fn(Float) -> Result, - arg1: Argument, - expected_result: Float, -) -> bool { - let result = tested_function(arg1.def_val).unwrap(); - assert_approx_eq!(Float, result, expected_result, epsilon = 0.01); - - let results = vec![tested_function(0.0)]; - - for result in results { - if result.is_ok() { - assert!(result.unwrap().is_finite()); - } - } - - for arg1_itr in 0..=100 { - let arg1_tmp = - (((arg1.range[1] - arg1.range[0]) / 100.0) * arg1_itr as Float) + arg1.range[0]; - - let result = tested_function(arg1_tmp); - - if result.is_err() { - assert!( - discriminant(&InputError::IncorrectArgumentSet(String::from(""))) - == discriminant(&result.unwrap_err()) - ); - } else { - assert!(result.unwrap().is_finite()); - } - } - - let expected = InputError::OutOfRange(String::from(arg1.name)); - let result = tested_function(arg1.range[0] - 0.1).unwrap_err(); - assert_eq!(result, expected); - let result = tested_function(arg1.range[1] + 0.1).unwrap_err(); - assert_eq!(result, expected); - - true -} - -#[allow(dead_code)] -pub fn test_with_3args( - tested_function: &dyn Fn(Float, Float, Float) -> Result, - arg1: Argument, - arg2: Argument, - arg3: Argument, - expected_result: Float, -) -> bool { - let result = tested_function(arg1.def_val, arg2.def_val, arg3.def_val).unwrap(); - assert_approx_eq!(Float, result, expected_result, epsilon = 0.01); - - let results = vec![ - tested_function(0.0, arg2.def_val, arg3.def_val), - tested_function(arg1.def_val, 0.0, arg3.def_val), - tested_function(arg1.def_val, arg2.def_val, 0.0), - tested_function(0.0, 0.0, 0.0), - ]; - - for result in results { - if result.is_ok() { - assert!(result.unwrap().is_finite()); - } - } - - for arg1_itr in 0..=100 { - for arg2_itr in 0..=100 { - for arg3_itr in 0..=100 { - let arg1_tmp = - (((arg1.range[1] - arg1.range[0]) / 100.0) * arg1_itr as Float) + arg1.range[0]; - let arg2_tmp = - (((arg2.range[1] - arg2.range[0]) / 100.0) * arg2_itr as Float) + arg2.range[0]; - let arg3_tmp = - (((arg3.range[1] - arg3.range[0]) / 100.0) * arg3_itr as Float) + arg3.range[0]; - - let result = tested_function(arg1_tmp, arg2_tmp, arg3_tmp); - - if result.is_err() { - assert!( - discriminant(&InputError::IncorrectArgumentSet(String::from(""))) - == discriminant(&result.unwrap_err()) - ); - } else { - println!("{} {} {} {:?}", arg1_tmp, arg2_tmp, arg3_tmp, result); - assert!(result.unwrap().is_finite()); - } - } - } - } - - let expected = InputError::OutOfRange(String::from(arg1.name)); - let result = tested_function(arg1.range[0] - 0.1, arg2.def_val, arg3.def_val).unwrap_err(); - assert_eq!(result, expected); - let result = tested_function(arg1.range[1] + 0.1, arg2.def_val, arg3.def_val).unwrap_err(); - assert_eq!(result, expected); - - let expected = InputError::OutOfRange(String::from(arg2.name)); - let result = tested_function(arg1.def_val, arg2.range[0] - 0.1, arg3.def_val).unwrap_err(); - assert_eq!(result, expected); - let result = tested_function(arg1.def_val, arg2.range[1] + 0.1, arg3.def_val).unwrap_err(); - assert_eq!(result, expected); - - let expected = InputError::OutOfRange(String::from(arg3.name)); - let result = tested_function(arg1.def_val, arg2.def_val, arg3.range[0] - 0.1).unwrap_err(); - assert_eq!(result, expected); - let result = tested_function(arg1.def_val, arg2.def_val, arg3.range[1] + 0.1).unwrap_err(); - assert_eq!(result, expected); - - true -} diff --git a/src/vapour_pressure.rs b/src/vapour_pressure.rs deleted file mode 100644 index ab98289..0000000 --- a/src/vapour_pressure.rs +++ /dev/null @@ -1,589 +0,0 @@ -//!Functions to calculate partial vapour pressure of the unsaturated air in Pa. -//! -//!To compute saturation vapour pressure input dry-bulb temperature in place of dewpoint temperature. - -use crate::Float; -use crate::{ - constants::{EPSILON, ZERO_CELSIUS}, - errors::InputError, -}; - -#[cfg(feature="debug")] -use floccus_proc::logerr; - -///Formula for computing vapour pressure from specific humidity and pressure. -///This function is theoretical not empirical. -/// -///Provided by [Rogers & Yau (1989)](https://www.elsevier.com/books/a-short-course-in-cloud-physics/yau/978-0-08-057094-5). -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `specific_humidity` range: 0.00001 - 2.0\ -///Valid `pressure` range: 100Pa - 150000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn general1(specific_humidity: Float, pressure: Float) -> Result { - //validate inputs - if !(0.00001..=2.0).contains(&specific_humidity) { - return Err(InputError::OutOfRange(String::from("specific_humidity"))); - } - - if !(100.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - let result = - -((pressure * specific_humidity) / ((specific_humidity * (EPSILON - 1.0)) - EPSILON)); - - Ok(result) -} - -///Formula for computing vapour pressure from dewpoint temperature and pressure. -///Should be used for air over water when accuracy is desired. -/// -///Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `dewpoint` range: 232K - 324K\ -///Valid `pressure` range: 100Pa - 150000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn buck1(dewpoint: Float, pressure: Float) -> Result { - //validate inputs - if !(232.0..=324.0).contains(&dewpoint) { - return Err(InputError::OutOfRange(String::from("dewpoint"))); - } - - if !(100.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - let dewpoint = dewpoint - ZERO_CELSIUS; //convert to C - let pressure = pressure / 100.0; //convert to hPa - - let lower_a = 6.1121; - let lower_b = 18.729; - let lower_c = 257.87; - let lower_d = 227.3; - - let upper_a = 0.000_72; - let upper_b = 0.000_003_2; - let upper_c = 0.000_000_000_59; - - let lower_e = - lower_a * (((lower_b - (dewpoint / lower_d)) * dewpoint) / (dewpoint + lower_c)).exp(); - let lower_f = 1.0 + upper_a + (pressure * (upper_b + (upper_c * dewpoint * dewpoint))); - - Ok((lower_e * lower_f) * 100.0) //return in Pa -} - -///Formula for computing vapour pressure from dewpoint temperature and pressure. -///Should be used for air over ice when accuracy is desired. -/// -///Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `dewpoint` range: 193K - 274K\ -///Valid `pressure` range: 100Pa - 150000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn buck2(dewpoint: Float, pressure: Float) -> Result { - //validate inputs - if !(193.0..=274.0).contains(&dewpoint) { - return Err(InputError::OutOfRange(String::from("dewpoint"))); - } - - if !(100.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - let dewpoint = dewpoint - ZERO_CELSIUS; //convert to C - let pressure = pressure / 100.0; //convert to hPa - - let lower_a = 6.1115; - let lower_b = 23.036; - let lower_c = 279.82; - let lower_d = 333.7; - - let upper_a = 0.000_22; - let upper_b = 0.000_003_83; - let upper_c = 0.000_000_000_64; - - let lower_e = - lower_a * (((lower_b - (dewpoint / lower_d)) * dewpoint) / (dewpoint + lower_c)).exp(); - let lower_f = 1.0 + upper_a + (pressure * (upper_b + (upper_c * dewpoint * dewpoint))); - - Ok((lower_e * lower_f) * 100.0) //return in Pa -} - -///Formula for computing vapour pressure from dewpoint temperature and pressure. -///Should be used for air over water for general use. -/// -///Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `dewpoint` range: 253K - 324K\ -///Valid `pressure` range: 100Pa - 150000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn buck3(dewpoint: Float, pressure: Float) -> Result { - //validate inputs - if !(253.0..=324.0).contains(&dewpoint) { - return Err(InputError::OutOfRange(String::from("dewpoint"))); - } - - if !(100.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - let dewpoint = dewpoint - ZERO_CELSIUS; //convert to C - let pressure = pressure / 100.0; //convert to hPa - - let lower_a = 6.1121; - let lower_b = 17.502; - let lower_c = 240.97; - - let upper_a = 0.000_7; - let upper_b = 0.000_003_46; - - let lower_e = lower_a * ((lower_b * dewpoint) / (dewpoint + lower_c)).exp(); - let lower_f = 1.0 + upper_a + (pressure * upper_b); - - Ok((lower_e * lower_f) * 100.0) //return in Pa -} - -///Formula for computing vapour pressure from dewpoint temperature. -///Simplified version of [`buck3`]. Very popular in meteorological sources. -/// -///Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `dewpoint` range: 253K - 324K -#[cfg_attr(feature = "debug", logerr)] -pub fn buck3_simplified(dewpoint: Float) -> Result { - //validate inputs - if !(253.0..=324.0).contains(&dewpoint) { - return Err(InputError::OutOfRange(String::from("dewpoint"))); - } - - let dewpoint = dewpoint - ZERO_CELSIUS; //convert to C - - let lower_a = 6.1121; - let lower_b = 17.502; - let lower_c = 240.97; - - let lower_e = lower_a * ((lower_b * dewpoint) / (dewpoint + lower_c)).exp(); - - Ok(lower_e * 100.0) //return in Pa -} - -///Formula for computing vapour pressure from dewpoint temperature and pressure. -///Should be used for air over ice for general use. -/// -///Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `dewpoint` range: 223K - 274K\ -///Valid `pressure` range: 100Pa - 150000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn buck4(dewpoint: Float, pressure: Float) -> Result { - //validate inputs - if !(223.0..=274.0).contains(&dewpoint) { - return Err(InputError::OutOfRange(String::from("dewpoint"))); - } - - if !(100.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - let dewpoint = dewpoint - ZERO_CELSIUS; //convert to C - let pressure = pressure / 100.0; //convert to hPa - - let lower_a = 6.1115; - let lower_b = 22.452; - let lower_c = 272.55; - - let upper_a = 0.000_3; - let upper_b = 0.000_004_18; - - let lower_e = lower_a * ((lower_b * dewpoint) / (dewpoint + lower_c)).exp(); - let lower_f = 1.0 + upper_a + (pressure * upper_b); - - Ok((lower_e * lower_f) * 100.0) //return in Pa -} - -///Formula for computing vapour pressure from dewpoint temperature. -///Simplified version of [`buck4`], analogical to [`buck3_simplified`]. -/// -///Derived by A. L. Buck (1981) [(doi: 10.1175/1520-0450(1981)020<1527:nefcvp>2.0.co;2)](https://doi.org/10.1175/1520-0450(1981)020%3C1527:NEFCVP%3E2.0.CO;2). -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `dewpoint` range: 223K - 274K -#[cfg_attr(feature = "debug", logerr)] -pub fn buck4_simplified(dewpoint: Float) -> Result { - //validate inputs - if !(223.0..=274.0).contains(&dewpoint) { - return Err(InputError::OutOfRange(String::from("dewpoint"))); - } - - let dewpoint = dewpoint - ZERO_CELSIUS; //convert to C - - let lower_a = 6.1115; - let lower_b = 22.452; - let lower_c = 272.55; - - let lower_e = lower_a * ((lower_b * dewpoint) / (dewpoint + lower_c)).exp(); - - Ok(lower_e * 100.0) //return in Pa -} - -///Formula for computing vapour pressure over water from dewpoint temperature. -///Should be used for temperatures above 273K. -/// -///Derived by O. Tetens (1930). -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when input is out of range.\ -///Valid `dewpoint` range: 273K - 353K -#[cfg_attr(feature = "debug", logerr)] -pub fn tetens1(dewpoint: Float) -> Result { - //validate inputs - if !(273.0..=353.0).contains(&dewpoint) { - return Err(InputError::OutOfRange(String::from("dewpoint"))); - } - - let dewpoint = dewpoint - ZERO_CELSIUS; //convert to C - - let lower_a = 0.61078; - let lower_b = 17.27; - let lower_c = 237.3; - - let result = lower_a * ((lower_b * dewpoint) / (dewpoint + lower_c)).exp(); - - Ok(result * 1000.0) //return in Pa -} - -///Formula for computing **ONLY** vapour pressure from saturation vapour pressure and relative humidity. -///For saturation vapour pressure use [`saturation_specific2`] -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when input is out of range.\ -///Valid `saturation_vapour_pressure` range: 0Pa - 10000Pa\ -///Valid `relative_humidity` range: 0.0 - 1.0 -#[cfg_attr(feature = "debug", logerr)] -pub fn saturation_specific1( - saturation_vapour_pressure: Float, - relative_humidity: Float, -) -> Result { - if !(0.0..=2.0).contains(&relative_humidity) { - return Err(InputError::OutOfRange(String::from("relative_humidity"))); - } - - if !(0.0..=50_000.0).contains(&saturation_vapour_pressure) { - return Err(InputError::OutOfRange(String::from( - "saturation_vapour_pressure", - ))); - } - - Ok(saturation_vapour_pressure * relative_humidity) -} - -///Formula for computing **ONLY** saturation vapour pressure from vapour pressure and relative humidity. -///For vapour pressure use [`saturation_specific1`] -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when input is out of range.\ -///Valid `vapour_pressure` range: 0Pa - 10000Pa\ -///Valid `relative_humidity` range: 0.00001 - 1.0 -#[cfg_attr(feature = "debug", logerr)] -pub fn saturation_specific2( - vapour_pressure: Float, - relative_humidity: Float, -) -> Result { - if !(0.00001..=2.0).contains(&relative_humidity) { - return Err(InputError::OutOfRange(String::from("relative_humidity"))); - } - - if !(0.0..=10_000.0).contains(&vapour_pressure) { - return Err(InputError::OutOfRange(String::from("vapour_pressure"))); - } - - Ok(vapour_pressure / relative_humidity) -} - -///Formula for computing vapour pressure over water from dewpoint temperature. -///Should be used when accuracy is required as it is -///computationally expensive. -/// -///Derived by A. Wexler (1976) [(doi: 10.6028/jres.080A.071)](https://dx.doi.org/10.6028%2Fjres.080A.071). -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `dewpoint` range: 273K - 374K -#[cfg_attr(feature = "debug", logerr)] -pub fn wexler1(dewpoint: Float) -> Result { - if !(273.0..=374.0).contains(&dewpoint) { - return Err(InputError::OutOfRange(String::from("dewpoint"))); - } - - // constants from the paper - let g: [Float; 8] = [ - -2991.2729, - -6017.0128, - 18.87643854, - -0.028354721, - 0.0000178383, - -0.00000000084150417, - 0.00000000000044412543, - 2.858487, - ]; - - let mut ln_p = g[7] * dewpoint.ln(); - - for i in 0..=6 { - ln_p += g[i] * dewpoint.powi(i as i32 - 2); - } - - Ok(ln_p.exp()) -} - -///Formula for computing vapour over ice pressure from dewpoint temperature. -///Should be used when accuracy is required as it is -///computationally expensive. -/// -///Derived by A. Wexler (1977) [(doi: 10.6028/jres.081A.003)](https://dx.doi.org/10.6028%2Fjres.081A.003). -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `dewpoint` range: 173K - 274K -#[cfg_attr(feature = "debug", logerr)] -pub fn wexler2(dewpoint: Float) -> Result { - if !(173.0..=274.0).contains(&dewpoint) { - return Err(InputError::OutOfRange(String::from("dewpoint"))); - } - - // constants from the paper - let big_k: [Float; 6] = [ - -5865.3696, - 22.241033, - 0.013749042, - -0.00003403177, - 0.000000026967687, - 0.6918651, - ]; - - let mut ln_p = big_k[5] * dewpoint.ln(); - - for j in 0..=4 { - ln_p += big_k[j] * dewpoint.powi(j as i32 - 1); - } - - Ok(ln_p.exp()) -} - -#[cfg(test)] -mod tests { - use crate::{ - tests_framework::{self, Argument}, - vapour_pressure, - }; - - #[test] - fn general1() { - assert!(tests_framework::test_with_2args( - &vapour_pressure::general1, - Argument { - name: "specific_humidity", - def_val: 0.022, - range: [0.00001, 2.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [100.0, 150_000.0] - }, - 3536.6680935251343 - )); - } - - #[test] - fn buck1() { - assert!(tests_framework::test_with_2args( - &vapour_pressure::buck1, - Argument { - name: "dewpoint", - def_val: 300.0, - range: [232.0, 324.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [100.0, 150_000.0] - }, - 3550.6603579471303 - )); - } - - #[test] - fn buck2() { - assert!(tests_framework::test_with_2args( - &vapour_pressure::buck2, - Argument { - name: "dewpoint", - def_val: 250.0, - range: [193.0, 274.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [100.0, 150_000.0] - }, - 76.38781790372722 - )); - } - - #[test] - fn buck3() { - assert!(tests_framework::test_with_2args( - &vapour_pressure::buck3, - Argument { - name: "dewpoint", - def_val: 300.0, - range: [253.0, 324.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [100.0, 150_000.0] - }, - 3548.5041048035896 - )); - } - - #[test] - fn buck4() { - assert!(tests_framework::test_with_2args( - &vapour_pressure::buck4, - Argument { - name: "dewpoint", - def_val: 250.0, - range: [223.0, 274.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [100.0, 150_000.0] - }, - 76.38685471836712 - )); - } - - #[test] - fn buck3_simplified() { - assert!(tests_framework::test_with_1arg( - &vapour_pressure::buck3_simplified, - Argument { - name: "dewpoint", - def_val: 300.0, - range: [253.0, 324.0] - }, - 3533.6421536199978 - )); - } - - #[test] - fn buck4_simplified() { - assert!(tests_framework::test_with_1arg( - &vapour_pressure::buck4_simplified, - Argument { - name: "dewpoint", - def_val: 250.0, - range: [223.0, 274.0] - }, - 76.04197508519536 - )); - } - - #[test] - fn tetens1() { - assert!(tests_framework::test_with_1arg( - &vapour_pressure::tetens1, - Argument { - name: "dewpoint", - def_val: 300.0, - range: [273.0, 353.0] - }, - 3533.969137160892 - )); - } - - #[test] - fn saturation_specific1() { - assert!(tests_framework::test_with_2args( - &vapour_pressure::saturation_specific1, - Argument { - name: "saturation_vapour_pressure", - def_val: 3550.0, - range: [0.0, 50_000.0] - }, - Argument { - name: "relative_humidity", - def_val: 0.5, - range: [0.0, 2.0] - }, - 1775.0 - )); - } - - #[test] - fn saturation_specific2() { - assert!(tests_framework::test_with_2args( - &vapour_pressure::saturation_specific2, - Argument { - name: "vapour_pressure", - def_val: 3000.0, - range: [0.0, 10_000.0] - }, - Argument { - name: "relative_humidity", - def_val: 0.5, - range: [0.00001, 2.0] - }, - 6000.0 - )); - } - - #[test] - fn wexler1() { - assert!(tests_framework::test_with_1arg( - &vapour_pressure::wexler1, - Argument { - name: "dewpoint", - def_val: 300.0, - range: [273.0, 374.0] - }, - 3535.4235919263083 - )); - } - - #[test] - fn wexler2() { - assert!(tests_framework::test_with_1arg( - &vapour_pressure::wexler2, - Argument { - name: "dewpoint", - def_val: 250.0, - range: [173.0, 274.0] - }, - 76.04351136780438 - )); - } -} diff --git a/src/vapour_pressure_deficit.rs b/src/vapour_pressure_deficit.rs deleted file mode 100644 index f42d17a..0000000 --- a/src/vapour_pressure_deficit.rs +++ /dev/null @@ -1,170 +0,0 @@ -//!Functions to calculate vapour pressure deficit in Pa. -//! -//!Vapour-pressure deficit, is the difference (deficit) between -//!the amount of moisture in the air and how much moisture the air can hold -//!when it is saturated ([Wikipedia](https://en.wikipedia.org/wiki/Vapour-pressure_deficit)). - -use crate::{errors::InputError, vapour_pressure}; -use crate::Float; - -#[cfg(feature="debug")] -use floccus_proc::logerr; - -///Formula for computing vapour pressure deficit from vapour pressure and saturation vapour pressure -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `vapour_pressure` range: 0Pa - 10000Pa -///Valid `saturation_vapour_pressure` range: 0Pa - 10000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn general1(vapour_pressure: Float, saturation_vapour_pressure: Float) -> Result { - if !(0.0..=50_000.0).contains(&vapour_pressure) { - return Err(InputError::OutOfRange(String::from("vapour_pressure"))); - } - - if !(0.0..=50_000.0).contains(&saturation_vapour_pressure) { - return Err(InputError::OutOfRange(String::from( - "saturation_vapour_pressure", - ))); - } - - Ok(saturation_vapour_pressure - vapour_pressure) -} - -///Formula for computing vapour pressure deficit from temperature, dewpoint and pressure -///using [`buck3`](vapour_pressure::buck3) function for vapour pressure calculation -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `vapour_pressure` range: 0Pa - 10000Pa -///Valid `saturation_vapour_pressure` range: 0Pa - 10000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn general2(temperature: Float, dewpoint: Float, pressure: Float) -> Result { - if !(253.0..=324.0).contains(&temperature) { - return Err(InputError::OutOfRange(String::from("temperature"))); - } - - if !(253.0..=324.0).contains(&dewpoint) { - return Err(InputError::OutOfRange(String::from("dewpoint"))); - } - - if !(100.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - let vapour_pressure = vapour_pressure::buck3(dewpoint, pressure)?; - let saturation_vapour_pressure = vapour_pressure::buck3(temperature, pressure)?; - - let result = general1(vapour_pressure, saturation_vapour_pressure)?; - - Ok(result) -} - -///Formula for computing vapour pressure deficit from temperature, relative humidity and pressure -///using [`buck3`](vapour_pressure::buck3) function for vapour pressure calculation -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `vapour_pressure` range: 0Pa - 10000Pa -///Valid `saturation_vapour_pressure` range: 0Pa - 10000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn general3( - temperature: Float, - relative_humidity: Float, - pressure: Float, -) -> Result { - if !(253.0..=319.0).contains(&temperature) { - return Err(InputError::OutOfRange(String::from("temperature"))); - } - - if !(0.05..=1.0).contains(&relative_humidity) { - return Err(InputError::OutOfRange(String::from("relative_humidity"))); - } - - if !(10000.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - let saturation_vapour_pressure = vapour_pressure::buck3(temperature, pressure)?; - let vapour_pressure = - vapour_pressure::saturation_specific1(saturation_vapour_pressure, relative_humidity)?; - - let result = general1(vapour_pressure, saturation_vapour_pressure)?; - - Ok(result) -} - -#[cfg(test)] -mod tests { - use crate::{ - tests_framework::{self, Argument}, - vapour_pressure_deficit, - }; - - #[test] - fn general1() { - assert!(tests_framework::test_with_2args( - &vapour_pressure_deficit::general1, - Argument { - name: "vapour_pressure", - def_val: 3000.0, - range: [0.0, 50_000.0] - }, - Argument { - name: "saturation_vapour_pressure", - def_val: 3550.0, - range: [0.0, 50_000.0] - }, - 550.0 - )); - } - - #[test] - fn general2() { - assert!(tests_framework::test_with_3args( - &vapour_pressure_deficit::general2, - Argument { - name: "temperature", - def_val: 300.0, - range: [253.0, 324.0] - }, - Argument { - name: "dewpoint", - def_val: 290.0, - range: [253.0, 324.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [100.0, 150_000.0] - }, - 1621.9415403325527 - )); - } - - #[test] - fn general3() { - assert!(tests_framework::test_with_3args( - &vapour_pressure_deficit::general3, - Argument { - name: "temperature", - def_val: 300.0, - range: [253.0, 319.0] - }, - Argument { - name: "relative_humidity", - def_val: 0.5, - range: [0.05, 1.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [10000.0, 150_000.0] - }, - 1774.2520524017948 - )); - } -} diff --git a/src/virtual_temperature.rs b/src/virtual_temperature.rs deleted file mode 100644 index 80c26ed..0000000 --- a/src/virtual_temperature.rs +++ /dev/null @@ -1,149 +0,0 @@ -//!Functions to calculate virtual temperature of air in K. -//! -//!In atmospheric thermodynamics, the virtual temperature of a moist air parcel is the temperature -//!at which a theoretical dry air parcel would have a total pressure and density equal -//!to the moist parcel of air ([Wikipedia](https://en.wikipedia.org/wiki/Virtual_temperature)). - -use crate::{constants::EPSILON, errors::InputError}; -use crate::Float; - -#[cfg(feature="debug")] -use floccus_proc::logerr; - -///Formula for computing virtual temperature from temperature and mixing ratio. -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `temperature` range: 173K - 373K\ -///Valid `mixing_ratio` range: 0.0000000001 - 0.5 -#[cfg_attr(feature = "debug", logerr)] -pub fn general1(temperature: Float, mixing_ratio: Float) -> Result { - if !(173.0..=354.0).contains(&temperature) { - return Err(InputError::OutOfRange(String::from("temperature"))); - } - - if !(0.000_000_000_1..=0.5).contains(&mixing_ratio) { - return Err(InputError::OutOfRange(String::from("mixing_ratio"))); - } - - let result = temperature * ((mixing_ratio + EPSILON) / (EPSILON * (1.0 + mixing_ratio))); - - Ok(result) -} - -///Formula for computing virtual temperature from air temperature, pressure and vapour pressure. -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `temperature` range: 173K - 373K\ -///Valid `pressure` range: 100Pa - 150000Pa\ -///Valid `vapour_pressure` range: 0Pa - 10000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn general2(temperature: Float, pressure: Float, vapour_pressure: Float) -> Result { - if !(173.0..=354.0).contains(&temperature) { - return Err(InputError::OutOfRange(String::from("temperature"))); - } - - if !(100.0..=150_000.0).contains(&pressure) { - return Err(InputError::OutOfRange(String::from("pressure"))); - } - - if !(0.0..=10_000.0).contains(&vapour_pressure) { - return Err(InputError::OutOfRange(String::from("vapour_pressure"))); - } - - let result = temperature / (1.0 - ((vapour_pressure / pressure) * (1.0 - EPSILON))); - - Ok(result) -} - -///Formula for computing virtual temperature from air temperature and specific humidity. -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `temperature` range: 173K - 373K\ -///Valid `specific_humidity` range: 100Pa - 150000Pa -#[cfg_attr(feature = "debug", logerr)] -pub fn general3(temperature: Float, specific_humidity: Float) -> Result { - if !(173.0..=354.0).contains(&temperature) { - return Err(InputError::OutOfRange(String::from("temperature"))); - } - - if !(0.000000001..=2.0).contains(&specific_humidity) { - return Err(InputError::OutOfRange(String::from("specific_humidity"))); - } - - let result = temperature * (1.0 + (specific_humidity * ((1.0 / EPSILON) - 1.0))); - - Ok(result) -} - -#[cfg(test)] -mod tests { - use crate::{ - tests_framework::{self, Argument}, - virtual_temperature, - }; - - #[test] - fn general1() { - assert!(tests_framework::test_with_2args( - &virtual_temperature::general1, - Argument { - name: "temperature", - def_val: 300.0, - range: [173.0, 354.0] - }, - Argument { - name: "mixing_ratio", - def_val: 0.022, - range: [0.000_000_000_1, 0.5] - }, - 303.9249219815806 - )); - } - - #[test] - fn general2() { - assert!(tests_framework::test_with_3args( - &virtual_temperature::general2, - Argument { - name: "temperature", - def_val: 300.0, - range: [173.0, 354.0] - }, - Argument { - name: "pressure", - def_val: 101325.0, - range: [100.0, 150_000.0] - }, - Argument { - name: "vapour_pressure", - def_val: 3550.0, - range: [0.0, 10_000.0] - }, - 304.0265941965307 - )); - } - - #[test] - fn general3() { - assert!(tests_framework::test_with_2args( - &virtual_temperature::general3, - Argument { - name: "temperature", - def_val: 300.0, - range: [173.0, 354.0] - }, - Argument { - name: "specific_humidity", - def_val: 0.022, - range: [0.000000001, 2.0] - }, - 304.0112702651753 - )); - } -} diff --git a/src/wet_bulb_potential_temperature.rs b/src/wet_bulb_potential_temperature.rs deleted file mode 100644 index 192e303..0000000 --- a/src/wet_bulb_potential_temperature.rs +++ /dev/null @@ -1,53 +0,0 @@ -//!Functions to calculate dry bulb potential temperature of unsaturated air in K. - -use crate::Float; -use crate::{ - constants::{C_P, R_D, ZERO_CELSIUS}, - errors::InputError, -}; - -#[cfg(feature="debug")] -use floccus_proc::logerr; - -///Formula for computing wet bulb potential temperature from equivalent potential temperature. -/// -///Derived by R. Davies-Jones (2008) [(doi:10.1175/2007MWR2224.1)](https://doi.org/10.1175/2007MWR2224.1) -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `temperature` range: 257K - 377K\ -#[cfg_attr(feature = "debug", logerr)] -pub fn davies_jones1(equivalent_potential_temperature: Float) -> Result { - if !(257.0..=377.0).contains(&equivalent_potential_temperature) { - return Err(InputError::OutOfRange(String::from( - "equivalent_potential_temperature", - ))); - } - - let lambda = C_P / R_D; - - let result = 45.114 - 51.489 * (ZERO_CELSIUS / equivalent_potential_temperature).powf(lambda); - Ok(result + ZERO_CELSIUS) -} - -#[cfg(test)] -mod tests { - use crate::{ - tests_framework::{self, Argument}, - wet_bulb_potential_temperature, - }; - - #[test] - fn davies_jones1() { - assert!(tests_framework::test_with_1arg( - &wet_bulb_potential_temperature::davies_jones1, - Argument { - name: "equivalent_potential_temperature", - def_val: 300.0, - range: [257.0, 377.0] - }, - 281.17941447108467 - )); - } -} diff --git a/src/wet_bulb_temperature.rs b/src/wet_bulb_temperature.rs deleted file mode 100644 index 1801500..0000000 --- a/src/wet_bulb_temperature.rs +++ /dev/null @@ -1,67 +0,0 @@ -//!Functions to calculate wet bulb temperature of unsaturated air in K. - -use crate::{constants::ZERO_CELSIUS, errors::InputError}; -use crate::Float; - -#[cfg(feature="debug")] -use floccus_proc::logerr; - -///Formula for computing wet bulb temperature pressure from dry bulb temperature and relative humidity. -/// -///Derived by R. Stull (2011) [(doi:10.1175/JAMC-D-11-0143.1)](https://doi.org/10.1175/JAMC-D-11-0143.1) -///Created with use of gene-expression programming.\ -///Result error is within −1K to +0.65K, with mean absolute error of 0.28K -/// -///# Errors -/// -///Returns [`InputError::OutOfRange`] when one of inputs is out of range.\ -///Valid `temperature` range: 253K - 324K\ -///Valid `relative_humidity` range: 0.05 - 0.99 -#[cfg_attr(feature = "debug", logerr)] -pub fn stull1(temperature: Float, relative_humidity: Float) -> Result { - if !(253.0..=324.0).contains(&temperature) { - return Err(InputError::OutOfRange(String::from("temperature"))); - } - - if !(0.05..=0.99).contains(&relative_humidity) { - return Err(InputError::OutOfRange(String::from("relative_humidity"))); - } - - //convert units - let temperature = temperature - ZERO_CELSIUS; - let relative_humidity = relative_humidity * 100.0; - - let result = (temperature * (0.151_977 * (relative_humidity + 8.313_659).sqrt()).atan()) - + (temperature + relative_humidity).atan() - - (relative_humidity - 1.676_331).atan() - + (0.003_918_38 * relative_humidity.powf(1.5) * (0.023_101 * relative_humidity).atan()) - - 4.686_035; - - Ok(result + ZERO_CELSIUS) -} - -#[cfg(test)] -mod tests { - use crate::{ - tests_framework::{self, Argument}, - wet_bulb_temperature, - }; - - #[test] - fn stull1() { - assert!(tests_framework::test_with_2args( - &wet_bulb_temperature::stull1, - Argument { - name: "temperature", - def_val: 300.0, - range: [253.0, 324.0] - }, - Argument { - name: "relative_humidity", - def_val: 0.5, - range: [0.05, 0.99] - }, - 292.73867410526674 - )); - } -}