diff --git a/src/dipole/force.rs b/src/dipole/force.rs
index 233c0ca2..f58e12b7 100644
--- a/src/dipole/force.rs
+++ b/src/dipole/force.rs
@@ -151,7 +151,7 @@ pub mod tests {
         assert_approx_eq!(-3.11151847e-23, sim_result_force[2], 2e-24_f64);
     }
 
-    #[test]
+    #[test] // Modified by Yijun on 2024/07/16
     fn test_apply_dipole_force_and_gradient_system() {
         let mut test_world = World::new();
 
@@ -171,7 +171,7 @@ pub mod tests {
             intersection: Vector3::new(0.0, 0.0, 0.0),
             e_radius,
             power,
-            direction: Vector3::x(),
+            direction: Vector3::y(),
             rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range(&1064.0e-9, &e_radius),
             ellipticity: 0.0,
         };
@@ -186,15 +186,15 @@ pub mod tests {
                 initiated: true,
             })
             .with(laser::frame::Frame {
-                x_vector: Vector3::y(),
-                y_vector: Vector3::z(),
+                x_vector: Vector3::z(),
+                y_vector: Vector3::x(),
             })
             .build();
         let gaussian_beam = GaussianBeam {
             intersection: Vector3::new(0.0, 0.0, 0.0),
             e_radius,
             power,
-            direction: Vector3::y(),
+            direction: Vector3::x(),
             rayleigh_range: crate::laser::gaussian::calculate_rayleigh_range(&1064.0e-9, &e_radius),
             ellipticity: 0.0,
         };
@@ -209,7 +209,7 @@ pub mod tests {
                 initiated: true,
             })
             .with(laser::frame::Frame {
-                x_vector: Vector3::x(),
+                x_vector: Vector3::y(),
                 y_vector: Vector3::z(),
             })
             .build();
@@ -218,7 +218,7 @@ pub mod tests {
         let atom1 = test_world
             .create_entity()
             .with(crate::atom::Position {
-                pos: Vector3::new(-1.0e-4, -1.0e-4, -2.0e-4),
+                pos: Vector3::new(-1.0e-4, -1.0e-4, 2.0e-4),
             })
             .with(Force {
                 force: Vector3::new(0.0, 0.0, 0.0),
@@ -248,21 +248,23 @@ pub mod tests {
         //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 x is: {}", sim_result_force[0]);
+        println!("force y is: {}", sim_result_force[1]);
+        println!("force z is: {}", sim_result_force[2]);
         assert_approx_eq!(
-            0.000000000000000000000000000000000127913190642808,
+            0.00000000000000000000000000000000012749214960,
             sim_result_force[0],
-            3e-46_f64
+            3e-40_f64
         );
         assert_approx_eq!(
-            0.000000000000000000000000000000000127913190642808,
+            0.00000000000000000000000000000000012749214960,
             sim_result_force[1],
-            2e-46_f64
+            2e-40_f64
         );
         assert_approx_eq!(
-            0.000000000000000000000000000000000511875188257342,
+            -0.0000000000000000000000000000000005101862262458,
             sim_result_force[2],
-            2e-46_f64
+            2e-40_f64
         );
     }
 }
diff --git a/src/laser/gaussian.rs b/src/laser/gaussian.rs
index 4da1d11f..e0589024 100644
--- a/src/laser/gaussian.rs
+++ b/src/laser/gaussian.rs
@@ -232,7 +232,8 @@ pub fn get_gaussian_beam_intensity_gradient(
 
     let spot_size_squared =
         2.0 * beam.e_radius.powf(2.0) * (1. + (z / beam.rayleigh_range).powf(2.0));
-    let vector = -4. * (reference_frame.x_vector * x + reference_frame.y_vector * y)
+    let vector = -4. * (reference_frame.x_vector / semi_major_axis.powf(0.5) * x + reference_frame.y_vector * y * semi_major_axis.powf(0.5)) 
+    // /semi_major_axis.powf(0.5) AND *semi_major_axis.powf(0.5) here accounts for the chain rule when calculation gradient force
         + 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
@@ -274,7 +275,31 @@ pub mod tests {
         let gradient = get_gaussian_beam_intensity_gradient(&beam, &pos1, &grf);
         assert_approx_eq!(gradient[0], -2.49605032e+13, 1e+8_f64);
         assert_approx_eq!(gradient[1], 0.0, 1e+9_f64);
-        assert_approx_eq!(gradient[2], -2.06143366e+08, 1e+6_f64);
+        assert_approx_eq!(gradient[2], -4.20876029e+08, 1e+6_f64); // Modified by Yijun on 2024/07/16
+    }
+
+    #[test] // added by Yijun on 2024/07/16 to test the gradient calculation for a beam with ellipticity
+    fn test_get_elliptic_gaussian_beam_intensity_gradient() {
+        let beam = GaussianBeam {
+            direction: Vector3::z(),
+            intersection: Vector3::new(0.0, 0.0, 0.0),
+            e_radius: 70.71067812e-6,
+            power: 100.0,
+            rayleigh_range: calculate_rayleigh_range(&1064.0e-9, &70.71067812e-6),
+            ellipticity: 0.5,
+        };
+        let pos1 = Position {
+            pos: Vector3::new(10.0e-6, 10.0e-6, 30.0e-6),
+        };
+        let grf = Frame {
+            x_vector: Vector3::x(),
+            y_vector: Vector3::y(),
+        };
+
+        let gradient = get_gaussian_beam_intensity_gradient(&beam, &pos1, &grf);
+        assert_approx_eq!(gradient[0], -2.11796152e+13, 1e+8_f64);
+        assert_approx_eq!(gradient[1], -2.82394869e+13, 1e+9_f64);
+        assert_approx_eq!(gradient[2], -4.03781010e+08, 1e+6_f64); 
     }
 
     #[test]
diff --git a/src/laser/intensity_gradient.rs b/src/laser/intensity_gradient.rs
index 395206b7..8a795fd5 100644
--- a/src/laser/intensity_gradient.rs
+++ b/src/laser/intensity_gradient.rs
@@ -173,7 +173,7 @@ pub mod tests {
             1e+5_f64
         );
     }
-    #[test]
+    #[test] // modified by Yijun on 2024/07/16
     fn test_sample_laser_intensity_gradient_again_system() {
         let mut test_world = World::new();
 
@@ -204,8 +204,8 @@ pub mod tests {
             })
             .with(beam)
             .with(Frame {
-                x_vector: Vector3::y(),
-                y_vector: Vector3::z(),
+                x_vector: Vector3::z(),
+                y_vector: Vector3::y(),
             })
             .with(DipoleLight {
                 wavelength: 1064.0e-9,
@@ -215,7 +215,7 @@ pub mod tests {
         let atom1 = test_world
             .create_entity()
             .with(Position {
-                pos: Vector3::new(20.0e-6, 20.0e-6, 20.0e-6),
+                pos: Vector3::new(20.0e-6, 25.0e-6, 20.0e-6),
             })
             .with(LaserIntensityGradientSamplers {
                 contents: [LaserIntensityGradientSampler::default();
@@ -232,9 +232,9 @@ pub mod tests {
             .expect("Entity not found!")
             .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!("Gradient: x = {}, y = {}, z = {}", sim_result_gradient[0], sim_result_gradient[1], sim_result_gradient[2]);
+        assert_approx_eq!(-1.89173404e+08, sim_result_gradient[0], 1e+5_f64);
+        assert_approx_eq!(-5.18620162e+13, sim_result_gradient[1], 1e+8_f64);
+        assert_approx_eq!(-4.14896130e+13, sim_result_gradient[2], 1e+8_f64);
     }
 }