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

use crate::kernel::Kernel;
use na::RealField;

/// The cubic spline smoothing kernel.
///
/// See https://pysph.readthedocs.io/en/latest/reference/kernels.html
#[derive(Copy, Clone, Debug)]
pub struct CubicSplineKernel;

impl Kernel for CubicSplineKernel {
    fn scalar_apply<N: RealField>(r: N, h: N) -> N {
        assert!(r >= N::zero());

        #[cfg(feature = "dim2")]
        let normalizer = na::convert::<_, N>(40.0 / 7.0) / (N::pi() * h * h);
        #[cfg(feature = "dim3")]
        let normalizer = na::convert::<_, N>(8.0) / (N::pi() * h * h * h);

        let _2: N = na::convert(2.0);
        let q = r / h;

        let rhs = if q <= na::convert(0.5) {
            let q2 = q * q;
            N::one() + (q2 * q - q2) * na::convert(6.0)
        } else if q <= N::one() {
            (N::one() - q).powi(3) * _2
        } else {
            N::zero()
        };

        normalizer * rhs

        /*
        let q = r / h;
        #[cfg(feature = "dim2")]
            let normalizer = na::convert::<_, N>(10.0 / 7.0) / (N::pi() * h * h);
        #[cfg(feature = "dim3")]
            let normalizer = N::one() / (N::pi() * h * h * h);

        let _2: N = na::convert(2.0);
        let _3: N = na::convert(3.0);
        let rhs = if q <= N::one() {
            N::one() - _3 / _2 * q * q * (N::one() - q / _2)
        } else if q <= _2 {
            (_2 - q).powi(3) / na::convert(4.0)
        } else {
            N::zero()
        };

        normalizer * rhs
        */
    }

    fn scalar_apply_diff<N: RealField>(r: N, h: N) -> N {
        assert!(r >= N::zero());

        #[cfg(feature = "dim2")]
        let normalizer = na::convert::<_, N>(40.0 / 7.0) / (N::pi() * h * h);
        #[cfg(feature = "dim3")]
        let normalizer = na::convert::<_, N>(8.0) / (N::pi() * h * h * h);

        let _2: N = na::convert(2.0);
        let _3: N = na::convert(3.0);
        let q = r / h;

        let rhs = if q <= na::convert(0.5) {
            (q * _3 - _2) * q * na::convert(6.0)
        } else if q <= N::one() {
            let one_q = N::one() - q;
            -one_q * one_q * na::convert(6.0)
        } else {
            N::zero()
        };

        normalizer * rhs / h

        /*
        let q = r / h;
        #[cfg(feature = "dim2")]
            let normalizer = na::convert::<_, N>(10.0 / 7.0) / (N::pi() * h * h);
        #[cfg(feature = "dim3")]
            let normalizer = N::one() / (N::pi() * h * h * h);

        let _2: N = na::convert(2.0);
        let _3: N = na::convert(3.0);
        let rhs = if q <= N::one() {
            -_3 * q * (N::one() - q * na::convert(3.0 / 4.0))
        } else if q <= _2 {
            -(_2 - q).powi(2) * na::convert(3.0 / 4.0)
        } else {
            N::zero()
        };

        normalizer * rhs
        */
    }
}