192 lines
7.3 KiB
Rust
192 lines
7.3 KiB
Rust
use crate::mat::Material;
|
|
use crate::real::Real;
|
|
use crate::vec::Vec3;
|
|
|
|
#[cfg(feature = "serde")]
|
|
use serde::{Serialize, Deserialize};
|
|
|
|
#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))]
|
|
#[cfg_attr(feature = "fmt", derive(Debug))]
|
|
#[derive(Copy, Clone, Default, PartialEq)]
|
|
pub struct MHPgram<R> {
|
|
/// optimized form of mu_0^-1 * (1-mu_r^-1)
|
|
b_mult: R,
|
|
/// optimized form of h_intercept * (1-mu_r^-1)
|
|
m_offset: R,
|
|
/// Vertical range of the graph
|
|
pub max_m: R,
|
|
}
|
|
|
|
impl<R: Real> MHPgram<R> {
|
|
/// h_intercept: X coordinate at which M is always zero.
|
|
/// mu_r: relative mu value along the non-flat edges of the parallelogram.
|
|
pub fn new(h_intercept: R, mu_r: R, max_m: R) -> Self {
|
|
let one_minus_mu_r_inv = R::one() - R::one()/mu_r;
|
|
Self {
|
|
b_mult: R::mu0_inv() * one_minus_mu_r_inv,
|
|
m_offset: h_intercept * one_minus_mu_r_inv,
|
|
max_m,
|
|
}
|
|
}
|
|
|
|
pub fn h_intercept(&self) -> R {
|
|
R::mu0_inv() * self.m_offset / self.b_mult
|
|
}
|
|
pub fn mu_r(&self) -> R {
|
|
let mu_r_inv = R::one() - R::mu0()*self.b_mult;
|
|
R::one() / mu_r_inv
|
|
}
|
|
|
|
/// Return the new `M`
|
|
pub fn move_b(self, m: R, target_b: R) -> R {
|
|
// The point may exist outside the parallelogram.
|
|
// The right slope is defined by:
|
|
// B(H)/mu0 = h_intercept + (h-h_intercept)*mu_r
|
|
// "unoptimized" solution:
|
|
// let target_mh = target_b*MU0_INV;
|
|
// let h_on_right = (target_mh-self.h_intercept)/self.mu_r + self.h_intercept;
|
|
// let m_on_right = target_mh - h_on_right;
|
|
// Left:
|
|
// B(H)/mu0 = -h_intercept + (h+h_intercept)*mu_r
|
|
// "unoptimized" solution:
|
|
// let h_on_left = (target_mh+self.h_intercept)/self.mu_r - self.h_intercept;
|
|
// let m_on_left = target_mh - h_on_left;
|
|
|
|
let m_on_right = target_b*self.b_mult - self.m_offset;
|
|
let m_on_left = target_b*self.b_mult + self.m_offset;
|
|
|
|
if m_on_right >= m {
|
|
// if m_on_right <= self.max_m {
|
|
// // rightward edge movement
|
|
// m_on_right
|
|
// } else {
|
|
// // right of saturation
|
|
// self.max_m
|
|
// }
|
|
m_on_right.min_or_undefined(self.max_m)
|
|
} else if m_on_left <= m {
|
|
// if m_on_left >= -self.max_m {
|
|
// // leftward edge movement
|
|
// m_on_left
|
|
// } else {
|
|
// // left of saturation
|
|
// -self.max_m
|
|
// }
|
|
m_on_left.max_or_undefined(-self.max_m)
|
|
} else {
|
|
// interior movement
|
|
m
|
|
}
|
|
// this boils down to:
|
|
// m.clamp(m_on_left, m_on_right).clamp(-max_m, max_m)
|
|
}
|
|
}
|
|
|
|
impl<R: Real> Material<R> for MHPgram<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()),
|
|
)
|
|
}
|
|
}
|
|
|
|
/// MHPgram that's vaguely similar to Ferroxcube's 3R1 material
|
|
#[cfg_attr(feature = "serde", derive(Deserialize, Serialize))]
|
|
#[cfg_attr(feature = "fmt", derive(Debug))]
|
|
#[derive(Copy, Clone, Default, PartialEq)]
|
|
pub struct Ferroxcube3R1MH;
|
|
|
|
impl Ferroxcube3R1MH {
|
|
pub fn new() -> Self {
|
|
Self
|
|
}
|
|
}
|
|
|
|
impl<R: Real> Into<MHPgram<R>> for Ferroxcube3R1MH {
|
|
fn into(self) -> MHPgram<R> {
|
|
// TODO: how much (if any) penalty do we pay for this not being `const`?
|
|
MHPgram::new(25.0f32.cast(), 881.33f32.cast(), 44000.0f32.cast())
|
|
}
|
|
}
|
|
|
|
impl<R: Real> Material<R> for Ferroxcube3R1MH {
|
|
fn move_b_vec(&self, m: Vec3<R>, target_b: Vec3<R>) -> Vec3<R> {
|
|
let curve: MHPgram<R> = (*self).into();
|
|
curve.move_b_vec(m, target_b)
|
|
}
|
|
}
|
|
|
|
|
|
#[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);
|
|
}
|
|
|
|
#[test]
|
|
fn mh_curve_parameters() {
|
|
let curve = MHPgram::new(50.0, 101.0, 500.0);
|
|
assert_eq_approx(curve.h_intercept(), 50.0, 1e-3);
|
|
assert_eq_approx(curve.mu_r(), 101.0, 1e-3);
|
|
}
|
|
|
|
#[test]
|
|
fn mh_curve_edge_travel() {
|
|
let mu0 = f32::mu0();
|
|
let curve = MHPgram::new(50.0, 101.0, 500.0);
|
|
|
|
assert_eq_approx(curve.move_b(0.0, 151.0*mu0), 100.0, 0.1); // rightward travel along edge
|
|
assert_eq_approx(curve.move_b(100.0, 252.0*mu0), 200.0, 0.1);
|
|
assert_eq_approx(curve.move_b(200.0, 555.0*mu0), 500.0, 0.1);
|
|
assert_eq_approx(curve.move_b(500.0, 500.0*mu0), 500.0, 0.1); // back (leftward) travel to H=0
|
|
assert_eq_approx(curve.move_b(500.0, 455.0*mu0), 500.0, 0.1); // back (leftward) travel to H=-45
|
|
assert_eq_approx(curve.move_b(500.0, 354.0*mu0), 400.0, 0.1); // back (leftward) travel to H=-46
|
|
assert_eq_approx(curve.move_b(400.0, 253.0*mu0), 300.0, 0.1); // back (leftward) travel to H=-47
|
|
assert_eq_approx(curve.move_b(300.0, -50.0*mu0), 0.0, 0.1); // back (leftward) travel to H=-50
|
|
assert_eq_approx(curve.move_b(0.0, -151.0*mu0), -100.0, 0.1); // back (leftward) travel to H=-51
|
|
assert_eq_approx(curve.move_b(-100.0, -555.0*mu0), -500.0, 0.1); // back (leftward) travel to H=-55
|
|
assert_eq_approx(curve.move_b(-500.0, -456.0*mu0), -500.0, 0.1); // interior (rightward) travel to H=44
|
|
assert_eq_approx(curve.move_b(-500.0, -354.0*mu0), -400.0, 0.1); // interior (rightward) travel to H=46
|
|
assert_eq_approx(curve.move_b(-400.0, -253.0*mu0), -300.0, 0.1); // interior (rightward) travel to H=47
|
|
assert_eq_approx(curve.move_b(-300.0, 50.0*mu0), 0.0, 0.1); // interior (rightward) travel to H=50
|
|
}
|
|
|
|
#[test]
|
|
fn mh_curve_interior() {
|
|
let mu0 = f32::mu0();
|
|
let curve = MHPgram::new(50.0, 101.0, 500.0);
|
|
// max M (right) happens at H=55; B = 555*mu0;
|
|
// min M (right) happens at H=45; B = -455*mu0;
|
|
// max M (left) happens at H=-45; B = 455*mu0;
|
|
// min M (left) happens at H=-55; B = -555*mu0;
|
|
|
|
assert_eq_approx(curve.move_b(0.0, 40.0*mu0), 0.0, 0.1);
|
|
assert_eq_approx(curve.move_b(0.0, 50.0*mu0), 0.0, 0.1);
|
|
assert_eq_approx(curve.move_b(0.0, -40.0*mu0), 0.0, 0.1);
|
|
assert_eq_approx(curve.move_b(0.0, -50.0*mu0), 0.0, 0.1);
|
|
|
|
assert_eq_approx(curve.move_b(-400.0, -400.0*mu0), -400.0, 0.1); // rightward travel from H=-54 to NOT H=-53.5
|
|
assert_eq_approx(curve.move_b(-400.0, -355.0*mu0), -400.0, 0.1); // rightward travel from H=-54 to H=45
|
|
assert_eq_approx(curve.move_b(-400.0, -354.0*mu0), -400.0, 0.1); // rightward travel from H=-54 to H=46
|
|
assert_eq_approx(curve.move_b(-400.0, 5000.0*mu0), 500.0, 0.1); // rightward travel from H=-54 to H>>55
|
|
}
|
|
|
|
#[test]
|
|
fn mh_curve_exterior() {
|
|
let mu0 = f32::mu0();
|
|
let curve = MHPgram::new(50.0, 101.0, 500.0);
|
|
|
|
assert_eq_approx(curve.move_b(500.0, 556.0*mu0), 500.0, 0.1); // exterior travel to H=56
|
|
assert_eq_approx(curve.move_b(500.0, 5000.0*mu0), 500.0, 0.1); // exterior travel to H>>55
|
|
|
|
assert_eq_approx(curve.move_b(500.0, -5000.0*mu0), -500.0, 0.1); // exterior travel from H>>55 to H << -55
|
|
assert_eq_approx(curve.move_b(-501.0, 454.0*mu0), 400.0, 0.1); // exterior travel from H<<-55 (rounding error) to H=54
|
|
|
|
}
|
|
}
|
|
|