stim: break apart into smaller modules
This commit is contained in:
parent
9301734fcf
commit
7f089bad45
|
|
@ -1,756 +0,0 @@
|
|||
use crate::real::{self, Real};
|
||||
use crate::cross::vec::Vec3;
|
||||
use crate::geom::{Coord as _, HasCrossSection, Index, Meters, Region};
|
||||
use coremem_cross::dim::DimSlice;
|
||||
use coremem_cross::vec::Vec3u;
|
||||
|
||||
use std::borrow::Cow;
|
||||
use std::ops::Deref;
|
||||
|
||||
use rand;
|
||||
|
||||
/// field densities
|
||||
#[derive(Clone, Copy, Debug, Default, PartialEq)]
|
||||
pub struct Fields<R> {
|
||||
pub e: Vec3<R>,
|
||||
pub h: Vec3<R>,
|
||||
}
|
||||
|
||||
impl<R> Fields<R> {
|
||||
pub fn new_eh(e: Vec3<R>, h: Vec3<R>) -> Self {
|
||||
Self { e, h}
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> Fields<R> {
|
||||
pub fn e(&self) -> Vec3<R> {
|
||||
self.e
|
||||
}
|
||||
pub fn h(&self) -> Vec3<R> {
|
||||
self.h
|
||||
}
|
||||
pub fn new_e(e: Vec3<R>) -> Self {
|
||||
Self::new_eh(e, Vec3::zero())
|
||||
}
|
||||
pub fn new_h(h: Vec3<R>) -> Self {
|
||||
Self::new_eh(Vec3::zero(), h)
|
||||
}
|
||||
pub fn elem_mul(self, other: FieldMags<R>) -> Fields<R> {
|
||||
Fields {
|
||||
e: self.e * other.e,
|
||||
h: self.h * other.h,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> std::ops::AddAssign for Fields<R> {
|
||||
fn add_assign(&mut self, other: Self) {
|
||||
self.e += other.e;
|
||||
self.h += other.h;
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> std::ops::Add for Fields<R> {
|
||||
type Output = Self;
|
||||
fn add(mut self, other: Self) -> Self::Output {
|
||||
self += other;
|
||||
self
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> std::ops::Mul<R> for Fields<R> {
|
||||
type Output = Self;
|
||||
fn mul(self, scale: R) -> Self::Output {
|
||||
Fields {
|
||||
e: self.e * scale,
|
||||
h: self.h * scale,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// field magnitude densities (really, signed magnitude)
|
||||
#[derive(Clone, Copy, Debug, Default, PartialEq)]
|
||||
pub struct FieldMags<R> {
|
||||
pub e: R,
|
||||
pub h: R,
|
||||
}
|
||||
|
||||
impl<R: Real> std::ops::AddAssign for FieldMags<R> {
|
||||
fn add_assign(&mut self, other: Self) {
|
||||
self.e += other.e;
|
||||
self.h += other.h;
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> std::ops::Add for FieldMags<R> {
|
||||
type Output = Self;
|
||||
fn add(mut self, other: Self) -> Self::Output {
|
||||
self += other;
|
||||
self
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> std::ops::Mul<R> for FieldMags<R> {
|
||||
type Output = Self;
|
||||
fn mul(self, scale: R) -> Self::Output {
|
||||
FieldMags {
|
||||
e: self.e * scale,
|
||||
h: self.h * scale,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl<R> FieldMags<R> {
|
||||
pub fn new_eh(e: R, h: R) -> Self {
|
||||
Self { e, h }
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> FieldMags<R> {
|
||||
pub fn e(&self) -> R {
|
||||
self.e
|
||||
}
|
||||
pub fn h(&self) -> R {
|
||||
self.h
|
||||
}
|
||||
pub fn new_e(e: R) -> Self {
|
||||
Self::new_eh(e, R::zero())
|
||||
}
|
||||
pub fn new_h(h: R) -> Self {
|
||||
Self::new_eh(R::zero(), h)
|
||||
}
|
||||
pub fn elem_mul(self, other: Self) -> Self {
|
||||
FieldMags {
|
||||
e: self.e * other.e,
|
||||
h: self.h * other.h,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pub trait Stimulus<R: Real>: Sync {
|
||||
/// Return the (E, H) field which should be added PER-SECOND to the provided position/time.
|
||||
fn at(&self, t_sec: R, feat_size: R, loc: Index) -> Fields<R>;
|
||||
|
||||
/// compute the value of this stimulus across all the simulation space
|
||||
fn rendered<'a>(
|
||||
&'a self, scale: R, t_sec: R, feature_size: R, dim: Vec3u
|
||||
) -> Cow<'a, RenderedStimulus<R>> {
|
||||
Cow::Owned(render_stim(self, scale, t_sec, feature_size, dim))
|
||||
}
|
||||
}
|
||||
fn render_stim<R: Real, S: Stimulus<R> + ?Sized>(
|
||||
stim: &S, scale: R, t_sec: R, feature_size: R, dim: Vec3u
|
||||
) -> RenderedStimulus<R> {
|
||||
let dim_len = dim.product_sum_usize();
|
||||
let mut e = Vec::new();
|
||||
e.resize_with(dim_len, Default::default);
|
||||
let mut h = Vec::new();
|
||||
h.resize_with(dim_len, Default::default);
|
||||
|
||||
rayon::scope(|s| {
|
||||
let mut undispatched_e = &mut e[..];
|
||||
let mut undispatched_h = &mut h[..];
|
||||
for z in 0..dim.z() {
|
||||
for y in 0..dim.y() {
|
||||
let (this_e, this_h);
|
||||
(this_e, undispatched_e) = undispatched_e.split_at_mut(dim.x() as usize);
|
||||
(this_h, undispatched_h) = undispatched_h.split_at_mut(dim.x() as usize);
|
||||
s.spawn(move |_| {
|
||||
for (x, (out_e, out_h)) in this_e.iter_mut().zip(this_h.iter_mut()).enumerate() {
|
||||
|
||||
let Fields { e, h } = stim.at(t_sec, feature_size, Index::new(x as u32, y, z));
|
||||
*out_e = e * scale;
|
||||
*out_h = h * scale;
|
||||
}
|
||||
});
|
||||
}
|
||||
}
|
||||
});
|
||||
|
||||
let field_e = DimSlice::new(dim, e);
|
||||
let field_h = DimSlice::new(dim, h);
|
||||
RenderedStimulus::new(
|
||||
field_e, field_h, scale, feature_size, t_sec
|
||||
)
|
||||
}
|
||||
|
||||
/// a static vector field. different value at each location, but constant in time.
|
||||
/// often used as a building block by wrapping it in something which modulates the fields over
|
||||
/// time.
|
||||
pub trait VectorField<R> {
|
||||
fn at(&self, feat_size: R, loc: Index) -> Fields<R>;
|
||||
}
|
||||
|
||||
// uniform vector field
|
||||
impl<R: Real> VectorField<R> for Fields<R> {
|
||||
fn at(&self, _feat_size: R, _loc: Index) -> Fields<R> {
|
||||
*self
|
||||
}
|
||||
}
|
||||
|
||||
// could broaden this and implement directly on T, but blanket impls
|
||||
// are unwieldy
|
||||
impl<R: Real, T> VectorField<R> for DimSlice<T>
|
||||
where
|
||||
DimSlice<T>: core::ops::Index<Vec3u, Output=Fields<R>>
|
||||
{
|
||||
fn at(&self, _feat_size: R, loc: Index) -> Fields<R> {
|
||||
self[loc.into()]
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Clone)]
|
||||
pub struct RenderedStimulus<R> {
|
||||
e: DimSlice<Vec<Vec3<R>>>,
|
||||
h: DimSlice<Vec<Vec3<R>>>,
|
||||
scale: R,
|
||||
feature_size: R,
|
||||
t_sec: R,
|
||||
}
|
||||
|
||||
impl<R> RenderedStimulus<R> {
|
||||
pub fn new(
|
||||
e: DimSlice<Vec<Vec3<R>>>,
|
||||
h: DimSlice<Vec<Vec3<R>>>,
|
||||
scale: R,
|
||||
feature_size: R,
|
||||
t_sec: R,
|
||||
) -> Self {
|
||||
Self { e, h, scale, feature_size, t_sec }
|
||||
}
|
||||
pub fn e<'a>(&'a self) -> DimSlice<&'a [Vec3<R>]> {
|
||||
self.e.as_ref()
|
||||
}
|
||||
pub fn h<'a>(&'a self) -> DimSlice<&'a [Vec3<R>]> {
|
||||
self.h.as_ref()
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> RenderedStimulus<R> {
|
||||
pub fn scale(&self) -> R {
|
||||
self.scale
|
||||
}
|
||||
pub fn feature_size(&self) -> R {
|
||||
self.feature_size
|
||||
}
|
||||
pub fn time(&self) -> R {
|
||||
self.t_sec
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> VectorField<R> for RenderedStimulus<R> {
|
||||
fn at(&self, _feat_size: R, loc: Index) -> Fields<R> {
|
||||
Fields::new_eh(self.e[loc.into()], self.h[loc.into()])
|
||||
}
|
||||
}
|
||||
impl<R: Real> Stimulus<R> for RenderedStimulus<R> {
|
||||
fn at(&self, _t_sec: R, _feat_size: R, loc: Index) -> Fields<R> {
|
||||
Fields::new_eh(self.e[loc.into()], self.h[loc.into()])
|
||||
}
|
||||
fn rendered<'a>(
|
||||
&'a self, scale: R, t_sec: R, feature_size: R, dim: Vec3u
|
||||
) -> Cow<'a, RenderedStimulus<R>> {
|
||||
if (self.scale, self.t_sec, self.feature_size, self.e.dim()) == (scale, t_sec, feature_size, dim) {
|
||||
Cow::Borrowed(self)
|
||||
} else {
|
||||
Cow::Owned(render_stim(self, scale, t_sec, feature_size, dim))
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> Stimulus<R> for Fields<R> {
|
||||
fn at(&self, _t_sec: R, _feat_size: R, _loc: Index) -> Fields<R> {
|
||||
*self
|
||||
}
|
||||
}
|
||||
|
||||
pub struct ModulatedVectorField<V, T> {
|
||||
fields: V,
|
||||
modulation: T,
|
||||
}
|
||||
|
||||
impl<V, T> ModulatedVectorField<V, T> {
|
||||
pub fn new(fields: V, modulation: T) -> Self {
|
||||
Self { fields, modulation }
|
||||
}
|
||||
pub fn into_inner(self) -> (V, T) {
|
||||
(self.fields, self.modulation)
|
||||
}
|
||||
pub fn fields(&self) -> &V {
|
||||
&self.fields
|
||||
}
|
||||
pub fn modulation(&self) -> &T {
|
||||
&self.modulation
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real, V: VectorField<R> + Sync, T: TimeVarying<R> + Sync> Stimulus<R> for ModulatedVectorField<V, T> {
|
||||
fn at(&self, t_sec: R, feat_size: R, loc: Index) -> Fields<R> {
|
||||
self.fields.at(feat_size, loc).elem_mul(self.modulation.at(t_sec))
|
||||
}
|
||||
}
|
||||
|
||||
/// used as a MapVisitor in order to evaluate each Stimulus in a List at a specific time/place.
|
||||
// struct StimulusEvaluator {
|
||||
// fields: Fields,
|
||||
// t_sec: f32,
|
||||
// feat_size: f32,
|
||||
// loc: Index,
|
||||
// }
|
||||
//
|
||||
// impl<S: Stimulus> Visitor<&S> for &mut StimulusEvaluator {
|
||||
// fn visit(&mut self, next: &S) {
|
||||
// self.fields += next.at(self.t_sec, self.feat_size, self.loc);
|
||||
// }
|
||||
// }
|
||||
//
|
||||
// impl<L: Sync> Stimulus for L
|
||||
// where
|
||||
// for<'a, 'b> &'a L: Visit<&'b mut StimulusEvaluator>,
|
||||
// {
|
||||
// fn at(&self, t_sec: f32, pos: Meters) -> Fields {
|
||||
// let mut ev = StimulusEvaluator { t_sec, pos, fields: Fields::default()};
|
||||
// self.visit(&mut ev);
|
||||
// ev.fields
|
||||
// }
|
||||
// }
|
||||
|
||||
// conflicts with List implementation
|
||||
// impl<T: Stimulus> Stimulus for &T {
|
||||
// fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
|
||||
// (*self).at(t_sec, feat_size, loc)
|
||||
// }
|
||||
// }
|
||||
|
||||
|
||||
pub struct StimuliVec<S>(Vec<S>);
|
||||
|
||||
pub type DynStimuli<R> = StimuliVec<Box<dyn Stimulus<R> + Send>>;
|
||||
|
||||
impl<S> Default for StimuliVec<S> {
|
||||
fn default() -> Self {
|
||||
Self(Vec::new())
|
||||
}
|
||||
}
|
||||
|
||||
impl<S> Deref for StimuliVec<S> {
|
||||
type Target = Vec<S>;
|
||||
fn deref(&self) -> &Self::Target {
|
||||
&self.0
|
||||
}
|
||||
}
|
||||
|
||||
impl<S> StimuliVec<S> {
|
||||
pub fn new() -> Self {
|
||||
Self::default()
|
||||
}
|
||||
pub fn from_vec(stim: Vec<S>) -> Self {
|
||||
Self(stim)
|
||||
}
|
||||
pub fn push(&mut self, a: S) {
|
||||
self.0.push(a)
|
||||
}
|
||||
}
|
||||
impl<R: Real, S: Stimulus<R>> Stimulus<R> for StimuliVec<S> {
|
||||
fn at(&self, t_sec: R, feat_size: R, loc: Index) -> Fields<R> {
|
||||
self.0.iter().map(|i| i.at(t_sec, feat_size, loc))
|
||||
.fold(Fields::default(), core::ops::Add::add)
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> Stimulus<R> for Box<dyn Stimulus<R> + Send> {
|
||||
fn at(&self, t_sec: R, feat_size: R, loc: Index) -> Fields<R> {
|
||||
(**self).at(t_sec, feat_size, loc)
|
||||
}
|
||||
}
|
||||
|
||||
pub struct NoopStimulus;
|
||||
impl<R: Real> Stimulus<R> for NoopStimulus {
|
||||
fn at(&self, _t_sec: R, _feat_size: R, _loc: Index) -> Fields<R> {
|
||||
Fields::default()
|
||||
}
|
||||
}
|
||||
|
||||
pub struct RngStimulus {
|
||||
seed: u64,
|
||||
e_scale: f32,
|
||||
h_scale: f32,
|
||||
}
|
||||
|
||||
impl RngStimulus {
|
||||
pub fn new(seed: u64) -> Self {
|
||||
Self { seed, e_scale: 1e15, h_scale: 1e15 }
|
||||
}
|
||||
pub fn new_e(seed: u64) -> Self {
|
||||
Self { seed, e_scale: 1e15, h_scale: 0.0 }
|
||||
}
|
||||
fn gen(&self, t_sec: f32, pos: Meters, scale: f32, salt: u64) -> Vec3<f32> {
|
||||
use rand::{Rng as _, SeedableRng as _};
|
||||
let seed = self.seed
|
||||
^ (t_sec.to_bits() as u64)
|
||||
^ ((pos.x().to_bits() as u64) << 8)
|
||||
^ ((pos.y().to_bits() as u64) << 16)
|
||||
^ ((pos.z().to_bits() as u64) << 24)
|
||||
^ (salt << 32);
|
||||
let mut rng = rand::rngs::StdRng::seed_from_u64(seed);
|
||||
Vec3::new(
|
||||
rng.gen_range(-scale..=scale),
|
||||
rng.gen_range(-scale..=scale),
|
||||
rng.gen_range(-scale..=scale),
|
||||
)
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> Stimulus<R> for RngStimulus {
|
||||
fn at(&self, t_sec: R, feat_size: R, loc: Index) -> Fields<R> {
|
||||
Fields {
|
||||
e: self.gen(t_sec.cast(), loc.to_meters(feat_size.cast()), self.e_scale, 0).cast(),
|
||||
h: self.gen(t_sec.cast(), loc.to_meters(feat_size.cast()), self.h_scale, 0x7de3).cast(),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Apply a time-varying stimulus uniformly across some region
|
||||
#[derive(Clone)]
|
||||
pub struct RegionGated<G, S> {
|
||||
region: G,
|
||||
stim: S,
|
||||
}
|
||||
|
||||
impl<G, S> RegionGated<G, S> {
|
||||
pub fn new(region: G, stim: S) -> Self {
|
||||
Self {
|
||||
region, stim
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real, G: Region + Sync, S: Stimulus<R> + Sync> Stimulus<R> for RegionGated<G, S> {
|
||||
fn at(&self, t_sec: R, feat_size: R, loc: Index) -> Fields<R> {
|
||||
if self.region.contains(loc.to_meters(feat_size.cast())) {
|
||||
self.stim.at(t_sec, feat_size, loc)
|
||||
} else {
|
||||
Fields::default()
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real, G: Region + Sync, S: VectorField<R>> VectorField<R> for RegionGated<G, S> {
|
||||
fn at(&self, feat_size: R, loc: Index) -> Fields<R> {
|
||||
if self.region.contains(loc.to_meters(feat_size.cast())) {
|
||||
self.stim.at(feat_size, loc)
|
||||
} else {
|
||||
Fields::default()
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// apply a stimulus across some region.
|
||||
/// the stimulus seen at each point is based on its angle about the specified ray.
|
||||
/// the stimulus has equal E and H vectors. if you want just one, filter it one with `Scaled`.
|
||||
#[derive(Clone)]
|
||||
pub struct CurlVectorField<G> {
|
||||
region: G,
|
||||
}
|
||||
|
||||
impl<G> CurlVectorField<G> {
|
||||
pub fn new(region: G) -> Self {
|
||||
Self { region }
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real, G: Region + HasCrossSection> VectorField<R> for CurlVectorField<G> {
|
||||
fn at(&self, feat_size: R, loc: Index) -> Fields<R> {
|
||||
let pos = loc.to_meters(feat_size.cast());
|
||||
if self.region.contains(pos) {
|
||||
// TODO: do we *want* this to be normalized?
|
||||
let rotational = self.region.cross_section_normal(pos).norm().cast();
|
||||
Fields::new_eh(rotational, rotational)
|
||||
} else {
|
||||
Fields::default()
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
pub trait StimExt<R>: Sized {
|
||||
fn shifted(self, new_start: R) -> Shifted<R, Self> {
|
||||
Shifted::new(self, new_start)
|
||||
}
|
||||
fn gated(self, from: R, to: R) -> Gated<R, Self> {
|
||||
Gated::new(Pulse::new(from, to), self)
|
||||
}
|
||||
fn scaled<T: TimeVarying<R>>(self, scale: T) -> Scaled<Self, T> {
|
||||
Scaled::new(self, scale)
|
||||
}
|
||||
fn summed<T: TimeVarying<R>>(self, with: T) -> Summed<Self, T> {
|
||||
Summed::new(self, with)
|
||||
}
|
||||
}
|
||||
|
||||
impl<R, T> StimExt<R> for T {}
|
||||
|
||||
pub trait TimeVarying<R> {
|
||||
fn at(&self, t_sec: R) -> FieldMags<R>;
|
||||
}
|
||||
|
||||
impl<R: Real> TimeVarying<R> for FieldMags<R> {
|
||||
fn at(&self, _t_sec: R) -> FieldMags<R> {
|
||||
*self
|
||||
}
|
||||
}
|
||||
|
||||
// assumed to represent the E field
|
||||
impl<R: Real> TimeVarying<real::Finite<R>> for real::Finite<R> {
|
||||
fn at(&self, _t_sec: real::Finite<R>) -> FieldMags<real::Finite<R>> {
|
||||
FieldMags::new_e(*self)
|
||||
}
|
||||
}
|
||||
impl TimeVarying<f32> for f32 {
|
||||
fn at(&self, _t_sec: f32) -> FieldMags<f32> {
|
||||
FieldMags::new_e(*self)
|
||||
}
|
||||
}
|
||||
impl TimeVarying<f64> for f64 {
|
||||
fn at(&self, _t_sec: f64) -> FieldMags<f64> {
|
||||
FieldMags::new_e(*self)
|
||||
}
|
||||
}
|
||||
|
||||
// Vec<T> at any `t_sec` behaves as the sum of all its components at that time.
|
||||
impl<R: Real, T: TimeVarying<R>> TimeVarying<R> for Vec<T> {
|
||||
fn at(&self, t_sec: R) -> FieldMags<R> {
|
||||
self.iter().fold(FieldMags::default(), |acc, i| acc + i.at(t_sec))
|
||||
}
|
||||
}
|
||||
|
||||
pub struct UnitEH;
|
||||
impl<R: Real> TimeVarying<R> for UnitEH {
|
||||
fn at(&self, _t_sec: R) -> FieldMags<R> {
|
||||
FieldMags::new_eh(R::one(), R::one())
|
||||
}
|
||||
}
|
||||
|
||||
/// E field which changes magnitude sinusoidally as a function of t
|
||||
#[derive(Clone)]
|
||||
pub struct Sinusoid<R> {
|
||||
amp: R,
|
||||
omega: R,
|
||||
}
|
||||
|
||||
impl<R: Real> Sinusoid<R> {
|
||||
pub fn new(amp: R, freq: R) -> Self {
|
||||
Self {
|
||||
amp,
|
||||
omega: freq * R::two_pi(),
|
||||
}
|
||||
}
|
||||
pub fn from_wavelength(amp: R, lambda: R) -> Self {
|
||||
Self::new(amp, lambda.inv())
|
||||
}
|
||||
pub fn freq(&self) -> R {
|
||||
self.omega / R::two_pi()
|
||||
}
|
||||
pub fn wavelength(&self) -> R {
|
||||
self.freq().inv()
|
||||
}
|
||||
pub fn one_cycle(self) -> Gated<R, Self> {
|
||||
let wl = self.wavelength();
|
||||
self.gated(R::zero(), wl)
|
||||
}
|
||||
pub fn half_cycle(self) -> Gated<R, Self> {
|
||||
let wl = self.wavelength();
|
||||
self.gated(R::zero(), R::half() * wl)
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> TimeVarying<R> for Sinusoid<R> {
|
||||
fn at(&self, t_sec: R) -> FieldMags<R> {
|
||||
FieldMags {
|
||||
e: self.amp * (t_sec * self.omega).sin(),
|
||||
h: R::zero(),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// E field that decays exponentially over t.
|
||||
#[derive(Clone)]
|
||||
pub struct Exp<R> {
|
||||
tau: R,
|
||||
}
|
||||
|
||||
impl<R: Real + TimeVarying<R>> Exp<R> {
|
||||
pub fn new(half_life: R) -> Self {
|
||||
let tau = R::ln2()/half_life;
|
||||
Self { tau }
|
||||
}
|
||||
|
||||
pub fn new_at(amp: R, start: R, half_life: R) -> Shifted<R, Gated<R, Scaled<Self, R>>> {
|
||||
Self::new(half_life).scaled(amp).gated(R::zero(), half_life*100f32.cast::<R>()).shifted(start)
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> TimeVarying<R> for Exp<R> {
|
||||
fn at(&self, t_sec: R) -> FieldMags<R> {
|
||||
let a = (t_sec * -self.tau).exp();
|
||||
FieldMags::new_eh(a, a)
|
||||
}
|
||||
}
|
||||
|
||||
/// pulses E=1.0 and H=1.0 over the provided duration.
|
||||
/// this is used as a building block to gate some VectorField over a specific time.
|
||||
#[derive(Clone)]
|
||||
pub struct Pulse<R> {
|
||||
start: R,
|
||||
end: R,
|
||||
}
|
||||
|
||||
impl<R> Pulse<R> {
|
||||
pub fn new(start: R, end: R) -> Self {
|
||||
Self { start, end }
|
||||
}
|
||||
}
|
||||
impl<R: Real> Pulse<R> {
|
||||
fn contains(&self, t: R) -> bool {
|
||||
t >= self.start && t < self.end
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> TimeVarying<R> for Pulse<R> {
|
||||
fn at(&self, t: R) -> FieldMags<R> {
|
||||
if self.contains(t) {
|
||||
FieldMags::new_eh(R::one(), R::one())
|
||||
} else {
|
||||
FieldMags::new_eh(R::zero(), R::zero())
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
pub type Gated<R, T> = Scaled<Pulse<R>, T>;
|
||||
|
||||
#[derive(Clone)]
|
||||
pub struct Shifted<R, T> {
|
||||
start_at: R,
|
||||
inner: T,
|
||||
}
|
||||
|
||||
impl<R, T> Shifted<R, T> {
|
||||
pub fn new(inner: T, start_at: R) -> Self {
|
||||
Self { inner, start_at }
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real, T: TimeVarying<R>> TimeVarying<R> for Shifted<R, T> {
|
||||
fn at(&self, t_sec: R) -> FieldMags<R> {
|
||||
self.inner.at(t_sec - self.start_at)
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real, T: Stimulus<R>> Stimulus<R> for Shifted<R, T> {
|
||||
fn at(&self, t_sec: R, feat_size: R, loc: Index) -> Fields<R> {
|
||||
self.inner.at(t_sec - self.start_at, feat_size, loc)
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#[derive(Clone)]
|
||||
pub struct Scaled<A, B>(A, B);
|
||||
|
||||
impl<A, B> Scaled<A, B> {
|
||||
pub fn new(a: A, b: B) -> Self {
|
||||
Self(a, b)
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real, A: TimeVarying<R>, B: TimeVarying<R>> TimeVarying<R> for Scaled<A, B> {
|
||||
fn at(&self, t_sec: R) -> FieldMags<R> {
|
||||
self.0.at(t_sec).elem_mul(self.1.at(t_sec))
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Clone)]
|
||||
pub struct Summed<A, B>(A, B);
|
||||
|
||||
impl<A, B> Summed<A, B> {
|
||||
pub fn new(a: A, b: B) -> Self {
|
||||
Self(a, b)
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real, A: TimeVarying<R>, B: TimeVarying<R>> TimeVarying<R> for Summed<A, B> {
|
||||
fn at(&self, t_sec: R) -> FieldMags<R> {
|
||||
self.0.at(t_sec) + self.1.at(t_sec)
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod test {
|
||||
use super::*;
|
||||
|
||||
macro_rules! assert_approx_eq {
|
||||
($x:expr, $e:expr, $h:expr) => {
|
||||
let x = $x;
|
||||
let e = $e;
|
||||
let h = $h;
|
||||
let diff_e = (x.e - e).abs();
|
||||
assert!(diff_e <= 0.001, "{:?} != {:?}", x, e);
|
||||
let diff_h = (x.h - h).abs();
|
||||
assert!(diff_h <= 0.001, "{:?} != {:?}", x, h);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn sinusoid() {
|
||||
let s = Sinusoid::new(10.0, 1000.0);
|
||||
assert_eq!(s.at(0.0), FieldMags::default());
|
||||
assert_approx_eq!(s.at(0.00025), 10.0, 0.0);
|
||||
assert_approx_eq!(s.at(0.00050), 0.0, 0.0);
|
||||
assert_approx_eq!(s.at(0.00075), -10.0, 0.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn sinusoid3_from_wavelength() {
|
||||
let s = Sinusoid::from_wavelength(10.0, 0.001);
|
||||
assert_eq!(s.at(0.0), FieldMags::default());
|
||||
assert_approx_eq!(s.at(0.00025), 10.0, 0.0);
|
||||
assert_approx_eq!(s.at(0.00050), 0.0, 0.0);
|
||||
assert_approx_eq!(s.at(0.00075), -10.0, 0.0);
|
||||
}
|
||||
|
||||
struct MockRegion {
|
||||
normal: Vec3<f32>,
|
||||
}
|
||||
impl HasCrossSection for MockRegion {
|
||||
fn cross_section_normal(&self, _p: Meters) -> Vec3<f32> {
|
||||
self.normal
|
||||
}
|
||||
}
|
||||
impl Region for MockRegion {
|
||||
fn contains(&self, _p: Meters) -> bool {
|
||||
true
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn curl_stimulus_trivial() {
|
||||
let region = MockRegion {
|
||||
normal: Vec3::new(1.0, 0.0, 0.0)
|
||||
};
|
||||
let stim = CurlVectorField::new(region);
|
||||
assert_eq!(stim.at(1.0, Index::new(0, 0, 0)), Fields {
|
||||
e: Vec3::new(1.0, 0.0, 0.0),
|
||||
h: Vec3::new(1.0, 0.0, 0.0),
|
||||
});
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn curl_stimulus_multi_axis() {
|
||||
let region = MockRegion {
|
||||
normal: Vec3::new(0.0, -1.0, 1.0)
|
||||
};
|
||||
let stim = CurlVectorField::new(region);
|
||||
let Fields { e, h } = stim.at(1.0, Index::new(0, 0, 0));
|
||||
assert_eq!(e, h);
|
||||
assert!(e.distance(Vec3::new(0.0, -1.0, 1.0).norm()) < 1e-6);
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,123 @@
|
|||
//! supporting types for basic Stimulus trait/impls
|
||||
use crate::real::Real;
|
||||
use crate::cross::vec::Vec3;
|
||||
|
||||
|
||||
/// field densities
|
||||
#[derive(Clone, Copy, Debug, Default, PartialEq)]
|
||||
pub struct Fields<R> {
|
||||
pub e: Vec3<R>,
|
||||
pub h: Vec3<R>,
|
||||
}
|
||||
|
||||
impl<R> Fields<R> {
|
||||
pub fn new_eh(e: Vec3<R>, h: Vec3<R>) -> Self {
|
||||
Self { e, h}
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> Fields<R> {
|
||||
pub fn e(&self) -> Vec3<R> {
|
||||
self.e
|
||||
}
|
||||
pub fn h(&self) -> Vec3<R> {
|
||||
self.h
|
||||
}
|
||||
pub fn new_e(e: Vec3<R>) -> Self {
|
||||
Self::new_eh(e, Vec3::zero())
|
||||
}
|
||||
pub fn new_h(h: Vec3<R>) -> Self {
|
||||
Self::new_eh(Vec3::zero(), h)
|
||||
}
|
||||
pub fn elem_mul(self, other: FieldMags<R>) -> Fields<R> {
|
||||
Fields {
|
||||
e: self.e * other.e,
|
||||
h: self.h * other.h,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> std::ops::AddAssign for Fields<R> {
|
||||
fn add_assign(&mut self, other: Self) {
|
||||
self.e += other.e;
|
||||
self.h += other.h;
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> std::ops::Add for Fields<R> {
|
||||
type Output = Self;
|
||||
fn add(mut self, other: Self) -> Self::Output {
|
||||
self += other;
|
||||
self
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> std::ops::Mul<R> for Fields<R> {
|
||||
type Output = Self;
|
||||
fn mul(self, scale: R) -> Self::Output {
|
||||
Fields {
|
||||
e: self.e * scale,
|
||||
h: self.h * scale,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// field magnitude densities (really, signed magnitude)
|
||||
#[derive(Clone, Copy, Debug, Default, PartialEq)]
|
||||
pub struct FieldMags<R> {
|
||||
pub e: R,
|
||||
pub h: R,
|
||||
}
|
||||
|
||||
impl<R: Real> std::ops::AddAssign for FieldMags<R> {
|
||||
fn add_assign(&mut self, other: Self) {
|
||||
self.e += other.e;
|
||||
self.h += other.h;
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> std::ops::Add for FieldMags<R> {
|
||||
type Output = Self;
|
||||
fn add(mut self, other: Self) -> Self::Output {
|
||||
self += other;
|
||||
self
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> std::ops::Mul<R> for FieldMags<R> {
|
||||
type Output = Self;
|
||||
fn mul(self, scale: R) -> Self::Output {
|
||||
FieldMags {
|
||||
e: self.e * scale,
|
||||
h: self.h * scale,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl<R> FieldMags<R> {
|
||||
pub fn new_eh(e: R, h: R) -> Self {
|
||||
Self { e, h }
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> FieldMags<R> {
|
||||
pub fn e(&self) -> R {
|
||||
self.e
|
||||
}
|
||||
pub fn h(&self) -> R {
|
||||
self.h
|
||||
}
|
||||
pub fn new_e(e: R) -> Self {
|
||||
Self::new_eh(e, R::zero())
|
||||
}
|
||||
pub fn new_h(h: R) -> Self {
|
||||
Self::new_eh(R::zero(), h)
|
||||
}
|
||||
pub fn elem_mul(self, other: Self) -> Self {
|
||||
FieldMags {
|
||||
e: self.e * other.e,
|
||||
h: self.h * other.h,
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -0,0 +1,295 @@
|
|||
use crate::cross::vec::Vec3;
|
||||
use crate::geom::{Coord as _, Index, Meters};
|
||||
use crate::real::Real;
|
||||
use coremem_cross::dim::DimSlice;
|
||||
use coremem_cross::vec::Vec3u;
|
||||
|
||||
use std::borrow::Cow;
|
||||
use std::ops::Deref;
|
||||
|
||||
use rand;
|
||||
|
||||
mod fields;
|
||||
mod time_varying;
|
||||
mod vector_field;
|
||||
pub use fields::{Fields, FieldMags};
|
||||
pub use time_varying::{
|
||||
Exp,
|
||||
Gated,
|
||||
Pulse,
|
||||
Scaled,
|
||||
Shifted,
|
||||
Sinusoid,
|
||||
StimExt,
|
||||
Summed,
|
||||
TimeVarying,
|
||||
UnitEH,
|
||||
};
|
||||
pub use vector_field::{
|
||||
CurlVectorField,
|
||||
RegionGated,
|
||||
VectorField,
|
||||
};
|
||||
|
||||
|
||||
pub trait Stimulus<R: Real>: Sync {
|
||||
/// Return the (E, H) field which should be added PER-SECOND to the provided position/time.
|
||||
fn at(&self, t_sec: R, feat_size: R, loc: Index) -> Fields<R>;
|
||||
|
||||
/// compute the value of this stimulus across all the simulation space
|
||||
fn rendered<'a>(
|
||||
&'a self, scale: R, t_sec: R, feature_size: R, dim: Vec3u
|
||||
) -> Cow<'a, RenderedStimulus<R>> {
|
||||
Cow::Owned(render_stim(self, scale, t_sec, feature_size, dim))
|
||||
}
|
||||
}
|
||||
fn render_stim<R: Real, S: Stimulus<R> + ?Sized>(
|
||||
stim: &S, scale: R, t_sec: R, feature_size: R, dim: Vec3u
|
||||
) -> RenderedStimulus<R> {
|
||||
let dim_len = dim.product_sum_usize();
|
||||
let mut e = Vec::new();
|
||||
e.resize_with(dim_len, Default::default);
|
||||
let mut h = Vec::new();
|
||||
h.resize_with(dim_len, Default::default);
|
||||
|
||||
rayon::scope(|s| {
|
||||
let mut undispatched_e = &mut e[..];
|
||||
let mut undispatched_h = &mut h[..];
|
||||
for z in 0..dim.z() {
|
||||
for y in 0..dim.y() {
|
||||
let (this_e, this_h);
|
||||
(this_e, undispatched_e) = undispatched_e.split_at_mut(dim.x() as usize);
|
||||
(this_h, undispatched_h) = undispatched_h.split_at_mut(dim.x() as usize);
|
||||
s.spawn(move |_| {
|
||||
for (x, (out_e, out_h)) in this_e.iter_mut().zip(this_h.iter_mut()).enumerate() {
|
||||
|
||||
let Fields { e, h } = stim.at(t_sec, feature_size, Index::new(x as u32, y, z));
|
||||
*out_e = e * scale;
|
||||
*out_h = h * scale;
|
||||
}
|
||||
});
|
||||
}
|
||||
}
|
||||
});
|
||||
|
||||
let field_e = DimSlice::new(dim, e);
|
||||
let field_h = DimSlice::new(dim, h);
|
||||
RenderedStimulus::new(
|
||||
field_e, field_h, scale, feature_size, t_sec
|
||||
)
|
||||
}
|
||||
|
||||
|
||||
#[derive(Clone)]
|
||||
pub struct RenderedStimulus<R> {
|
||||
e: DimSlice<Vec<Vec3<R>>>,
|
||||
h: DimSlice<Vec<Vec3<R>>>,
|
||||
scale: R,
|
||||
feature_size: R,
|
||||
t_sec: R,
|
||||
}
|
||||
|
||||
impl<R> RenderedStimulus<R> {
|
||||
pub fn new(
|
||||
e: DimSlice<Vec<Vec3<R>>>,
|
||||
h: DimSlice<Vec<Vec3<R>>>,
|
||||
scale: R,
|
||||
feature_size: R,
|
||||
t_sec: R,
|
||||
) -> Self {
|
||||
Self { e, h, scale, feature_size, t_sec }
|
||||
}
|
||||
pub fn e<'a>(&'a self) -> DimSlice<&'a [Vec3<R>]> {
|
||||
self.e.as_ref()
|
||||
}
|
||||
pub fn h<'a>(&'a self) -> DimSlice<&'a [Vec3<R>]> {
|
||||
self.h.as_ref()
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> RenderedStimulus<R> {
|
||||
pub fn scale(&self) -> R {
|
||||
self.scale
|
||||
}
|
||||
pub fn feature_size(&self) -> R {
|
||||
self.feature_size
|
||||
}
|
||||
pub fn time(&self) -> R {
|
||||
self.t_sec
|
||||
}
|
||||
}
|
||||
|
||||
// TODO: is this necessary?
|
||||
impl<R: Real> VectorField<R> for RenderedStimulus<R> {
|
||||
fn at(&self, _feat_size: R, loc: Index) -> Fields<R> {
|
||||
Fields::new_eh(self.e[loc.into()], self.h[loc.into()])
|
||||
}
|
||||
}
|
||||
impl<R: Real> Stimulus<R> for RenderedStimulus<R> {
|
||||
fn at(&self, _t_sec: R, _feat_size: R, loc: Index) -> Fields<R> {
|
||||
Fields::new_eh(self.e[loc.into()], self.h[loc.into()])
|
||||
}
|
||||
fn rendered<'a>(
|
||||
&'a self, scale: R, t_sec: R, feature_size: R, dim: Vec3u
|
||||
) -> Cow<'a, RenderedStimulus<R>> {
|
||||
if (self.scale, self.t_sec, self.feature_size, self.e.dim()) == (scale, t_sec, feature_size, dim) {
|
||||
Cow::Borrowed(self)
|
||||
} else {
|
||||
Cow::Owned(render_stim(self, scale, t_sec, feature_size, dim))
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> Stimulus<R> for Fields<R> {
|
||||
fn at(&self, _t_sec: R, _feat_size: R, _loc: Index) -> Fields<R> {
|
||||
*self
|
||||
}
|
||||
}
|
||||
|
||||
/// a VectorField type whose amplitude is modulated by a TimeVarying component.
|
||||
/// users will almost always use this as their stimulus implementation
|
||||
pub struct ModulatedVectorField<V, T> {
|
||||
fields: V,
|
||||
modulation: T,
|
||||
}
|
||||
|
||||
impl<V, T> ModulatedVectorField<V, T> {
|
||||
pub fn new(fields: V, modulation: T) -> Self {
|
||||
Self { fields, modulation }
|
||||
}
|
||||
pub fn into_inner(self) -> (V, T) {
|
||||
(self.fields, self.modulation)
|
||||
}
|
||||
pub fn fields(&self) -> &V {
|
||||
&self.fields
|
||||
}
|
||||
pub fn modulation(&self) -> &T {
|
||||
&self.modulation
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real, V: VectorField<R> + Sync, T: TimeVarying<R> + Sync> Stimulus<R> for ModulatedVectorField<V, T> {
|
||||
fn at(&self, t_sec: R, feat_size: R, loc: Index) -> Fields<R> {
|
||||
self.fields.at(feat_size, loc).elem_mul(self.modulation.at(t_sec))
|
||||
}
|
||||
}
|
||||
|
||||
/// used as a MapVisitor in order to evaluate each Stimulus in a List at a specific time/place.
|
||||
// struct StimulusEvaluator {
|
||||
// fields: Fields,
|
||||
// t_sec: f32,
|
||||
// feat_size: f32,
|
||||
// loc: Index,
|
||||
// }
|
||||
//
|
||||
// impl<S: Stimulus> Visitor<&S> for &mut StimulusEvaluator {
|
||||
// fn visit(&mut self, next: &S) {
|
||||
// self.fields += next.at(self.t_sec, self.feat_size, self.loc);
|
||||
// }
|
||||
// }
|
||||
//
|
||||
// impl<L: Sync> Stimulus for L
|
||||
// where
|
||||
// for<'a, 'b> &'a L: Visit<&'b mut StimulusEvaluator>,
|
||||
// {
|
||||
// fn at(&self, t_sec: f32, pos: Meters) -> Fields {
|
||||
// let mut ev = StimulusEvaluator { t_sec, pos, fields: Fields::default()};
|
||||
// self.visit(&mut ev);
|
||||
// ev.fields
|
||||
// }
|
||||
// }
|
||||
|
||||
// conflicts with List implementation
|
||||
// impl<T: Stimulus> Stimulus for &T {
|
||||
// fn at(&self, t_sec: f32, feat_size: f32, loc: Index) -> Fields {
|
||||
// (*self).at(t_sec, feat_size, loc)
|
||||
// }
|
||||
// }
|
||||
|
||||
|
||||
pub struct StimuliVec<S>(Vec<S>);
|
||||
|
||||
pub type DynStimuli<R> = StimuliVec<Box<dyn Stimulus<R> + Send>>;
|
||||
|
||||
impl<S> Default for StimuliVec<S> {
|
||||
fn default() -> Self {
|
||||
Self(Vec::new())
|
||||
}
|
||||
}
|
||||
|
||||
impl<S> Deref for StimuliVec<S> {
|
||||
type Target = Vec<S>;
|
||||
fn deref(&self) -> &Self::Target {
|
||||
&self.0
|
||||
}
|
||||
}
|
||||
|
||||
impl<S> StimuliVec<S> {
|
||||
pub fn new() -> Self {
|
||||
Self::default()
|
||||
}
|
||||
pub fn from_vec(stim: Vec<S>) -> Self {
|
||||
Self(stim)
|
||||
}
|
||||
pub fn push(&mut self, a: S) {
|
||||
self.0.push(a)
|
||||
}
|
||||
}
|
||||
impl<R: Real, S: Stimulus<R>> Stimulus<R> for StimuliVec<S> {
|
||||
fn at(&self, t_sec: R, feat_size: R, loc: Index) -> Fields<R> {
|
||||
self.0.iter().map(|i| i.at(t_sec, feat_size, loc))
|
||||
.fold(Fields::default(), core::ops::Add::add)
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> Stimulus<R> for Box<dyn Stimulus<R> + Send> {
|
||||
fn at(&self, t_sec: R, feat_size: R, loc: Index) -> Fields<R> {
|
||||
(**self).at(t_sec, feat_size, loc)
|
||||
}
|
||||
}
|
||||
|
||||
pub struct NoopStimulus;
|
||||
impl<R: Real> Stimulus<R> for NoopStimulus {
|
||||
fn at(&self, _t_sec: R, _feat_size: R, _loc: Index) -> Fields<R> {
|
||||
Fields::default()
|
||||
}
|
||||
}
|
||||
|
||||
pub struct RngStimulus {
|
||||
seed: u64,
|
||||
e_scale: f32,
|
||||
h_scale: f32,
|
||||
}
|
||||
|
||||
impl RngStimulus {
|
||||
pub fn new(seed: u64) -> Self {
|
||||
Self { seed, e_scale: 1e15, h_scale: 1e15 }
|
||||
}
|
||||
pub fn new_e(seed: u64) -> Self {
|
||||
Self { seed, e_scale: 1e15, h_scale: 0.0 }
|
||||
}
|
||||
fn gen(&self, t_sec: f32, pos: Meters, scale: f32, salt: u64) -> Vec3<f32> {
|
||||
use rand::{Rng as _, SeedableRng as _};
|
||||
let seed = self.seed
|
||||
^ (t_sec.to_bits() as u64)
|
||||
^ ((pos.x().to_bits() as u64) << 8)
|
||||
^ ((pos.y().to_bits() as u64) << 16)
|
||||
^ ((pos.z().to_bits() as u64) << 24)
|
||||
^ (salt << 32);
|
||||
let mut rng = rand::rngs::StdRng::seed_from_u64(seed);
|
||||
Vec3::new(
|
||||
rng.gen_range(-scale..=scale),
|
||||
rng.gen_range(-scale..=scale),
|
||||
rng.gen_range(-scale..=scale),
|
||||
)
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> Stimulus<R> for RngStimulus {
|
||||
fn at(&self, t_sec: R, feat_size: R, loc: Index) -> Fields<R> {
|
||||
Fields {
|
||||
e: self.gen(t_sec.cast(), loc.to_meters(feat_size.cast()), self.e_scale, 0).cast(),
|
||||
h: self.gen(t_sec.cast(), loc.to_meters(feat_size.cast()), self.h_scale, 0x7de3).cast(),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,254 @@
|
|||
//! time-varying portions of a Stimulus
|
||||
|
||||
use crate::real::{self, Real};
|
||||
use crate::stim::{Fields, FieldMags, Stimulus};
|
||||
use crate::geom::Index;
|
||||
|
||||
pub trait TimeVarying<R> {
|
||||
fn at(&self, t_sec: R) -> FieldMags<R>;
|
||||
}
|
||||
|
||||
// TODO: rename to TimeVaryingExt?
|
||||
pub trait StimExt<R>: Sized {
|
||||
fn shifted(self, new_start: R) -> Shifted<R, Self> {
|
||||
Shifted::new(self, new_start)
|
||||
}
|
||||
fn gated(self, from: R, to: R) -> Gated<R, Self> {
|
||||
Gated::new(Pulse::new(from, to), self)
|
||||
}
|
||||
fn scaled<T: TimeVarying<R>>(self, scale: T) -> Scaled<Self, T> {
|
||||
Scaled::new(self, scale)
|
||||
}
|
||||
fn summed<T: TimeVarying<R>>(self, with: T) -> Summed<Self, T> {
|
||||
Summed::new(self, with)
|
||||
}
|
||||
}
|
||||
|
||||
impl<R, T> StimExt<R> for T {}
|
||||
|
||||
impl<R: Real> TimeVarying<R> for FieldMags<R> {
|
||||
fn at(&self, _t_sec: R) -> FieldMags<R> {
|
||||
*self
|
||||
}
|
||||
}
|
||||
|
||||
// assumed to represent the E field
|
||||
impl<R: Real> TimeVarying<real::Finite<R>> for real::Finite<R> {
|
||||
fn at(&self, _t_sec: real::Finite<R>) -> FieldMags<real::Finite<R>> {
|
||||
FieldMags::new_e(*self)
|
||||
}
|
||||
}
|
||||
impl TimeVarying<f32> for f32 {
|
||||
fn at(&self, _t_sec: f32) -> FieldMags<f32> {
|
||||
FieldMags::new_e(*self)
|
||||
}
|
||||
}
|
||||
impl TimeVarying<f64> for f64 {
|
||||
fn at(&self, _t_sec: f64) -> FieldMags<f64> {
|
||||
FieldMags::new_e(*self)
|
||||
}
|
||||
}
|
||||
|
||||
// Vec<T> at any `t_sec` behaves as the sum of all its components at that time.
|
||||
impl<R: Real, T: TimeVarying<R>> TimeVarying<R> for Vec<T> {
|
||||
fn at(&self, t_sec: R) -> FieldMags<R> {
|
||||
self.iter().fold(FieldMags::default(), |acc, i| acc + i.at(t_sec))
|
||||
}
|
||||
}
|
||||
|
||||
pub struct UnitEH;
|
||||
impl<R: Real> TimeVarying<R> for UnitEH {
|
||||
fn at(&self, _t_sec: R) -> FieldMags<R> {
|
||||
FieldMags::new_eh(R::one(), R::one())
|
||||
}
|
||||
}
|
||||
|
||||
/// E field which changes magnitude sinusoidally as a function of t
|
||||
#[derive(Clone)]
|
||||
pub struct Sinusoid<R> {
|
||||
amp: R,
|
||||
omega: R,
|
||||
}
|
||||
|
||||
impl<R: Real> Sinusoid<R> {
|
||||
pub fn new(amp: R, freq: R) -> Self {
|
||||
Self {
|
||||
amp,
|
||||
omega: freq * R::two_pi(),
|
||||
}
|
||||
}
|
||||
pub fn from_wavelength(amp: R, lambda: R) -> Self {
|
||||
Self::new(amp, lambda.inv())
|
||||
}
|
||||
pub fn freq(&self) -> R {
|
||||
self.omega / R::two_pi()
|
||||
}
|
||||
pub fn wavelength(&self) -> R {
|
||||
self.freq().inv()
|
||||
}
|
||||
pub fn one_cycle(self) -> Gated<R, Self> {
|
||||
let wl = self.wavelength();
|
||||
self.gated(R::zero(), wl)
|
||||
}
|
||||
pub fn half_cycle(self) -> Gated<R, Self> {
|
||||
let wl = self.wavelength();
|
||||
self.gated(R::zero(), R::half() * wl)
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> TimeVarying<R> for Sinusoid<R> {
|
||||
fn at(&self, t_sec: R) -> FieldMags<R> {
|
||||
FieldMags {
|
||||
e: self.amp * (t_sec * self.omega).sin(),
|
||||
h: R::zero(),
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// E field that decays exponentially over t.
|
||||
#[derive(Clone)]
|
||||
pub struct Exp<R> {
|
||||
tau: R,
|
||||
}
|
||||
|
||||
impl<R: Real + TimeVarying<R>> Exp<R> {
|
||||
pub fn new(half_life: R) -> Self {
|
||||
let tau = R::ln2()/half_life;
|
||||
Self { tau }
|
||||
}
|
||||
|
||||
pub fn new_at(amp: R, start: R, half_life: R) -> Shifted<R, Gated<R, Scaled<Self, R>>> {
|
||||
Self::new(half_life).scaled(amp).gated(R::zero(), half_life*100f32.cast::<R>()).shifted(start)
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> TimeVarying<R> for Exp<R> {
|
||||
fn at(&self, t_sec: R) -> FieldMags<R> {
|
||||
let a = (t_sec * -self.tau).exp();
|
||||
FieldMags::new_eh(a, a)
|
||||
}
|
||||
}
|
||||
|
||||
/// pulses E=1.0 and H=1.0 over the provided duration.
|
||||
/// this is used as a building block to gate some VectorField over a specific time.
|
||||
#[derive(Clone)]
|
||||
pub struct Pulse<R> {
|
||||
start: R,
|
||||
end: R,
|
||||
}
|
||||
|
||||
impl<R> Pulse<R> {
|
||||
pub fn new(start: R, end: R) -> Self {
|
||||
Self { start, end }
|
||||
}
|
||||
}
|
||||
impl<R: Real> Pulse<R> {
|
||||
fn contains(&self, t: R) -> bool {
|
||||
t >= self.start && t < self.end
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real> TimeVarying<R> for Pulse<R> {
|
||||
fn at(&self, t: R) -> FieldMags<R> {
|
||||
if self.contains(t) {
|
||||
FieldMags::new_eh(R::one(), R::one())
|
||||
} else {
|
||||
FieldMags::new_eh(R::zero(), R::zero())
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
pub type Gated<R, T> = Scaled<Pulse<R>, T>;
|
||||
|
||||
#[derive(Clone)]
|
||||
pub struct Shifted<R, T> {
|
||||
start_at: R,
|
||||
inner: T,
|
||||
}
|
||||
|
||||
impl<R, T> Shifted<R, T> {
|
||||
pub fn new(inner: T, start_at: R) -> Self {
|
||||
Self { inner, start_at }
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real, T: TimeVarying<R>> TimeVarying<R> for Shifted<R, T> {
|
||||
fn at(&self, t_sec: R) -> FieldMags<R> {
|
||||
self.inner.at(t_sec - self.start_at)
|
||||
}
|
||||
}
|
||||
|
||||
// TODO: is this necessary?
|
||||
impl<R: Real, T: Stimulus<R>> Stimulus<R> for Shifted<R, T> {
|
||||
fn at(&self, t_sec: R, feat_size: R, loc: Index) -> Fields<R> {
|
||||
self.inner.at(t_sec - self.start_at, feat_size, loc)
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#[derive(Clone)]
|
||||
pub struct Scaled<A, B>(A, B);
|
||||
|
||||
impl<A, B> Scaled<A, B> {
|
||||
pub fn new(a: A, b: B) -> Self {
|
||||
Self(a, b)
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real, A: TimeVarying<R>, B: TimeVarying<R>> TimeVarying<R> for Scaled<A, B> {
|
||||
fn at(&self, t_sec: R) -> FieldMags<R> {
|
||||
self.0.at(t_sec).elem_mul(self.1.at(t_sec))
|
||||
}
|
||||
}
|
||||
|
||||
#[derive(Clone)]
|
||||
pub struct Summed<A, B>(A, B);
|
||||
|
||||
impl<A, B> Summed<A, B> {
|
||||
pub fn new(a: A, b: B) -> Self {
|
||||
Self(a, b)
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real, A: TimeVarying<R>, B: TimeVarying<R>> TimeVarying<R> for Summed<A, B> {
|
||||
fn at(&self, t_sec: R) -> FieldMags<R> {
|
||||
self.0.at(t_sec) + self.1.at(t_sec)
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
#[cfg(test)]
|
||||
mod test {
|
||||
use super::*;
|
||||
|
||||
macro_rules! assert_approx_eq {
|
||||
($x:expr, $e:expr, $h:expr) => {
|
||||
let x = $x;
|
||||
let e = $e;
|
||||
let h = $h;
|
||||
let diff_e = (x.e - e).abs();
|
||||
assert!(diff_e <= 0.001, "{:?} != {:?}", x, e);
|
||||
let diff_h = (x.h - h).abs();
|
||||
assert!(diff_h <= 0.001, "{:?} != {:?}", x, h);
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn sinusoid() {
|
||||
let s = Sinusoid::new(10.0, 1000.0);
|
||||
assert_eq!(s.at(0.0), FieldMags::default());
|
||||
assert_approx_eq!(s.at(0.00025), 10.0, 0.0);
|
||||
assert_approx_eq!(s.at(0.00050), 0.0, 0.0);
|
||||
assert_approx_eq!(s.at(0.00075), -10.0, 0.0);
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn sinusoid_from_wavelength() {
|
||||
let s = Sinusoid::from_wavelength(10.0, 0.001);
|
||||
assert_eq!(s.at(0.0), FieldMags::default());
|
||||
assert_approx_eq!(s.at(0.00025), 10.0, 0.0);
|
||||
assert_approx_eq!(s.at(0.00050), 0.0, 0.0);
|
||||
assert_approx_eq!(s.at(0.00075), -10.0, 0.0);
|
||||
}
|
||||
}
|
||||
|
|
@ -0,0 +1,137 @@
|
|||
use crate::geom::{Coord as _, HasCrossSection, Index, Region};
|
||||
use crate::real::Real;
|
||||
use crate::stim::{Fields, Stimulus};
|
||||
use coremem_cross::dim::DimSlice;
|
||||
use coremem_cross::vec::Vec3u;
|
||||
|
||||
/// a static vector field. different value at each location, but constant in time.
|
||||
/// often used as a building block by wrapping it in something which modulates the fields over
|
||||
/// time.
|
||||
pub trait VectorField<R> {
|
||||
fn at(&self, feat_size: R, loc: Index) -> Fields<R>;
|
||||
}
|
||||
|
||||
// uniform vector field
|
||||
impl<R: Real> VectorField<R> for Fields<R> {
|
||||
fn at(&self, _feat_size: R, _loc: Index) -> Fields<R> {
|
||||
*self
|
||||
}
|
||||
}
|
||||
|
||||
// could broaden this and implement directly on T, but blanket impls
|
||||
// are unwieldy
|
||||
impl<R: Real, T> VectorField<R> for DimSlice<T>
|
||||
where
|
||||
DimSlice<T>: core::ops::Index<Vec3u, Output=Fields<R>>
|
||||
{
|
||||
fn at(&self, _feat_size: R, loc: Index) -> Fields<R> {
|
||||
self[loc.into()]
|
||||
}
|
||||
}
|
||||
|
||||
/// Apply a time-varying stimulus uniformly across some region
|
||||
#[derive(Clone)]
|
||||
pub struct RegionGated<G, S> {
|
||||
region: G,
|
||||
stim: S,
|
||||
}
|
||||
|
||||
impl<G, S> RegionGated<G, S> {
|
||||
pub fn new(region: G, stim: S) -> Self {
|
||||
Self {
|
||||
region, stim
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// TODO: is this necessary?
|
||||
impl<R: Real, G: Region + Sync, S: Stimulus<R> + Sync> Stimulus<R> for RegionGated<G, S> {
|
||||
fn at(&self, t_sec: R, feat_size: R, loc: Index) -> Fields<R> {
|
||||
if self.region.contains(loc.to_meters(feat_size.cast())) {
|
||||
self.stim.at(t_sec, feat_size, loc)
|
||||
} else {
|
||||
Fields::default()
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real, G: Region + Sync, S: VectorField<R>> VectorField<R> for RegionGated<G, S> {
|
||||
fn at(&self, feat_size: R, loc: Index) -> Fields<R> {
|
||||
if self.region.contains(loc.to_meters(feat_size.cast())) {
|
||||
self.stim.at(feat_size, loc)
|
||||
} else {
|
||||
Fields::default()
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// apply a stimulus across some region.
|
||||
/// the stimulus seen at each point is based on its angle about the specified ray.
|
||||
/// the stimulus has equal E and H vectors. if you want just one, filter it one with `Scaled`.
|
||||
#[derive(Clone)]
|
||||
pub struct CurlVectorField<G> {
|
||||
region: G,
|
||||
}
|
||||
|
||||
impl<G> CurlVectorField<G> {
|
||||
pub fn new(region: G) -> Self {
|
||||
Self { region }
|
||||
}
|
||||
}
|
||||
|
||||
impl<R: Real, G: Region + HasCrossSection> VectorField<R> for CurlVectorField<G> {
|
||||
fn at(&self, feat_size: R, loc: Index) -> Fields<R> {
|
||||
let pos = loc.to_meters(feat_size.cast());
|
||||
if self.region.contains(pos) {
|
||||
// TODO: do we *want* this to be normalized?
|
||||
let rotational = self.region.cross_section_normal(pos).norm().cast();
|
||||
Fields::new_eh(rotational, rotational)
|
||||
} else {
|
||||
Fields::default()
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#[cfg(test)]
|
||||
mod test {
|
||||
use super::*;
|
||||
use crate::cross::vec::Vec3;
|
||||
use crate::geom::Meters;
|
||||
|
||||
struct MockRegion {
|
||||
normal: Vec3<f32>,
|
||||
}
|
||||
impl HasCrossSection for MockRegion {
|
||||
fn cross_section_normal(&self, _p: Meters) -> Vec3<f32> {
|
||||
self.normal
|
||||
}
|
||||
}
|
||||
impl Region for MockRegion {
|
||||
fn contains(&self, _p: Meters) -> bool {
|
||||
true
|
||||
}
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn curl_stimulus_trivial() {
|
||||
let region = MockRegion {
|
||||
normal: Vec3::new(1.0, 0.0, 0.0)
|
||||
};
|
||||
let stim = CurlVectorField::new(region);
|
||||
assert_eq!(stim.at(1.0, Index::new(0, 0, 0)), Fields {
|
||||
e: Vec3::new(1.0, 0.0, 0.0),
|
||||
h: Vec3::new(1.0, 0.0, 0.0),
|
||||
});
|
||||
}
|
||||
|
||||
#[test]
|
||||
fn curl_stimulus_multi_axis() {
|
||||
let region = MockRegion {
|
||||
normal: Vec3::new(0.0, -1.0, 1.0)
|
||||
};
|
||||
let stim = CurlVectorField::new(region);
|
||||
let Fields { e, h } = stim.at(1.0, Index::new(0, 0, 0));
|
||||
assert_eq!(e, h);
|
||||
assert!(e.distance(Vec3::new(0.0, -1.0, 1.0).norm()) < 1e-6);
|
||||
}
|
||||
}
|
||||
Loading…
Reference in New Issue