%include lhs2TeX.fmt
%align 33

\begin{code}
module Boson (Boson(Boson),CBoson(CBoson),ncVop) where
import Tcsa
import qualified List
import Complex

-- DECREASING POSITIVE modes (with possible multiplicities)
type Mode                       = Int

-- p=n/r +- qr/2 (vertex op. momentum)
-- the modes are the multiplicities, starting from mode a_(-1)
newtype CBoson                  =  CBoson (Rational,[Mode]) deriving (Eq,Show)
newtype Boson                   =  Boson (CBoson,CBoson) deriving (Eq,Show)

-- put (m<=) instead of (m<) to get the fermionic version
partitions                      :: [[[Int]]]
partitions                      =  [[]]:[[n]:concat [map (m:) $ dropWhile ((m<).head) pars
                                                     | (m,pars) <- zip [n-1,n-2..1] (tail partitions)]
                                         | n <- [1..]]

incrementalBase                 :: [[([Mode],[Mode])]]
incrementalBase                 =  map (concat . map pairs) (diagSquare countedParts)
    where pairs (cls,crs) = [(cl,cr) | cl <- cls, cr <- crs]
          countedParts    = map (map (counted 0 . reverse)) partitions
          diagSquare cs   = [zip (reverse $ take n cs) cs | n <- [1..]]

inModule                        :: [[([Mode],[Mode])]] -> (Rational,Rational) -> [[Boson]]
levels `inModule` (pl,pr)       =  map (map attach) levels
    where attach (cl,cr) = Boson (CBoson (pl,cl),CBoson (pr,cr))

-- call with prev=0
counted                         :: Mode -> [Mode] -> [Mode]
counted prev []                 =  []
counted prev (m:ms)             =  replicate (m-prev-1) 0 ++ 1+length same:counted m rest
    where (same,rest) = span (m==) ms

allTowers                       :: Int -> Rational -> [[[Boson]]]
allTowers q r                   =  incrementalBase `inModule` (p2,-p2):[bothBases (n/r) | n <- [1..]]
    where bothBases p1 = zipWith (++) (incrementalBase `inModule` (p1+p2,p1-p2))
                                      (incrementalBase `inModule` (-p1+p2,-p1-p2))
          p2           = fromIntegral q * r/2

-- preconditions: 1. on the next list the element with the same index is not less than this
--                2. we have infinitely many infinite lists
mergeUp                         :: [[[Boson]]] -> [[Boson]]
mergeUp towers                  =  List.sortBy scaleCmp heads ++ mergeUp (dropFirst (length heads) towers)
    where heads                = takeWhile listTest (map head towers)
          scaleCmp (a:_) (b:_) = compare (scaleDimension a) (scaleDimension b)
          scaleCmp _     _     = EQ
          listTest []          = True
          listTest (a:_)       = (scaleDimension $ head $ head towers !! 1) > scaleDimension a
          dropFirst n list     = map tail (take n list) ++ (drop n list)

instance Chiral CBoson where
    weight (CBoson (p,c))       =  p^2/2 + fromIntegral (sum $ zipWith (*) c [1..])
    norm' (CBoson (_,ms))       =  sqrt $ fromIntegral $ product $ [n^exp*factorial!!exp | (n,exp) <- zip [1..] ms]

instance BaseVector Boson where
    scaleDimension (Boson (l,r))
                                =  weight l + weight r
    spinOf (Boson (l,r))        =  truncate (weight l - weight r)
    norm (Boson (l,r))          =  norm' l * norm' r
    baseLevels params           =  mergeUp $ allTowers (topCharge params) (radius params)
    energy s                    =  fromRational $ scaleDimension s - 1/12
    showParams _ params         =  unlines ["cut: " ++ show (fromRational $ cut params),
                                            "radius: " ++ show (fromRational $ radius params),
                                            "topcharge: " ++ show (topCharge params),
                                            "spin: " ++ show (spin params),
                                            "delta: " ++ show (delta params)]

factorial                       :: [Int]
factorial                       =  1:1:zipWith (*) [2..13] (tail factorial) -- avoid overflow

vopMode                         :: Rational -> Int -> Int -> Int -> Rational
vopMode p n l r                 =  sum [block (-p) r k * block p l (l-r+k) *
                                        fromIntegral (n^(r-k) * factorial!!(r-k)) | k <- [max 0 (r-l)..r]]
    where block p top bottom = p^bottom * fromIntegral (binomial top bottom)
          binomial n k       = round $ fromIntegral (factorial!!n) / fromIntegral (factorial!!k * factorial!!(n-k))

vop                             :: Rational -> CBoson -> CBoson -> Rational
vop p (CBoson (pOut,cOut)) (CBoson (pIn,cIn))
    | pOut - pIn == p           =  product $ zipWith3 (vopMode p) [1..nMax] (cOut++repeat 0) (cIn++repeat 0)
    | otherwise                 =  0
    where nMax = max (length cOut) (length cIn)

ncVop'                          :: Rational -> Boson -> Boson -> Rational
ncVop' p (Boson (lOut,rOut)) (Boson (lIn,rIn))
                                =  vop p lOut lIn * vop p rOut rIn

ncVop                           :: Params -> Operator Boson
ncVop params bra ket            =  prefac params * fromRational (ncVop' p bra ket + ncVop' (-p) bra ket) / (2 * norm bra * norm ket) :+ 0
    where p = 1/radius params
\end{code}

% Local Variables:
% mode: latex
% mode: auto-fill
% eval: (local-set-key "\C-C\C-c" 'compile)
% TeX-master: "Fermion"
% mmm-classes: literate-haskell
% End: