diff --git a/examples/1d_mot.rs b/examples/1d_mot.rs index ca052b21..74d1e16b 100644 --- a/examples/1d_mot.rs +++ b/examples/1d_mot.rs @@ -43,11 +43,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 0.01, + w0_x: 0.01, + w0_y: 0.01, power, direction: -Vector3::z(), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, @@ -58,11 +59,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 0.01, + w0_x: 0.01, + w0_y: 0.01, power, direction: Vector3::z(), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, diff --git a/examples/2d_plus_mot_from_oven.rs b/examples/2d_plus_mot_from_oven.rs index d07374e8..c50c703c 100644 --- a/examples/2d_plus_mot_from_oven.rs +++ b/examples/2d_plus_mot_from_oven.rs @@ -53,11 +53,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: push_beam_radius, + w0_x: push_beam_radius, + w0_y: push_beam_radius, power: push_beam_power, direction: Vector3::z(), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( push_beam_detuning, @@ -73,11 +74,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: radius, + w0_x: radius, + w0_y: radius, power, direction: Vector3::new(1.0, 1.0, 0.0).normalize(), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, @@ -88,11 +90,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: radius, + w0_x: radius, + w0_y: radius, power, direction: Vector3::new(1.0, -1.0, 0.0).normalize(), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, @@ -103,11 +106,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: radius, + w0_x: radius, + w0_y: radius, power, direction: Vector3::new(-1.0, 1.0, 0.0).normalize(), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, @@ -118,11 +122,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: radius, + w0_x: radius, + w0_y: radius, power, direction: Vector3::new(-1.0, -1.0, 0.0).normalize(), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, diff --git a/examples/benchmark.rs b/examples/benchmark.rs index f74d73cb..d623ccba 100644 --- a/examples/benchmark.rs +++ b/examples/benchmark.rs @@ -86,11 +86,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: beam_centre, - e_radius: radius, + w0_x: radius, + w0_y: radius, power, direction: Vector3::new(0.0, 0.0, 1.0), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, @@ -101,11 +102,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: beam_centre, - e_radius: radius, + w0_x: radius, + w0_y: radius, power, direction: Vector3::new(0.0, 0.0, -1.0), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, @@ -116,11 +118,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: beam_centre, - e_radius: radius, + w0_x: radius, + w0_y: radius, power, direction: Vector3::new(-1.0, 0.0, 0.0), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, @@ -131,11 +134,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: beam_centre, - e_radius: radius, + w0_x: radius, + w0_y: radius, power, direction: Vector3::new(1.0, 0.0, 0.0), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, @@ -146,11 +150,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: beam_centre, - e_radius: radius, + w0_x: radius, + w0_y: radius, power, direction: Vector3::new(0.0, 1.0, 0.0), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, @@ -161,11 +166,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: beam_centre, - e_radius: radius, + w0_x: radius, + w0_y: radius, power, direction: Vector3::new(0.0, -1.0, 0.0), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, diff --git a/examples/cross_section.rs b/examples/cross_section.rs index a6001ea3..608e2aae 100644 --- a/examples/cross_section.rs +++ b/examples/cross_section.rs @@ -41,11 +41,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: radius, + w0_x: radius, + w0_y: radius, power, direction: Vector3::x(), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, diff --git a/examples/dipole_trap.rs b/examples/dipole_trap.rs index 2d615906..2b44f53f 100644 --- a/examples/dipole_trap.rs +++ b/examples/dipole_trap.rs @@ -30,16 +30,18 @@ fn main() { // Create dipole laser. let power = 10.0; - let e_radius = 60.0e-6 / (2.0_f64.sqrt()); + let w0_x = 60.0e-6; + let w0_y = 60.0e-6; let wavelength = 1064.0e-9; let gaussian_beam = GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius, + w0_x, + w0_y, power, direction: Vector3::x(), - rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range(&wavelength, &e_radius), - ellipticity: 0.0, + rayleigh_range_x: crate::laser::gaussian::calculate_rayleigh_range(&wavelength, &w0_x), + rayleigh_range_y: crate::laser::gaussian::calculate_rayleigh_range(&wavelength, &w0_y), }; sim.world .create_entity() @@ -53,11 +55,12 @@ fn main() { let gaussian_beam = GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius, + w0_x, + w0_y, power, direction: Vector3::y(), - rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range(&wavelength, &e_radius), - ellipticity: 0.0, + rayleigh_range_x: crate::laser::gaussian::calculate_rayleigh_range(&wavelength, &w0_x), + rayleigh_range_y: crate::laser::gaussian::calculate_rayleigh_range(&wavelength, &w0_y), }; sim.world .create_entity() diff --git a/examples/doppler_limit.rs b/examples/doppler_limit.rs index 10a2d34f..3a711c4d 100644 --- a/examples/doppler_limit.rs +++ b/examples/doppler_limit.rs @@ -79,18 +79,19 @@ fn main() { // Create cooling lasers. let detuning = configuration.detuning; let power = 0.02; - let radius = 66.7e-3 / (2.0_f64.sqrt()); + let radius = 66.7e-3 ; let beam_centre = Vector3::new(0.0, 0.0, 0.0); sim.world .create_entity() .with(GaussianBeam { intersection: beam_centre, - e_radius: radius, + w0_x: radius, + w0_y: radius, power, direction: Vector3::new(0.0, 0.0, 1.0), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, @@ -101,11 +102,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: beam_centre, - e_radius: radius, + w0_x: radius, + w0_y: radius, power, direction: Vector3::new(0.0, 0.0, -1.0), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, @@ -116,11 +118,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: beam_centre, - e_radius: radius, + w0_x: radius, + w0_y: radius, power, - direction: Vector3::new(-1.0, 0.0, 0.0), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + direction: Vector3::new(-1.0, 0.0, 1.0), + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, @@ -131,11 +134,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: beam_centre, - e_radius: radius, + w0_x: radius, + w0_y: radius, power, - direction: Vector3::new(1.0, 0.0, 0.0), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + direction: Vector3::new(1.0, 0.0, 1.0), + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, @@ -146,11 +150,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: beam_centre, - e_radius: radius, + w0_x: radius, + w0_y: radius, power, direction: Vector3::new(0.0, 1.0, 0.0), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, @@ -161,11 +166,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: beam_centre, - e_radius: radius, + w0_x: radius, + w0_y: radius, power, direction: Vector3::new(0.0, -1.0, 0.0), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( detuning, diff --git a/examples/molasses_1d.rs b/examples/molasses_1d.rs index 89fdb1b7..ba761e28 100644 --- a/examples/molasses_1d.rs +++ b/examples/molasses_1d.rs @@ -48,11 +48,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 0.01, + w0_x: 0.01, + w0_y: 0.01, power: 0.01, direction: -Vector3::z(), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( -6.0, @@ -63,11 +64,12 @@ fn main() { .create_entity() .with(GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 0.01, + w0_x: 0.01, + w0_y: 0.01, power: 0.01, direction: Vector3::z(), - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, }) .with(CoolingLight::for_transition::( -6.0, diff --git a/src/atom_sources/emit.rs b/src/atom_sources/emit.rs index 79b671bc..92fdec2a 100644 --- a/src/atom_sources/emit.rs +++ b/src/atom_sources/emit.rs @@ -49,7 +49,7 @@ impl<'a> System<'a> for EmitNumberPerFrameSystem { ); fn run(&mut self, (emit_numbers, mut numbers_to_emit): Self::SystemData) { - for (emit_number, mut number_to_emit) in (&emit_numbers, &mut numbers_to_emit).join() { + for (emit_number, mut number_to_emit) in (&emit_numbers, &mut numbers_to_emit).join() { number_to_emit.number = emit_number.number; } } diff --git a/src/dipole/force.rs b/src/dipole/force.rs index 233c0ca2..b3b7e874 100644 --- a/src/dipole/force.rs +++ b/src/dipole/force.rs @@ -22,6 +22,9 @@ impl<'a, const N: usize> System<'a> for ApplyDipoleForceSystem { WriteStorage<'a, Force>, ); + /// Here we are applying dipole force, it extract polarizability factor and gassuain intensity gradient, + /// their product is the dipole force + fn run( &mut self, (dipole_light, dipole_index, polarizability, gradient_sampler, mut force): Self::SystemData, @@ -165,15 +168,17 @@ pub mod tests { test_world.register::(); let power = 10.0; - let e_radius = 60.0e-6 / (2.0_f64.sqrt()); + let w0_x = 50.0e-6; + let w0_y = 60.0e-6; 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(&1064.0e-9, &e_radius), - ellipticity: 0.0, + w0_x: w0_x, + w0_y: w0_y, + power: power, + direction: Vector3::z(), + rayleigh_range_x: crate::laser::gaussian::calculate_rayleigh_range(&1064.0e-9, &w0_x), + rayleigh_range_y: crate::laser::gaussian::calculate_rayleigh_range(&1064.0e-9, &w0_y), }; test_world .create_entity() @@ -186,17 +191,18 @@ pub mod tests { initiated: true, }) .with(laser::frame::Frame { - x_vector: Vector3::y(), - y_vector: Vector3::z(), + x_vector: Vector3::x(), + y_vector: Vector3::y(), }) .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(&1064.0e-9, &e_radius), - ellipticity: 0.0, + w0_x: w0_x, + w0_y: w0_y, + power: power, + direction: Vector3::x(), + rayleigh_range_x: crate::laser::gaussian::calculate_rayleigh_range(&1064.0e-9, &w0_x), + rayleigh_range_y: crate::laser::gaussian::calculate_rayleigh_range(&1064.0e-9, &w0_y), }; test_world .create_entity() @@ -209,7 +215,7 @@ pub mod tests { initiated: true, }) .with(laser::frame::Frame { - x_vector: Vector3::x(), + x_vector: Vector3::y(), y_vector: Vector3::z(), }) .build(); @@ -241,28 +247,37 @@ pub mod tests { let grad_sampler_storage = test_world.read_storage::>(); let sim_result_force = sampler_storage.get(atom1).expect("Entity not found!").force; - let _sim_result_grad = grad_sampler_storage + let sim_result_grad = grad_sampler_storage .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[0].gradient); + println!("gradient 2 is: {}", sim_result_grad[1].gradient); + + let pol = Polarizability::calculate_for( + 1064e-9, // dipole wavelength + 461e-9, // transition wavelength + 32e6, // linewidth (Hz) + ); + + println!("Polarizability prefactor = {}", pol.prefactor); assert_approx_eq!( - 0.000000000000000000000000000000000127913190642808, + 3.3243137806595777e-28, sim_result_force[0], - 3e-46_f64 + 3e-34_f64 ); assert_approx_eq!( - 0.000000000000000000000000000000000127913190642808, + 2.3094285637507477e-28, sim_result_force[1], - 2e-46_f64 + 2e-34_f64 ); assert_approx_eq!( - 0.000000000000000000000000000000000511875188257342, + -1.5146534391703284e-31, sim_result_force[2], - 2e-46_f64 + 2e-37_f64 ); } } + diff --git a/src/integration_tests/rate_equation.rs b/src/integration_tests/rate_equation.rs index edcf481b..1b3a7dad 100644 --- a/src/integration_tests/rate_equation.rs +++ b/src/integration_tests/rate_equation.rs @@ -61,11 +61,12 @@ pub mod tests { 1, )) .with(LaserIndex::default()) - .with(GaussianBeam::from_peak_intensity_with_rayleigh_range( + .with(GaussianBeam::from_peak_intensity_ellipticity_with_rayleigh_range( Vector3::new(0.0, 0.0, 0.0), Vector3::new(-1.0, 0.0, 0.0), intensity, 0.01, + 0.01, 780.0e-9, )) .build(); diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index 4da1d11f..1f2cd3f6 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -26,26 +26,22 @@ pub struct GaussianBeam { /// A point that the laser beam intersects pub intersection: Vector3, - /// Direction the beam propagates with respect to cartesian `x,y,z` axes. + /// Direction the beam propagates, the axis direction of the gaussian beam. + /// the transverse direction is not defined here, they will be defined by frame vectors when extracting + /// gaussian intensity and gaussian intensity gradient pub direction: Vector3, - - /// Radius of the beam at which the intensity is 1/e of the peak value, SI units of m. - /// - /// Since in the literature the e^2_radius (where intensity is 1/e^2 of peak value) is used - /// very often as well, it is useful to note the following relation: - /// - /// e_radius = e^2_radius / sqrt(2) - pub e_radius: f64, - + /// Power of the laser in W pub power: f64, - /// The distance along the propagation direction of a beam from the - /// waist to the place where the area of the cross section is doubled in units of metres - pub rayleigh_range: f64, + /// standard gaussian beam waist (1/e^2 radius) in units of m. + pub w0_x: f64, + pub w0_y: f64, + + /// Rayleigh range of the elliptical beam in units of m. + pub rayleigh_range_x: f64, + pub rayleigh_range_y: f64, - /// ellipticity - pub ellipticity: f64, } impl Component for GaussianBeam { type Storage = HashMapStorage; @@ -61,22 +57,27 @@ impl GaussianBeam { /// /// `peak_intensity`: peak intensity in units of W/m^2. /// - /// `e_radius`: radius of beam in units of m. + /// `w0_x, w0_y`: radius of beam in units of m. + + /// Here we are defining and unphysical beam, kinda "collimated" beam, + /// this bit is used in laser cooling modules. pub fn from_peak_intensity( intersection: Vector3, direction: Vector3, peak_intensity: f64, - e_radius: f64, + w0_x: f64, + w0_y: 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 = 1.0/2.0 * std::f64::consts::PI * w0_x * w0_y * peak_intensity; GaussianBeam { intersection, direction, power, - e_radius, - rayleigh_range: f64::INFINITY, - ellipticity: 0.0, + w0_x, + w0_y, + rayleigh_range_x: f64::INFINITY, + rayleigh_range_y: f64::INFINITY, + } } } @@ -92,25 +93,26 @@ impl GaussianBeam { /// /// `peak_intensity`: peak intensity in units of W/m^2. /// - /// `e_radius`: radius of beam in units of m. + /// `w0_x, w0_y`: radius of beam in units of m. /// /// `wavelength`: wavelength of the electromagnetic light - pub fn from_peak_intensity_with_rayleigh_range( + pub fn from_peak_intensity_ellipticity_with_rayleigh_range( intersection: Vector3, direction: Vector3, peak_intensity: f64, - e_radius: f64, + w0_x: f64, + w0_y: 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 * w0_x * w0_y * peak_intensity / 2.0; GaussianBeam { intersection, direction, power, - e_radius, - rayleigh_range: calculate_rayleigh_range(&wavelength, &e_radius), - ellipticity: 0.0, + w0_x, + w0_y, + rayleigh_range_x: calculate_rayleigh_range(&wavelength, &w0_x), + rayleigh_range_y: calculate_rayleigh_range(&wavelength, &w0_y), } } /// Create a GaussianBeam component by specifying the peak intensity, rather than power. @@ -132,17 +134,18 @@ impl GaussianBeam { intersection: Vector3, direction: Vector3, power: f64, - e_radius: f64, + w0_x: f64, + w0_y: f64, wavelength: f64, - ellipiticity: f64, ) -> Self { GaussianBeam { intersection, direction: direction.normalize(), power, - e_radius, - rayleigh_range: calculate_rayleigh_range(&wavelength, &e_radius), - ellipticity: ellipiticity, + w0_x, + w0_y, + rayleigh_range_x: calculate_rayleigh_range(&wavelength, &w0_x), + rayleigh_range_y: calculate_rayleigh_range(&wavelength, &w0_y), } } } @@ -160,60 +163,93 @@ impl Component for CircularMask { } /// Returns the intensity of a gaussian laser beam at the specified position. + pub fn get_gaussian_beam_intensity( beam: &GaussianBeam, pos: &Position, - mask: Option<&CircularMask>, + _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( - &pos.pos, - &beam.intersection, - &beam.direction, - ); - (z, distance * distance) - } - }; - let power = match mask { - Some(mask) => { - if distance_squared.powf(0.5) < mask.radius { - 0.0 - } else { - beam.power - } - } - None => beam.power, + /// Frame vector specify transverses direction of the beam, then convert pos (which is cartesian standard) + /// into beam frame, then fed into standard gaussian intensity formula to extract intensity. + let binding = Frame{ + x_vector: Vector3::x(), + y_vector: Vector3::y(), }; - power / PI / beam.e_radius.powf(2.0) / (1.0 + (z / beam.rayleigh_range).powf(2.0)) + let frame = frame.unwrap_or(&binding); + + let (x, y, z) = maths::get_relative_coordinates_line_point( + &pos.pos, + &beam.intersection, + &beam.direction, + &frame, + ); + + 2.0 * beam.power / PI / beam.w0_x / beam.w0_y / (1.0 + (z / beam.rayleigh_range_x).powf(2.0)).powf(0.5) + / (1.0 + (z / beam.rayleigh_range_y).powf(2.0)).powf(0.5) * EXP.powf( - -distance_squared - / (beam.e_radius.powf(2.0) * (1. + (z / beam.rayleigh_range).powf(2.0))), + -2.0 * x.powf(2.0) / (beam.w0_x.powf(2.0) * (1. + (z / beam.rayleigh_range_x).powf(2.0))) + -2.0 * y.powf(2.0) / (beam.w0_y.powf(2.0) * (1. + (z / beam.rayleigh_range_y).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 +/// Computes the rayleigh range for a given beam and wavelength, w0 is standard 1/e^2 gaussian beam waist. +pub fn calculate_rayleigh_range(wavelength: &f64, w0: &f64) -> f64 { + PI * w0.powf(2.0) / wavelength } + +/// Returns the intensity of a gaussian laser beam at the specified position. +// pub fn get_gaussian_beam_intensity( +// beam: &GaussianBeam, +// pos: &Position, +// 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( +// &pos.pos, +// &beam.intersection, +// &beam.direction, +// ); +// (z, distance * distance) +// } +// }; +// let power = match mask { +// Some(mask) => { +// if distance_squared.powf(0.5) < mask.radius { +// 0.0 +// } else { +// beam.power +// } +// } +// 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))), +// ) +// } + + /// Computes the intensity gradient of a given beam and returns it as /// a three-dimensional vector pub fn get_gaussian_beam_intensity_gradient( @@ -221,24 +257,56 @@ pub fn get_gaussian_beam_intensity_gradient( pos: &Position, reference_frame: &Frame, ) -> Vector3 { + // Relative coordinate from beam intersection 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); + // Ellipticity treatment: waists along long and short axes + let wx = beam.w0_x; // long axis waist squared + let wy = beam.w0_y; // short axis waist squared - 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); + // Rayleigh ranges for each axis (waist squared) + let zx = beam.rayleigh_range_x; + let zy = beam.rayleigh_range_y; // convert for long axis + + // Local coordinates along beam axes + let x = rela_coord.dot(&reference_frame.x_vector); + let y = rela_coord.dot(&reference_frame.y_vector); 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); + // z-dependent spot sizes + let wx_z = wx * ( (1.0 + (z / zx).powi(2))).sqrt(); + let wy_z = wy * ( (1.0 + (z / zy).powi(2))).sqrt(); + + // Intensity prefactor + let intensity_prefactor = 2.0 * beam.power / std::f64::consts::PI / (wx_z * wy_z); + + // Exponential factor + let exp_factor = (-2.0 * (x.powi(2) / wx_z.powi(2) + y.powi(2) / wy_z.powi(2))).exp(); + + // Gradient along x and y + let grad_x = -4.0 * x / wx_z.powi(2) * intensity_prefactor * exp_factor; + let grad_y = -4.0 * y / wy_z.powi(2) * intensity_prefactor * exp_factor; + + // Gradient along z + // Term 1: derivative of prefactor 1/(wx_z wy_z) + let term1 = -intensity_prefactor / (wx_z * wy_z) + * ( + z * wx.powi(2) * wy_z/ (zx.powi(2) * wx_z) + + z * wy.powi(2) * wx_z / (zy.powi(2) * wy_z) + ); + + // Term 2: derivative of exponential + let term2 = intensity_prefactor * 4.0 * ( + x.powi(2) * z * wx.powi(2) / (wx_z.powi(4) * zx.powi(2)) + + y.powi(2) * z * wy.powi(2) / (wy_z.powi(4) * zy.powi(2)) + ); - intensity / spot_size_squared * vector + let grad_z = (term1 + term2) * exp_factor; + + // Combine into Vector3 + reference_frame.x_vector * grad_x + + reference_frame.y_vector * grad_y + + beam.direction * grad_z } #[cfg(test)] @@ -258,10 +326,11 @@ pub mod tests { let beam = GaussianBeam { direction: Vector3::z(), intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 70.71067812e-6, + w0_x: 70.71067812e-6, + w0_y: 70.71067812e-6, power: 100.0, - rayleigh_range: calculate_rayleigh_range(&1064.0e-9, &70.71067812e-6), - ellipticity: 0.0, + rayleigh_range_x: calculate_rayleigh_range(&1064.0e-9, &70.71067812e-6), + rayleigh_range_y: calculate_rayleigh_range(&1064.0e-9, &70.71067812e-6), }; let pos1 = Position { pos: Vector3::new(10.0e-6, 0.0, 30.0e-6), @@ -272,9 +341,13 @@ pub mod tests { }; let gradient = get_gaussian_beam_intensity_gradient(&beam, &pos1, &grf); - assert_approx_eq!(gradient[0], -2.49605032e+13, 1e+8_f64); + + println!("Gradient: {:?}", gradient); + + assert_approx_eq!(gradient[0], -97864416562438.33, 1e+8_f64); assert_approx_eq!(gradient[1], 0.0, 1e+9_f64); - assert_approx_eq!(gradient[2], -2.06143366e+08, 1e+6_f64); + assert_approx_eq!(gradient[2], -3232964116.6507816, 1e+6_f64); + } #[test] @@ -282,16 +355,17 @@ pub mod tests { let beam = GaussianBeam { direction: Vector3::x(), intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 2.0, + w0_x: 2.0, + w0_y: 2.0, power: 1.0, - rayleigh_range: calculate_rayleigh_range(&1064.0e-9, &2.0), - ellipticity: 0.0, + rayleigh_range_x: calculate_rayleigh_range(&1064.0e-9, &2.0), + rayleigh_range_y: calculate_rayleigh_range(&1064.0e-9, &2.0), }; let pos1 = Position { pos: Vector3::x() }; assert_approx_eq!( beam.power - / (PI.powf(0.5) * beam.e_radius).powf(2.0) + / (PI.powf(0.5) * beam.w0_x.sqrt() * PI.powf(0.5) * beam.w0_y.sqrt()) / (1.0 + 1.0 / calculate_rayleigh_range(&1064.0e-9, &2.0).powf(2.0)), get_gaussian_beam_intensity(&beam, &pos1, None, None), 1e-6_f64 @@ -299,96 +373,46 @@ pub mod tests { let pos2 = Position { pos: Vector3::y() }; assert_approx_eq!( - 1.0 / (PI.powf(0.5) * beam.e_radius).powf(2.0) - * (-pos2.pos[1] / beam.e_radius.powf(2.0)).exp(), + 1.0 / (PI.powf(0.5) * beam.w0_x.sqrt() * PI.powf(0.5) * beam.w0_y.sqrt()) + * (-pos2.pos[1] / beam.w0_x.sqrt() / beam.w0_y.sqrt()).exp(), get_gaussian_beam_intensity(&beam, &pos2, None, None), 1e-6_f64 ); assert_approx_eq!( beam.power - / (PI.powf(0.5) * beam.e_radius).powf(2.0) - / (1.0 + 1.0 / calculate_rayleigh_range(&1064.0e-9, &2.0).powf(2.0)), + / (PI.powf(0.5) * beam.w0_x/(2.0f64).sqrt() * PI.powf(0.5) * beam.w0_y/(2.0f64).sqrt()) + / (1.0 + 1.0 / calculate_rayleigh_range(&1064.0e-9, &(2.0*(2.0f64).sqrt())).powf(2.0)), get_gaussian_beam_intensity(&beam, &pos1, None, None), 1e-6_f64 ); assert_approx_eq!( - 1.0 / (PI.powf(0.5) * beam.e_radius).powf(2.0) - * (-pos2.pos[1] / beam.e_radius.powf(2.0)).exp(), + 1.0 / (PI.powf(0.5) * beam.w0_x.sqrt() * PI.powf(0.5) * beam.w0_y.sqrt()) + * (-pos2.pos[1] / beam.w0_x.sqrt() / beam.w0_y.sqrt()).exp(), get_gaussian_beam_intensity(&beam, &pos2, None, None), 1e-6_f64 ); - let rayleigh_range_2 = calculate_rayleigh_range(&1064.0e-6, &beam.e_radius); - let pos3 = Position { - pos: Vector3::x() * rayleigh_range_2, - }; - // Test with a frame but ellipticity = 0 - let frame = Frame::from_direction(beam.direction, Vector3::new(0.0, 1.0, 0.0)); - assert_approx_eq!( - beam.power / (PI.powf(0.5) * beam.e_radius).powf(2.0), - get_gaussian_beam_intensity(&beam, &pos3, None, Some(&frame)), - 1e-6_f64 - ); - // Position along the focused axis - let pos4 = Position { - pos: Vector3::x() + Vector3::y(), - }; - // Now with an ellipticity, that implies a/b = 2 - let beam = GaussianBeam { - direction: Vector3::x(), + let beam2 = GaussianBeam { + direction: Vector3::z(), intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 2.0, + w0_x: 2.0, + w0_y: 2.0, power: 1.0, - rayleigh_range: calculate_rayleigh_range(&1064.0e-9, &2.0), - ellipticity: (3.0 / 4.0_f64).powf(0.5), + rayleigh_range_x: calculate_rayleigh_range(&1064.0e-9, &2.0), + rayleigh_range_y: calculate_rayleigh_range(&1064.0e-9, &2.0), }; - // checking if value on x-axis stays the same (as without ellipticity and frame) - assert_approx_eq!( - beam.power / (PI.powf(0.5) * beam.e_radius).powf(2.0), - get_gaussian_beam_intensity(&beam, &pos3, None, Some(&frame)), - 1e-6_f64 - ); - - // manual calculation to get beam intensity - let intensity_0 = beam.power / (PI * beam.e_radius.powf(2.0)); - let broadening = 1.0 / (1.0 + (1.0 / beam.rayleigh_range).powf(2.0)); - // factor of 0.5 in exponent because of rescaling the axis by a/b = 2 - let intensity = - intensity_0 * broadening * EXP.powf(-0.5 * broadening / beam.e_radius.powf(2.0)); - + let pos4 = Position { pos: Vector3::new(1.0, 2.0, 3.0) }; + let intensity_from_rust = get_gaussian_beam_intensity(&beam2, &pos4, None, None); + println!("Intensity from rust: {}", intensity_from_rust); assert_approx_eq!( - intensity, - get_gaussian_beam_intensity(&beam, &pos4, None, Some(&frame)), - 1e-6_f64 + 0.013064233284686179, + intensity_from_rust, + 1e-12_f64 ); - // now the ration is a/b = 4 - let beam = GaussianBeam { - direction: Vector3::x(), - intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 2.0, - power: 1.0, - rayleigh_range: calculate_rayleigh_range(&1064.0e-9, &2.0), - ellipticity: (15.0 / 16.0_f64).powf(0.5), - }; - // but we check along the de-focused axis (so intensity is lower than in symmetrical case) - let intensity = - intensity_0 * broadening * EXP.powf(-4.0 * broadening / beam.e_radius.powf(2.0)); - assert_approx_eq!( - intensity, - get_gaussian_beam_intensity( - &beam, - &Position { - pos: Vector3::x() + Vector3::z(), - }, - None, - Some(&frame) - ), - 1e-6_f64 - ); } } diff --git a/src/laser/intensity.rs b/src/laser/intensity.rs index e3d8b590..fbc9ad36 100644 --- a/src/laser/intensity.rs +++ b/src/laser/intensity.rs @@ -162,10 +162,11 @@ pub mod tests { .with(GaussianBeam { direction: Vector3::new(1.0, 0.0, 0.0), intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 2.0, + w0_x: 2.0, + w0_y: 2.0, power: 1.0, - rayleigh_range: gaussian::calculate_rayleigh_range(&461.0e-9, &2.0), - ellipticity: 0.0, + rayleigh_range_x: gaussian::calculate_rayleigh_range(&461.0e-9, &(2.0)), + rayleigh_range_y: gaussian::calculate_rayleigh_range(&461.0e-9, &(2.0)), }) .build(); @@ -187,16 +188,19 @@ pub mod tests { &GaussianBeam { direction: Vector3::new(1.0, 0.0, 0.0), intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 2.0, + w0_x: 2.0, + w0_y: 2.0, power: 1.0, - rayleigh_range: gaussian::calculate_rayleigh_range(&461.0e-9, &2.0), - ellipticity: 0.0, + rayleigh_range_x: gaussian::calculate_rayleigh_range(&461.0e-9, &(2.0)), + rayleigh_range_y: gaussian::calculate_rayleigh_range(&461.0e-9, &(2.0)), }, &Position { pos: Vector3::y() }, None, None, ); + println!("{}", 5.0); + assert_approx_eq!( sampler_storage .get(atom1) diff --git a/src/laser/intensity_gradient.rs b/src/laser/intensity_gradient.rs index 395206b7..50fd7025 100644 --- a/src/laser/intensity_gradient.rs +++ b/src/laser/intensity_gradient.rs @@ -99,13 +99,17 @@ pub mod tests { let beam = GaussianBeam { direction: Vector3::z(), intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 70.71067812e-6, + w0_x: 70.71067812e-6 * f64::sqrt(2.0), + w0_y: 70.71067812e-6 * f64::sqrt(2.0), power: 100.0, - rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range( + rayleigh_range_x: crate::laser::gaussian::calculate_rayleigh_range( &1064.0e-9, - &70.71067812e-6, + &(70.71067812e-6 * f64::sqrt(2.0)), + ), + rayleigh_range_y: crate::laser::gaussian::calculate_rayleigh_range( + &1064.0e-9, + &(70.71067812e-6 * f64::sqrt(2.0)), ), - ellipticity: 0.0, }; test_world @@ -185,15 +189,19 @@ pub mod tests { test_world.register::(); let beam = GaussianBeam { - direction: Vector3::x(), + direction: Vector3::z(), intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 70.71067812e-6, + w0_x: 70.71067812e-6, + w0_y: 50.71067812e-6, power: 100.0, - rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range( + rayleigh_range_x: crate::laser::gaussian::calculate_rayleigh_range( + &1064.0e-9, + &(70.71067812e-6), + ), + rayleigh_range_y: crate::laser::gaussian::calculate_rayleigh_range( &1064.0e-9, - &70.71067812e-6, + &(50.71067812e-6), ), - ellipticity: 0.0, }; test_world @@ -204,8 +212,8 @@ pub mod tests { }) .with(beam) .with(Frame { - x_vector: Vector3::y(), - y_vector: Vector3::z(), + x_vector: Vector3::x(), + y_vector: Vector3::y(), }) .with(DipoleLight { wavelength: 1064.0e-9, @@ -233,8 +241,10 @@ 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); + println!("Simulated Gradient: {:?}", sim_result_gradient); + + assert_approx_eq!(-177345704576703.28, sim_result_gradient[0], 1e+5_f64); + assert_approx_eq!(-344817759803007.9, sim_result_gradient[1], 1e+8_f64); + assert_approx_eq!(-2144411873.4435744, sim_result_gradient[2], 1e+8_f64); } } diff --git a/src/laser_cooling/doppler.rs b/src/laser_cooling/doppler.rs index dfe18452..4c0508ea 100644 --- a/src/laser_cooling/doppler.rs +++ b/src/laser_cooling/doppler.rs @@ -139,10 +139,11 @@ pub mod tests { .with(GaussianBeam { direction: Vector3::new(1.0, 0.0, 0.0), intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 2.0, + w0_x: 2.0, + w0_y: 2.0, power: 1.0, - rayleigh_range: gaussian::calculate_rayleigh_range(&wavelength, &2.0), - ellipticity: 0.0, + rayleigh_range_x: gaussian::calculate_rayleigh_range(&wavelength, &2.0), + rayleigh_range_y: gaussian::calculate_rayleigh_range(&wavelength, &2.0), }) .build(); diff --git a/src/laser_cooling/force.rs b/src/laser_cooling/force.rs index 6e16d593..4a78fed1 100644 --- a/src/laser_cooling/force.rs +++ b/src/laser_cooling/force.rs @@ -229,10 +229,11 @@ pub mod tests { .with(GaussianBeam { direction: Vector3::new(1.0, 0.0, 0.0), intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 2.0, + w0_x: 2.0, + w0_y: 2.0, power: 1.0, - rayleigh_range: gaussian::calculate_rayleigh_range(&wavelength, &2.0), - ellipticity: 0.0, + rayleigh_range_x: gaussian::calculate_rayleigh_range(&wavelength, &2.0), + rayleigh_range_y: gaussian::calculate_rayleigh_range(&wavelength, &2.0), }) .build(); diff --git a/src/laser_cooling/rate.rs b/src/laser_cooling/rate.rs index 751ff9c3..77a7ea07 100644 --- a/src/laser_cooling/rate.rs +++ b/src/laser_cooling/rate.rs @@ -188,10 +188,11 @@ pub mod tests { .with(GaussianBeam { direction: Vector3::new(1.0, 0.0, 0.0), intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius: 2.0, + w0_x: 2.0, + w0_y: 2.0, power: 1.0, - rayleigh_range: 1.0, - ellipticity: 0.0, + rayleigh_range_x: crate::laser::gaussian::calculate_rayleigh_range(&wavelength, &2.0), + rayleigh_range_y: crate::laser::gaussian::calculate_rayleigh_range(&wavelength, &2.0), }) .build();