{-# LANGUAGE FunctionalDependencies
           , FlexibleInstances      #-}
---------------------------------------------------------------------------------------------------
-- |
-- Module      : SU3
-- Description : Showcase
-- Copyright   : (c) Felix Springer, 2019
-- License     : BSD3
-- Maintainer  : felixspringer149@gmail.com
-- Stability   : experimental
-- Portability : POSIX
--
-- This is a definition of the SU(3), useful to load into GHCi.
--
---------------------------------------------------------------------------------------------------

module LieExample.SU3 ( CD
                      , Matrix (..)
                      , (|||)
                      , (|-|)
                      ) where

import Lie.LieAlgebra

import Data.Complex
import Data.List (transpose)

-- | Complex Number (real and imaginary part have Real coefficients)
type CD = Complex Double

-- | Matrix in ℝ^(3⨯3)
data Matrix = M CD CD CD
                CD CD CD
                CD CD CD
                deriving Eq

instance Show Matrix where
  show (M a b c
          d e f
          g h i) = "\n" ++ (foldl (++) "" $ map showRow $ eqLen entries)
    where entries = [ [ show a , show b , show c ]
                    , [ show d , show e , show f ]
                    , [ show g , show h , show i ]
                    ]
          showRow :: [String] -> String
          showRow ls = (foldl conn "" ls) ++ "\n"
            where conn a b = a ++ " " ++ b
          eqLen :: [[String]] -> [[String]]
          eqLen mat = transpose $ map toMaxLen matT
            where matT = transpose mat
          toMaxLen :: [String] -> [String]
          toMaxLen ls = map (toLen n) ls
            where n = maximum $ map length ls
          toLen :: Int -> String -> String
          toLen n s = [' ' | i <- [(length s)..(n)]] ++ s

-- | Matrix Multiplication
(|||) :: Matrix -> Matrix -> Matrix
(|||) (M a b c
         d e f
         g h i)
                (M j k l
                   m n o
                   p q r) = M (a*j + b*m + c*p) (a*k + b*n + c*q) (a*l + b*o + c*r)
                              (d*j + e*m + f*p) (d*k + e*n + f*q) (d*l + e*o + f*r)
                              (g*j + h*m + i*p) (g*k + h*n + i*q) (g*l + h*o + i*r)

gellMannRepr :: [Matrix]
gellMannRepr = map ((|*|) (1/2))
               [ M (  0  :+ 0) (  1  :+ 0) (  0  :+ 0)
                   (  1  :+ 0) (  0  :+ 0) (  0  :+ 0)
                   (  0  :+ 0) (  0  :+ 0) (  0  :+ 0)
               , M (  0  :+ 0) (0 :+ (-1)) (  0  :+ 0)
                   (  0  :+ 1) (  0  :+ 0) (  0  :+ 0)
                   (  0  :+ 0) (  0  :+ 0) (  0  :+ 0)
               , M (  1  :+ 0) (  0  :+ 0) (  0  :+ 0)
                   (  0  :+ 0) ((-1) :+ 0) (  0  :+ 0)
                   (  0  :+ 0) (  0  :+ 0) (  0  :+ 0)
               , M (  0  :+ 0) (  0  :+ 0) (  1  :+ 0)
                   (  0  :+ 0) (  0  :+ 0) (  0  :+ 0)
                   (  1  :+ 0) (  0  :+ 0) (  0  :+ 0)
               , M (  0  :+ 0) (  0  :+ 0) (0 :+ (-1))
                   (  0  :+ 1) (  0  :+ 0) (  0  :+ 0)
                   (  0  :+ 0) (  0  :+ 0) (  0  :+ 0)
               , M (  0  :+ 0) (  0  :+ 0) (  0  :+ 0)
                   (  0  :+ 0) (  0  :+ 0) (  1  :+ 0)
                   (  0  :+ 0) (  1  :+ 0) (  0  :+ 0)
               , M (  0  :+ 0) (  0  :+ 0) (  0  :+ 0)
                   (  0  :+ 0) (  0  :+ 0) (0 :+ (-1))
                   (  0  :+ 0) (  0  :+ 1) (  0  :+ 0)
               , (1/(sqrt 3)) |*|
                 M (  1  :+ 0) (  0  :+ 0) (  0  :+ 0)
                   (  0  :+ 0) (  1  :+ 0) (  0  :+ 0)
                   (  0  :+ 0) (  0  :+ 0) ((-2) :+ 0)
               ]

instance LieAlgebra Matrix CD where
  (|+|) (M a b c
           d e f
           g h i)
                  (M j k l
                     m n o
                     p q r) = M (a+j) (b+k) (c+l)
                                (d+m) (e+n) (f+o)
                                (g+p) (h+q) (i+r)

  (|*|) alpha (M a b c
                 d e f
                 g h i) = M (alpha*a) (alpha*b) (alpha*c)
                            (alpha*d) (alpha*e) (alpha*f)
                            (alpha*g) (alpha*h) (alpha*i)

  (|.|) a b = (a ||| b) |+| (((-1) :+ 0) |*| (b ||| a))

  basis = gellMannRepr

  trace f = sumDiag $ f one
    where one :: Matrix
          one = M (1 :+ 0) (0 :+ 0) (0 :+ 0)
                  (0 :+ 0) (1 :+ 0) (0 :+ 0)
                  (0 :+ 0) (0 :+ 0) (1 :+ 0)
          sumDiag :: Matrix -> CD
          sumDiag (M a _ _
                     _ b _
                     _ _ c) = a + b + c


-- | Subtraction, which is Inverse to Vectorspace Addition
(|-|) :: Matrix -> Matrix -> Matrix
(|-|) a b = a |+| (inv b)
  where inv :: Matrix -> Matrix
        inv a = zero |+| (((-1) :+ 0) |*| a)
          where zero = M (0 :+ 0) (0 :+ 0) (0 :+ 0)
                         (0 :+ 0) (0 :+ 0) (0 :+ 0)
                         (0 :+ 0) (0 :+ 0) (0 :+ 0)