coremem_types: split MBPgram to its own file

This commit is contained in:
2022-07-18 18:20:27 -07:00
parent bb35a12c69
commit 6be7026c95
2 changed files with 148 additions and 129 deletions

View File

@@ -0,0 +1,145 @@
use crate::mat::Material;
use crate::real::Real;
use crate::vec::Vec3;
#[cfg(feature = "serde")]
use serde::{Serialize, Deserialize};
/// M(B) parallelogram
///
///```ignore
/// ____________
/// / /
/// / . /
/// / /
/// /___________/
/// ```
///
/// The `.` depicts (0, 0). X axis is B; y axis is M.
/// As B increases, M remains constant until it hits an edge.
/// Then M rises up to its max.
/// Same thing happens on the left edge, as B decreases and M falls to its min.
#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))]
#[cfg_attr(feature = "fmt", derive(Debug))]
#[derive(Copy, Clone, Default, PartialEq)]
pub struct MBPgram<R> {
/// X coordinate at which the upward slope starts
pub b_start: R,
/// X coordinate at which the upward slope ends
pub b_end: R,
/// Vertical range of the graph
pub max_m: R,
}
/// f(x0) = -ymax
/// f(x1) = ymax
fn eval_bounded_line<R: Real>(x: R, x0: R, x1: R, ymax: R) -> R {
let x = x - x0;
let x1 = x1 - x0;
// Now, f(0) = -ymax, f(x1) = ymax
let x = x.max_or_undefined(R::zero()).min_or_undefined(x1);
ymax * (R::two()*x/x1 - R::one())
}
impl<R> MBPgram<R> {
pub fn new(b_start: R, b_end: R, max_m: R) -> Self {
Self { b_start, b_end, max_m }
}
}
impl<R: Real> MBPgram<R> {
/// Return the new `M`
pub fn move_b(self, m: R, target_b: R) -> R {
// Evaluate the curve on the right:
let y_right = eval_bounded_line(target_b, self.b_start, self.b_end, self.max_m);
let y_left = eval_bounded_line(target_b, -self.b_end, -self.b_start, self.max_m);
m.max_or_undefined(y_right).min_or_undefined(y_left)
}
}
impl<R: Real> Material<R> for MBPgram<R> {
fn move_b_vec(&self, m: Vec3<R>, target_b: Vec3<R>) -> Vec3<R> {
Vec3::new(
self.move_b(m.x(), target_b.x()),
self.move_b(m.y(), target_b.y()),
self.move_b(m.z(), target_b.z()),
)
}
}
#[cfg(test)]
mod test {
use super::*;
fn assert_eq_approx(actual: f32, expected: f32, max_error: f32) {
assert!((actual - expected).abs() < max_error, "{} != {}", actual, expected);
}
fn assert_eq_approx_tight(actual: f32, expected: f32) {
assert_eq_approx(actual, expected, 1e-5);
}
#[test]
fn mb_curve_interior() {
let curve = MBPgram::new(4.0, 6.0, 20.0f32);
assert_eq_approx_tight(curve.move_b(0.0, 2.0), 0.0);
assert_eq_approx_tight(curve.move_b(0.0, 5.0), 0.0);
assert_eq_approx_tight(curve.move_b(1.0, 5.0), 1.0);
assert_eq_approx_tight(curve.move_b(20.0, 5.0), 20.0);
assert_eq_approx_tight(curve.move_b(-20.0, 4.0), -20.0);
assert_eq_approx_tight(curve.move_b(-20.0, -6.0), -20.0);
assert_eq_approx_tight(curve.move_b(20.0, -4.0), 20.0);
assert_eq_approx_tight(curve.move_b(10.0, -2.0), 10.0);
}
#[test]
fn mb_curve_exterior() {
let curve = MBPgram::new(4.0, 6.0, 20.0f32);
assert_eq_approx_tight(curve.move_b(0.0, 6.0), 20.0);
assert_eq_approx_tight(curve.move_b(0.0, 7.0), 20.0);
assert_eq_approx_tight(curve.move_b(0.0, -6.0), -20.0);
assert_eq_approx_tight(curve.move_b(0.0, -7.0), -20.0);
assert_eq_approx_tight(curve.move_b(2.0, -6.0), -20.0);
assert_eq_approx_tight(curve.move_b(20.0, -6.0), -20.0);
assert_eq_approx_tight(curve.move_b(20.0, -5.0), 0.0);
assert_eq_approx_tight(curve.move_b(20.0, -4.5), 10.0);
assert_eq_approx_tight(curve.move_b(-15.0, 4.5), -10.0);
assert_eq_approx_tight(curve.move_b(-15.0, 5.0), 0.0);
assert_eq_approx_tight(curve.move_b(-15.0, 5.5), 10.0);
assert_eq_approx_tight(curve.move_b(-15.0, 7.5), 20.0);
}
#[test]
fn mb_curve_3r1() {
// slope of 3r1 is about M=793210*B
// This is almost identical to iron (795615)!
let curve = MBPgram::new(-0.3899, 0.3900, 310000f32);
// magnetizing:
// v.s. 198893 in B(H) curve
assert_eq_approx(curve.move_b(0.0, 0.250), 198703.0, 1.0);
// v.s. 278321 in B(H) curve
assert_eq_approx(curve.move_b(198703.0, 0.350), 278201.0, 1.0);
assert_eq_approx(curve.move_b(278201.0, 0.390), 310000.0, 1.0);
// de-magnetizing:
// From saturation, decreasing B causes NO decrease in M: instead, it causes a decrease in
// H. This is probably BAD: in the B(H) curve, a large change in H always causes a large
// change in B. The movement of H here is likely to induce current, whereas it SHOULDN'T.
assert_eq_approx(curve.move_b(310000.0, 0.38995), 310000.0, 1.0);
// v.s. 258626 in B(H); H = 220
assert_eq_approx(curve.move_b(310000.0, 0.325), 258406.0, 1.0);
// here's where H crosses 0 (v.s. B=0.325 in the B(H) curve... quite a difference)
assert_eq_approx(curve.move_b(310000.0, 0.050), 39788.438, 1.0);
// v.s. 35.0 in B(H)
assert_eq_approx(curve.move_b(310000.0, 0.0), 39.75, 1.0);
// negative magnetization:
// erase the magnetization: H = -40
assert_eq_approx(curve.move_b(310000.0, -0.00005), 0.0, 1.0);
// the magnetization has been completely erased:
assert_eq_approx(curve.move_b(310000.0, -0.25), -198703.0, 1.0);
}
}

View File

@@ -4,6 +4,9 @@ use crate::vec::Vec3;
#[cfg(feature = "serde")] #[cfg(feature = "serde")]
use serde::{Serialize, Deserialize}; use serde::{Serialize, Deserialize};
mod mb_pgram;
pub use mb_pgram::MBPgram;
pub trait Material<R: Real>: Sized { pub trait Material<R: Real>: Sized {
fn conductivity(&self) -> Vec3<R> { fn conductivity(&self) -> Vec3<R> {
Default::default() Default::default()
@@ -36,67 +39,6 @@ impl<R: Real> Material<R> for IsomorphicConductor<R> {
} }
} }
/// M(B) parallelogram
///
///```ignore
/// ____________
/// / /
/// / . /
/// / /
/// /___________/
/// ```
///
/// The `.` depicts (0, 0). X axis is B; y axis is M.
/// As B increases, M remains constant until it hits an edge.
/// Then M rises up to its max.
/// Same thing happens on the left edge, as B decreases and M falls to its min.
#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))]
#[cfg_attr(feature = "fmt", derive(Debug))]
#[derive(Copy, Clone, Default, PartialEq)]
pub struct MBPgram<R> {
/// X coordinate at which the upward slope starts
pub b_start: R,
/// X coordinate at which the upward slope ends
pub b_end: R,
/// Vertical range of the graph
pub max_m: R,
}
/// f(x0) = -ymax
/// f(x1) = ymax
fn eval_bounded_line<R: Real>(x: R, x0: R, x1: R, ymax: R) -> R {
let x = x - x0;
let x1 = x1 - x0;
// Now, f(0) = -ymax, f(x1) = ymax
let x = x.max_or_undefined(R::zero()).min_or_undefined(x1);
ymax * (R::two()*x/x1 - R::one())
}
impl<R> MBPgram<R> {
pub fn new(b_start: R, b_end: R, max_m: R) -> Self {
Self { b_start, b_end, max_m }
}
}
impl<R: Real> MBPgram<R> {
/// Return the new `M`
pub fn move_b(self, m: R, target_b: R) -> R {
// Evaluate the curve on the right:
let y_right = eval_bounded_line(target_b, self.b_start, self.b_end, self.max_m);
let y_left = eval_bounded_line(target_b, -self.b_end, -self.b_start, self.max_m);
m.max_or_undefined(y_right).min_or_undefined(y_left)
}
}
impl<R: Real> Material<R> for MBPgram<R> {
fn move_b_vec(&self, m: Vec3<R>, target_b: Vec3<R>) -> Vec3<R> {
Vec3::new(
self.move_b(m.x(), target_b.x()),
self.move_b(m.y(), target_b.y()),
self.move_b(m.z(), target_b.z()),
)
}
}
#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))] #[cfg_attr(feature = "serde", derive(Deserialize, Serialize))]
#[cfg_attr(feature = "fmt", derive(Debug))] #[cfg_attr(feature = "fmt", derive(Debug))]
@@ -220,74 +162,6 @@ mod test {
assert!((actual - expected).abs() < max_error, "{} != {}", actual, expected); assert!((actual - expected).abs() < max_error, "{} != {}", actual, expected);
} }
fn assert_eq_approx_tight(actual: f32, expected: f32) {
assert_eq_approx(actual, expected, 1e-5);
}
#[test]
fn mb_curve_interior() {
let curve = MBPgram::new(4.0, 6.0, 20.0f32);
assert_eq_approx_tight(curve.move_b(0.0, 2.0), 0.0);
assert_eq_approx_tight(curve.move_b(0.0, 5.0), 0.0);
assert_eq_approx_tight(curve.move_b(1.0, 5.0), 1.0);
assert_eq_approx_tight(curve.move_b(20.0, 5.0), 20.0);
assert_eq_approx_tight(curve.move_b(-20.0, 4.0), -20.0);
assert_eq_approx_tight(curve.move_b(-20.0, -6.0), -20.0);
assert_eq_approx_tight(curve.move_b(20.0, -4.0), 20.0);
assert_eq_approx_tight(curve.move_b(10.0, -2.0), 10.0);
}
#[test]
fn mb_curve_exterior() {
let curve = MBPgram::new(4.0, 6.0, 20.0f32);
assert_eq_approx_tight(curve.move_b(0.0, 6.0), 20.0);
assert_eq_approx_tight(curve.move_b(0.0, 7.0), 20.0);
assert_eq_approx_tight(curve.move_b(0.0, -6.0), -20.0);
assert_eq_approx_tight(curve.move_b(0.0, -7.0), -20.0);
assert_eq_approx_tight(curve.move_b(2.0, -6.0), -20.0);
assert_eq_approx_tight(curve.move_b(20.0, -6.0), -20.0);
assert_eq_approx_tight(curve.move_b(20.0, -5.0), 0.0);
assert_eq_approx_tight(curve.move_b(20.0, -4.5), 10.0);
assert_eq_approx_tight(curve.move_b(-15.0, 4.5), -10.0);
assert_eq_approx_tight(curve.move_b(-15.0, 5.0), 0.0);
assert_eq_approx_tight(curve.move_b(-15.0, 5.5), 10.0);
assert_eq_approx_tight(curve.move_b(-15.0, 7.5), 20.0);
}
#[test]
fn mb_curve_3r1() {
// slope of 3r1 is about M=793210*B
// This is almost identical to iron (795615)!
let curve = MBPgram::new(-0.3899, 0.3900, 310000f32);
// magnetizing:
// v.s. 198893 in B(H) curve
assert_eq_approx(curve.move_b(0.0, 0.250), 198703.0, 1.0);
// v.s. 278321 in B(H) curve
assert_eq_approx(curve.move_b(198703.0, 0.350), 278201.0, 1.0);
assert_eq_approx(curve.move_b(278201.0, 0.390), 310000.0, 1.0);
// de-magnetizing:
// From saturation, decreasing B causes NO decrease in M: instead, it causes a decrease in
// H. This is probably BAD: in the B(H) curve, a large change in H always causes a large
// change in B. The movement of H here is likely to induce current, whereas it SHOULDN'T.
assert_eq_approx(curve.move_b(310000.0, 0.38995), 310000.0, 1.0);
// v.s. 258626 in B(H); H = 220
assert_eq_approx(curve.move_b(310000.0, 0.325), 258406.0, 1.0);
// here's where H crosses 0 (v.s. B=0.325 in the B(H) curve... quite a difference)
assert_eq_approx(curve.move_b(310000.0, 0.050), 39788.438, 1.0);
// v.s. 35.0 in B(H)
assert_eq_approx(curve.move_b(310000.0, 0.0), 39.75, 1.0);
// negative magnetization:
// erase the magnetization: H = -40
assert_eq_approx(curve.move_b(310000.0, -0.00005), 0.0, 1.0);
// the magnetization has been completely erased:
assert_eq_approx(curve.move_b(310000.0, -0.25), -198703.0, 1.0);
}
#[test] #[test]
fn mh_curve_parameters() { fn mh_curve_parameters() {
let curve = MHPgram::new(50.0, 101.0, 500.0); let curve = MHPgram::new(50.0, 101.0, 500.0);