From 39f3f72c88ca75432334748d58b0b537d114743d Mon Sep 17 00:00:00 2001 From: Yijun Tang Date: Fri, 21 Nov 2025 17:17:01 +0000 Subject: [PATCH 01/15] modify gaussian.rs to be inputting 2 seperate waists and seperate rayleigh ranges --- src/laser/gaussian.rs | 393 ++++++++++++++++++++++++------------------ 1 file changed, 228 insertions(+), 165 deletions(-) diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index 4da1d11f..b40115d6 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -34,18 +34,18 @@ pub struct GaussianBeam { /// 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; @@ -66,17 +66,18 @@ impl GaussianBeam { 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, } } } @@ -99,18 +100,19 @@ impl GaussianBeam { 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; 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, 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), } } } @@ -164,56 +167,82 @@ pub fn get_gaussian_beam_intensity( beam: &GaussianBeam, pos: &Position, mask: Option<&CircularMask>, - frame: Option<&Frame>, + frame: &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)) + + + 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 +250,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 + + // Rayleigh ranges for each axis (waist squared) + let zx = beam.rayleigh_range_x; + let zy = beam.rayleigh_range_y; // convert for long axis - 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); + // 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 / spot_size_squared * vector + // 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)) + ); + + 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 +319,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 * (2.0f64).sqrt(), + w0_y: 70.71067812e-6 * (2.0f64).sqrt(), 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), @@ -282,16 +344,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 * (2.0f64).sqrt(), + w0_y: 2.0 * (2.0f64).sqrt(), 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/(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).powf(2.0)), get_gaussian_beam_intensity(&beam, &pos1, None, None), 1e-6_f64 @@ -299,96 +362,96 @@ 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/(2.0f64).sqrt() * PI.powf(0.5) * beam.w0_y/(2.0f64).sqrt()) + * (-pos2.pos[1] / beam.w0_x/(2.0f64).sqrt() / beam.w0_y/(2.0f64).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/(2.0f64).sqrt() * PI.powf(0.5) * beam.w0_y/(2.0f64).sqrt()) + * (-pos2.pos[1] / beam.w0_x/(2.0f64).sqrt() / beam.w0_y/(2.0f64).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(), - 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: (3.0 / 4.0_f64).powf(0.5), - }; - - // 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)); - - assert_approx_eq!( - intensity, - get_gaussian_beam_intensity(&beam, &pos4, None, Some(&frame)), - 1e-6_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 - ); + // 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(), + // 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: (3.0 / 4.0_f64).powf(0.5), + // }; + + // // 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)); + + // assert_approx_eq!( + // intensity, + // get_gaussian_beam_intensity(&beam, &pos4, None, Some(&frame)), + // 1e-6_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 + // ); } } From dc8a946ffcd92647c0225a6ba15459d381e9615a Mon Sep 17 00:00:00 2001 From: Yijun Tang Date: Fri, 21 Nov 2025 17:30:40 +0000 Subject: [PATCH 02/15] myre ference frame is &frame --- src/dipole/force.rs | 21 ++++++++++++--------- src/integration_tests/rate_equation.rs | 3 ++- src/laser/gaussian.rs | 2 +- src/laser/intensity.rs | 2 +- 4 files changed, 16 insertions(+), 12 deletions(-) diff --git a/src/dipole/force.rs b/src/dipole/force.rs index 233c0ca2..73541fef 100644 --- a/src/dipole/force.rs +++ b/src/dipole/force.rs @@ -165,15 +165,17 @@ pub mod tests { test_world.register::(); 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 gaussian_beam = GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius, - power, + w0_x: w0_x, + w0_y: w0_y, + power: power, direction: Vector3::x(), - rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range(&1064.0e-9, &e_radius), - ellipticity: 0.0, + 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() @@ -192,11 +194,12 @@ pub mod tests { .build(); let gaussian_beam = GaussianBeam { intersection: Vector3::new(0.0, 0.0, 0.0), - e_radius, - power, + w0_x: w0_x, + w0_y: w0_y, + power: power, direction: Vector3::y(), - rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range(&1064.0e-9, &e_radius), - ellipticity: 0.0, + 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() diff --git a/src/integration_tests/rate_equation.rs b/src/integration_tests/rate_equation.rs index edcf481b..4015e8d2 100644 --- a/src/integration_tests/rate_equation.rs +++ b/src/integration_tests/rate_equation.rs @@ -65,7 +65,8 @@ pub mod tests { Vector3::new(0.0, 0.0, 0.0), Vector3::new(-1.0, 0.0, 0.0), intensity, - 0.01, + 0.01 * (2.0_f64).sqrt(), + 0.01 * (2.0_f64).sqrt(), 780.0e-9, )) .build(); diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index b40115d6..b36da1df 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -135,7 +135,7 @@ impl GaussianBeam { direction: Vector3, power: f64, w0_x: f64 - w0_y, + w0_y: f64 wavelength: f64, ) -> Self { GaussianBeam { diff --git a/src/laser/intensity.rs b/src/laser/intensity.rs index e3d8b590..87a078b6 100644 --- a/src/laser/intensity.rs +++ b/src/laser/intensity.rs @@ -123,7 +123,7 @@ impl<'a, const N: usize> System<'a> for SampleLaserIntensitySystem { gaussian, pos, mask.as_ref(), - frame.as_ref(), + frame.as_ref().expect("REASON"), ); } }); From e8e1ce9ab86e618bc3c3f9d8ae5a154d51e3686e Mon Sep 17 00:00:00 2001 From: Yijun Tang Date: Fri, 21 Nov 2025 19:57:29 +0000 Subject: [PATCH 03/15] commit changes of gaussian beam in all related usage --- examples/1d_mot.rs | 14 +++++---- examples/2d_plus_mot_from_oven.rs | 35 ++++++++++++---------- examples/benchmark.rs | 42 +++++++++++++++------------ examples/cross_section.rs | 7 +++-- examples/dipole_trap.rs | 17 ++++++----- examples/doppler_limit.rs | 48 +++++++++++++++++-------------- examples/molasses_1d.rs | 14 +++++---- src/laser/gaussian.rs | 8 +++--- src/laser/intensity.rs | 16 ++++++----- src/laser/intensity_gradient.rs | 26 +++++++++++------ src/laser_cooling/doppler.rs | 7 +++-- src/laser_cooling/force.rs | 7 +++-- src/laser_cooling/rate.rs | 7 +++-- 13 files changed, 143 insertions(+), 105 deletions(-) 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/laser/gaussian.rs b/src/laser/gaussian.rs index b36da1df..fb9ab1ac 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -134,8 +134,8 @@ impl GaussianBeam { intersection: Vector3, direction: Vector3, power: f64, - w0_x: f64 - w0_y: f64 + w0_x: f64, + w0_y: f64, wavelength: f64, ) -> Self { GaussianBeam { @@ -167,7 +167,7 @@ pub fn get_gaussian_beam_intensity( beam: &GaussianBeam, pos: &Position, mask: Option<&CircularMask>, - frame: &Frame, + frame: Option<&Frame>, ) -> f64 { @@ -175,7 +175,7 @@ pub fn get_gaussian_beam_intensity( &pos.pos, &beam.intersection, &beam.direction, - frame, + frame.expect("REASON"), ); 2.0 * beam.power / PI / beam.w0_x / beam.w0_y / (1.0 + (z / beam.rayleigh_range_x).powf(2.0)).powf(0.5) diff --git a/src/laser/intensity.rs b/src/laser/intensity.rs index 87a078b6..0276f210 100644 --- a/src/laser/intensity.rs +++ b/src/laser/intensity.rs @@ -123,7 +123,7 @@ impl<'a, const N: usize> System<'a> for SampleLaserIntensitySystem { gaussian, pos, mask.as_ref(), - frame.as_ref().expect("REASON"), + Some(frame.as_ref().expect("REASON")), ); } }); @@ -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 * f64::sqrt(2.0), + w0_y: 2.0 * f64::sqrt(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 * f64::sqrt(2.0))), + rayleigh_range_y: gaussian::calculate_rayleigh_range(&461.0e-9, &(2.0 * f64::sqrt(2.0))), }) .build(); @@ -187,10 +188,11 @@ 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 * f64::sqrt(2.0), + w0_y: 2.0 * f64::sqrt(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 * f64::sqrt(2.0))), + rayleigh_range_y: gaussian::calculate_rayleigh_range(&461.0e-9, &(2.0 * f64::sqrt(2.0))), }, &Position { pos: Vector3::y() }, None, diff --git a/src/laser/intensity_gradient.rs b/src/laser/intensity_gradient.rs index 395206b7..f5ce6cd3 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 * 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 * f64::sqrt(2.0)), + ), + rayleigh_range_y: crate::laser::gaussian::calculate_rayleigh_range( &1064.0e-9, - &70.71067812e-6, + &(70.71067812e-6 * f64::sqrt(2.0)), ), - ellipticity: 0.0, }; test_world 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..35a6ae94 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: gaussian::calculate_rayleigh_range(&wavelength, &2.0), + rayleigh_range_y: gaussian::calculate_rayleigh_range(&wavelength, &2.0), }) .build(); From 0869e499530c7814f9ea06a91325c86efd359e5e Mon Sep 17 00:00:00 2001 From: Yijun Tang Date: Fri, 21 Nov 2025 20:02:50 +0000 Subject: [PATCH 04/15] commit changes to further place all old gaussian beam with new ones --- examples/top_trap_with_collisions.rs | 2 +- src/atom_sources/emit.rs | 6 +++--- src/atom_sources/mass.rs | 2 +- src/collisions.rs | 2 +- src/integrator.rs | 2 +- src/laser/gaussian.rs | 2 +- src/laser/index.rs | 2 +- src/laser/intensity.rs | 2 +- src/laser_cooling/doppler.rs | 2 +- src/laser_cooling/photons_scattered.rs | 2 +- src/laser_cooling/rate.rs | 3 ++- src/laser_cooling/sampler.rs | 2 +- src/magnetic/mod.rs | 6 +++--- 13 files changed, 18 insertions(+), 17 deletions(-) diff --git a/examples/top_trap_with_collisions.rs b/examples/top_trap_with_collisions.rs index 5e5174a8..c246ad90 100644 --- a/examples/top_trap_with_collisions.rs +++ b/examples/top_trap_with_collisions.rs @@ -7,7 +7,7 @@ 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::{MagneticDipole}; use lib::magnetic::quadrupole::QuadrupoleField3D; use lib::magnetic::top::TimeOrbitingPotential; use lib::magnetic::MagneticTrapPlugin; diff --git a/src/atom_sources/emit.rs b/src/atom_sources/emit.rs index 79b671bc..dee46d73 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, number_to_emit) in (&emit_numbers, &mut numbers_to_emit).join() { number_to_emit.number = emit_number.number; } } @@ -69,7 +69,7 @@ impl<'a> System<'a> for EmitFixedRateSystem { fn run(&mut self, (rates, timestep, mut emit_numbers): Self::SystemData) { let mut rng = rand::thread_rng(); - for (rate, mut emit_numbers) in (&rates, &mut emit_numbers).join() { + for (rate, emit_numbers) in (&rates, &mut emit_numbers).join() { let avg_number_to_emit = rate.rate * timestep.delta; let guaranteed_number = avg_number_to_emit.floor(); let number: i32; @@ -92,7 +92,7 @@ impl<'a> System<'a> for EmitOnceSystem { ); fn run(&mut self, (emit_onces, mut emit_numbers): Self::SystemData) { - for (_, mut emit_numbers) in (&emit_onces, &mut emit_numbers).join() { + for (_, emit_numbers) in (&emit_onces, &mut emit_numbers).join() { emit_numbers.number = 0; } } diff --git a/src/atom_sources/mass.rs b/src/atom_sources/mass.rs index 6923aec8..85d70f66 100644 --- a/src/atom_sources/mass.rs +++ b/src/atom_sources/mass.rs @@ -48,7 +48,7 @@ impl MassDistribution { total += mr.ratio; } - for mut mr in &mut self.distribution { + for mr in &mut self.distribution { mr.ratio /= total; } self.normalised = true diff --git a/src/collisions.rs b/src/collisions.rs index 498bdfb5..b6882d4c 100644 --- a/src/collisions.rs +++ b/src/collisions.rs @@ -198,7 +198,7 @@ impl<'a> System<'a> for ApplyCollisionsSystem { // build list of ids for each atom (&positions, &mut boxids) .par_join() - .for_each(|(position, mut boxid)| { + .for_each(|(position, boxid)| { boxid.id = pos_to_id(position.pos, n, params.box_width); }); diff --git a/src/integrator.rs b/src/integrator.rs index d632dd01..e056c909 100644 --- a/src/integrator.rs +++ b/src/integrator.rs @@ -86,7 +86,7 @@ impl<'a> System<'a> for VelocityVerletIntegratePositionSystem { (&mut pos, &vel, &mut old_force, &force, &mass) .par_join() - .for_each(|(mut pos, vel, mut old_force, force, mass)| { + .for_each(|(pos, vel, old_force, force, mass)| { pos.pos = pos.pos + vel.vel * dt + force.force / (constant::AMU * mass.value) / 2.0 * dt * dt; diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index fb9ab1ac..816f9de1 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -166,7 +166,7 @@ impl Component for CircularMask { pub fn get_gaussian_beam_intensity( beam: &GaussianBeam, pos: &Position, - mask: Option<&CircularMask>, + _mask: Option<&CircularMask>, frame: Option<&Frame>, ) -> f64 { diff --git a/src/laser/index.rs b/src/laser/index.rs index 7d25098b..7fe7653d 100644 --- a/src/laser/index.rs +++ b/src/laser/index.rs @@ -31,7 +31,7 @@ impl<'a> System<'a> for IndexLasersSystem { } } if need_to_assign_indices { - for mut index in (&mut indices).join() { + for index in (&mut indices).join() { index.index = iter; index.initiated = true; iter += 1; diff --git a/src/laser/intensity.rs b/src/laser/intensity.rs index 0276f210..42d65165 100644 --- a/src/laser/intensity.rs +++ b/src/laser/intensity.rs @@ -53,7 +53,7 @@ impl<'a, const N: usize> System<'a> for InitialiseLaserIntensitySamplersSystem System<'a> for InitialiseDopplerShiftSamplersSystem fn run(&mut self, (mut samplers,): Self::SystemData) { use rayon::prelude::*; - (&mut samplers).par_join().for_each(|mut sampler| { + (&mut samplers).par_join().for_each(|sampler| { sampler.contents = [DopplerShiftSampler::default(); N]; }); } diff --git a/src/laser_cooling/photons_scattered.rs b/src/laser_cooling/photons_scattered.rs index bab58164..167d0e77 100644 --- a/src/laser_cooling/photons_scattered.rs +++ b/src/laser_cooling/photons_scattered.rs @@ -123,7 +123,7 @@ impl<'a, T, const N: usize> System<'a> for InitialiseExpectedPhotonsScatteredVec fn run(&mut self, (mut expected_photons,): Self::SystemData) { use rayon::prelude::*; - (&mut expected_photons).par_join().for_each(|mut expected| { + (&mut expected_photons).par_join().for_each(|expected| { expected.contents = [ExpectedPhotonsScattered::default(); N]; }); } diff --git a/src/laser_cooling/rate.rs b/src/laser_cooling/rate.rs index 35a6ae94..90fd64a2 100644 --- a/src/laser_cooling/rate.rs +++ b/src/laser_cooling/rate.rs @@ -6,6 +6,7 @@ use std::marker::PhantomData; use super::CoolingLight; use super::transition::{TransitionComponent}; +use crate::laser::gaussian; use crate::laser::gaussian::GaussianBeam; use crate::laser::index::LaserIndex; use crate::laser::intensity::LaserIntensitySamplers; @@ -57,7 +58,7 @@ impl<'a, T, const N: usize> System<'a> for InitialiseRateCoefficientsSystem System<'a> for InitialiseLaserDetuningSamplersSystem fn run(&mut self, (mut samplers,): Self::SystemData) { use rayon::prelude::*; - (&mut samplers).par_join().for_each(|mut sampler| { + (&mut samplers).par_join().for_each(|sampler| { sampler.contents = [LaserDetuningSampler::default(); N]; }); } diff --git a/src/magnetic/mod.rs b/src/magnetic/mod.rs index dce97c51..5a1eba3a 100644 --- a/src/magnetic/mod.rs +++ b/src/magnetic/mod.rs @@ -76,7 +76,7 @@ impl<'a> System<'a> for ClearMagneticFieldSamplerSystem { fn run(&mut self, mut sampler: Self::SystemData) { use rayon::prelude::*; - (&mut sampler).par_join().for_each(|mut sampler| { + (&mut sampler).par_join().for_each(|sampler| { sampler.magnitude = 0.; sampler.field = Vector3::new(0.0, 0.0, 0.0); sampler.gradient = Vector3::new(0.0, 0.0, 0.0); @@ -96,7 +96,7 @@ impl<'a> System<'a> for CalculateMagneticFieldMagnitudeSystem { fn run(&mut self, mut sampler: Self::SystemData) { use rayon::prelude::*; - (&mut sampler).par_join().for_each(|mut sampler| { + (&mut sampler).par_join().for_each(|sampler| { sampler.magnitude = sampler.field.norm(); if sampler.magnitude.is_nan() { sampler.magnitude = 0.0; @@ -115,7 +115,7 @@ impl<'a> System<'a> for CalculateMagneticMagnitudeGradientSystem { fn run(&mut self, mut sampler: Self::SystemData) { use rayon::prelude::*; - (&mut sampler).par_join().for_each(|mut sampler| { + (&mut sampler).par_join().for_each(|sampler| { let mut gradient = Vector3::new(0.0, 0.0, 0.0); for i in 0..3 { gradient[i] = From 8841fd66fb64dec3ac98b3e15a1bf10564d5d5e2 Mon Sep 17 00:00:00 2001 From: Yijun Tang Date: Fri, 21 Nov 2025 20:03:30 +0000 Subject: [PATCH 05/15] cargo fix lib --- src/laser_cooling/rate.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/laser_cooling/rate.rs b/src/laser_cooling/rate.rs index 90fd64a2..b7329667 100644 --- a/src/laser_cooling/rate.rs +++ b/src/laser_cooling/rate.rs @@ -6,7 +6,7 @@ use std::marker::PhantomData; use super::CoolingLight; use super::transition::{TransitionComponent}; -use crate::laser::gaussian; + use crate::laser::gaussian::GaussianBeam; use crate::laser::index::LaserIndex; use crate::laser::intensity::LaserIntensitySamplers; From 5287f45f5b955bb3b9e80f149149dfb36070e0ed Mon Sep 17 00:00:00 2001 From: Yijun Tang Date: Fri, 21 Nov 2025 21:20:56 +0000 Subject: [PATCH 06/15] correct rate.rs --- src/laser_cooling/rate.rs | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/laser_cooling/rate.rs b/src/laser_cooling/rate.rs index b7329667..9f96dca7 100644 --- a/src/laser_cooling/rate.rs +++ b/src/laser_cooling/rate.rs @@ -6,7 +6,6 @@ use std::marker::PhantomData; use super::CoolingLight; use super::transition::{TransitionComponent}; - use crate::laser::gaussian::GaussianBeam; use crate::laser::index::LaserIndex; use crate::laser::intensity::LaserIntensitySamplers; @@ -192,8 +191,8 @@ pub mod tests { w0_x: 2.0, w0_y: 2.0, power: 1.0, - rayleigh_range_x: gaussian::calculate_rayleigh_range(&wavelength, &2.0), - rayleigh_range_y: gaussian::calculate_rayleigh_range(&wavelength, &2.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(); From a3f16ddf1ca271684559128b1896c5e772a7c970 Mon Sep 17 00:00:00 2001 From: Yijun Tang Date: Fri, 21 Nov 2025 23:47:33 +0000 Subject: [PATCH 07/15] Finish first unit test --- src/dipole/force.rs | 43 ++++++++++++++++++++++++++----------------- 1 file changed, 26 insertions(+), 17 deletions(-) diff --git a/src/dipole/force.rs b/src/dipole/force.rs index 73541fef..3b9f8c17 100644 --- a/src/dipole/force.rs +++ b/src/dipole/force.rs @@ -165,7 +165,7 @@ pub mod tests { test_world.register::(); let power = 10.0; - let w0_x = 60.0e-6; + let w0_x = 50.0e-6; let w0_y = 60.0e-6; let gaussian_beam = GaussianBeam { @@ -173,7 +173,7 @@ pub mod tests { w0_x: w0_x, w0_y: w0_y, power: power, - direction: Vector3::x(), + 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), }; @@ -188,16 +188,16 @@ 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), w0_x: w0_x, w0_y: w0_y, - power: power, - direction: Vector3::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), }; @@ -212,7 +212,7 @@ pub mod tests { initiated: true, }) .with(laser::frame::Frame { - x_vector: Vector3::x(), + x_vector: Vector3::y(), y_vector: Vector3::z(), }) .build(); @@ -244,28 +244,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 ); } } + From f07fe844c9ada0eb6db16b59b465ed67d4768451 Mon Sep 17 00:00:00 2001 From: Yijun Tang Date: Mon, 24 Nov 2025 15:33:22 +0000 Subject: [PATCH 08/15] intensity gradient unit test fixed --- src/integration_tests/rate_equation.rs | 4 ++-- src/laser/gaussian.rs | 12 ++++++++---- src/laser/intensity.rs | 18 ++++++++++-------- src/laser/intensity_gradient.rs | 20 +++++++++++--------- 4 files changed, 31 insertions(+), 23 deletions(-) diff --git a/src/integration_tests/rate_equation.rs b/src/integration_tests/rate_equation.rs index 4015e8d2..d23398d5 100644 --- a/src/integration_tests/rate_equation.rs +++ b/src/integration_tests/rate_equation.rs @@ -65,8 +65,8 @@ pub mod tests { Vector3::new(0.0, 0.0, 0.0), Vector3::new(-1.0, 0.0, 0.0), intensity, - 0.01 * (2.0_f64).sqrt(), - 0.01 * (2.0_f64).sqrt(), + 0.01, + 0.01, 780.0e-9, )) .build(); diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index 816f9de1..887eed1c 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -319,8 +319,8 @@ pub mod tests { let beam = GaussianBeam { direction: Vector3::z(), intersection: Vector3::new(0.0, 0.0, 0.0), - w0_x: 70.71067812e-6 * (2.0f64).sqrt(), - w0_y: 70.71067812e-6 * (2.0f64).sqrt(), + w0_x: 70.71067812e-6, + w0_y: 70.71067812e-6, power: 100.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), @@ -334,9 +334,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] diff --git a/src/laser/intensity.rs b/src/laser/intensity.rs index 42d65165..e92c3b97 100644 --- a/src/laser/intensity.rs +++ b/src/laser/intensity.rs @@ -162,11 +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), - w0_x: 2.0 * f64::sqrt(2.0), - w0_y: 2.0 * f64::sqrt(2.0), + w0_x: 2.0, + w0_y: 2.0, power: 1.0, - rayleigh_range_x: gaussian::calculate_rayleigh_range(&461.0e-9, &(2.0 * f64::sqrt(2.0))), - rayleigh_range_y: gaussian::calculate_rayleigh_range(&461.0e-9, &(2.0 * f64::sqrt(2.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(); @@ -188,17 +188,19 @@ pub mod tests { &GaussianBeam { direction: Vector3::new(1.0, 0.0, 0.0), intersection: Vector3::new(0.0, 0.0, 0.0), - w0_x: 2.0 * f64::sqrt(2.0), - w0_y: 2.0 * f64::sqrt(2.0), + w0_x: 2.0, + w0_y: 2.0, power: 1.0, - rayleigh_range_x: gaussian::calculate_rayleigh_range(&461.0e-9, &(2.0 * f64::sqrt(2.0))), - rayleigh_range_y: gaussian::calculate_rayleigh_range(&461.0e-9, &(2.0 * f64::sqrt(2.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 f5ce6cd3..50fd7025 100644 --- a/src/laser/intensity_gradient.rs +++ b/src/laser/intensity_gradient.rs @@ -191,16 +191,16 @@ pub mod tests { let beam = GaussianBeam { direction: Vector3::z(), intersection: Vector3::new(0.0, 0.0, 0.0), - w0_x: 70.71067812e-6 * f64::sqrt(2.0), - w0_y: 70.71067812e-6 * f64::sqrt(2.0), + w0_x: 70.71067812e-6, + w0_y: 50.71067812e-6, power: 100.0, rayleigh_range_x: crate::laser::gaussian::calculate_rayleigh_range( &1064.0e-9, - &(70.71067812e-6 * f64::sqrt(2.0)), + &(70.71067812e-6), ), rayleigh_range_y: crate::laser::gaussian::calculate_rayleigh_range( &1064.0e-9, - &(70.71067812e-6 * f64::sqrt(2.0)), + &(50.71067812e-6), ), }; @@ -212,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, @@ -241,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); } } From 7ae6a85648f69c8fd55a37410e52117cfee3d3bd Mon Sep 17 00:00:00 2001 From: Yijun Tang Date: Mon, 24 Nov 2025 16:11:27 +0000 Subject: [PATCH 09/15] fixing gaussian intensity unit test --- src/laser/gaussian.rs | 113 +++++++++++++---------------------------- src/laser/intensity.rs | 2 +- 2 files changed, 35 insertions(+), 80 deletions(-) diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index 887eed1c..7b7bbad7 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -170,12 +170,17 @@ pub fn get_gaussian_beam_intensity( frame: Option<&Frame>, ) -> f64 { + let binding = Frame{ + x_vector: Vector3::x(), + y_vector: Vector3::y(), + }; + let frame = frame.unwrap_or(&binding); let (x, y, z) = maths::get_relative_coordinates_line_point( &pos.pos, &beam.intersection, &beam.direction, - frame.expect("REASON"), + &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) @@ -348,8 +353,8 @@ pub mod tests { let beam = GaussianBeam { direction: Vector3::x(), intersection: Vector3::new(0.0, 0.0, 0.0), - w0_x: 2.0 * (2.0f64).sqrt(), - w0_y: 2.0 * (2.0f64).sqrt(), + w0_x: 2.0, + w0_y: 2.0, power: 1.0, rayleigh_range_x: calculate_rayleigh_range(&1064.0e-9, &2.0), rayleigh_range_y: calculate_rayleigh_range(&1064.0e-9, &2.0), @@ -358,7 +363,7 @@ pub mod tests { let pos1 = Position { pos: Vector3::x() }; assert_approx_eq!( beam.power - / (PI.powf(0.5) * beam.w0_x/(2.0f64).sqrt() * PI.powf(0.5) * beam.w0_y/(2.0f64).sqrt()) + / (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 @@ -366,8 +371,8 @@ pub mod tests { let pos2 = Position { pos: Vector3::y() }; assert_approx_eq!( - 1.0 / (PI.powf(0.5) * beam.w0_x/(2.0f64).sqrt() * PI.powf(0.5) * beam.w0_y/(2.0f64).sqrt()) - * (-pos2.pos[1] / beam.w0_x/(2.0f64).sqrt() / beam.w0_y/(2.0f64).sqrt()).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 ); @@ -381,81 +386,31 @@ pub mod tests { ); assert_approx_eq!( - 1.0 / (PI.powf(0.5) * beam.w0_x/(2.0f64).sqrt() * PI.powf(0.5) * beam.w0_y/(2.0f64).sqrt()) - * (-pos2.pos[1] / beam.w0_x/(2.0f64).sqrt() / beam.w0_y/(2.0f64).sqrt()).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(), - // 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: (3.0 / 4.0_f64).powf(0.5), - // }; - - // // 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)); - - // assert_approx_eq!( - // intensity, - // get_gaussian_beam_intensity(&beam, &pos4, None, Some(&frame)), - // 1e-6_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 - // ); + + + let beam2 = GaussianBeam { + direction: Vector3::z(), + intersection: Vector3::new(0.0, 0.0, 0.0), + w0_x: 2.0, + w0_y: 2.0, + power: 1.0, + rayleigh_range_x: calculate_rayleigh_range(&1064.0e-9, &2.0), + rayleigh_range_y: calculate_rayleigh_range(&1064.0e-9, &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!( + 0.013064233284686179, + intensity_from_rust, + 1e-12_f64 + ); + } } diff --git a/src/laser/intensity.rs b/src/laser/intensity.rs index e92c3b97..e738205f 100644 --- a/src/laser/intensity.rs +++ b/src/laser/intensity.rs @@ -123,7 +123,7 @@ impl<'a, const N: usize> System<'a> for SampleLaserIntensitySystem { gaussian, pos, mask.as_ref(), - Some(frame.as_ref().expect("REASON")), + frame.as_ref(), ); } }); From 0c6929f78aaf98c9cd1cefc8afe600149da4a61c Mon Sep 17 00:00:00 2001 From: Yijun Tang Date: Mon, 24 Nov 2025 16:15:32 +0000 Subject: [PATCH 10/15] forget /2 in gaussain beam from peak intensity (create gaussian beram_) --- 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 7b7bbad7..c3f4cd5e 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -104,7 +104,7 @@ impl GaussianBeam { w0_y: f64, wavelength: f64, ) -> Self { - let power = std::f64::consts::PI * w0_x * w0_y * peak_intensity; + let power = std::f64::consts::PI * w0_x * w0_y * peak_intensity / 2.0; GaussianBeam { intersection, direction, From cd4344057a41bd6e821e2423156dea5d560d22bf Mon Sep 17 00:00:00 2001 From: Yijun Tang Date: Thu, 27 Nov 2025 11:31:02 +0000 Subject: [PATCH 11/15] comment added about direction and frame, also made some renaing to avoid confusion --- src/dipole/force.rs | 3 +++ src/integration_tests/rate_equation.rs | 2 +- src/laser/gaussian.rs | 24 +++++++++++++----------- 3 files changed, 17 insertions(+), 12 deletions(-) diff --git a/src/dipole/force.rs b/src/dipole/force.rs index 3b9f8c17..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, diff --git a/src/integration_tests/rate_equation.rs b/src/integration_tests/rate_equation.rs index d23398d5..1b3a7dad 100644 --- a/src/integration_tests/rate_equation.rs +++ b/src/integration_tests/rate_equation.rs @@ -61,7 +61,7 @@ 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, diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs index c3f4cd5e..1f2cd3f6 100644 --- a/src/laser/gaussian.rs +++ b/src/laser/gaussian.rs @@ -26,14 +26,10 @@ 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: - /// /// Power of the laser in W pub power: f64, @@ -61,7 +57,10 @@ 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, @@ -78,6 +77,7 @@ impl GaussianBeam { w0_y, rayleigh_range_x: f64::INFINITY, rayleigh_range_y: f64::INFINITY, + } } } @@ -93,10 +93,10 @@ 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, @@ -163,13 +163,15 @@ 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>, frame: Option<&Frame>, ) -> f64 { - + /// 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(), From 71b5c770cda9bdd49e319ae526a2eff4078064aa Mon Sep 17 00:00:00 2001 From: Yijun Tang Date: Thu, 27 Nov 2025 16:37:02 +0000 Subject: [PATCH 12/15] revert wierd changes in mut, don't know why a lot of mut are deleted in this comment --- examples/top_trap_with_collisions.rs | 2 +- src/atom_sources/emit.rs | 6 +++--- src/atom_sources/mass.rs | 2 +- src/collisions.rs | 2 +- src/integrator.rs | 2 +- src/laser/index.rs | 2 +- src/laser/intensity.rs | 2 +- src/laser_cooling/doppler.rs | 2 +- src/laser_cooling/photons_scattered.rs | 2 +- src/laser_cooling/rate.rs | 2 +- src/laser_cooling/sampler.rs | 2 +- src/magnetic/mod.rs | 2 +- 12 files changed, 14 insertions(+), 14 deletions(-) diff --git a/examples/top_trap_with_collisions.rs b/examples/top_trap_with_collisions.rs index c246ad90..5e5174a8 100644 --- a/examples/top_trap_with_collisions.rs +++ b/examples/top_trap_with_collisions.rs @@ -7,7 +7,7 @@ use lib::collisions::CollisionPlugin; use lib::collisions::{ApplyCollisionsOption, CollisionParameters, CollisionsTracker}; use lib::initiate::NewlyCreated; use lib::integrator::Timestep; -use lib::magnetic::force::{MagneticDipole}; +use lib::magnetic::force::{ApplyMagneticForceSystem, MagneticDipole}; use lib::magnetic::quadrupole::QuadrupoleField3D; use lib::magnetic::top::TimeOrbitingPotential; use lib::magnetic::MagneticTrapPlugin; diff --git a/src/atom_sources/emit.rs b/src/atom_sources/emit.rs index dee46d73..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, 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; } } @@ -69,7 +69,7 @@ impl<'a> System<'a> for EmitFixedRateSystem { fn run(&mut self, (rates, timestep, mut emit_numbers): Self::SystemData) { let mut rng = rand::thread_rng(); - for (rate, emit_numbers) in (&rates, &mut emit_numbers).join() { + for (rate, mut emit_numbers) in (&rates, &mut emit_numbers).join() { let avg_number_to_emit = rate.rate * timestep.delta; let guaranteed_number = avg_number_to_emit.floor(); let number: i32; @@ -92,7 +92,7 @@ impl<'a> System<'a> for EmitOnceSystem { ); fn run(&mut self, (emit_onces, mut emit_numbers): Self::SystemData) { - for (_, emit_numbers) in (&emit_onces, &mut emit_numbers).join() { + for (_, mut emit_numbers) in (&emit_onces, &mut emit_numbers).join() { emit_numbers.number = 0; } } diff --git a/src/atom_sources/mass.rs b/src/atom_sources/mass.rs index 85d70f66..6923aec8 100644 --- a/src/atom_sources/mass.rs +++ b/src/atom_sources/mass.rs @@ -48,7 +48,7 @@ impl MassDistribution { total += mr.ratio; } - for mr in &mut self.distribution { + for mut mr in &mut self.distribution { mr.ratio /= total; } self.normalised = true diff --git a/src/collisions.rs b/src/collisions.rs index b6882d4c..498bdfb5 100644 --- a/src/collisions.rs +++ b/src/collisions.rs @@ -198,7 +198,7 @@ impl<'a> System<'a> for ApplyCollisionsSystem { // build list of ids for each atom (&positions, &mut boxids) .par_join() - .for_each(|(position, boxid)| { + .for_each(|(position, mut boxid)| { boxid.id = pos_to_id(position.pos, n, params.box_width); }); diff --git a/src/integrator.rs b/src/integrator.rs index e056c909..d632dd01 100644 --- a/src/integrator.rs +++ b/src/integrator.rs @@ -86,7 +86,7 @@ impl<'a> System<'a> for VelocityVerletIntegratePositionSystem { (&mut pos, &vel, &mut old_force, &force, &mass) .par_join() - .for_each(|(pos, vel, old_force, force, mass)| { + .for_each(|(mut pos, vel, mut old_force, force, mass)| { pos.pos = pos.pos + vel.vel * dt + force.force / (constant::AMU * mass.value) / 2.0 * dt * dt; diff --git a/src/laser/index.rs b/src/laser/index.rs index 7fe7653d..7d25098b 100644 --- a/src/laser/index.rs +++ b/src/laser/index.rs @@ -31,7 +31,7 @@ impl<'a> System<'a> for IndexLasersSystem { } } if need_to_assign_indices { - for index in (&mut indices).join() { + for mut index in (&mut indices).join() { index.index = iter; index.initiated = true; iter += 1; diff --git a/src/laser/intensity.rs b/src/laser/intensity.rs index e738205f..fbc9ad36 100644 --- a/src/laser/intensity.rs +++ b/src/laser/intensity.rs @@ -53,7 +53,7 @@ impl<'a, const N: usize> System<'a> for InitialiseLaserIntensitySamplersSystem System<'a> for InitialiseDopplerShiftSamplersSystem fn run(&mut self, (mut samplers,): Self::SystemData) { use rayon::prelude::*; - (&mut samplers).par_join().for_each(|sampler| { + (&mut samplers).par_join().for_each(|mut sampler| { sampler.contents = [DopplerShiftSampler::default(); N]; }); } diff --git a/src/laser_cooling/photons_scattered.rs b/src/laser_cooling/photons_scattered.rs index 167d0e77..bab58164 100644 --- a/src/laser_cooling/photons_scattered.rs +++ b/src/laser_cooling/photons_scattered.rs @@ -123,7 +123,7 @@ impl<'a, T, const N: usize> System<'a> for InitialiseExpectedPhotonsScatteredVec fn run(&mut self, (mut expected_photons,): Self::SystemData) { use rayon::prelude::*; - (&mut expected_photons).par_join().for_each(|expected| { + (&mut expected_photons).par_join().for_each(|mut expected| { expected.contents = [ExpectedPhotonsScattered::default(); N]; }); } diff --git a/src/laser_cooling/rate.rs b/src/laser_cooling/rate.rs index 9f96dca7..77a7ea07 100644 --- a/src/laser_cooling/rate.rs +++ b/src/laser_cooling/rate.rs @@ -57,7 +57,7 @@ impl<'a, T, const N: usize> System<'a> for InitialiseRateCoefficientsSystem System<'a> for InitialiseLaserDetuningSamplersSystem fn run(&mut self, (mut samplers,): Self::SystemData) { use rayon::prelude::*; - (&mut samplers).par_join().for_each(|sampler| { + (&mut samplers).par_join().for_each(|mut sampler| { sampler.contents = [LaserDetuningSampler::default(); N]; }); } diff --git a/src/magnetic/mod.rs b/src/magnetic/mod.rs index 5a1eba3a..710d5c16 100644 --- a/src/magnetic/mod.rs +++ b/src/magnetic/mod.rs @@ -76,7 +76,7 @@ impl<'a> System<'a> for ClearMagneticFieldSamplerSystem { fn run(&mut self, mut sampler: Self::SystemData) { use rayon::prelude::*; - (&mut sampler).par_join().for_each(|sampler| { + (&mut sampler).par_join().for_each(|mut sampler| { sampler.magnitude = 0.; sampler.field = Vector3::new(0.0, 0.0, 0.0); sampler.gradient = Vector3::new(0.0, 0.0, 0.0); From b3567b71a93c58b3c479a30fcaf7566e1765dc53 Mon Sep 17 00:00:00 2001 From: Yijun Tang Date: Thu, 27 Nov 2025 16:42:19 +0000 Subject: [PATCH 13/15] some dummy change in the storgae --- src/dipole/force.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dipole/force.rs b/src/dipole/force.rs index b3b7e874..3de298f7 100644 --- a/src/dipole/force.rs +++ b/src/dipole/force.rs @@ -247,7 +247,7 @@ 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; From cc303e506ca4a85fd65f8faa0cdd3b60bcd05b6c Mon Sep 17 00:00:00 2001 From: Yijun Tang Date: Thu, 27 Nov 2025 16:43:10 +0000 Subject: [PATCH 14/15] correct it --- src/dipole/force.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dipole/force.rs b/src/dipole/force.rs index 3de298f7..b3b7e874 100644 --- a/src/dipole/force.rs +++ b/src/dipole/force.rs @@ -247,7 +247,7 @@ 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; From 512d022c9272c786992302fdf05a67596140663c Mon Sep 17 00:00:00 2001 From: Yijun Tang Date: Mon, 1 Dec 2025 15:21:08 +0000 Subject: [PATCH 15/15] mut in magnetic mod module --- src/magnetic/mod.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/magnetic/mod.rs b/src/magnetic/mod.rs index 710d5c16..dce97c51 100644 --- a/src/magnetic/mod.rs +++ b/src/magnetic/mod.rs @@ -96,7 +96,7 @@ impl<'a> System<'a> for CalculateMagneticFieldMagnitudeSystem { fn run(&mut self, mut sampler: Self::SystemData) { use rayon::prelude::*; - (&mut sampler).par_join().for_each(|sampler| { + (&mut sampler).par_join().for_each(|mut sampler| { sampler.magnitude = sampler.field.norm(); if sampler.magnitude.is_nan() { sampler.magnitude = 0.0; @@ -115,7 +115,7 @@ impl<'a> System<'a> for CalculateMagneticMagnitudeGradientSystem { fn run(&mut self, mut sampler: Self::SystemData) { use rayon::prelude::*; - (&mut sampler).par_join().for_each(|sampler| { + (&mut sampler).par_join().for_each(|mut sampler| { let mut gradient = Vector3::new(0.0, 0.0, 0.0); for i in 0..3 { gradient[i] =