atm_refraction/air/
refractive.rs

1// calculation of the refractive index of air based on
2// https://emtoolbox.nist.gov/wavelength/Documentation.asp#ComparisonCiddorandEdlenEquations
3// Uses the modified Edlen equation
4
5use super::{dp_sv, p_sv};
6
7const A: f64 = 8342.54;
8const B: f64 = 2406147.0;
9const C: f64 = 15998.0;
10const D: f64 = 96095.43;
11const E: f64 = 0.601;
12const F: f64 = 0.00972;
13const G: f64 = 0.003661;
14
15/// Returns the air refractive index for the given wavelength (`lambda`), at the given pressure
16/// (`p`), temperature (`t`) and relative humidity (`rh`)
17pub fn air_index(lambda: f64, p: f64, t: f64, rh: f64) -> f64 {
18    let lambda_um = lambda * 1e6;
19    let s = 1.0 / lambda_um / lambda_um;
20    let t1 = t - 273.15;
21
22    let alpha = 1e-8 * (A + B / (130.0 - s) + C / (38.9 - s));
23    let beta = 1e-8 * E;
24    let gamma = -1e-8 * F;
25    let delta = D;
26    let epsilon = D * G;
27    let zeta = (3.7345 - s * 0.0401) * 1e-10;
28
29    let pv = rh / 100.0 * p_sv(t);
30
31    1.0 + alpha * p * (1.0 + beta * p + gamma * t1 * p) / (delta + epsilon * t1)
32        - (292.75 / t) * zeta * pv
33}
34
35/// Returns the derivative of the air refractive index for the given wavelength (`lambda`) as a
36/// function of pressure (`p`), temperature (`t`), relative humidity (`rh`) and their derivatives
37/// (`dp`, `dt`, `drh`)
38pub fn d_air_index(lambda: f64, p: f64, t: f64, rh: f64, dp: f64, dt: f64, drh: f64) -> f64 {
39    let lambda_um = lambda * 1e6;
40    let s = 1.0 / lambda_um / lambda_um;
41    let t1 = t - 273.15;
42
43    let alpha = 1e-8 * (A + B / (130.0 - s) + C / (38.9 - s));
44    let beta = 1e-8 * E;
45    let gamma = -1e-8 * F;
46    let delta = D;
47    let epsilon = D * G;
48    let zeta = (3.7345 - s * 0.0401) * 1e-10;
49
50    let pv = rh / 100.0 * p_sv(t);
51    let dpv = drh / 100.0 * p_sv(t) + rh / 100.0 * dp_sv(t) * dt;
52
53    alpha * dp * (1.0 + beta * p + gamma * t1 * p) / (delta + epsilon * t1)
54        + alpha
55            * p
56            * ((beta * dp + gamma * t1 * dp + gamma * p * dt) * (delta + epsilon * t1)
57                - epsilon * dt * (1.0 + beta * p + gamma * t1 * p))
58            / (delta + epsilon * t1)
59            / (delta + epsilon * t1)
60        + 292.75 / t / t * dt * zeta * pv
61        - 292.75 / t * zeta * dpv
62}