From c82f4dc87ff7a224c1d5d0749644b68bbf0c045e Mon Sep 17 00:00:00 2001
From: Elias <elias.isak.borjesson@gmail.com>
Date: Tue, 7 Mar 2023 09:25:28 +0100
Subject: [PATCH] add orthogonal isparallel function

---
 src/Tensors.jl   |  2 +-
 src/math_ops.jl  | 29 +++++++++++++++++++++++++++++
 test/test_ops.jl | 19 +++++++++++++++++++
 3 files changed, 49 insertions(+), 1 deletion(-)

diff --git a/src/Tensors.jl b/src/Tensors.jl
index 89bfb8b4..14d5488b 100644
--- a/src/Tensors.jl
+++ b/src/Tensors.jl
@@ -20,7 +20,7 @@ export tdot, dott, dotdot
 export hessian, gradient, curl, divergence, laplace
 export @implement_gradient
 export basevec, eᵢ
-export rotate, rotation_tensor
+export rotate, rotation_tensor, orthogonal
 export tovoigt, tovoigt!, fromvoigt, tomandel, tomandel!, frommandel
 #########
 # Types #
diff --git a/src/math_ops.jl b/src/math_ops.jl
index d9b2925a..58d87e26 100644
--- a/src/math_ops.jl
+++ b/src/math_ops.jl
@@ -378,3 +378,32 @@ function rotate(x::SymmetricTensor{4}, args...)
     R = rotation_tensor(args...)
     return unsafe_symmetric(otimesu(R, R) ⊡ x ⊡ otimesu(R', R'))
 end
+
+isparallel(v1::Vec{3}, v2::Vec{3}) = isapprox(Tensors.cross(v1,v2), zero(v1), atol = 1e-14)
+isparallel(v1::Vec{2}, v2::Vec{2}) = isapprox(v1[2]*v2[1], v1[1]*v2[2], atol = 1e-14)
+
+"""
+    orthogonal(u::Vec)
+
+Return a orthogonal vector `v` to u (v ⋅ u = 0.0)
+"""
+function orthogonal(u::Vec{3}) 
+    iszero(u) && error("Cannot construct a orthogonal vector to a vector with zero length")
+    r = rand(u)
+    while isparallel(r,u) || iszero(r) #Is this needed?
+        r = rand(u)
+    end
+    q = r - u*(r⋅u)/(u⋅u)
+    return q
+end
+
+
+function orthogonal(u::Vec{2}) 
+    iszero(u) && error("Cannot construct a orthogonal vector to a vector with zero length")
+    return Vec((u[2], -u[1]))
+end
+
+function orthogonal(u::Vec{1}) 
+    error("Cannot construct a orthogonal vector for a one dimensional vector")
+    #return Vec((0.0,))
+end
\ No newline at end of file
diff --git a/test/test_ops.jl b/test/test_ops.jl
index c457fdce..cb97ab27 100644
--- a/test/test_ops.jl
+++ b/test/test_ops.jl
@@ -319,5 +319,24 @@ end
         @test isa(fromvoigt(SymmetricTensor{4,dim}, rand(num_components,num_components) .> 0.5), SymmetricTensor{4,dim,Bool})
     end
 end
+
+@testsection "isparallel" begin
+    if dim == 3 || dim == 2
+        u = rand(Vec{dim})
+        v = 2u
+        @test Tensors.isparallel(u,v) == true
+        v = orthogonal(u)
+        @test Tensors.isparallel(u,v) == false
+    end
+end
+
+@testsection "orthogonal" begin
+    if dim == 3 || dim == 2
+        u = rand(Vec{dim})
+        v = orthogonal(u)
+        @test u⋅v ≈ 0.0 atol = 1e-14
+    end
+end
+
 end # of testsection
 end # of testsection