roche/planck.rs
1use pyo3::prelude::*;
2use crate::constants::{C, H, K};
3///
4/// Computes the Planck function Bnu = (2 h \nu^3/c^2)/(exp(h \nu/kT) - 1)
5/// as a function of wavelength and temperature. Output units are W/m**2/Hz/sr.
6///
7/// Arguments:
8///
9/// * `wave`: wavelength in nanometres
10/// * `temp`: temperature in K
11///
12#[pyfunction]
13pub fn planck(wave: f64, temp: f64) -> f64 {
14 let fac1: f64 = 2.0e27 * H * C;
15 let fac2: f64 = 1.0e9 * H * C / K;
16
17 let exponent: f64 = fac2 / (wave * temp);
18 if exponent > 40.0 {
19 fac1 * ((-exponent).exp()) / (wave * wave * wave)
20 } else {
21 fac1 / exponent.exp_m1() / (wave * wave * wave)
22 }
23}
24
25///
26/// Computes the logarithmic derivative of the Planck function Bnu wrt
27/// wavelength (i.e. d ln(Bnu) / d ln(lambda)) as a function of wavelength and temperature
28///
29/// Arguments:
30///
31/// * `wave`: wavelength in nanometres
32/// * `temp`: temperature in K
33///
34#[pyfunction]
35pub fn dplanck(wave: f64, temp: f64) -> f64 {
36 let fac2: f64 = 1.0e9 * H * C / K;
37
38 let exponent: f64 = fac2 / (wave * temp);
39 exponent / (1.0 - -exponent.exp()) - 3.0
40}
41
42///
43/// Computes the logarithmic derivative of the Planck function Bnu wrt
44/// T (i.e. d ln(Bnu) / d ln(T)) as a function of wavelength and temperature
45///
46/// Arguments:
47///
48/// * `wave`: wavelength in nanometres
49/// * `temp`: temperature in K
50#[pyfunction]
51pub fn dlpdlt(wave: f64, temp: f64) -> f64 {
52 let fac2: f64 = 1.0e9 * H * C / K;
53
54 let exponent: f64 = fac2 / (wave * temp);
55 exponent / (1.0 - (-exponent).exp())
56}