comp_flow/
mach_from.rs

1//! Collection of functions for isentropic compressible flow.
2//! TODO: finish documentation
3
4use crate::{
5    der_mach_to_f_mcpt0, der_normal_mach2, der_normal_p02_p01, mach_to_a_ac, mach_to_f_mcpt,
6    mach_to_pm_angle, normal_mach2, normal_p02_p01,
7};
8use eqsolver::single_variable::{FDNewton, Newton};
9use num::Float;
10
11/// Mach number for a given Prandtl-Meyer angle in radians.
12///
13/// <div class="warning">
14///
15/// This function uses Newton's method to solve for the Mach number. If this
16/// function must be called many times, it may be preferable to make a look up
17/// table with `mach_to_pm_angle` and interpolate those values.
18///
19/// </div>
20///
21/// # Examples
22///
23/// ```
24/// use comp_flow::mach_from_pm_angle;
25///
26/// assert_eq!(mach_from_pm_angle(0.4604136818474_f32, 1.4_f32), 2.0);
27/// assert_eq!(mach_from_pm_angle(0.0_f64, 1.4_f64),  1.00000022981460310);
28/// ```
29pub fn mach_from_pm_angle<F: Float>(pm_angle: F, gamma: F) -> F {
30    let f = |m| mach_to_pm_angle(m, gamma) - pm_angle;
31    let x0 = F::from(2.).unwrap();
32    match FDNewton::new(f).solve(x0) {
33        Ok(x) => x,
34        Err(_) => F::nan(),
35    }
36}
37
38/// Mach number for a given mach angle in radians.
39///
40/// # Examples
41///
42/// ```
43/// use comp_flow::mach_from_mach_angle;
44///
45/// assert_eq!(mach_from_mach_angle(0.5235988_f32), 2.0);
46/// assert_eq!(mach_from_mach_angle(1.5707963267948966_f64), 1.0);
47/// ```
48pub fn mach_from_mach_angle<F: Float>(mach_angle: F) -> F {
49    // TODO: check for invalid input i.e. mach_angle > 90 deg
50    (F::one()) / mach_angle.sin()
51}
52
53/// Mach number for a given total temperature ratio.
54///
55/// # Examples
56///
57/// ```
58/// use comp_flow::mach_from_t_t0;
59///
60/// assert_eq!(mach_from_t_t0(1.0, 1.4), 0.0);
61/// assert_eq!(mach_from_t_t0(0.8333333333333334, 1.4), 1.0);
62/// assert_eq!(mach_from_t_t0(0.55555556_f32, 1.4), 2.0);
63/// ```
64pub fn mach_from_t_t0<F: Float>(t_t0: F, gamma: F) -> F {
65    let two = F::from(2.0).unwrap();
66    (two / (gamma - F::one()) * (F::one() / t_t0 - F::one())).sqrt()
67}
68
69/// Mach number for a given total temperature ratio.
70///
71pub fn mach_from_t0_t<F: Float>(t0_t: F, gamma: F) -> F {
72    let two = F::from(2.0).unwrap();
73    (two / (gamma - F::one()) * (t0_t - F::one())).sqrt()
74}
75
76/// Mach number for a given total pressure ratio.
77///
78/// # Examples
79///
80/// ```
81/// use comp_flow::mach_from_p_p0;
82///
83/// assert_eq!(mach_from_p_p0(1.0, 1.4), 0.0);
84/// assert_eq!(mach_from_p_p0(0.5282817877171742, 1.4), 1.0);
85/// assert_eq!(mach_from_p_p0(0.1278045254629509, 1.4), 2.0);
86/// ```
87pub fn mach_from_p_p0<F: Float>(p_p0: F, gamma: F) -> F {
88    let two = F::from(2.0).unwrap();
89    (two / (gamma - F::one()) * (p_p0.powf((F::one() - gamma) / gamma) - F::one())).sqrt()
90}
91
92pub fn mach_from_p0_p<F: Float>(p0_p: F, gamma: F) -> F {
93    let two = F::from(2.0).unwrap();
94    (two / (gamma - F::one()) * (p0_p.powf((gamma - F::one()) / gamma) - F::one())).sqrt()
95}
96
97/// Mach number for a given stagnation density ratio.
98///
99/// # Examples
100///
101/// ```
102/// use comp_flow::mach_from_rho_rho0;
103///
104/// assert_eq!(mach_from_rho_rho0(1.0, 1.4), 0.0);
105/// assert_eq!(mach_from_rho_rho0(0.633938145260609, 1.4), 1.0);
106/// assert_eq!(mach_from_rho_rho0(0.2300481458333117, 1.4), 2.0);
107/// ```
108pub fn mach_from_rho_rho0<F: Float>(rho_rho0: F, gamma: F) -> F {
109    let two = F::from(2.0).unwrap();
110    (two / (gamma - F::one()) * (rho_rho0.powf(F::one() - gamma) - F::one())).sqrt()
111}
112
113pub fn mach_from_rho0_rho<F: Float>(rho0_rho: F, gamma: F) -> F {
114    let two = F::from(2.0).unwrap();
115    (two / (gamma - F::one()) * (rho0_rho.powf(gamma - F::one()) - F::one())).sqrt()
116}
117
118pub fn mach_from_v_cpt0<F: Float>(v_cpt0: F, gamma: F) -> F {
119    let half = F::from(0.5).unwrap();
120    (F::one() / (gamma - F::one()) * v_cpt0.powi(2) / (F::one() - half * v_cpt0.powi(2))).sqrt()
121}
122
123pub fn mach_from_f_mcpt0<F: Float>(f_mcpt0: F, gamma: F, supersonic: bool) -> F {
124    let f = |m| mach_to_f_mcpt(m, gamma) - f_mcpt0;
125    let df = |m| der_mach_to_f_mcpt0(m, gamma) - f_mcpt0;
126    let x0 = match supersonic {
127        true => F::from(1.5).unwrap(),
128        false => F::from(0.5).unwrap(),
129    };
130    match Newton::new(f, df).with_itermax(100).solve(x0) {
131        Ok(x) => x,
132        Err(_) => F::nan(),
133    }
134}
135
136pub fn mach_from_mcpt0_ap0<F: Float>(mcpt0_ap0: F, gamma: F, supersonic: bool) -> F {
137    let half = F::from(0.5).unwrap();
138    let gm1 = gamma - F::one();
139    let gp1 = gamma + F::one();
140    let gm1_2 = gm1 * half;
141    let gp1_2 = gp1 * half;
142    let m_gp1_gm1_2 = gp1 / gm1 * -half;
143    let g_sq_gm1 = gamma / gm1.sqrt();
144    let fcrit = g_sq_gm1 * (F::one() + gm1_2).powf(m_gp1_gm1_2);
145    if mcpt0_ap0 > fcrit {
146        return F::nan();
147    }
148    let mut m_try: F;
149    if supersonic {
150        m_try = F::from(1.5).unwrap();
151    } else {
152        m_try = half;
153    }
154    let mut k = 0;
155    let mut err = F::infinity();
156    let tol = F::from(1e-6).unwrap();
157    while err > tol && k < 100 {
158        k = k + 1;
159        // USE NEWTON'S METHOD FOR NEW GUESS OF MA
160        let to_t = F::one() + gm1_2 * m_try.powi(2);
161        let f = g_sq_gm1 * m_try * to_t.powf(m_gp1_gm1_2);
162        let df = f * (F::one() - gp1_2 * m_try.powi(2) / to_t) / m_try;
163        let m_new = m_try - (f - mcpt0_ap0) / df;
164        err = (m_new - m_try).abs();
165        m_try = m_new;
166    }
167    let mach: F;
168    if err < tol {
169        mach = m_try;
170    } else {
171        mach = F::nan();
172    }
173    mach
174}
175
176pub fn mach_from_mcpt0_ap<F: Float>(mcpt0_ap: F, gamma: F) -> F {
177    let half = F::from(0.5).unwrap();
178    let gm1 = gamma - F::one();
179    let gm1_2 = gm1 * half;
180    let g_sq_gm1 = gamma / gm1.sqrt();
181    let tol = F::from(1e-6).unwrap();
182    let mut m_try = F::from(0.5).unwrap();
183    let mut err = F::infinity();
184    let mut i: i32 = 0;
185    while err > tol && i < 100 {
186        i = i + 1;
187        // use newton's method
188        let to_t = F::one() + gm1_2 * m_try.powi(2);
189        let f = g_sq_gm1 * m_try * to_t.sqrt();
190        let df = g_sq_gm1 * (to_t.powi(2) - gm1_2 * m_try.powi(2) / to_t.sqrt());
191        let m_new = m_try - (f - mcpt0_ap) / df;
192        err = (m_new - m_try).abs();
193        m_try = m_new;
194    }
195    let mach: F;
196    if err < tol {
197        mach = m_try;
198    } else {
199        mach = F::nan();
200    }
201    mach
202}
203
204/// Mach number for a given critical area ratio.
205///
206/// <div class="warning">
207///
208/// This function uses a bisection algorythm to solve for the Mach number. If this
209/// function must be called many times, it may be preferable to make a look up
210/// table with `mach_to_a_ac` and interpolate those values.
211///
212/// </div>
213///
214/// # Examples
215///
216/// ```
217/// use comp_flow::mach_from_a_ac;
218///
219/// assert_eq!(mach_from_a_ac(1.6875000000000002_f64, 1.4, true), 2.0);
220/// assert_eq!(mach_from_a_ac(5.821828750000001_f64, 1.4, false), 0.1);
221/// assert_eq!(mach_from_a_ac(6.25, 1.4, false), 0.09307469911759117);
222/// assert_eq!(mach_from_a_ac(1.0, 1.4, false), 1.0);
223/// assert_eq!(mach_from_a_ac(1.0, 1.4, true), 1.0);
224/// ```
225pub fn mach_from_a_ac<F: Float>(a_ac: F, gamma: F, supersonic: bool) -> F {
226    if a_ac.is_one() {
227        return F::one();
228    }
229    let mut m_max: F;
230    let mut m_min: F;
231    if supersonic {
232        m_max = F::max_value();
233        m_min = F::one();
234    } else {
235        m_max = F::one();
236        m_min = F::zero();
237    }
238    let mut m_try = (m_max + m_min) / F::from(2.).unwrap();
239    let mut val: F;
240    let mut i: usize = 0;
241    let max_iter: usize = 10000;
242    while i < max_iter {
243        val = mach_to_a_ac(m_try, gamma);
244        if val == a_ac {
245            return m_try;
246        } else if val < a_ac {
247            if supersonic {
248                m_min = m_try;
249            } else {
250                m_max = m_try;
251            }
252        } else if val > a_ac {
253            if supersonic {
254                m_max = m_try;
255            } else {
256                m_min = m_try;
257            }
258        }
259        m_try = (m_max + m_min) / F::from(2.).unwrap();
260        i += 1;
261    }
262    return m_try;
263}
264
265pub fn mach_from_normal_mach2<F: Float>(mach2: F, gamma: F) -> F {
266    if mach2 > F::one() {
267        return F::nan();
268    }
269    let f = |m| normal_mach2(m, gamma) - mach2;
270    let df = |m| der_normal_mach2(m, gamma) - mach2;
271    let x0 = F::from(0.5).unwrap();
272    match Newton::new(f, df).with_itermax(100).solve(x0) {
273        Ok(x) => x,
274        Err(_) => F::nan(),
275    }
276}
277
278pub fn mach_from_normal_p02_p01<F: Float>(p02_p01: F, gamma: F) -> F {
279    let f = |m| normal_p02_p01(m, gamma) - p02_p01;
280    let df = |m| der_normal_p02_p01(m, gamma) - p02_p01;
281    let x0 = F::from(1.5).unwrap();
282    match Newton::new(f, df).with_itermax(100).solve(x0) {
283        Ok(x) => x,
284        Err(_) => F::nan(),
285    }
286}
comp_flow/mach_from.rs

comp_flow/
mach_from.rs
