darkmatter 0.0.2

Dark matter simulation engine — gravitational fields, particle dynamics, halo stability, and cosmological constants
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104

use std::f64::consts::PI;

pub fn erf_approx(x: f64) -> f64 {
    let t = 1.0 / (1.0 + 0.3275911 * x.abs());
    let poly = t
        * (0.254829592
            + t * (-0.284496736 + t * (1.421413741 + t * (-1.453152027 + t * 1.061405429))));
    let sign = if x < 0.0 { -1.0 } else { 1.0 };
    sign * (1.0 - poly * (-x * x).exp())
}

pub fn bessel_i0(x: f64) -> f64 {
    let ax = x.abs();
    if ax < 3.75 {
        let t = (x / 3.75).powi(2);
        1.0 + t
            * (3.5156229
                + t * (3.0899424
                    + t * (1.2067492 + t * (0.2659732 + t * (0.0360768 + t * 0.0045813)))))
    } else {
        let t = 3.75 / ax;
        (ax.exp() / ax.sqrt())
            * (0.39894228
                + t * (0.01328592
                    + t * (0.00225319
                        + t * (-0.00157565
                            + t * (0.00916281
                                + t * (-0.02057706
                                    + t * (0.02635537 + t * (-0.01647633 + t * 0.00392377))))))))
    }
}

pub fn bessel_i1(x: f64) -> f64 {
    let ax = x.abs();
    if ax < 3.75 {
        let t = (x / 3.75).powi(2);
        x * (0.5
            + t * (0.87890594
                + t * (0.51498869
                    + t * (0.15084934 + t * (0.02658733 + t * (0.00301532 + t * 0.00032411))))))
    } else {
        let t = 3.75 / ax;
        let result = (ax.exp() / ax.sqrt())
            * (0.39894228
                + t * (-0.03988024
                    + t * (-0.00362018
                        + t * (0.00163801
                            + t * (-0.01031555
                                + t * (0.02282967
                                    + t * (-0.02895312
                                        + t * (0.01787654 + t * (-0.00420059)))))))));
        if x < 0.0 { -result } else { result }
    }
}

pub fn bessel_k0(x: f64) -> f64 {
    if x <= 2.0 {
        let t = x * x / 4.0;
        -x.ln() * bessel_i0(x)
            + (-0.57721566
                + t * (0.42278420
                    + t * (0.23069756
                        + t * (0.03488590 + t * (0.00262698 + t * (0.00010750 + t * 0.00000740))))))
    } else {
        let t = 2.0 / x;
        ((-x).exp() / x.sqrt())
            * (1.25331414
                + t * (-0.07832358
                    + t * (0.02189568
                        + t * (-0.01062446
                            + t * (0.00587872 + t * (-0.00251540 + t * 0.00053208))))))
    }
}

pub fn bessel_k1(x: f64) -> f64 {
    if x <= 2.0 {
        let t = x * x / 4.0;
        x.ln() * bessel_i1(x)
            + (1.0 / x)
                * (1.0
                    + t * (0.15443144
                        + t * (-0.67278579
                            + t * (-0.18156897
                                + t * (-0.01919402 + t * (-0.00110404 + t * (-0.00004686)))))))
    } else {
        let t = 2.0 / x;
        ((-x).exp() / x.sqrt())
            * (1.25331414
                + t * (0.23498619
                    + t * (-0.03655620
                        + t * (0.01504268
                            + t * (-0.00780353 + t * (0.00325614 + t * (-0.00068245)))))))
    }
}

pub fn nfw_density(r: f64, rho_s: f64, r_s: f64) -> f64 {
    let x = r / r_s;
    rho_s / (x * (1.0 + x) * (1.0 + x))
}

pub fn nfw_enclosed_mass(r: f64, rho_s: f64, r_s: f64) -> f64 {
    let x = r / r_s;
    4.0 * PI * rho_s * r_s * r_s * r_s * ((1.0 + x).ln() - x / (1.0 + x))
}