Fixed trap frequency and updated unit tests

.gitignore

@@ -16,3 +16,4 @@
 *.blend
 *.blend1
 .idea/
+*/.ds_store

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,39 +23,21 @@ 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::<Position, Text, Atom>::new("pos.txt".to_string(), 100));
-    sim_builder.add_plugin(FileOutputPlugin::<Velocity, Text, Atom>::new("vel.txt".to_string(), 100));
-    sim_builder.add_plugin(FileOutputPlugin::<Position, XYZ, Atom>::new("position.xyz".to_string(), 100));
+    sim_builder.add_plugin(FileOutputPlugin::<Position, Text, Atom>::new("pos.txt".to_string(), 1));
+    sim_builder.add_plugin(FileOutputPlugin::<Velocity, Text, Atom>::new("vel.txt".to_string(), 1));
+    sim_builder.add_plugin(FileOutputPlugin::<Position, XYZ, Atom>::new("position.xyz".to_string(), 1));
     let mut sim = sim_builder.build();
 
     // Create dipole laser.
-    let power = 10.0;
-    let e_radius = 60.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(),
-        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::y(),
-            y_vector: Vector3::z(),
-        })
-        .build();
-
-    let gaussian_beam = GaussianBeam {
-        intersection: Vector3::new(0.0, 0.0, 0.0),
-        e_radius,
-        power,
-        direction: Vector3::y(),
+        direction: Vector3::z(),
         rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range(&wavelength, &e_radius),
         ellipticity: 0.0,
     };
@@ -65,20 +47,17 @@ fn main() {
         .with(dipole::DipoleLight { wavelength })
         .with(laser::frame::Frame {
             x_vector: Vector3::x(),
-            y_vector: Vector3::z(),
+            y_vector: Vector3::y(),
         })
         .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, 1.0e-7, 1.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-7 });
+
     // Run the simulation for a number of steps.
-    for _i in 0..100_000 {
+    for _i in 0..200_000 {
         sim.step();
     }
 

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;

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
@@ -22,7 +22,10 @@ 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;
+
+/// Vacuum permittivity [F/m]
+pub const EPSILON0: f64 = 8.8541878128e-12;

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,10 @@ impl<'a, const N: usize> System<'a> for ApplyDipoleForceSystem<N> {
     ) {
         (&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 += polarizability.prefactor*sampler.contents[index.index].gradient/
+                                   ( 2.0 * constant::EPSILON0 * constant::C );
                 }
             });
     }
@@ -50,6 +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;
 
     #[test]
     fn test_apply_dipole_force_system() {
@@ -110,46 +114,46 @@ pub mod tests {
     fn test_apply_dipole_force_again_system() {
         let mut test_world = World::new();
 
-        test_world.register::<LaserIndex>();
-        test_world.register::<DipoleLight>();
-        test_world.register::<Force>();
-        test_world.register::<LaserIntensityGradientSamplers<{ DEFAULT_BEAM_LIMIT }>>();
-        test_world.register::<Polarizability>();
+         test_world.register::<LaserIndex>();
+         test_world.register::<DipoleLight>();
+         test_world.register::<Force>();
+         test_world.register::<LaserIntensityGradientSamplers<{ DEFAULT_BEAM_LIMIT }>>();
+         test_world.register::<Polarizability>();
 
-        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::<Force>();
-        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::<Force>();
+         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!(-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]
     fn test_apply_dipole_force_and_gradient_system() {
@@ -175,6 +179,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)
@@ -190,6 +195,7 @@ pub mod tests {
                 y_vector: Vector3::z(),
             })
             .build();
+
         let gaussian_beam = GaussianBeam {
             intersection: Vector3::new(0.0, 0.0, 0.0),
             e_radius,
@@ -198,6 +204,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)
@@ -215,6 +222,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 {
@@ -245,24 +253,21 @@ 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);
 
         assert_approx_eq!(
-            0.000000000000000000000000000000000127913190642808,
+            3.17e-32,
             sim_result_force[0],
-            3e-46_f64
+            5e-33_f64
         );
         assert_approx_eq!(
-            0.000000000000000000000000000000000127913190642808,
+            3.17e-32,
             sim_result_force[1],
-            2e-46_f64
+            5e-33_f64
         );
         assert_approx_eq!(
-            0.000000000000000000000000000000000511875188257342,
+            1.28e-31,
             sim_result_force[2],
-            2e-46_f64
+            5e-33_f64
         );
     }
 }

src/dipole/mod.rs

@@ -46,7 +46,7 @@ impl Component for Polarizability {
     type Storage = VecStorage<Self>;
 }
 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.
@@ -57,10 +57,11 @@ impl Polarizability {
     ) -> 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))
-            * optical_transition_linewidth
-            * -(1. / (transition_f - dipole_f) + 1. / (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)
+            );
         Polarizability { prefactor }
     }
 }
@@ -83,11 +84,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<const N : usize>;
 impl<const N: usize> Plugin for DipolePlugin<N> {