-- |
-- Module      : Main
-- Copyright   : (c) Justus Sagemüller 2018
-- License     : GPL v3
-- 
-- Maintainer  : (@) sagemueller $ geo.uni-koeln.de
-- Stability   : experimental
-- Portability : portable
-- 

{-# LANGUAGE OverloadedLists, TypeFamilies, FlexibleContexts, UndecidableInstances #-}
{-# LANGUAGE TypeOperators #-}

module Main where

import qualified Text.Show.Pragmatic as SP
import Math.Rotations.Class

import Data.Manifold

import Test.Tasty
import Test.Tasty.HUnit
import qualified Test.Tasty.QuickCheck as QC
import Test.Tasty.QuickCheck ((==>))


main = defaultMain tests

tests :: TestTree
tests = testGroup "Tests"
 [ testGroup "Rotation matrices"
    [ testCase "180° x-axis" $ rotmatrixForAxis xAxis (rad π)
                      @?≈ [ [1,0,0], [0,-1,0], [0,0,-1] ]
    , testCase " 90° x-axis" $ rotmatrixForAxis xAxis (rad π'₂)
                      @?≈ [ [1,0,0], [0,0,-1], [0,1,0] ]
    , testCase "-90° x-axis" $ rotmatrixForAxis xAxis (rad $ -π'₂)
                      @?≈ [ [1,0,0], [0,0,1], [0,-1,0] ]
    , testCase " 45° x-axis" $ rotmatrixForAxis xAxis (rad π'₄)
                      @?≈ [ [1,0,0], [0,sqrt 2/2,-sqrt 2/2], [0,sqrt 2/2,sqrt 2/2] ]
    , testCase "180° y-axis" $ rotmatrixForAxis yAxis (rad π)
                      @?≈ [ [-1,0,0], [0,1,0], [0,0,-1] ]
    , testCase " 90° y-axis" $ rotmatrixForAxis yAxis (rad π'₂)
                      @?≈ [ [0,0,1], [0,1,0], [-1,0,0] ]
    , testCase "-90° y-axis" $ rotmatrixForAxis yAxis (rad $ -π'₂)
                      @?≈ [ [0,0,-1], [0,1,0], [1,0,0] ]
    , testCase " 45° y-axis" $ rotmatrixForAxis yAxis (rad π'₄)
                      @?≈ [ [sqrt 2/2,0,sqrt 2/2], [0,1,0], [-sqrt 2/2,0,sqrt 2/2] ]
    , testCase "180° z-axis" $ rotmatrixForAxis zAxis (rad π)
                      @?≈ [ [-1,0,0], [0,-1,0], [0,0,1] ]
    , testCase " 90° z-axis" $ rotmatrixForAxis zAxis (rad π'₂)
                      @?≈ [ [0,-1,0], [1,0,0], [0,0,1] ]
    , testCase "-90° z-axis" $ rotmatrixForAxis zAxis (rad $ -π'₂)
                      @?≈ [ [0,1,0], [-1,0,0], [0,0,1] ]
    , testCase " 45° z-axis" $ rotmatrixForAxis zAxis (rad π'₄)
                      @?≈ [ [sqrt 2/2,-sqrt 2/2,0], [sqrt 2/2,sqrt 2/2,0], [0,0,1] ]
    ]
 , testGroup "Euler angles"
    [ testCase "180° x-axis" $ eulerAnglesZYZForMatrix (rotmatrixForAxis xAxis (rad π))
                      @?≈ [ π³₄,  π ,-π'₄]
    , testCase " 90° x-axis" $ eulerAnglesZYZForMatrix (rotmatrixForAxis xAxis (rad π'₂))
                      @?≈ [ π'₂, π'₂,-π'₂]
    , testCase "-90° x-axis" $ eulerAnglesZYZForMatrix (rotmatrixForAxis xAxis (rad $ -π'₂))
                      @?≈ [-π'₂, π'₂, π'₂]
    , testCase " 45° x-axis" $ eulerAnglesZYZForMatrix (rotmatrixForAxis xAxis (rad π'₄))
                      @?≈ [ π'₂, π'₄,-π'₂]
    , testCase "180° y-axis" $ eulerAnglesZYZForMatrix (rotmatrixForAxis yAxis (rad π))
                      @?≈ [ π'₄, π  , π'₄]
    , testCase " 90° y-axis" $ eulerAnglesZYZForMatrix (rotmatrixForAxis yAxis (rad π'₂))
                      @?≈ [ 0  , π'₂, 0  ]
    , testCase "-90° y-axis" $ eulerAnglesZYZForMatrix (rotmatrixForAxis yAxis (rad $ -π'₂))
                      @?≈ [-π  , π'₂, π  ]
    , testCase " 45° y-axis" $ eulerAnglesZYZForMatrix (rotmatrixForAxis yAxis (rad π'₄))
                      @?≈ [ 0  , π'₄, 0  ]
    , testCase "180° z-axis" $ eulerAnglesZYZForMatrix (rotmatrixForAxis zAxis (rad π))
                      @?≈ [ π  ,  0 , 0  ]
    , testCase " 90° z-axis" $ eulerAnglesZYZForMatrix (rotmatrixForAxis zAxis (rad π'₂))
                      @?≈ [ π  ,  0 ,-π'₂]
    , testCase "-90° z-axis" $ eulerAnglesZYZForMatrix (rotmatrixForAxis zAxis (rad $ -π'₂))
                      @?≈ [ π  ,  0 , π'₂]
    , testCase " 45° z-axis" $ eulerAnglesZYZForMatrix (rotmatrixForAxis zAxis (rad π'₄))
                      @?≈ [ π  , 0  ,-π³₄]
    , QC.testProperty "Arbitrary rotation-matrix reconstruction"
         $ \ax α -> let mat = rotmatrixForAxis ax α
                    in rotmatrixForEulerAnglesZYZ (eulerAnglesZYZForMatrix mat) ≈ mat
    ]
 , testGroup "Concrete 180° rotations"
    [ testCase "x around x" $ (180° xAxis) (π'₂ ⦺ 0  ) @?≈ π'₂ ⦺  0
    , testCase "..full rot" $ (180° xAxis)                (π'₂ ⦺  0  )@?≈ π'₂ ⦺ 0
    , testCase "x around y" $ (180° yAxis) (π'₂ ⦺ 0  ) @?≈ π'₂ ⦺  π
    , testCase "..full rot" $ (180° yAxis)                (π'₂ ⦺  π  )@?≈ π'₂ ⦺ 0
    , testCase "x around z" $ (180° zAxis) (π'₂ ⦺ 0  ) @?≈ π'₂ ⦺  π
    , testCase "..full rot" $ (180° zAxis)                (π'₂ ⦺  π  )@?≈ π'₂ ⦺ 0
    , testCase "y around x" $ (180° xAxis) (π'₂ ⦺ π'₂) @?≈ π'₂ ⦺ -π'₂
    , testCase "..full rot" $ (180° xAxis)                (π'₂ ⦺ -π'₂)@?≈ π'₂ ⦺ π'₂
    , testCase "y around y" $ (180° yAxis) (π'₂ ⦺ π'₂) @?≈ π'₂ ⦺  π'₂
    , testCase "..full rot" $ (180° yAxis)                (π'₂ ⦺  π'₂)@?≈ π'₂ ⦺ π'₂
    , testCase "y around z" $ (180° zAxis) (π'₂ ⦺ π'₂) @?≈ π'₂ ⦺ -π'₂
    , testCase "..full rot" $ (180° zAxis)                (π'₂ ⦺ -π'₂)@?≈ π'₂ ⦺ π'₂
    , testCase "z around x" $ (180° xAxis) (0   ⦺ 0  ) @?≈ π   ⦺  0
    , testCase "..full rot" $ (180° xAxis)                (π   ⦺  0  )@?≈ 0 ⦺ 0
    , testCase "z around y" $ (180° yAxis) (0   ⦺ 0  ) @?≈ π   ⦺  0
    , testCase "..full rot" $ (180° yAxis)                (π   ⦺  0  )@?≈ 0 ⦺ 0
    , testCase "z around z" $ (180° zAxis) (0   ⦺ 0  ) @?≈ 0   ⦺  0
    , testCase "..full rot" $ (180° zAxis)                (0   ⦺  0  )@?≈ 0 ⦺ 0
    ]
 , testGroup "Concrete 90° rotations"
    [ testCase "x around x" $ (90° xAxis) (π'₂ ⦺ 0  ) @?≈ π'₂ ⦺ 0
    , testCase "x around y" $ (90° yAxis) (π'₂ ⦺ 0  ) @?≈ π   ⦺ 0
    , testCase "x around z" $ (90° zAxis) (π'₂ ⦺ 0  ) @?≈ π'₂ ⦺ π'₂
    , testCase "y around x" $ (90° xAxis) (π'₂ ⦺ π'₂) @?≈ 0   ⦺ 0
    , testCase "y around y" $ (90° yAxis) (π'₂ ⦺ π'₂) @?≈ π'₂ ⦺ π'₂
    , testCase "y around z" $ (90° zAxis) (π'₂ ⦺ π'₂) @?≈ π'₂ ⦺ π
    , testCase "z around x" $ (90° xAxis) (0   ⦺ 0  ) @?≈ π'₂ ⦺ (-π'₂)
    , testCase "z around y" $ (90° yAxis) (0   ⦺ 0  ) @?≈ π'₂ ⦺ 0
    , testCase "z around z" $ (90° zAxis) (0   ⦺ 0  ) @?≈ 0   ⦺ π'₂
    
    , testCase "around y, over x" $ (90° yAxis) (π'₄ ⦺  0  ) @?≈ π³₄ ⦺ 0
    , testCase "around z, over x" $ (90° zAxis) (π'₂ ⦺ -π'₄) @?≈ π'₂ ⦺ π'₄
    , testCase "around x, over y" $ (90° xAxis) (π³₄ ⦺  π'₂) @?≈ π'₄ ⦺ π'₂
    , testCase "around z, over y" $ (90° zAxis) (π'₂ ⦺  π'₄) @?≈ π'₂ ⦺ π³₄
    , testCase "around x, over z" $ (90° xAxis) (π'₄ ⦺  π'₂) @?≈ π'₄ ⦺ -π'₂
    , testCase "around y, over z" $ (90° yAxis) (π'₄ ⦺  π  ) @?≈ π'₄ ⦺ 0
    ]
 , testGroup "Concrete 45° rotations"
    [ testCase "x around x" $ (45° xAxis) (π'₂ ⦺ 0  ) @?≈ π'₂ ⦺  0
    , testCase "x around y" $ (45° yAxis) (π'₂ ⦺ 0  ) @?≈ π³₄ ⦺  0
    , testCase "x around z" $ (45° zAxis) (π'₂ ⦺ 0  ) @?≈ π'₂ ⦺  π'₄
    , testCase "y around x" $ (45° xAxis) (π'₂ ⦺ π'₂) @?≈ π'₄ ⦺  π'₂
    , testCase "y around y" $ (45° yAxis) (π'₂ ⦺ π'₂) @?≈ π'₂ ⦺  π'₂
    , testCase "y around z" $ (45° zAxis) (π'₂ ⦺ π'₂) @?≈ π'₂ ⦺  π³₄
    , testCase "z around x" $ (45° xAxis) (0   ⦺ 0  ) @?≈ π'₄ ⦺ -π'₂
    , testCase "z around y" $ (45° yAxis) (0   ⦺ 0  ) @?≈ π'₄ ⦺  0
    , testCase "z around z" $ (45° zAxis) (0   ⦺ 0  ) @?≈ 0   ⦺  π'₂
    ]
 , testGroup "Reversability"
    [ QC.testProperty "Arbitrary axis, angles and points"
           $ \ax ψ p -> rotateAboutThenUndo ax ψ p ≈ p
    ]
 ]
 where π = pi
       π'₂ = pi/2
       π'₄ = pi/4
       π³₄ = 3*pi/4



infix 4 ≈
class AEq e where
  (≈) :: e -> e -> Bool

instance AEq Double where
  x ≈ y = abs (x-y) < 1e-9
instance AEq S¹ where
  S¹Polar φ ≈ S¹Polar ϕ
   | φ > pi/2, ϕ < -pi/2  = S¹Polar (φ - 2*pi) ≈ S¹Polar ϕ
   | ϕ > pi/2, φ < -pi/2  = S¹Polar φ ≈ S¹Polar (ϕ - 2*pi)
   | otherwise            = φ ≈ ϕ
instance AEq S² where
  S²Polar θ φ ≈ S²Polar ϑ ϕ
   | φ > pi/2, ϕ < -pi/2  = S²Polar θ (φ - 2*pi) ≈ S²Polar ϑ ϕ
   | ϕ > pi/2, φ < -pi/2  = S²Polar θ φ ≈ S²Polar ϑ (ϕ - 2*pi)
   | otherwise            = θ ≈ ϑ && abs (φ - ϕ) * sin θ < 1e-9
instance AEq a => AEq [a] where
  [] ≈ [] = True
  x:xs ≈ y:ys = x ≈ y && xs ≈ ys
  _ ≈ _ = False

                  
infix 1 @?≈       
(@?≈) :: (AEq e, Show e) => e -> e -> Assertion
a@?≈b
 | a≈b        = return ()
 | otherwise  = assertFailure $ "Expected "++show b++", but got "++show a

rad :: Double -> S¹
rad = S¹Polar

infix 5 ⦺
(⦺) :: Double -> Double -> S²
(⦺) = S²Polar

rotateAboutThenUndo :: ℝP² -> S¹ -> S² -> S²
rotateAboutThenUndo ax ψ@(S¹Polar w) p = rotateAbout ax (S¹Polar $ -w) $ rotateAbout ax ψ p