From f446c7fbc3c2134ea30d0d92f95099144ae35aba Mon Sep 17 00:00:00 2001
From: Mo8it <76752051+Mo8it@users.noreply.github.com>
Date: Fri, 7 Jan 2022 00:28:40 +0100
Subject: [PATCH] Added elliptical distance

---
 src/Shape.jl | 32 +++++++++++++++++++++++++++++---
 1 file changed, 29 insertions(+), 3 deletions(-)

diff --git a/src/Shape.jl b/src/Shape.jl
index 1ede4af..854fe8c 100644
--- a/src/Shape.jl
+++ b/src/Shape.jl
@@ -1,20 +1,23 @@
 module Shape
 
-export center_of_mass, gyration_tensor_eigvals_ratio
+export center_of_mass,
+    gyration_tensor_eigvals_ratio, gyration_tensor_eigvecs, elliptical_distance
 
 using StaticArrays: SVector, SMatrix
-using LinearAlgebra: eigvals, Hermitian
+using LinearAlgebra: eigvals, eigvecs, Hermitian, dot
 
 using ..ReCo: Particle, restrict_coordinate, restrict_coordinates
 
 function project_to_unit_circle(x::Float64, half_box_len::Float64)
     φ = (x + half_box_len) * π / half_box_len
     si, co = sincos(φ)
+
     return SVector(co, si)
 end
 
 function project_back_from_unit_circle(θ::T, half_box_len::Float64) where {T<:Real}
     x = θ * half_box_len / π - half_box_len
+
     return restrict_coordinate(x, half_box_len)
 end
 
@@ -87,8 +90,31 @@ function gyration_tensor(particles::Vector{Particle}, half_box_len::Float64)
 end
 
 function gyration_tensor_eigvals_ratio(particles::Vector{Particle}, half_box_len::Float64)
-    ev = eigvals(gyration_tensor(particles, half_box_len)) # Eigenvalues are sorted
+    g_tensor = gyration_tensor(particles, half_box_len)
+    ev = eigvals(g_tensor) # Eigenvalues are sorted
     return ev[1] / ev[2]
 end
 
+function gyration_tensor_eigvecs(particles::Vector{Particle}, half_box_len::Float64)
+    g_tensor = gyration_tensor(particles, half_box_len)
+    eig_vecs = eigvecs(g_tensor)
+
+    v1 = eig_vecs[:, 1]
+    v2 = eig_vecs[:, 2]
+
+    return (v1, v2)
+end
+
+function elliptical_distance(
+    particle::Particle,
+    gyration_tensor_eigvec_to_smaller_eigval::SVector{2,Float64},
+    gyration_tensor_eigvec_to_bigger_eigval::SVector{2,Float64},
+    goal_gyration_tensor_eigvals_ratio::Float64,
+)
+    cx′ = dot(particle.c, gyration_tensor_eigvec_to_bigger_eigval)
+    cy′ = dot(particle.c, gyration_tensor_eigvec_to_smaller_eigval)
+
+    return cx′^2 + (cy′ / goal_gyration_tensor_eigvals_ratio)^2
+end
+
 end # module
\ No newline at end of file