cfsem 5.1.0

Quasi-steady electromagnetics including filamentized approximations, Biot-Savart, and Grad-Shafranov.
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
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155

//! Test utilities

use core::f64::consts::{E, PI};

use itertools::{self, Itertools};

/// Divide-by-zero-resistant approximate comparison
pub(crate) fn approx(truth: f64, val: f64, rtol: f64, atol: f64) -> bool {
    let abs_err = (val - truth).abs();
    let lim = rtol * truth.abs() + atol;
    abs_err < lim
}

/// Evenly spaced values from start to end
pub(crate) fn linspace(start: f64, end: f64, n: usize) -> Vec<f64> {
    (0..n)
        .map(|i| start + (i as f64 / (n - 1) as f64) * (end - start))
        .collect::<Vec<f64>>()
}

/// Dense N-dimensional meshgrid; cartesian product of
/// input grids, in canonical array order
pub(crate) fn meshgrid(grids: &[&[f64]]) -> Vec<Vec<f64>> {
    let ngrids = grids.len();

    // Interleaved cartesian product
    let interleaved = grids
        .iter()
        .map(|&x| x.iter())
        .multi_cartesian_product()
        .flatten()
        .cloned()
        .collect_vec();

    // Deinterleave
    let mut meshes = Vec::with_capacity(ngrids);
    for i in 0..ngrids {
        meshes.push(
            interleaved[i..]
                .iter()
                .step_by(ngrids)
                .cloned()
                .collect_vec(),
        )
    }

    meshes
}

/// First-order forward difference; returns n-1 sized output
pub(crate) fn diff(v: &[f64]) -> Vec<f64> {
    v[1..]
        .iter()
        .zip(v[0..v.len() - 1].iter())
        .map(|(&b, &a)| b - a)
        .collect::<Vec<f64>>()
}

/// Convert circular filament to ndiscr-1 piecewise linear segments
/// (ndiscr points)
pub(crate) fn discretize_circular_filament(
    r: f64,
    z: f64,
    ndiscr: usize,
) -> (Vec<f64>, Vec<f64>, Vec<f64>) {
    let x: Vec<f64> = linspace(0.0, 2.0 * PI, ndiscr)
        .iter()
        .map(|v| r * libm::cos(*v))
        .collect();
    let y: Vec<f64> = linspace(0.0, 2.0 * PI, ndiscr)
        .iter()
        .map(|v| r * libm::sin(*v))
        .collect();
    let z: Vec<f64> = (0..ndiscr).map(|_| z).collect();
    (x, y, z)
}

/// Make an example helical path for testing
pub(crate) fn example_helix() -> (Vec<f64>, Vec<f64>, Vec<f64>) {
    // Make a slightly tilted helical piecewise-linear filament
    let n = 10_000;
    let xc = [0.1, -0.1]; // Start and end of centerline path
    let yc = [-0.05, 0.2];
    let zc = [-2.0, 2.0];

    let xc: Vec<f64> = linspace(xc[0], xc[1], n);
    let yc: Vec<f64> = linspace(yc[0], yc[1], n);
    let zc: Vec<f64> = linspace(zc[0], zc[1], n);

    let mut x = xc.clone();
    let mut y = xc.clone();
    let mut z = xc.clone();
    crate::mesh::filament_helix_path(
        (&xc, &yc, &zc),
        (2.0 * E / 3.0, 0.0, 0.0),
        0.5,
        0.0,
        (&mut x, &mut y, &mut z),
    )
    .unwrap();

    (x, y, z)
}

/// Example set of circular filaments with (r, z, n_turns) values
pub(crate) fn example_circular_filaments() -> (Vec<f64>, Vec<f64>, Vec<f64>) {
    // Make some circular filaments
    let r = 1.0 / PI; // [m] some number
    let z = 1.0 / E; // [m] some number

    let rfil = [r, r + E / 4.0, r + E / 2.0].to_vec();
    let zfil = [z, -z, 0.0].to_vec();
    let nfil = [PI, E, E / 2.0].to_vec();

    (rfil, zfil, nfil)
}

// Midpoints between each pair of points.
// Output is 1 index shorter than input.
pub(crate) fn midpoints(x: &[f64]) -> Vec<f64> {
    x[..]
        .iter()
        .zip(x[1..].iter())
        .map(|(this, next)| (this + next) / 2.0)
        .collect()
}

#[cfg(test)]
mod test {
    use super::*;

    /// Check that meshgrid returns the correct shape and values for
    /// a 2-dimensional input.
    #[test]
    fn test_meshgrid_2d() {
        let x = [0.0, 1.0, 2.0];
        let y = [3.0, 4.0, 5.0, 6.0];
        let m = meshgrid(&[&x, &y]);

        // Check that the meshgrid has the correct shape
        assert_eq!(2, m.len()); // should be number of dimensions
        for mi in m.iter() {
            assert_eq!(x.len() * y.len(), mi.len());
        }

        // Check that the meshgrid contains the correct values
        for (i, &xi) in x.iter().enumerate() {
            for (j, &yj) in y.iter().enumerate() {
                let k = i * y.len() + j;
                assert_eq!(xi, m[0][k]);
                assert_eq!(yj, m[1][k]);
            }
        }
    }
}