comp_flow/
mach_from.rs

1//! Collection of functions for isentropic compressible flow.
2
3use crate::{mach_to_a_ac, mach_to_pm_angle};
4use eqsolver::single_variable::FDNewton;
5use num::Float;
6
7/// Mach number for a given Prandtl-Meyer angle in radians.
8///
9/// <div class="warning">
10///
11/// This function uses Newton's method to solve for the Mach number. If this
12/// function must be called many times, it may be preferable to make a look up
13/// table with `mach_to_pm_angle` and interpolate those values.
14///
15/// </div>
16///
17/// # Examples
18///
19/// ```
20/// use comp_flow::mach_from_pm_angle;
21///
22/// assert_eq!(mach_from_pm_angle(0.4604136818474_f32, 1.4_f32), 2.0);
23/// assert_eq!(mach_from_pm_angle(0.0_f64, 1.4_f64),  1.00000022981460310);
24/// ```
25pub fn mach_from_pm_angle<F: Float>(pm_angle: F, gamma: F) -> F {
26    let f = |m| mach_to_pm_angle(m, gamma) - pm_angle;
27    let x0 = F::from(2.).unwrap();
28    FDNewton::new(f).solve(x0).unwrap()
29}
30
31/// Mach number for a given mach angle in radians.
32///
33/// # Examples
34///
35/// ```
36/// use comp_flow::mach_from_mach_angle;
37///
38/// assert_eq!(mach_from_mach_angle(0.5235988_f32), 2.0);
39/// assert_eq!(mach_from_mach_angle(1.5707963267948966_f64), 1.0);
40/// ```
41pub fn mach_from_mach_angle<F: Float>(mach_angle: F) -> F {
42    // TODO check for invalid input i.e. mach_angle > 90 deg
43    (F::one()) / mach_angle.sin()
44}
45
46/// Mach number for a given total temperature ratio.
47///
48/// # Examples
49///
50/// ```
51/// use comp_flow::mach_from_t_t0;
52///
53/// assert_eq!(mach_from_t_t0(1.0, 1.4), 0.0);
54/// assert_eq!(mach_from_t_t0(0.8333333333333334, 1.4), 1.0);
55/// assert_eq!(mach_from_t_t0(0.55555556_f32, 1.4), 2.0);
56/// ```
57pub fn mach_from_t_t0<F: Float>(t_t0: F, gamma: F) -> F {
58    let two = F::from(2.0).unwrap();
59    (two / (gamma - F::one()) * (F::one() / t_t0 - F::one())).sqrt()
60}
61
62/// Mach number for a given total pressure ratio.
63///
64/// # Examples
65///
66/// ```
67/// use comp_flow::mach_from_p_p0;
68///
69/// assert_eq!(mach_from_p_p0(1.0, 1.4), 0.0);
70/// assert_eq!(mach_from_p_p0(0.5282817877171742, 1.4), 1.0);
71/// assert_eq!(mach_from_p_p0(0.1278045254629509, 1.4), 2.0);
72/// ```
73pub fn mach_from_p_p0<F: Float>(p_p0: F, gamma: F) -> F {
74    let two = F::from(2.0).unwrap();
75    (two / (gamma - F::one()) * (p_p0.powf((F::one() - gamma) / gamma) - F::one())).sqrt()
76}
77
78/// Mach number for a given stagnation density ratio.
79///
80/// # Examples
81///
82/// ```
83/// use comp_flow::mach_from_rho_rho0;
84///
85/// assert_eq!(mach_from_rho_rho0(1.0, 1.4), 0.0);
86/// assert_eq!(mach_from_rho_rho0(0.633938145260609, 1.4), 1.0);
87/// assert_eq!(mach_from_rho_rho0(0.2300481458333117, 1.4), 2.0);
88/// ```
89pub fn mach_from_rho_rho0<F: Float>(rho_rho0: F, gamma: F) -> F {
90    let two = F::from(2.0).unwrap();
91    (two / (gamma - F::one()) * (rho_rho0.powf(F::one() - gamma) - F::one())).sqrt()
92}
93
94/// Mach number for a given critical area ratio.
95///
96/// <div class="warning">
97///
98/// This function uses a bisection algorythm to solve for the Mach number. If this
99/// function must be called many times, it may be preferable to make a look up
100/// table with `mach_to_a_ac` and interpolate those values.
101///
102/// </div>
103///
104/// # Examples
105///
106/// ```
107/// use comp_flow::mach_from_a_ac;
108///
109/// assert_eq!(mach_from_a_ac(1.6875000000000002_f64, 1.4, true), 2.0);
110/// assert_eq!(mach_from_a_ac(5.821828750000001_f64, 1.4, false), 0.1);
111/// assert_eq!(mach_from_a_ac(6.25, 1.4, false), 0.09307469911759117);
112/// assert_eq!(mach_from_a_ac(1.0, 1.4, false), 1.0);
113/// assert_eq!(mach_from_a_ac(1.0, 1.4, true), 1.0);
114/// ```
115pub fn mach_from_a_ac<F: Float>(a_ac: F, gamma: F, supersonic: bool) -> F {
116    if a_ac.is_one() {
117        return F::one();
118    }
119    let mut m_max: F;
120    let mut m_min: F;
121    if supersonic {
122        m_max = F::max_value();
123        m_min = F::one();
124    } else {
125        m_max = F::one();
126        m_min = F::zero();
127    }
128    let mut m_try = (m_max + m_min) / F::from(2.).unwrap();
129    let mut val: F;
130    let mut i: usize = 0;
131    let max_iter: usize = 10000;
132    while i < max_iter {
133        val = mach_to_a_ac(m_try, gamma);
134        if val == a_ac {
135            return m_try;
136        } else if val < a_ac {
137            if supersonic {
138                m_min = m_try;
139            } else {
140                m_max = m_try;
141            }
142        } else if val > a_ac {
143            if supersonic {
144                m_max = m_try;
145            } else {
146                m_min = m_try;
147            }
148        }
149        m_try = (m_max + m_min) / F::from(2.).unwrap();
150        i += 1;
151    }
152    return m_try;
153}
comp_flow/mach_from.rs

comp_flow/
mach_from.rs
