From 4ce8d5389e7413afe46fb37a2428e75f62254b23 Mon Sep 17 00:00:00 2001 From: BrianBostwick <42308966+BrianBostwick@users.noreply.github.com> Date: Thu, 25 Aug 2022 18:24:01 +0100 Subject: [PATCH 01/19] Fixing beam intensity gradient A factor of two missing from the intensity gradient vector caused incorrect dipole trap frequencies. --- src/laser/gaussian.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index b3c94cc3..4da1d11f 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -234,7 +234,7 @@ pub fn get_gaussian_beam_intensity_gradient( 2.0 * beam.e_radius.powf(2.0) * (1. + (z / beam.rayleigh_range).powf(2.0)); let vector = -4. * (reference_frame.x_vector * x + reference_frame.y_vector * y) + beam.direction * z / (beam.rayleigh_range.powf(2.0) + z.powf(2.0)) - * (-spot_size_squared + 4. * (x.powf(2.0) + y.powf(2.0))); + * (- 2.0 * spot_size_squared + 4. * (x.powf(2.0) + y.powf(2.0))); let intensity = 2. * beam.power / PI / spot_size_squared * EXP.powf(-2. * (x.powf(2.0) + y.powf(2.0)) / spot_size_squared); From 99dd38c7d8bd624a576004d93e0f692bfd61f133 Mon Sep 17 00:00:00 2001 From: Brian Date: Fri, 9 Sep 2022 14:58:52 +0100 Subject: [PATCH 02/19] changed to focused beam trap --- examples/dipole_trap.rs | 36 +++++++++--------------------------- 1 file changed, 9 insertions(+), 27 deletions(-) diff --git a/examples/dipole_trap.rs b/examples/dipole_trap.rs index 2d615906..4942fac2 100644 --- a/examples/dipole_trap.rs +++ b/examples/dipole_trap.rs @@ -14,7 +14,7 @@ use nalgebra::Vector3; use specs::prelude::*; use std::time::Instant; -const BEAM_NUMBER: usize = 2; +const BEAM_NUMBER: usize = 1; fn main() { let now = Instant::now(); @@ -23,9 +23,9 @@ fn main() { let mut sim_builder = SimulationBuilder::default(); sim_builder.add_plugin(LaserPlugin::<{BEAM_NUMBER}>); sim_builder.add_plugin(DipolePlugin::<{BEAM_NUMBER}>); - sim_builder.add_plugin(FileOutputPlugin::::new("pos.txt".to_string(), 100)); - sim_builder.add_plugin(FileOutputPlugin::::new("vel.txt".to_string(), 100)); - sim_builder.add_plugin(FileOutputPlugin::::new("position.xyz".to_string(), 100)); + sim_builder.add_plugin(FileOutputPlugin::::new("pos.txt".to_string(), 1)); + sim_builder.add_plugin(FileOutputPlugin::::new("vel.txt".to_string(), 1)); + sim_builder.add_plugin(FileOutputPlugin::::new("position.xyz".to_string(), 1)); let mut sim = sim_builder.build(); // Create dipole laser. @@ -51,34 +51,13 @@ fn main() { }) .build(); - let gaussian_beam = GaussianBeam { - intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius, - power, - direction: Vector3::y(), - rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range(&wavelength, &e_radius), - ellipticity: 0.0, - }; - sim.world - .create_entity() - .with(gaussian_beam) - .with(dipole::DipoleLight { wavelength }) - .with(laser::frame::Frame { - x_vector: Vector3::x(), - y_vector: Vector3::z(), - }) - .build(); - - // Define timestep - sim.world.insert(Timestep { delta: 1.0e-5 }); - // Create a single test atom sim.world .create_entity() .with(atom::Mass { value: 87.0 }) .with(atom::Force::new()) .with(atom::Position { - pos: Vector3::new(-5.0e-6, 5.0e-6, 5.0e-6), + pos: Vector3::new(1.0e-7, 0.0e-7, 0.0e-7), }) .with(atom::Velocity { vel: Vector3::new(0.0, 0.0, 0.0), @@ -90,8 +69,11 @@ fn main() { .with(lib::initiate::NewlyCreated) .build(); + // Define timestep + sim.world.insert(Timestep { delta: 1.0e-5 }); + // Run the simulation for a number of steps. - for _i in 0..100_000 { + for _i in 0..300_000 { sim.step(); } From 1fa446153c79cad3df5b0efa6445bf9d68ebd2e7 Mon Sep 17 00:00:00 2001 From: Brian Date: Fri, 9 Sep 2022 15:04:49 +0100 Subject: [PATCH 03/19] editing force to fix trap freq --- .DS_Store | Bin 0 -> 6148 bytes src/constant.rs | 7 ++- src/dipole/force.rs | 96 ++++++++++++++++++-------------- src/dipole/mod.rs | 18 +++--- src/laser/gaussian.rs | 41 +++++++++----- src/laser/intensity_gradient.rs | 7 ++- 6 files changed, 99 insertions(+), 70 deletions(-) create mode 100644 .DS_Store diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000000000000000000000000000000000000..5f246ddb48f3aac0e500c7205417caf653dee7c1 GIT binary patch literal 6148 zcmeHLze@u#6#illJzEgRg3F=OQE+xSipLn7 zY&SbuIxIx`4zGY$;6EzB?{0uCoZu8A@&49ZH?A+ulB7Mz+SC`;HkUW5htIpx-<|0n z_MN?r{;4GAN62dF-$8`8an*XEs|{Gw5rtfLFjP@T~yv4-pl{*kW!_Zyi+WD*$0Y-P)+F z_e#QuI|yToxk2h7Qzn(rq$+#GP$r%FBO4c6%nh1!C_OWcV`o?T75D;5`nkCP literal 0 HcmV?d00001 diff --git a/src/constant.rs b/src/constant.rs index 4ce4d969..bc3991ff 100644 --- a/src/constant.rs +++ b/src/constant.rs @@ -9,10 +9,10 @@ pub const GC: f64 = 9.80665; /// Mathematical constant exp(1) pub const EXP: f64 = std::f64::consts::E; -/// Mathematica constant pi +/// Mathematical constant pi pub const PI: f64 = std::f64::consts::PI; -// The Bohr magneton, defined in SI units of Joules/Tesla. +/// The Bohr magneton, defined in SI units of Joules/Tesla. pub const BOHRMAG: f64 = 9.274e-24; /// Boltzmann constant in SI units @@ -26,3 +26,6 @@ pub const C: f64 = 299297458.0; /// Sqrt of 2 pub const SQRT2: f64 = std::f64::consts::SQRT_2; + +/// Vacuum permittivity [F/m] +pub const EPSILON0: f64 = 8.8541878128e-12; diff --git a/src/dipole/force.rs b/src/dipole/force.rs index 233c0ca2..aeb9c27e 100644 --- a/src/dipole/force.rs +++ b/src/dipole/force.rs @@ -6,6 +6,7 @@ use crate::atom::Force; use crate::dipole::DipoleLight; use crate::dipole::Polarizability; use crate::laser::index::LaserIndex; +use crate::constant; /// Calculates forces exerted onto the atoms by dipole laser beams. /// @@ -28,10 +29,13 @@ impl<'a, const N: usize> System<'a> for ApplyDipoleForceSystem { ) { (&mut force, &polarizability, &gradient_sampler) .par_join() - .for_each(|(force, polarizability, sampler)| { + .for_each( |(force, polarizability, sampler)| { for (index, _dipole) in (&dipole_index, &dipole_light).join() { - force.force += - polarizability.prefactor * sampler.contents[index.index].gradient; + + force.force += ( 1.0 / ( 2.0 * constant::EPSILON0 * constant::C ) ) + * polarizability.prefactor + * sampler.contents[index.index].gradient; + } }); } @@ -50,6 +54,9 @@ pub mod tests { use crate::laser::gaussian::GaussianBeam; use crate::laser::DEFAULT_BEAM_LIMIT; use nalgebra::Vector3; + use crate::atom::Position; + use crate::laser::frame::Frame; + use crate::laser::intensity_gradient::LaserIntensityGradientSampler; #[test] fn test_apply_dipole_force_system() { @@ -110,46 +117,46 @@ pub mod tests { fn test_apply_dipole_force_again_system() { let mut test_world = World::new(); - test_world.register::(); - test_world.register::(); - test_world.register::(); - test_world.register::>(); - test_world.register::(); + test_world.register::(); + test_world.register::(); + test_world.register::(); + test_world.register::>(); + test_world.register::(); - test_world - .create_entity() - .with(LaserIndex { - index: 0, - initiated: true, - }) - .with(DipoleLight { - wavelength: 1064.0e-9, - }) - .build(); + test_world + .create_entity() + .with(LaserIndex { + index: 0, + initiated: true, + }) + .with(DipoleLight { + wavelength: 1064.0e-9, + }) + .build(); - let transition = Polarizability::calculate_for(1064e-9, 461e-9, 32e6); - let atom1 = test_world - .create_entity() - .with(Force { - force: Vector3::new(0.0, 0.0, 0.0), - }) - .with(LaserIntensityGradientSamplers { - contents: [crate::laser::intensity_gradient::LaserIntensityGradientSampler { - gradient: Vector3::new(-8.4628e+7, -4.33992902e+13, -4.33992902e+13), - }; crate::laser::DEFAULT_BEAM_LIMIT], - }) - .with(transition) - .build(); - let mut system = ApplyDipoleForceSystem::<{ DEFAULT_BEAM_LIMIT }>; - system.run_now(&test_world); - test_world.maintain(); - let sampler_storage = test_world.read_storage::(); - let sim_result_force = sampler_storage.get(atom1).expect("Entity not found!").force; + let transition = Polarizability::calculate_for(1064e-9, 461e-9, 32e6); + let atom1 = test_world + .create_entity() + .with(Force { + force: Vector3::new(0.0, 0.0, 0.0), + }) + .with(LaserIntensityGradientSamplers { + contents: [crate::laser::intensity_gradient::LaserIntensityGradientSampler { + gradient: Vector3::new(-8.4628e+7, -4.33992902e+13, -4.33992902e+13), + }; crate::laser::DEFAULT_BEAM_LIMIT], + }) + .with(transition) + .build(); + let mut system = ApplyDipoleForceSystem::<{ DEFAULT_BEAM_LIMIT }>; + system.run_now(&test_world); + test_world.maintain(); + let sampler_storage = test_world.read_storage::(); + let sim_result_force = sampler_storage.get(atom1).expect("Entity not found!").force; - assert_approx_eq!(-6.386888332902177e-29, sim_result_force[0], 3e-30_f64); - assert_approx_eq!(-3.11151847e-23, sim_result_force[1], 2e-24_f64); - assert_approx_eq!(-3.11151847e-23, sim_result_force[2], 2e-24_f64); - } + assert_approx_eq!(-6.386888332902177e-29, sim_result_force[0], 3e-30_f64); + assert_approx_eq!(-3.11151847e-23, sim_result_force[1], 2e-24_f64); + assert_approx_eq!(-3.11151847e-23, sim_result_force[2], 2e-24_f64); + } #[test] fn test_apply_dipole_force_and_gradient_system() { @@ -245,9 +252,12 @@ pub mod tests { .get(atom1) .expect("Entity not found!") .contents; - //println!("force is: {}", sim_result_force); - //println!("gradient 1 is: {}", sim_result_grad[0].gradient); - //println!("gradient 2 is: {}", sim_result_grad[1].gradient); + + println!("force is: {}", sim_result_force); + // println!("gradient 1 is: {}", _sim_result_grad); + // println!("gradient 2 is: {}", _sim_result_grad[1].gradient); + + assert_approx_eq!( 0.000000000000000000000000000000000127913190642808, diff --git a/src/dipole/mod.rs b/src/dipole/mod.rs index 8c0d4b42..4c4c4d47 100644 --- a/src/dipole/mod.rs +++ b/src/dipole/mod.rs @@ -46,7 +46,7 @@ impl Component for Polarizability { type Storage = VecStorage; } impl Polarizability { - /// Calculate the polarizability of an atom in a dipole beam of given wavelength, detuned from a strong optical transition. + /// Calculate the real part of the polarizability of an atom in a dipole beam of given wavelength, detuned from a strong optical transition. /// /// The wavelengths of both transitions are in SI units of m. /// The linewidth of the optical transition is in SI units of Hz. @@ -55,12 +55,16 @@ impl Polarizability { optical_transition_wavelength: f64, optical_transition_linewidth: f64, ) -> Polarizability { + let transition_f = constant::C / optical_transition_wavelength; + let dipole_f = constant::C / dipole_beam_wavelength; - let prefactor = -3. * constant::PI * constant::C.powf(2.0) - / (2. * (2. * constant::PI * transition_f).powf(3.0)) + + let prefactor = 3.0 * constant::PI * constant::C.powf(3.0) + * constant::EPSILON0/(transition_f.powf(3.0)) * optical_transition_linewidth - * -(1. / (transition_f - dipole_f) + 1. / (transition_f + dipole_f)); + * (1.0 / (transition_f - dipole_f) + 1.0 / (transition_f + dipole_f)); + Polarizability { prefactor } } } @@ -83,11 +87,11 @@ impl<'a> System<'a> for AttachIndexToDipoleLightSystem { } /// This plugin implements a dipole force that can be used to confine cold atoms. -/// +/// /// See also [crate::dipole] -/// +/// /// # Generic Arguments -/// +/// /// * `N`: The maximum number of laser beams (must match the `LaserPlugin`). pub struct DipolePlugin; impl Plugin for DipolePlugin { diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index 4da1d11f..ecba1228 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -209,6 +209,7 @@ pub fn get_gaussian_beam_intensity( / (beam.e_radius.powf(2.0) * (1. + (z / beam.rayleigh_range).powf(2.0))), ) } + /// Computes the rayleigh range for a given beam and wavelength pub fn calculate_rayleigh_range(wavelength: &f64, e_radius: &f64) -> f64 { 2.0 * PI * e_radius.powf(2.0) / wavelength @@ -230,15 +231,24 @@ pub fn get_gaussian_beam_intensity_gradient( let y = rela_coord.dot(&reference_frame.y_vector) * semi_major_axis.powf(0.5); let z = rela_coord.dot(&beam.direction); - let spot_size_squared = - 2.0 * beam.e_radius.powf(2.0) * (1. + (z / beam.rayleigh_range).powf(2.0)); - let vector = -4. * (reference_frame.x_vector * x + reference_frame.y_vector * y) - + beam.direction * z / (beam.rayleigh_range.powf(2.0) + z.powf(2.0)) - * (- 2.0 * spot_size_squared + 4. * (x.powf(2.0) + y.powf(2.0))); - let intensity = 2. * beam.power / PI / spot_size_squared - * EXP.powf(-2. * (x.powf(2.0) + y.powf(2.0)) / spot_size_squared); + let spot_size_squared = 2.0 * beam.e_radius.powf(2.0) + * (1. + (z / beam.rayleigh_range).powf(2.0)); + + // let vector = -4.0 * (reference_frame.x_vector * x + reference_frame.y_vector * y) + // - 2.0 * beam.direction * z / ( beam.rayleigh_range.powf(2.0) + z.powf(2.0) ) * ( + // spot_size_squared * ( + // 1.0 + ( z / beam.rayleigh_range ).powf(2.0) + // ) - 2.0 * (x.powf(2.0) + y.powf(2.0)) + // ); + + let vector = -4.0 * (reference_frame.x_vector * x + reference_frame.y_vector * y) + -2.0 * beam.direction * z / ( beam.rayleigh_range.powf(2.0) + z.powf(2.0) ) + * ( spot_size_squared - 2.0 * ( x.powf(2.0) + y.powf(2.0) ) ); + + let intensity = 2.0 * beam.power / ( PI * spot_size_squared ) + * EXP.powf(-2.0 * (x.powf(2.0) + y.powf(2.0)) / spot_size_squared); - intensity / spot_size_squared * vector + - intensity / spot_size_squared * vector } #[cfg(test)] @@ -258,13 +268,13 @@ pub mod tests { let beam = GaussianBeam { direction: Vector3::z(), intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 70.71067812e-6, - power: 100.0, - rayleigh_range: calculate_rayleigh_range(&1064.0e-9, &70.71067812e-6), + e_radius: 60e-6, + power: 10.0, + rayleigh_range: calculate_rayleigh_range(&1064.0e-9, &60e-6), ellipticity: 0.0, }; let pos1 = Position { - pos: Vector3::new(10.0e-6, 0.0, 30.0e-6), + pos: Vector3::new(1.0e-7, 1.2e-6, 6.0e-7), }; let grf = Frame { x_vector: Vector3::x(), @@ -272,9 +282,10 @@ pub mod tests { }; let gradient = get_gaussian_beam_intensity_gradient(&beam, &pos1, &grf); - assert_approx_eq!(gradient[0], -2.49605032e+13, 1e+8_f64); - assert_approx_eq!(gradient[1], 0.0, 1e+9_f64); - assert_approx_eq!(gradient[2], -2.06143366e+08, 1e+6_f64); + println!("{}",gradient[0]); + assert_approx_eq!(gradient[0], -4.91021e+10, 1e+9_f64); + assert_approx_eq!(gradient[1], -5.89225e+11, 1e+9_f64); + assert_approx_eq!(gradient[2], -9.38334e+6, 1e+9_f64); } #[test] diff --git a/src/laser/intensity_gradient.rs b/src/laser/intensity_gradient.rs index 395206b7..23e974ba 100644 --- a/src/laser/intensity_gradient.rs +++ b/src/laser/intensity_gradient.rs @@ -173,6 +173,7 @@ pub mod tests { 1e+5_f64 ); } + #[test] fn test_sample_laser_intensity_gradient_again_system() { let mut test_world = World::new(); @@ -233,8 +234,8 @@ pub mod tests { .contents[0] .gradient; - assert_approx_eq!(-8.4628e+7, sim_result_gradient[0], 1e+5_f64); - assert_approx_eq!(-4.33992902e+13, sim_result_gradient[1], 1e+8_f64); - assert_approx_eq!(-4.33992902e+13, sim_result_gradient[2], 1e+8_f64); + assert_approx_eq!( -2.09081e+8, sim_result_gradient[0], 1e+5_f64); + assert_approx_eq!(-4.33993e+13, sim_result_gradient[1], 1e+8_f64); + assert_approx_eq!(-4.33993e+13, sim_result_gradient[2], 1e+8_f64); } } From 0d5412b62113f6e253835fa3aa64bc50cd2ebdef Mon Sep 17 00:00:00 2001 From: Brian Date: Sat, 10 Sep 2022 22:55:55 +0100 Subject: [PATCH 04/19] editing force to fix trap freq --- examples/dipole_trap.rs | 14 +++++++------- src/dipole/force.rs | 1 + src/laser/gaussian.rs | 32 ++++++++++++++++++-------------- 3 files changed, 26 insertions(+), 21 deletions(-) diff --git a/examples/dipole_trap.rs b/examples/dipole_trap.rs index 4942fac2..1a4ebb57 100644 --- a/examples/dipole_trap.rs +++ b/examples/dipole_trap.rs @@ -29,7 +29,7 @@ fn main() { let mut sim = sim_builder.build(); // Create dipole laser. - let power = 10.0; + let power = 1.0; let e_radius = 60.0e-6 / (2.0_f64.sqrt()); let wavelength = 1064.0e-9; @@ -37,7 +37,7 @@ fn main() { intersection: Vector3::new(0.0, 0.0, 0.0), e_radius, power, - direction: Vector3::x(), + direction: Vector3::z(), rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range(&wavelength, &e_radius), ellipticity: 0.0, }; @@ -46,8 +46,8 @@ fn main() { .with(gaussian_beam) .with(dipole::DipoleLight { wavelength }) .with(laser::frame::Frame { - x_vector: Vector3::y(), - y_vector: Vector3::z(), + x_vector: Vector3::x(), + y_vector: Vector3::y(), }) .build(); @@ -57,7 +57,7 @@ fn main() { .with(atom::Mass { value: 87.0 }) .with(atom::Force::new()) .with(atom::Position { - pos: Vector3::new(1.0e-7, 0.0e-7, 0.0e-7), + pos: Vector3::new(1.0e-7, 1.0e-7, 1.0e-7), }) .with(atom::Velocity { vel: Vector3::new(0.0, 0.0, 0.0), @@ -70,10 +70,10 @@ fn main() { .build(); // Define timestep - sim.world.insert(Timestep { delta: 1.0e-5 }); + sim.world.insert(Timestep { delta: 1.0e-6 }); // Run the simulation for a number of steps. - for _i in 0..300_000 { + for _i in 0..200_000 { sim.step(); } diff --git a/src/dipole/force.rs b/src/dipole/force.rs index aeb9c27e..a247fe39 100644 --- a/src/dipole/force.rs +++ b/src/dipole/force.rs @@ -36,6 +36,7 @@ impl<'a, const N: usize> System<'a> for ApplyDipoleForceSystem { * polarizability.prefactor * sampler.contents[index.index].gradient; + } }); } diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index ecba1228..df90feb8 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -222,33 +222,37 @@ pub fn get_gaussian_beam_intensity_gradient( pos: &Position, reference_frame: &Frame, ) -> Vector3 { + + let rela_coord = pos.pos - beam.intersection; // ellipticity treatment - let semi_major_axis = 1.0 / (1.0 - beam.ellipticity.powf(2.0)).powf(0.5); + //let semi_major_axis = 1.0 / (1.0 - beam.ellipticity.powf(2.0)).powf(0.5); - let x = rela_coord.dot(&reference_frame.x_vector) / semi_major_axis.powf(0.5); - let y = rela_coord.dot(&reference_frame.y_vector) * semi_major_axis.powf(0.5); + // println!("{}", beam.intersection); + let x = rela_coord.dot(&reference_frame.x_vector);// / semi_major_axis.powf(0.5); + let y = rela_coord.dot(&reference_frame.y_vector);// * semi_major_axis.powf(0.5); let z = rela_coord.dot(&beam.direction); - let spot_size_squared = 2.0 * beam.e_radius.powf(2.0) - * (1. + (z / beam.rayleigh_range).powf(2.0)); + let spot_size_squared = + 2.0 * beam.e_radius.powf(2.0) * (1. + (z / beam.rayleigh_range).powf(2.0)); + + let vector = -4.0 * (reference_frame.x_vector * x + reference_frame.y_vector * y) + - 2.0 * beam.direction * z / ( beam.rayleigh_range.powf(2.0) + z.powf(2.0) ) * ( + spot_size_squared * ( + 1.0 + ( z / beam.rayleigh_range ).powf(2.0) + ) - 2.0 * (x.powf(2.0) + y.powf(2.0)) + ); // let vector = -4.0 * (reference_frame.x_vector * x + reference_frame.y_vector * y) - // - 2.0 * beam.direction * z / ( beam.rayleigh_range.powf(2.0) + z.powf(2.0) ) * ( - // spot_size_squared * ( - // 1.0 + ( z / beam.rayleigh_range ).powf(2.0) - // ) - 2.0 * (x.powf(2.0) + y.powf(2.0)) - // ); + // - 2.0 * beam.direction * ( z / ( beam.rayleigh_range.powf(2.0) + z.powf(2.0) ) ) + // * ( spot_size_squared - 2.0 * ( x.powf(2.0) + y.powf(2.0) ) ); - let vector = -4.0 * (reference_frame.x_vector * x + reference_frame.y_vector * y) - -2.0 * beam.direction * z / ( beam.rayleigh_range.powf(2.0) + z.powf(2.0) ) - * ( spot_size_squared - 2.0 * ( x.powf(2.0) + y.powf(2.0) ) ); let intensity = 2.0 * beam.power / ( PI * spot_size_squared ) * EXP.powf(-2.0 * (x.powf(2.0) + y.powf(2.0)) / spot_size_squared); - - intensity / spot_size_squared * vector + intensity / spot_size_squared * vector } #[cfg(test)] From 764261ffff859c567e9ab9983640a2c5aa5f060f Mon Sep 17 00:00:00 2001 From: Brian Date: Sun, 11 Sep 2022 21:58:07 +0100 Subject: [PATCH 05/19] Changing grad(Intensity) --- examples/dipole_trap.rs | 6 +++--- src/laser/gaussian.rs | 19 ++++++++++--------- 2 files changed, 13 insertions(+), 12 deletions(-) diff --git a/examples/dipole_trap.rs b/examples/dipole_trap.rs index 1a4ebb57..16d31ec8 100644 --- a/examples/dipole_trap.rs +++ b/examples/dipole_trap.rs @@ -37,7 +37,7 @@ fn main() { intersection: Vector3::new(0.0, 0.0, 0.0), e_radius, power, - direction: Vector3::z(), + direction: Vector3::x(), rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range(&wavelength, &e_radius), ellipticity: 0.0, }; @@ -46,8 +46,8 @@ fn main() { .with(gaussian_beam) .with(dipole::DipoleLight { wavelength }) .with(laser::frame::Frame { - x_vector: Vector3::x(), - y_vector: Vector3::y(), + x_vector: Vector3::y(), + y_vector: Vector3::z(), }) .build(); diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index df90feb8..571ef364 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -237,22 +237,23 @@ pub fn get_gaussian_beam_intensity_gradient( let spot_size_squared = 2.0 * beam.e_radius.powf(2.0) * (1. + (z / beam.rayleigh_range).powf(2.0)); - let vector = -4.0 * (reference_frame.x_vector * x + reference_frame.y_vector * y) - - 2.0 * beam.direction * z / ( beam.rayleigh_range.powf(2.0) + z.powf(2.0) ) * ( - spot_size_squared * ( - 1.0 + ( z / beam.rayleigh_range ).powf(2.0) - ) - 2.0 * (x.powf(2.0) + y.powf(2.0)) - ); - // let vector = -4.0 * (reference_frame.x_vector * x + reference_frame.y_vector * y) - // - 2.0 * beam.direction * ( z / ( beam.rayleigh_range.powf(2.0) + z.powf(2.0) ) ) - // * ( spot_size_squared - 2.0 * ( x.powf(2.0) + y.powf(2.0) ) ); + // - 2.0 * beam.direction * z / ( beam.rayleigh_range.powf(2.0) + z.powf(2.0) ) * ( + // spot_size_squared * ( + // 1.0 + ( z / beam.rayleigh_range ).powf(2.0) + // ) - 2.0 * (x.powf(2.0) + y.powf(2.0)) + // ); + + let vector = -4.0 * ( reference_frame.x_vector * x + reference_frame.y_vector * y ) + - 2.0 * beam.direction * ( z / ( beam.rayleigh_range.powf(2.0) + z.powf(2.0) ) ) + * ( spot_size_squared - 2.0 * ( x.powf(2.0) + y.powf(2.0) ) ); let intensity = 2.0 * beam.power / ( PI * spot_size_squared ) * EXP.powf(-2.0 * (x.powf(2.0) + y.powf(2.0)) / spot_size_squared); intensity / spot_size_squared * vector + } #[cfg(test)] From daaa2a7709af4f4331a569868c204beb8c7ab72f Mon Sep 17 00:00:00 2001 From: Brian Date: Mon, 12 Sep 2022 12:30:43 +0100 Subject: [PATCH 06/19] fixing polarizability --- src/dipole/mod.rs | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/src/dipole/mod.rs b/src/dipole/mod.rs index 4c4c4d47..2196d087 100644 --- a/src/dipole/mod.rs +++ b/src/dipole/mod.rs @@ -60,10 +60,22 @@ impl Polarizability { let dipole_f = constant::C / dipole_beam_wavelength; - let prefactor = 3.0 * constant::PI * constant::C.powf(3.0) - * constant::EPSILON0/(transition_f.powf(3.0)) - * optical_transition_linewidth - * (1.0 / (transition_f - dipole_f) + 1.0 / (transition_f + dipole_f)); + let prefactor = ( 3.0 * constant::PI * constant::C.powf(3.0) * constant::EPSILON0 / + transition_f.powf(3.0) ) * + ( optical_transition_linewidth / (transition_f - dipole_f) + + optical_transition_linewidth / (transition_f + dipole_f) + ); + + + println!("optical transition linewidth {}", optical_transition_linewidth); + println!("polarizability {}", prefactor); + + println!("constant::PI = {}", constant::PI); + println!("constant::C.powf(3.0) = {}", constant::C.powf(3.0)); + println!("constant::C = {}", constant::C); + println!("constant::EPSILON0 = {}", constant::EPSILON0); + println!("transition_f = {}", transition_f); + println!("dipole_f = {}", dipole_f); Polarizability { prefactor } } From 4d47ecccd527d5fd36bce3c5622d49389eb1d0b1 Mon Sep 17 00:00:00 2001 From: Brian Date: Mon, 12 Sep 2022 15:39:28 +0100 Subject: [PATCH 07/19] updating --- src/constant.rs | 2 +- src/dipole/force.rs | 9 +++------ 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/src/constant.rs b/src/constant.rs index bc3991ff..4f36537f 100644 --- a/src/constant.rs +++ b/src/constant.rs @@ -22,7 +22,7 @@ pub const BOLTZCONST: f64 = 1.38e-23; pub const AMU: f64 = 1.6605e-27; /// Speed of light in SI units of m/s -pub const C: f64 = 299297458.0; +pub const C: f64 = 299792458.0; /// Sqrt of 2 pub const SQRT2: f64 = std::f64::consts::SQRT_2; diff --git a/src/dipole/force.rs b/src/dipole/force.rs index a247fe39..57133b81 100644 --- a/src/dipole/force.rs +++ b/src/dipole/force.rs @@ -31,12 +31,9 @@ impl<'a, const N: usize> System<'a> for ApplyDipoleForceSystem { .par_join() .for_each( |(force, polarizability, sampler)| { for (index, _dipole) in (&dipole_index, &dipole_light).join() { - - force.force += ( 1.0 / ( 2.0 * constant::EPSILON0 * constant::C ) ) - * polarizability.prefactor - * sampler.contents[index.index].gradient; - - + force.force += polarizability.prefactor*sampler.contents[index.index].gradient/ + ( 2.0 * constant::EPSILON0 * constant::C ); + println!("force = {}", force.force) } }); } From 8b837e24a7c5af6571b76ee3786f6c2f43d585e6 Mon Sep 17 00:00:00 2001 From: Brian Bostwick Date: Mon, 12 Sep 2022 17:10:05 +0100 Subject: [PATCH 08/19] update --- examples/dipole_trap.rs | 2 +- src/dipole/force.rs | 1 - src/dipole/mod.rs | 19 ++----------------- src/laser/gaussian.rs | 26 ++++++++------------------ src/laser/intensity_gradient.rs | 8 ++++---- 5 files changed, 15 insertions(+), 41 deletions(-) diff --git a/examples/dipole_trap.rs b/examples/dipole_trap.rs index 16d31ec8..471034e3 100644 --- a/examples/dipole_trap.rs +++ b/examples/dipole_trap.rs @@ -30,7 +30,7 @@ fn main() { // Create dipole laser. let power = 1.0; - let e_radius = 60.0e-6 / (2.0_f64.sqrt()); + let e_radius = 60.0e-6 / 2.0_f64.sqrt(); let wavelength = 1064.0e-9; let gaussian_beam = GaussianBeam { diff --git a/src/dipole/force.rs b/src/dipole/force.rs index 57133b81..cd02c433 100644 --- a/src/dipole/force.rs +++ b/src/dipole/force.rs @@ -33,7 +33,6 @@ impl<'a, const N: usize> System<'a> for ApplyDipoleForceSystem { for (index, _dipole) in (&dipole_index, &dipole_light).join() { force.force += polarizability.prefactor*sampler.contents[index.index].gradient/ ( 2.0 * constant::EPSILON0 * constant::C ); - println!("force = {}", force.force) } }); } diff --git a/src/dipole/mod.rs b/src/dipole/mod.rs index 2196d087..ed45da8d 100644 --- a/src/dipole/mod.rs +++ b/src/dipole/mod.rs @@ -55,28 +55,13 @@ impl Polarizability { optical_transition_wavelength: f64, optical_transition_linewidth: f64, ) -> Polarizability { - let transition_f = constant::C / optical_transition_wavelength; - let dipole_f = constant::C / dipole_beam_wavelength; - let prefactor = ( 3.0 * constant::PI * constant::C.powf(3.0) * constant::EPSILON0 / - transition_f.powf(3.0) ) * - ( optical_transition_linewidth / (transition_f - dipole_f) + + transition_f.powf(3.0) ) * ( + optical_transition_linewidth / (transition_f - dipole_f) + optical_transition_linewidth / (transition_f + dipole_f) ); - - - println!("optical transition linewidth {}", optical_transition_linewidth); - println!("polarizability {}", prefactor); - - println!("constant::PI = {}", constant::PI); - println!("constant::C.powf(3.0) = {}", constant::C.powf(3.0)); - println!("constant::C = {}", constant::C); - println!("constant::EPSILON0 = {}", constant::EPSILON0); - println!("transition_f = {}", transition_f); - println!("dipole_f = {}", dipole_f); - Polarizability { prefactor } } } diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index 571ef364..bfbf011f 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -178,10 +178,7 @@ pub fn get_gaussian_beam_intensity( let semi_major_axis = 1.0 / (1.0 - beam.ellipticity.powf(2.0)).powf(0.5); // the factor (1.0 / semi_major_axis) is necessary so the overall power of the beam is not changed. - ( - z, - (1.0 / semi_major_axis) * ((x).powf(2.0) + (y * semi_major_axis).powf(2.0)), - ) + ( z, (1.0 / semi_major_axis) * ((x).powf(2.0) + (y * semi_major_axis).powf(2.0)), ) } // ellipticity will be ignored (i.e. treated as zero) if no `Frame` is supplied. None => { @@ -193,6 +190,7 @@ pub fn get_gaussian_beam_intensity( (z, distance * distance) } }; + let power = match mask { Some(mask) => { if distance_squared.powf(0.5) < mask.radius { @@ -227,32 +225,24 @@ pub fn get_gaussian_beam_intensity_gradient( let rela_coord = pos.pos - beam.intersection; // ellipticity treatment - //let semi_major_axis = 1.0 / (1.0 - beam.ellipticity.powf(2.0)).powf(0.5); + let semi_major_axis = 1.0 / (1.0 - beam.ellipticity.powf(2.0)).powf(0.5); - // println!("{}", beam.intersection); - let x = rela_coord.dot(&reference_frame.x_vector);// / semi_major_axis.powf(0.5); - let y = rela_coord.dot(&reference_frame.y_vector);// * semi_major_axis.powf(0.5); + let x = rela_coord.dot(&reference_frame.x_vector) / semi_major_axis.powf(0.5); + let y = rela_coord.dot(&reference_frame.y_vector) * semi_major_axis.powf(0.5); let z = rela_coord.dot(&beam.direction); let spot_size_squared = 2.0 * beam.e_radius.powf(2.0) * (1. + (z / beam.rayleigh_range).powf(2.0)); - // let vector = -4.0 * (reference_frame.x_vector * x + reference_frame.y_vector * y) - // - 2.0 * beam.direction * z / ( beam.rayleigh_range.powf(2.0) + z.powf(2.0) ) * ( - // spot_size_squared * ( - // 1.0 + ( z / beam.rayleigh_range ).powf(2.0) - // ) - 2.0 * (x.powf(2.0) + y.powf(2.0)) - // ); - let vector = -4.0 * ( reference_frame.x_vector * x + reference_frame.y_vector * y ) - - 2.0 * beam.direction * ( z / ( beam.rayleigh_range.powf(2.0) + z.powf(2.0) ) ) - * ( spot_size_squared - 2.0 * ( x.powf(2.0) + y.powf(2.0) ) ); + + 2.0 * beam.direction * ( z / ( beam.rayleigh_range.powf(2.0) + z.powf(2.0) ) ) + * ( x.powf(2.0) + y.powf(2.0) - spot_size_squared/2.0 ); let intensity = 2.0 * beam.power / ( PI * spot_size_squared ) * EXP.powf(-2.0 * (x.powf(2.0) + y.powf(2.0)) / spot_size_squared); - intensity / spot_size_squared * vector + intensity / spot_size_squared * vector / 2.0 } diff --git a/src/laser/intensity_gradient.rs b/src/laser/intensity_gradient.rs index 23e974ba..0e2408fe 100644 --- a/src/laser/intensity_gradient.rs +++ b/src/laser/intensity_gradient.rs @@ -64,11 +64,11 @@ impl<'a, const N: usize> System<'a> for SampleGaussianLaserIntensityGradientSyst use rayon::prelude::*; for (_dipole, index, beam, reference) in - (&dipole, &index, &gaussian, &reference_frame).join() - { + (&dipole, &index, &gaussian, &reference_frame).join() { (&pos, &mut sampler).par_join().for_each(|(pos, sampler)| { - sampler.contents[index.index].gradient = - get_gaussian_beam_intensity_gradient(beam, pos, reference); + sampler.contents[index.index].gradient = get_gaussian_beam_intensity_gradient( + beam, pos, reference + ); }); } } From 1380dd8fb2873d067f96c75cbabbdc76e037ac4b Mon Sep 17 00:00:00 2001 From: Brian Bostwick Date: Mon, 12 Sep 2022 17:11:00 +0100 Subject: [PATCH 09/19] update --- src/laser/gaussian.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index bfbf011f..22ff37b6 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -242,7 +242,7 @@ pub fn get_gaussian_beam_intensity_gradient( let intensity = 2.0 * beam.power / ( PI * spot_size_squared ) * EXP.powf(-2.0 * (x.powf(2.0) + y.powf(2.0)) / spot_size_squared); - intensity / spot_size_squared * vector / 2.0 + intensity / spot_size_squared * vector } From 0819d14e3990f6aa753ae8bec85e143f1350ee40 Mon Sep 17 00:00:00 2001 From: Brian Bostwick Date: Tue, 13 Sep 2022 10:53:04 +0100 Subject: [PATCH 10/19] fixing sqrt(2) for beam intesity in gaussian beam --- src/laser/gaussian.rs | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index 22ff37b6..a0dab0a5 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -68,8 +68,8 @@ impl GaussianBeam { peak_intensity: f64, e_radius: f64, ) -> Self { - let std = e_radius / 2.0_f64.powf(0.5); - let power = 2.0 * std::f64::consts::PI * std.powi(2) * peak_intensity; + let std = e_radius * 2.0_f64.powf(0.5); + let power = 0.5 * std::f64::consts::PI * std.powi(2) * peak_intensity; GaussianBeam { intersection, direction, @@ -166,20 +166,25 @@ pub fn get_gaussian_beam_intensity( mask: Option<&CircularMask>, frame: Option<&Frame>, ) -> f64 { + let (z, distance_squared) = match frame { + // checking if frame is given (for calculating ellipticity) Some(frame) => { + let (x, y, z) = maths::get_relative_coordinates_line_point( &pos.pos, &beam.intersection, &beam.direction, frame, ); + let semi_major_axis = 1.0 / (1.0 - beam.ellipticity.powf(2.0)).powf(0.5); // the factor (1.0 / semi_major_axis) is necessary so the overall power of the beam is not changed. ( z, (1.0 / semi_major_axis) * ((x).powf(2.0) + (y * semi_major_axis).powf(2.0)), ) } + // ellipticity will be ignored (i.e. treated as zero) if no `Frame` is supplied. None => { let (distance, z) = maths::get_minimum_distance_line_point( @@ -201,13 +206,17 @@ pub fn get_gaussian_beam_intensity( } None => beam.power, }; - power / PI / beam.e_radius.powf(2.0) / (1.0 + (z / beam.rayleigh_range).powf(2.0)) + power/ PI / beam.e_radius.powf(2.0) / (1.0 + (z / beam.rayleigh_range).powf(2.0)) * EXP.powf( -distance_squared / (beam.e_radius.powf(2.0) * (1. + (z / beam.rayleigh_range).powf(2.0))), ) } + + + + /// Computes the rayleigh range for a given beam and wavelength pub fn calculate_rayleigh_range(wavelength: &f64, e_radius: &f64) -> f64 { 2.0 * PI * e_radius.powf(2.0) / wavelength From a7f20302596c32f930fb22963a51663983aa6d47 Mon Sep 17 00:00:00 2001 From: Brian Bostwick Date: Tue, 13 Sep 2022 16:41:13 +0100 Subject: [PATCH 11/19] Fixed inconsitant use of e^2 radiouse --- examples/top_trap_with_collisions.rs | 3 +- src/dipole/force.rs | 21 ++++++---- src/laser/gaussian.rs | 26 ++++++------- src/laser/intensity.rs | 57 ++++++++++++++++++++++------ src/laser/intensity_gradient.rs | 18 +++++++-- 5 files changed, 86 insertions(+), 39 deletions(-) diff --git a/examples/top_trap_with_collisions.rs b/examples/top_trap_with_collisions.rs index 5e5174a8..395d0d33 100644 --- a/examples/top_trap_with_collisions.rs +++ b/examples/top_trap_with_collisions.rs @@ -7,7 +7,8 @@ use lib::collisions::CollisionPlugin; use lib::collisions::{ApplyCollisionsOption, CollisionParameters, CollisionsTracker}; use lib::initiate::NewlyCreated; use lib::integrator::Timestep; -use lib::magnetic::force::{ApplyMagneticForceSystem, MagneticDipole}; +// use lib::magnetic::force::{ApplyMagneticForceSystem, MagneticDipole}; +use lib::magnetic::force::{MagneticDipole}; use lib::magnetic::quadrupole::QuadrupoleField3D; use lib::magnetic::top::TimeOrbitingPotential; use lib::magnetic::MagneticTrapPlugin; diff --git a/src/dipole/force.rs b/src/dipole/force.rs index cd02c433..ffa92386 100644 --- a/src/dipole/force.rs +++ b/src/dipole/force.rs @@ -51,9 +51,9 @@ pub mod tests { use crate::laser::gaussian::GaussianBeam; use crate::laser::DEFAULT_BEAM_LIMIT; use nalgebra::Vector3; - use crate::atom::Position; - use crate::laser::frame::Frame; - use crate::laser::intensity_gradient::LaserIntensityGradientSampler; + // use crate::atom::Position; + // use crate::laser::frame::Frame; + // use crate::laser::intensity_gradient::LaserIntensityGradientSampler; #[test] fn test_apply_dipole_force_system() { @@ -179,6 +179,9 @@ pub mod tests { rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range(&1064.0e-9, &e_radius), ellipticity: 0.0, }; + + println!("rayleigh_range {}", gaussian_beam.rayleigh_range); + test_world .create_entity() .with(gaussian_beam) @@ -194,6 +197,7 @@ pub mod tests { y_vector: Vector3::z(), }) .build(); + let gaussian_beam = GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), e_radius, @@ -202,6 +206,7 @@ pub mod tests { rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range(&1064.0e-9, &e_radius), ellipticity: 0.0, }; + test_world .create_entity() .with(gaussian_beam) @@ -219,6 +224,7 @@ pub mod tests { .build(); let transition = Polarizability::calculate_for(1064e-9, 460.7e-9, 32e6); + let atom1 = test_world .create_entity() .with(crate::atom::Position { @@ -251,23 +257,22 @@ pub mod tests { .contents; println!("force is: {}", sim_result_force); - // println!("gradient 1 is: {}", _sim_result_grad); - // println!("gradient 2 is: {}", _sim_result_grad[1].gradient); + println!("gradient 1 is: {}", _sim_result_grad[0].gradient); assert_approx_eq!( - 0.000000000000000000000000000000000127913190642808, + 3.17164e-32, sim_result_force[0], 3e-46_f64 ); assert_approx_eq!( - 0.000000000000000000000000000000000127913190642808, + 3.17164e-32, sim_result_force[1], 2e-46_f64 ); assert_approx_eq!( - 0.000000000000000000000000000000000511875188257342, + 2.2692e-31, sim_result_force[2], 2e-46_f64 ); diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index a0dab0a5..7ca6acd0 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -69,7 +69,7 @@ impl GaussianBeam { e_radius: f64, ) -> Self { let std = e_radius * 2.0_f64.powf(0.5); - let power = 0.5 * std::f64::consts::PI * std.powi(2) * peak_intensity; + let power = std::f64::consts::PI * std.powi(2) * peak_intensity; GaussianBeam { intersection, direction, @@ -102,8 +102,8 @@ impl GaussianBeam { e_radius: f64, wavelength: f64, ) -> Self { - let std = e_radius / 2.0_f64.powf(0.5); - let power = 2.0 * std::f64::consts::PI * std.powi(2) * peak_intensity; + let std = e_radius * 2.0_f64.powf(0.5); + let power = std::f64::consts::PI * std.powi(2) * peak_intensity; GaussianBeam { intersection, direction, @@ -206,16 +206,12 @@ pub fn get_gaussian_beam_intensity( } None => beam.power, }; - power/ PI / beam.e_radius.powf(2.0) / (1.0 + (z / beam.rayleigh_range).powf(2.0)) - * EXP.powf( - -distance_squared - / (beam.e_radius.powf(2.0) * (1. + (z / beam.rayleigh_range).powf(2.0))), - ) -} - - + let spot_size_squared = + 2.0 * beam.e_radius.powf(2.0) * (1. + (z / beam.rayleigh_range).powf(2.0)); + 2.0 * power / ( PI * spot_size_squared ) * EXP.powf(-2.0 * distance_squared / spot_size_squared) +} /// Computes the rayleigh range for a given beam and wavelength pub fn calculate_rayleigh_range(wavelength: &f64, e_radius: &f64) -> f64 { @@ -243,15 +239,17 @@ pub fn get_gaussian_beam_intensity_gradient( let spot_size_squared = 2.0 * beam.e_radius.powf(2.0) * (1. + (z / beam.rayleigh_range).powf(2.0)); - let vector = -4.0 * ( reference_frame.x_vector * x + reference_frame.y_vector * y ) + + // this will match mathematica if the ebam directionis mult py 2 + let vector = - 4.0 * ( reference_frame.x_vector * x + reference_frame.y_vector * y ) + 2.0 * beam.direction * ( z / ( beam.rayleigh_range.powf(2.0) + z.powf(2.0) ) ) - * ( x.powf(2.0) + y.powf(2.0) - spot_size_squared/2.0 ); + * ( 2.0 * (x.powf(2.0) + y.powf(2.0)) - spot_size_squared); let intensity = 2.0 * beam.power / ( PI * spot_size_squared ) * EXP.powf(-2.0 * (x.powf(2.0) + y.powf(2.0)) / spot_size_squared); - intensity / spot_size_squared * vector + vector * intensity / spot_size_squared } diff --git a/src/laser/intensity.rs b/src/laser/intensity.rs index e3d8b590..ae6dc46a 100644 --- a/src/laser/intensity.rs +++ b/src/laser/intensity.rs @@ -153,39 +153,45 @@ pub mod tests { test_world.register::(); test_world.register::>(); + let beam = GaussianBeam { + direction: Vector3::x(), + intersection: Vector3::new(0.0, 0.0, 0.0), + e_radius: 2.0, + power: 1.0, + rayleigh_range: gaussian::calculate_rayleigh_range(&461.0e-9, &2.0), + ellipticity: 0.0, + }; + test_world .create_entity() .with(LaserIndex { index: 0, initiated: true, }) - .with(GaussianBeam { - direction: Vector3::new(1.0, 0.0, 0.0), - intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 2.0, - power: 1.0, - rayleigh_range: gaussian::calculate_rayleigh_range(&461.0e-9, &2.0), - ellipticity: 0.0, - }) + .with(beam) .build(); let atom1 = test_world .create_entity() - .with(Position { pos: Vector3::y() }) + .with(Position { pos: Vector3::new( 0.0, 1.0, 0.0 ) }) .with(LaserIntensitySamplers { contents: [LaserIntensitySampler::default(); crate::laser::DEFAULT_BEAM_LIMIT], }) .build(); let mut system = SampleLaserIntensitySystem::<{ DEFAULT_BEAM_LIMIT }>; + system.run_now(&test_world); + test_world.maintain(); - let sampler_storage = - test_world.read_storage::>(); + + let sampler_storage = test_world.read_storage::>(); let actual_intensity = gaussian::get_gaussian_beam_intensity( &GaussianBeam { - direction: Vector3::new(1.0, 0.0, 0.0), + direction: Vector3::x(), intersection: Vector3::new(0.0, 0.0, 0.0), e_radius: 2.0, power: 1.0, @@ -196,6 +202,33 @@ pub mod tests { None, None, ); + // + // println!("rayleigh = {}", beam.rayleigh_range); + // + // println!( "sampler_storage = {},", + // sampler_storage + // .get(atom1) + // .expect("entity not found") + // .contents[0] + // .intensity + // ); + // + // println!( "actual_intensity = {},", + // actual_intensity + // ); + // + // assert_approx_eq!( + // sampler_storage + // .get(atom1) + // .expect("entity not found") + // .contents[0] + // .intensity, + // actual_intensity, + // 1e-6_f64 + // ); + // + // assert_approx_eq!(1.0_f64, 2.0_f64); + assert_approx_eq!( sampler_storage diff --git a/src/laser/intensity_gradient.rs b/src/laser/intensity_gradient.rs index 0e2408fe..75f712c3 100644 --- a/src/laser/intensity_gradient.rs +++ b/src/laser/intensity_gradient.rs @@ -185,18 +185,22 @@ pub mod tests { test_world.register::(); test_world.register::(); + let e_rad = 70.71067812e-6; + let beam = GaussianBeam { direction: Vector3::x(), intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 70.71067812e-6, + e_radius: e_rad, power: 100.0, rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range( &1064.0e-9, - &70.71067812e-6, + &e_rad, ), ellipticity: 0.0, }; + // println!("rayleigh_range: {}", beam.rayleigh_range); + test_world .create_entity() .with(LaserIndex { @@ -205,8 +209,8 @@ pub mod tests { }) .with(beam) .with(Frame { - x_vector: Vector3::y(), - y_vector: Vector3::z(), + x_vector: Vector3::z(), + y_vector: Vector3::y(), }) .with(DipoleLight { wavelength: 1064.0e-9, @@ -223,6 +227,9 @@ pub mod tests { crate::laser::DEFAULT_BEAM_LIMIT], }) .build(); + + + let mut system = SampleGaussianLaserIntensityGradientSystem::<{ DEFAULT_BEAM_LIMIT }>; system.run_now(&test_world); test_world.maintain(); @@ -234,6 +241,9 @@ pub mod tests { .contents[0] .gradient; + println!("the sin_resutl_gradient = {}",sim_result_gradient ); + + assert_approx_eq!( -2.09081e+8, sim_result_gradient[0], 1e+5_f64); assert_approx_eq!(-4.33993e+13, sim_result_gradient[1], 1e+8_f64); assert_approx_eq!(-4.33993e+13, sim_result_gradient[2], 1e+8_f64); From 085800a035c815d190e68410701fd90846c14665 Mon Sep 17 00:00:00 2001 From: Brian Bostwick Date: Wed, 14 Sep 2022 10:52:05 +0100 Subject: [PATCH 12/19] Fixed issue with trap freq and unit tests --- examples/dipole_trap.rs | 2 +- src/dipole/force.rs | 23 +++++++++++------------ src/laser/gaussian.rs | 18 +++++++++--------- src/laser/intensity.rs | 26 -------------------------- src/laser/intensity_gradient.rs | 7 ------- 5 files changed, 21 insertions(+), 55 deletions(-) diff --git a/examples/dipole_trap.rs b/examples/dipole_trap.rs index 471034e3..4122980d 100644 --- a/examples/dipole_trap.rs +++ b/examples/dipole_trap.rs @@ -30,7 +30,7 @@ fn main() { // Create dipole laser. let power = 1.0; - let e_radius = 60.0e-6 / 2.0_f64.sqrt(); + let e_radius = 80.0e-6 / 2.0_f64.sqrt(); let wavelength = 1064.0e-9; let gaussian_beam = GaussianBeam { diff --git a/src/dipole/force.rs b/src/dipole/force.rs index ffa92386..4900fa8c 100644 --- a/src/dipole/force.rs +++ b/src/dipole/force.rs @@ -150,9 +150,11 @@ pub mod tests { let sampler_storage = test_world.read_storage::(); let sim_result_force = sampler_storage.get(atom1).expect("Entity not found!").force; - assert_approx_eq!(-6.386888332902177e-29, sim_result_force[0], 3e-30_f64); - assert_approx_eq!(-3.11151847e-23, sim_result_force[1], 2e-24_f64); - assert_approx_eq!(-3.11151847e-23, sim_result_force[2], 2e-24_f64); + println!("force = {}", sim_result_force); + + assert_approx_eq!(-1.579e-26, sim_result_force[0], 3e-30_f64); + assert_approx_eq!(-8.097e-21, sim_result_force[1], 2e-24_f64); + assert_approx_eq!(-8.097e-21, sim_result_force[2], 2e-24_f64); } #[test] @@ -257,24 +259,21 @@ pub mod tests { .contents; println!("force is: {}", sim_result_force); - println!("gradient 1 is: {}", _sim_result_grad[0].gradient); - - assert_approx_eq!( - 3.17164e-32, + 3.17e-32, sim_result_force[0], - 3e-46_f64 + 5e-33_f64 ); assert_approx_eq!( - 3.17164e-32, + 3.17e-32, sim_result_force[1], - 2e-46_f64 + 5e-33_f64 ); assert_approx_eq!( - 2.2692e-31, + 1.28e-31, sim_result_force[2], - 2e-46_f64 + 5e-33_f64 ); } } diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index 7ca6acd0..7515d815 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -68,8 +68,8 @@ impl GaussianBeam { peak_intensity: f64, e_radius: f64, ) -> Self { - let std = e_radius * 2.0_f64.powf(0.5); - let power = std::f64::consts::PI * std.powi(2) * peak_intensity; + let std = e_radius / 2.0_f64.powf(0.5); + let power = 2.0 * std::f64::consts::PI * std.powi(2) * peak_intensity; GaussianBeam { intersection, direction, @@ -102,8 +102,8 @@ impl GaussianBeam { e_radius: f64, wavelength: f64, ) -> Self { - let std = e_radius * 2.0_f64.powf(0.5); - let power = std::f64::consts::PI * std.powi(2) * peak_intensity; + let std = e_radius / 2.0_f64.powf(0.5); + let power = 2.0 * std::f64::consts::PI * std.powi(2) * peak_intensity; GaussianBeam { intersection, direction, @@ -210,7 +210,7 @@ pub fn get_gaussian_beam_intensity( let spot_size_squared = 2.0 * beam.e_radius.powf(2.0) * (1. + (z / beam.rayleigh_range).powf(2.0)); - 2.0 * power / ( PI * spot_size_squared ) * EXP.powf(-2.0 * distance_squared / spot_size_squared) + 2.0 * power / ( PI * spot_size_squared ) * EXP.powf(-2.0 * distance_squared / spot_size_squared) } /// Computes the rayleigh range for a given beam and wavelength @@ -231,24 +231,24 @@ pub fn get_gaussian_beam_intensity_gradient( // ellipticity treatment let semi_major_axis = 1.0 / (1.0 - beam.ellipticity.powf(2.0)).powf(0.5); - let x = rela_coord.dot(&reference_frame.x_vector) / semi_major_axis.powf(0.5); let y = rela_coord.dot(&reference_frame.y_vector) * semi_major_axis.powf(0.5); let z = rela_coord.dot(&beam.direction); + //Calculate the spot_size_squared for convenience of later calculations. let spot_size_squared = 2.0 * beam.e_radius.powf(2.0) * (1. + (z / beam.rayleigh_range).powf(2.0)); - - // this will match mathematica if the ebam directionis mult py 2 + //Calculating the gradient vector for grad(I(r,z)). let vector = - 4.0 * ( reference_frame.x_vector * x + reference_frame.y_vector * y ) + 2.0 * beam.direction * ( z / ( beam.rayleigh_range.powf(2.0) + z.powf(2.0) ) ) * ( 2.0 * (x.powf(2.0) + y.powf(2.0)) - spot_size_squared); - + //Calculate the intensity at pos (x, y, z). let intensity = 2.0 * beam.power / ( PI * spot_size_squared ) * EXP.powf(-2.0 * (x.powf(2.0) + y.powf(2.0)) / spot_size_squared); + //Returns the full grad(I) with correct magnitude and direction. vector * intensity / spot_size_squared } diff --git a/src/laser/intensity.rs b/src/laser/intensity.rs index ae6dc46a..16928549 100644 --- a/src/laser/intensity.rs +++ b/src/laser/intensity.rs @@ -202,32 +202,6 @@ pub mod tests { None, None, ); - // - // println!("rayleigh = {}", beam.rayleigh_range); - // - // println!( "sampler_storage = {},", - // sampler_storage - // .get(atom1) - // .expect("entity not found") - // .contents[0] - // .intensity - // ); - // - // println!( "actual_intensity = {},", - // actual_intensity - // ); - // - // assert_approx_eq!( - // sampler_storage - // .get(atom1) - // .expect("entity not found") - // .contents[0] - // .intensity, - // actual_intensity, - // 1e-6_f64 - // ); - // - // assert_approx_eq!(1.0_f64, 2.0_f64); assert_approx_eq!( diff --git a/src/laser/intensity_gradient.rs b/src/laser/intensity_gradient.rs index 75f712c3..9d5a0a56 100644 --- a/src/laser/intensity_gradient.rs +++ b/src/laser/intensity_gradient.rs @@ -199,8 +199,6 @@ pub mod tests { ellipticity: 0.0, }; - // println!("rayleigh_range: {}", beam.rayleigh_range); - test_world .create_entity() .with(LaserIndex { @@ -228,8 +226,6 @@ pub mod tests { }) .build(); - - let mut system = SampleGaussianLaserIntensityGradientSystem::<{ DEFAULT_BEAM_LIMIT }>; system.run_now(&test_world); test_world.maintain(); @@ -241,9 +237,6 @@ pub mod tests { .contents[0] .gradient; - println!("the sin_resutl_gradient = {}",sim_result_gradient ); - - assert_approx_eq!( -2.09081e+8, sim_result_gradient[0], 1e+5_f64); assert_approx_eq!(-4.33993e+13, sim_result_gradient[1], 1e+8_f64); assert_approx_eq!(-4.33993e+13, sim_result_gradient[2], 1e+8_f64); From fb90e645410613452924c76c3c33e691344179bb Mon Sep 17 00:00:00 2001 From: Brian Bostwick Date: Wed, 14 Sep 2022 10:55:48 +0100 Subject: [PATCH 13/19] removed debuging print statments from unit tests --- src/dipole/force.rs | 6 ------ src/laser/gaussian.rs | 2 +- 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/src/dipole/force.rs b/src/dipole/force.rs index 4900fa8c..7d1d422b 100644 --- a/src/dipole/force.rs +++ b/src/dipole/force.rs @@ -150,8 +150,6 @@ pub mod tests { let sampler_storage = test_world.read_storage::(); let sim_result_force = sampler_storage.get(atom1).expect("Entity not found!").force; - println!("force = {}", sim_result_force); - assert_approx_eq!(-1.579e-26, sim_result_force[0], 3e-30_f64); assert_approx_eq!(-8.097e-21, sim_result_force[1], 2e-24_f64); assert_approx_eq!(-8.097e-21, sim_result_force[2], 2e-24_f64); @@ -182,8 +180,6 @@ pub mod tests { ellipticity: 0.0, }; - println!("rayleigh_range {}", gaussian_beam.rayleigh_range); - test_world .create_entity() .with(gaussian_beam) @@ -258,8 +254,6 @@ pub mod tests { .expect("Entity not found!") .contents; - println!("force is: {}", sim_result_force); - assert_approx_eq!( 3.17e-32, sim_result_force[0], diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index 7515d815..7f844c8c 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -284,7 +284,7 @@ pub mod tests { }; let gradient = get_gaussian_beam_intensity_gradient(&beam, &pos1, &grf); - println!("{}",gradient[0]); + assert_approx_eq!(gradient[0], -4.91021e+10, 1e+9_f64); assert_approx_eq!(gradient[1], -5.89225e+11, 1e+9_f64); assert_approx_eq!(gradient[2], -9.38334e+6, 1e+9_f64); From 22fe48da4b8fab0ba1429da0b106268d77fcdfbc Mon Sep 17 00:00:00 2001 From: Brian Bostwick Date: Wed, 14 Sep 2022 11:01:16 +0100 Subject: [PATCH 14/19] removed debuging print statments from unit tests --- src/dipole/force.rs | 6 ------ src/laser/gaussian.rs | 2 +- 2 files changed, 1 insertion(+), 7 deletions(-) diff --git a/src/dipole/force.rs b/src/dipole/force.rs index 4900fa8c..7d1d422b 100644 --- a/src/dipole/force.rs +++ b/src/dipole/force.rs @@ -150,8 +150,6 @@ pub mod tests { let sampler_storage = test_world.read_storage::(); let sim_result_force = sampler_storage.get(atom1).expect("Entity not found!").force; - println!("force = {}", sim_result_force); - assert_approx_eq!(-1.579e-26, sim_result_force[0], 3e-30_f64); assert_approx_eq!(-8.097e-21, sim_result_force[1], 2e-24_f64); assert_approx_eq!(-8.097e-21, sim_result_force[2], 2e-24_f64); @@ -182,8 +180,6 @@ pub mod tests { ellipticity: 0.0, }; - println!("rayleigh_range {}", gaussian_beam.rayleigh_range); - test_world .create_entity() .with(gaussian_beam) @@ -258,8 +254,6 @@ pub mod tests { .expect("Entity not found!") .contents; - println!("force is: {}", sim_result_force); - assert_approx_eq!( 3.17e-32, sim_result_force[0], diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index 7515d815..7f844c8c 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -284,7 +284,7 @@ pub mod tests { }; let gradient = get_gaussian_beam_intensity_gradient(&beam, &pos1, &grf); - println!("{}",gradient[0]); + assert_approx_eq!(gradient[0], -4.91021e+10, 1e+9_f64); assert_approx_eq!(gradient[1], -5.89225e+11, 1e+9_f64); assert_approx_eq!(gradient[2], -9.38334e+6, 1e+9_f64); From 703c78172f641ff2dbd18c337b2260de1d9bb8ee Mon Sep 17 00:00:00 2001 From: Brian Bostwick Date: Thu, 15 Sep 2022 16:50:44 +0100 Subject: [PATCH 15/19] editing example --- examples/dipole_trap.rs | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/examples/dipole_trap.rs b/examples/dipole_trap.rs index 4122980d..57e71d29 100644 --- a/examples/dipole_trap.rs +++ b/examples/dipole_trap.rs @@ -29,15 +29,15 @@ fn main() { let mut sim = sim_builder.build(); // Create dipole laser. - let power = 1.0; - let e_radius = 80.0e-6 / 2.0_f64.sqrt(); + let power = 100.0; + let e_radius = 60.0e-6 / 2.0_f64.sqrt(); let wavelength = 1064.0e-9; let gaussian_beam = GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), e_radius, power, - direction: Vector3::x(), + direction: Vector3::z(), rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range(&wavelength, &e_radius), ellipticity: 0.0, }; @@ -46,8 +46,8 @@ fn main() { .with(gaussian_beam) .with(dipole::DipoleLight { wavelength }) .with(laser::frame::Frame { - x_vector: Vector3::y(), - y_vector: Vector3::z(), + x_vector: Vector3::x(), + y_vector: Vector3::y(), }) .build(); @@ -70,7 +70,7 @@ fn main() { .build(); // Define timestep - sim.world.insert(Timestep { delta: 1.0e-6 }); + sim.world.insert(Timestep { delta: 1.0e-7 }); // Run the simulation for a number of steps. for _i in 0..200_000 { From 5888d4fc1ce1cfa717dc242c821670ca43aa5d68 Mon Sep 17 00:00:00 2001 From: BrianBostwick <42308966+BrianBostwick@users.noreply.github.com> Date: Fri, 16 Sep 2022 11:31:44 +0100 Subject: [PATCH 16/19] Delete .DS_Store --- .DS_Store | Bin 6148 -> 0 bytes 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 .DS_Store diff --git a/.DS_Store b/.DS_Store deleted file mode 100644 index 5f246ddb48f3aac0e500c7205417caf653dee7c1..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHLze@u#6#illJzEgRg3F=OQE+xSipLn7 zY&SbuIxIx`4zGY$;6EzB?{0uCoZu8A@&49ZH?A+ulB7Mz+SC`;HkUW5htIpx-<|0n z_MN?r{;4GAN62dF-$8`8an*XEs|{Gw5rtfLFjP@T~yv4-pl{*kW!_Zyi+WD*$0Y-P)+F z_e#QuI|yToxk2h7Qzn(rq$+#GP$r%FBO4c6%nh1!C_OWcV`o?T75D;5`nkCP From 24500ec1f0f9b842242868ebff3a5da656fa5701 Mon Sep 17 00:00:00 2001 From: BrianBostwick <42308966+BrianBostwick@users.noreply.github.com> Date: Fri, 16 Sep 2022 11:33:00 +0100 Subject: [PATCH 17/19] adding */.ds_store --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 99be97d9..a9fd4aa9 100644 --- a/.gitignore +++ b/.gitignore @@ -16,3 +16,4 @@ *.blend *.blend1 .idea/ +*/.ds_store From 8505b831c6207e94cdd738894656389b7788cce1 Mon Sep 17 00:00:00 2001 From: Brian Bostwick Date: Fri, 16 Sep 2022 18:00:46 +0100 Subject: [PATCH 18/19] redundant factor of 2 --- src/laser/gaussian.rs | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index 7f844c8c..b08a4b0f 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -68,8 +68,9 @@ impl GaussianBeam { peak_intensity: f64, e_radius: f64, ) -> Self { - let std = e_radius / 2.0_f64.powf(0.5); - let power = 2.0 * std::f64::consts::PI * std.powi(2) * peak_intensity; + + let power = std::f64::consts::PI * e_radius.powi(2) * peak_intensity; + GaussianBeam { intersection, direction, @@ -102,8 +103,9 @@ impl GaussianBeam { e_radius: f64, wavelength: f64, ) -> Self { - let std = e_radius / 2.0_f64.powf(0.5); - let power = 2.0 * std::f64::consts::PI * std.powi(2) * peak_intensity; + + let power = std::f64::consts::PI * e_radius.powi(2) * peak_intensity; + GaussianBeam { intersection, direction, From 8d99d5b0aa04c3383a2641e3a19695cae760478d Mon Sep 17 00:00:00 2001 From: Brian Bostwick Date: Tue, 20 Sep 2022 15:36:16 +0100 Subject: [PATCH 19/19] adding printing for intensity --- src/laser/intensity.rs | 15 +++++++++++++++ src/output/file.rs | 12 ++++++------ 2 files changed, 21 insertions(+), 6 deletions(-) diff --git a/src/laser/intensity.rs b/src/laser/intensity.rs index 16928549..6512e74a 100644 --- a/src/laser/intensity.rs +++ b/src/laser/intensity.rs @@ -12,6 +12,7 @@ use crate::atom::Position; use crate::laser::index::LaserIndex; use serde::Serialize; use specs::prelude::*; +use std::fmt; const LASER_CACHE_SIZE: usize = 16; @@ -31,6 +32,14 @@ impl Default for LaserIntensitySampler { } } +impl Component for LaserIntensitySampler { + type Storage = VecStorage; +} +impl fmt::Display for LaserIntensitySampler { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "{:?}", self.intensity) +}} + /// Component that holds a list of `LaserIntensitySamplers` #[derive(Copy, Clone, Serialize)] pub struct LaserIntensitySamplers { @@ -43,6 +52,12 @@ impl Component for LaserIntensitySamplers { type Storage = VecStorage; } +impl fmt::Display for LaserIntensitySamplers { + fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result { + write!(f, "({}, {}, 0.0)", self.contents[0], self.contents[1]) + } +} + /// This system initialises all `LaserIntensitySamplers` to a NAN value. /// /// It also ensures that the size of the `LaserIntensitySamplers` components match the number of CoolingLight entities in the world. diff --git a/src/output/file.rs b/src/output/file.rs index e776c1d2..a2fc601c 100644 --- a/src/output/file.rs +++ b/src/output/file.rs @@ -40,11 +40,11 @@ pub struct FileOutputPlugin phantom_f: PhantomData, phantom_a: PhantomData } -impl FileOutputPlugin - where +impl FileOutputPlugin + where C: Component + Clone, A: Component, - F: Format> + F: Format> { pub fn new(file_name: String, interval: u64) -> FileOutputPlugin { @@ -53,13 +53,13 @@ impl FileOutputPlugin interval, phantom_a: PhantomData, phantom_c: PhantomData, - phantom_f: PhantomData + phantom_f: PhantomData } } } -impl Plugin for FileOutputPlugin -where +impl Plugin for FileOutputPlugin +where C: Component + Clone + Sync + Send + 'static, A: Component + Sync + Send + 'static, F: Format> + Sync + Send + 'static