Skip to content

Commit 29ecabf

Browse files
drobnyjtdrobnyjt
authored
Merge pull request #154 from lcpp-org/input_distributions
New feature: particle distributions
2 parents 093a857 + b434775 commit 29ecabf

File tree

5 files changed

+91
-30
lines changed

5 files changed

+91
-30
lines changed

‎Cargo.toml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -14,6 +14,7 @@ crate-type = ["cdylib"]
1414

1515
[dependencies]
1616
rand = "0.8.3"
17+
rand_distr = "0.4.2"
1718
toml = "0.5.8"
1819
anyhow = "1.0.38"
1920
itertools = "0.10.0"

‎src/enums.rs‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,14 @@ pub enum MaterialType {
1111
TRIMESH(material::Material<parry::ParryTriMesh>)
1212
}
1313

14+
#[derive(Deserialize, PartialEq, Clone, Copy)]
15+
#[serde(untagged)]
16+
pub enum Distributions {
17+
UNIFORM{min: f64, max: f64},
18+
NORMAL{mean: f64, std: f64},
19+
POINT(f64),
20+
}
21+
1422
#[derive(Deserialize)]
1523
pub enum GeometryType {
1624
MESH0D,

‎src/input.rs‎

Lines changed: 76 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
use super::*;
2+
use rand_distr::{Normal, Distribution, Uniform};
23

34
pub trait InputFile: GeometryInput {
45
fn new(string: &str) -> Self;
@@ -332,25 +333,53 @@ where <T as Geometry>::InputFileFormat: Deserialize<'static> + 'static {
332333
let Ec = particle_parameters.Ec[particle_index];
333334
let Es = particle_parameters.Es[particle_index];
334335
let interaction_index = particle_parameters.interaction_index[particle_index];
336+
335337
let (x, y, z) = particle_parameters.pos[particle_index];
336338
let (cosx, cosy, cosz) = particle_parameters.dir[particle_index];
337-
assert!(cosx < 1.,
338-
"Input error: particle x-direction cannot be exactly equal to 1 to avoid numerical gimbal lock.");
339+
339340
for sub_particle_index in 0..N_ {
340341
//Add new particle to particle vector
341342
particle_input.push(
342343
particle::ParticleInput{
343344
m: m*mass_unit,
344345
Z: Z,
345-
E: E*energy_unit,
346+
E: match E {
347+
Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*energy_unit},
348+
Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*energy_unit},
349+
Distributions::POINT(x) => x*energy_unit,
350+
},
346351
Ec: Ec*energy_unit,
347352
Es: Es*energy_unit,
348-
x: x*length_unit,
349-
y: y*length_unit,
350-
z: z*length_unit,
351-
ux: cosx,
352-
uy: cosy,
353-
uz: cosz,
353+
x: match x {
354+
Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
355+
Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit},
356+
Distributions::POINT(x) => x*length_unit,
357+
},
358+
y: match y {
359+
Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
360+
Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit},
361+
Distributions::POINT(x) => x*length_unit,
362+
},
363+
z: match z {
364+
Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
365+
Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit},
366+
Distributions::POINT(x) => x*length_unit,
367+
},
368+
ux: match cosx {
369+
Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
370+
Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit},
371+
Distributions::POINT(x) => {assert!(x != 1.0, "ux cannot equal exactly 1.0 to prevent gimbal lock"); x*length_unit}
372+
},
373+
uy: match cosy {
374+
Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
375+
Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit},
376+
Distributions::POINT(x) => x*length_unit,
377+
},
378+
uz: match cosz {
379+
Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
380+
Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit},
381+
Distributions::POINT(x) => x*length_unit,
382+
},
354383
interaction_index: interaction_index,
355384
tag: 0,
356385
weight: 1.0,
@@ -379,27 +408,54 @@ where <T as Geometry>::InputFileFormat: Deserialize<'static> + 'static {
379408
let interaction_index = particle_parameters.interaction_index[particle_index];
380409
let (x, y, z) = particle_parameters.pos[particle_index];
381410
let (cosx, cosy, cosz) = particle_parameters.dir[particle_index];
382-
assert!(cosx < 1.,
383-
"Input error: particle x-direction cannot be exactly equal to 1 to avoid numerical gimbal lock.");
411+
384412
for sub_particle_index in 0..N_ {
385413

386414
//Add new particle to particle vector
387415
particle_input.push(
388416
particle::ParticleInput{
389417
m: m*mass_unit,
390418
Z: Z,
391-
E: E*energy_unit,
419+
E: match E {
420+
Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*energy_unit},
421+
Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*energy_unit},
422+
Distributions::POINT(x) => x*energy_unit,
423+
},
392424
Ec: Ec*energy_unit,
393425
Es: Es*energy_unit,
394-
x: x*length_unit,
395-
y: y*length_unit,
396-
z: z*length_unit,
397-
ux: cosx,
398-
uy: cosy,
399-
uz: cosz,
400-
interaction_index,
401-
weight: 1.0,
426+
x: match x {
427+
Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
428+
Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit},
429+
Distributions::POINT(x) => x*length_unit,
430+
},
431+
y: match y {
432+
Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
433+
Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit},
434+
Distributions::POINT(x) => x*length_unit,
435+
},
436+
z: match z {
437+
Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
438+
Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit},
439+
Distributions::POINT(x) => x*length_unit,
440+
},
441+
ux: match cosx {
442+
Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
443+
Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit},
444+
Distributions::POINT(x) => {assert!(x != 1.0, "ux cannot equal exactly 1.0 to prevent gimbal lock"); x*length_unit},
445+
},
446+
uy: match cosy {
447+
Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
448+
Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit},
449+
Distributions::POINT(x) => x*length_unit,
450+
},
451+
uz: match cosz {
452+
Distributions::NORMAL{mean, std} => {let normal = Normal::new(mean, std).unwrap(); normal.sample(&mut rand::thread_rng())*length_unit},
453+
Distributions::UNIFORM{min, max} => {let uniform = Uniform::from(min..max); uniform.sample(&mut rand::thread_rng())*length_unit},
454+
Distributions::POINT(x) => x*length_unit,
455+
},
456+
interaction_index: interaction_index,
402457
tag: 0,
458+
weight: 1.0,
403459
}
404460
);
405461
}

‎src/main.rs‎

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -43,9 +43,6 @@ use std::f64::consts::FRAC_2_SQRT_PI;
4343
use std::f64::consts::PI;
4444
use std::f64::consts::SQRT_2;
4545

46-
//rng
47-
//use rand::{Rng, thread_rng};
48-
4946
//Load internal modules
5047
pub mod material;
5148
pub mod particle;

‎src/particle.rs‎

Lines changed: 6 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -2,7 +2,6 @@ use super::*;
22

33
/// Rustbca's internal representation of the particle_parameters input.
44
5-
65
#[cfg(feature = "hdf5_input")]
76
#[derive(Deserialize, Clone)]
87
pub struct ParticleParameters {
@@ -13,11 +12,11 @@ pub struct ParticleParameters {
1312
pub N: Vec<usize>,
1413
pub m: Vec<f64>,
1514
pub Z: Vec<f64>,
16-
pub E: Vec<f64>,
15+
pub E: Vec<Distributions>,
1716
pub Ec: Vec<f64>,
1817
pub Es: Vec<f64>,
19-
pub pos: Vec<(f64, f64, f64)>,
20-
pub dir: Vec<(f64, f64, f64)>,
18+
pub pos: Vec<(Distributions, Distributions, Distributions)>,
19+
pub dir: Vec<(Distributions, Distributions, Distributions)>,
2120
pub interaction_index: Vec<usize>,
2221
}
2322

@@ -30,11 +29,11 @@ pub struct ParticleParameters {
3029
pub N: Vec<usize>,
3130
pub m: Vec<f64>,
3231
pub Z: Vec<f64>,
33-
pub E: Vec<f64>,
32+
pub E: Vec<Distributions>,
3433
pub Ec: Vec<f64>,
3534
pub Es: Vec<f64>,
36-
pub pos: Vec<(f64, f64, f64)>,
37-
pub dir: Vec<(f64, f64, f64)>,
35+
pub pos: Vec<(Distributions, Distributions, Distributions)>,
36+
pub dir: Vec<(Distributions, Distributions, Distributions)>,
3837
pub interaction_index: Vec<usize>,
3938
}
4039

0 commit comments

Comments
 (0)