From 946e646594396c6fb37e1af5d7570821a366ae63 Mon Sep 17 00:00:00 2001
From: Mo8it <mo8it@protonmail.com>
Date: Sun, 30 Jan 2022 14:45:22 +0100
Subject: [PATCH] Added eigvals_ratio for centers

---
 src/Shape.jl | 79 +++++++++++++++++++++++++++++++++++++---------------
 1 file changed, 56 insertions(+), 23 deletions(-)

diff --git a/src/Shape.jl b/src/Shape.jl
index 6068e0e..b01cba4 100644
--- a/src/Shape.jl
+++ b/src/Shape.jl
@@ -22,8 +22,8 @@ function project_back_from_unit_circle(θ::Real, half_box_len::Real)
 end
 
 function center_of_mass_from_proj_sums(
-    x_proj_sum::SVector{2,R}, y_proj_sum::SVector{2,R}, half_box_len::R
-) where {R<:Real}
+    x_proj_sum::SVector{2,<:Real}, y_proj_sum::SVector{2,<:Real}, half_box_len::Real
+)
     # Prevent for example atan(1e-16, 1e-15) != 0 with rounding
     digits = 5
 
@@ -47,9 +47,7 @@ function center_of_mass_from_proj_sums(
     return SVector(COM_x, COM_y)
 end
 
-function center_of_mass(
-    centers::AbstractVector{SVector{2,R}}, half_box_len::R
-) where {R<:Real}
+function center_of_mass(centers::AbstractVector{SVector{2,<:Real}}, half_box_len::Real)
     x_proj_sum = SVector(0.0, 0.0)
     y_proj_sum = SVector(0.0, 0.0)
 
@@ -74,8 +72,8 @@ function center_of_mass(particles::AbstractVector{Particle}, half_box_len::Real)
 end
 
 function gyration_tensor(
-    particles::AbstractVector{Particle}, half_box_len::R, COM::SVector{2,R}
-) where {R<:Real}
+    particles::AbstractVector{Particle}, half_box_len::Real, COM::SVector{2,<:Real}
+)
     S11 = 0.0
     S12 = 0.0
     S22 = 0.0
@@ -91,23 +89,58 @@ function gyration_tensor(
     return LA.Hermitian(SMatrix{2,2}(S11, S12, S12, S22))
 end
 
-function gyration_tensor(particles::AbstractVector{Particle}, half_box_len::Real)
-    COM = center_of_mass(particles, half_box_len)
+function gyration_tensor(
+    centers::AbstractVector{SVector{2,<:Real}}, half_box_len::Real, COM::SVector{2,<:Real}
+)
+    S11 = 0.0
+    S12 = 0.0
+    S22 = 0.0
+
+    for c in centers
+        shifted_c = ReCo.restrict_coordinates(c - COM, half_box_len)
+
+        S11 += shifted_c[1]^2
+        S12 += shifted_c[1] * shifted_c[2]
+        S22 += shifted_c[2]^2
+    end
+
+    return LA.Hermitian(SMatrix{2,2}(S11, S12, S12, S22))
+end
+
+function gyration_tensor(
+    particles_or_centers::Union{AbstractVector{Particle},AbstractVector{SVector{2,<:Real}}},
+    half_box_len::Real,
+)
+    COM = center_of_mass(particles_or_centers, half_box_len)
 
     return gyration_tensor(particles, half_box_len, COM)
 end
 
-function gyration_tensor_eigvals_ratio(
-    particles::AbstractVector{Particle}, half_box_len::Real
-)
-    g_tensor = gyration_tensor(particles, half_box_len)
-    ev = LA.eigvals(g_tensor) # Eigenvalues are sorted
+function eigvals_ratio(matrix)
+    ev = LA.eigvals(matrix) # Eigenvalues are sorted
     return abs(ev[1] / ev[2])
 end
 
+function gyration_tensor_eigvals_ratio(
+    particles_or_centers::Union{AbstractVector{Particle},AbstractVector{SVector{2,<:Real}}},
+    half_box_len::Real,
+)
+    g_tensor = gyration_tensor(particles_or_centers, half_box_len)
+    return eigvals_ratio(g_tensor)
+end
+
+function gyration_tensor_eigvals_ratio(
+    particles_or_centers::Union{AbstractVector{Particle},AbstractVector{SVector{2,<:Real}}},
+    half_box_len::Real,
+    COM::SVector{2,<:Real},
+)
+    g_tensor = gyration_tensor(particles_or_centers, half_box_len, COM)
+    return eigvals_ratio(g_tensor)
+end
+
 function gyration_tensor_eigvecs(
-    particles::AbstractVector{Particle}, half_box_len::R, COM::SVector{2,R}
-) where {R<:Real}
+    particles::AbstractVector{Particle}, half_box_len::Real, COM::SVector{2,<:Real}
+)
     g_tensor = gyration_tensor(particles, half_box_len, COM)
     eig_vecs = LA.eigvecs(g_tensor)
 
@@ -118,13 +151,13 @@ function gyration_tensor_eigvecs(
 end
 
 function elliptical_distance(
-    v::SVector{2,R},
-    COM::SVector{2,R},
-    gyration_tensor_eigvec_to_smaller_eigval::SVector{2,R},
-    gyration_tensor_eigvec_to_bigger_eigval::SVector{2,R},
-    elliptical_a_b_ratio::R,
-    half_box_len::R,
-) where {R<:Real}
+    v::SVector{2,<:Real},
+    COM::SVector{2,<:Real},
+    gyration_tensor_eigvec_to_smaller_eigval::SVector{2,<:Real},
+    gyration_tensor_eigvec_to_bigger_eigval::SVector{2,<:Real},
+    elliptical_a_b_ratio::Real,
+    half_box_len::Real,
+)
     v′ = ReCo.restrict_coordinates(v - COM, half_box_len)
 
     x = LA.dot(v′, gyration_tensor_eigvec_to_bigger_eigval)