diff --git a/src/shape.jl b/src/shape.jl
index 084a88b..daab1ae 100644
--- a/src/shape.jl
+++ b/src/shape.jl
@@ -1,4 +1,5 @@
 using StaticArrays: SVector, SMatrix
+using LinearAlgebra: eigvals, Hermitian
 
 function project_to_unit_circle(x, l)
     φ = (x + l) * π / l
@@ -24,15 +25,13 @@ function center_of_mass(args)
 
     # No need for 1/N with atan
     # If proj is (0, 0) then COM is 0 or L or -L. Here, 0 is choosen with θ = π
-    if (round(x_proj_sum[1]; digits=digits) == 0) &&
-        (round(x_proj_sum[2]; digits=digits) == 0)
+    if round(x_proj_sum[1]; digits=digits) == round(x_proj_sum[2]; digits=digits) == 0
         x_θ = π
     else
         x_θ = atan(x_proj_sum[2], x_proj_sum[1])
     end
 
-    if (round(y_proj_sum[1]; digits=digits) == 0) &&
-        (round(y_proj_sum[2]; digits=digits) == 0)
+    if round(y_proj_sum[1]; digits=digits) == round(y_proj_sum[2]; digits=digits) == 0
         y_θ = π
     else
         y_θ = atan(y_proj_sum[2], y_proj_sum[1])
@@ -53,13 +52,18 @@ function gyration_tensor(args)
 
     for p in args.particles
         c = p.c
-        x = c[1] - COM[1]
-        y = c[2] - COM[2]
+        x = restrict_coordinate(c[1] - COM[1]; l=args.l)
+        y = restrict_coordinate(c[2] - COM[2]; l=args.l)
 
         S11 += x^2
         S12 += x * y
         S22 += y^2
     end
 
-    return SMatrix{2,2}(S11, S12, S12, S22)
+    return Hermitian(SMatrix{2,2}(S11, S12, S12, S22))
+end
+
+function gyration_tensor_eigvals_ratio(args)
+    ev = eigvals(gyration_tensor(args)) # Eigenvalues are sorted
+    return ev[1] / ev[2]
 end
\ No newline at end of file