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

use std::ops::*;
use prelude::*;

pub struct Integration<V, F, N> {
    f:  F,
pub t:  N,
pub y:  V,
    h:  N,
    h_half: N,
    h_third: N,
    h_sixth: N,
    wrap: N
}
impl<V, F, N> Integration<V, F, N> where N: Real {
    pub fn new(f: F, s0: V, t0: N, dt: N, wrap: N) -> Integration<V, F, N> {
        Integration {
            f:  f,
            t:  t0,
            y:  s0,
            h:  dt,
            h_half: dt * N::frac(1, 2),
            h_third: dt * N::frac(1, 3),
            h_sixth: dt * N::frac(1, 6),
            wrap:   wrap
        }
    }
}
impl<V, F, N> Iterator for Integration<V, F, N> where
    V: Real<Scalar=N>,
    F: Fn(N, V) -> V,
    N: Real
{
    type Item = V;
    fn next(&mut self) -> Option<V> {
        let ref f = self.f;
        let t = self.t;
        let h = self.h;
        let h_half = self.h_half;
        let h_third = self.h_third;
        let h_sixth = self.h_sixth;
    
        let y = self.y;
        let k1 = f(t, y);
        let k2 = f(t + h_half, y + k1 * V::splat(h_half));
        let k3 = f(t + h_half, y + k2 * V::splat(h_half));
        let k4 = f(t + h, y + k3 * V::splat(h));
        
        self.y = y + (k1 + k4) * V::splat(h_sixth) + (k2 + k3) * V::splat(h_third);
        self.t = (t + h).wrap(self.wrap, self.wrap + self.wrap);
        
        Some(self.y)
    }
}