{-# LANGUAGE CPP #-}
module Data.Matrix.SymmetryOperationsSymbols.RotationCase (
nFoldRotationCase,
) where
import Data.Ratio
import Data.List (transpose)
import Data.Matrix hiding (transpose)
import Data.Matrix.SymmetryOperationsSymbols.Solve
import Data.Matrix.SymmetryOperationsSymbols.Common
import Data.Matrix.AsXYZ
import Data.Matrix.SymmetryOperationsSymbols.SymmetryOperation
#if MIN_VERSION_base(4,11,0)
import Control.Monad.Fail (MonadFail)
#endif
#if MIN_VERSION_base(4,11,0)
nFoldRotationCase :: (Monad m, MonadFail m, Integral a) => Matrix (Ratio a) -> m (SymmetryOperation a)
#else
nFoldRotationCase :: (Monad m, Integral a) => Matrix (Ratio a) -> m (SymmetryOperation a)
#endif
nFoldRotationCase m = arrange m <$> solvingEquation m
arrange :: Integral a => Matrix (Ratio a) -> [Ratio a] -> SymmetryOperation a
arrange m ans
| rt == 2 && noScrew = TwoFoldRotation { axis = location }
| rt == 2 = TwoFoldScrew { axis = location, vector = (sa,sb,sc) }
| noScrew = NFoldRotation { nFold = n, sense = sen, axis = location }
| otherwise = NFoldScrew { nFold = n, sense = sen, axis = location, vector = (sa,sb,sc) }
where
noScrew = all (==0) screwPart
rt = rotationType m
n | rt == 3 = ThreeFold
| rt == 4 = FourFold
| rt == 6 = SixFold
sen = case senseOf m of
"+" -> Positive
"-" -> Negative
_ -> error ""
s = if null . senseOf $ m then " " else senseOf m
sym = " " ++ show rt ++ s
screwPart@[sa,sb,sc] = toList (wg m)
location = locationOf m <|> fromList 3 1 ans
testMat = fromList 3 4 [ 0,0,-1,0, -1,0,0,1%2, 0,1,0,0 ]
rotationType :: (Num a,Eq a,Integral b) => Matrix a -> b
rotationType m
| tr == (-1) = 2
| tr == 0 = 3
| tr == 1 = 4
| tr == 2 = 6
where tr = trace $ rotPart m
pow :: (Num a, Integral b) => Matrix a -> b -> Matrix a
pow m 0 = identity 3
pow m n = foldl1 multStd $ replicate (fromIntegral n) (rotPart m)
t :: (Num b, Eq b) => Matrix b -> Matrix b
t mat = f [ pow mat n | n <- [0..(pred r)] ]
where
r = rotationType mat
w = transPart mat
g n m = multStd m n
f l = foldl1 (elementwise (+)) $ map (g w) l
wg :: (Fractional b, Eq b) => Matrix b -> Matrix b
wg mat = fmap (/ fromIntegral r) (t mat)
where
r = rotationType mat
wl :: (Fractional b, Eq b) => Matrix b -> Matrix b
wl mat = elementwise (-) (transPart mat) (wg mat)
solvingEquation' :: Integral a => Matrix (Ratio a) -> Maybe (Matrix (Ratio a))
solvingEquation' mat = solve (iw mat) (wl mat)
solvingEquation'' :: Integral a => Matrix (Ratio a) -> Maybe [Ratio a]
solvingEquation'' mat = do
sol <- solvingEquation' mat
adjustAnswerOnAxis mat (toList sol)
#if MIN_VERSION_base(4,11,0)
solvingEquation :: (Monad m, MonadFail m, Integral a) => Matrix (Ratio a) -> m [Ratio a]
#else
solvingEquation :: (Monad m, Integral a) => Matrix (Ratio a) -> m [Ratio a]
#endif
solvingEquation mat = case solvingEquation'' mat of
Nothing -> fail "<Rot> when calculate equation."
Just a -> return a
check mat sol = if toList (multStd (iw mat) (vec sol)) == toList (wl mat)
then sol
else error "General solution error."
where
vec l = fromList (length l) 1 l