--------------------------------------------------------------------------------
module AC where

import Circuit
import Array
import List
import Matrix.Complex.LU
import Complex 

---- Y node admittance matrix
---- v vector of unknown voltages
---- i equivalent current source vector

--- Build matrix Y of Yv = i
mk_Y :: [([Int], [Double])] -> [Double] -> Array (Int, Int) (Complex Double)       
mk_Y (e:es) param = array ((1, 1), (n, n)) 
          ([((i, j), z i j) | i <- [1..n], j <- [1..n]])
            where z i j | i <= m && j <= m 
                             = yy_1 i j (admList (e:es))
                        | i > m || j > m 
                             = yy_2 i j (srcList (e:es))
                        | otherwise = 0
                  
                  yy_2 _ _ [] = 0
                  yy_2 i j (e:es) | typeDev (fst e) == tVAC    
                                    = vsrc i j (e:es) + yy_2 i j es
                                  | otherwise = yy_2 i j es                              
                  yy_1 _ _ [] = 0
                  yy_1 i j (e:es) | typeDev (fst e) == tRES
                                    = adm i j (e:es) + yy_1 i j es
                                  | typeDev (fst e) == tCAP
                                    = cap i j (e:es) + yy_1 i j es 
                                  | typeDev (fst e) == tIND
                                    = ind i j (e:es) + yy_1 i j es 
                                  | otherwise = yy_1 i j es
                               
                  adm m n (e:es) | m == nodez 1 (fst e) && n == nodez 2 (fst e) 
                                      = -1 * (((snd e)!!0) :+ 0) 
                                 | n == nodez 1 (fst e) && m == nodez 2 (fst e) 
                                      = -1 * (((snd e)!!0) :+ 0)   
                                 | n == m && (n == nodez 1 (fst e) || 
                                   m == nodez 2 (fst e))   
                                      = (((snd e)!!0) :+ 0)
                                 | otherwise = 0
 
                  ind m n (e:es) | m == nodez 1 (fst e) && n == nodez 2 (fst e) 
                                      = -1.0 / (0 :+ (xL ((snd e)!!0))) 
                                 | n == nodez 1 (fst e) && m == nodez 2 (fst e) 
                                      = -1.0 / (0 :+ (xL ((snd e)!!0)))
                                 | n == m && (n == nodez 1 (fst e) || 
                                   m == nodez 2 (fst e))   
                                      = 1.0 / (0 :+ (xL ((snd e)!!0)))
                                 | otherwise = 0
 
                  cap m n (e:es) | m == nodez 1 (fst e) && n == nodez 2 (fst e) 
                                      = -1.0 / (0 :+ (xC ((snd e)!!0))) 
                                 | n == nodez 1 (fst e) && m == nodez 2 (fst e) 
                                      = -1.0 / (0 :+ (xC ((snd e)!!0))) 
                                 | n == m && (n == nodez 1 (fst e) || 
                                   m == nodez 2 (fst e))   
                                      = 1.0 / (0 :+ (xC ((snd e)!!0))) 
                                 | otherwise = 0
                               
                  vsrc i j (e:es) | length es == (u - i) && 
                                            j == nodez 1 (fst e) = 1
                                  | length es == (u - i) && 
                                            j == nodez 2 (fst e) = 1
                                  | length es == (u - j) && 
                                            i == nodez 1 (fst e) = 1
                                  | length es == (u - j) && 
                                            i == nodez 2 (fst e) = -1
                                  | otherwise = 0      
                  
                  m = ncol (e:es)
                  n = nrow (e:es) 
                  u = m + length (srcList (e:es))
                  freq = param!!pFREQ 
                  xC value = -1.0 / (value * 2 * pi  * freq) 
                  xL value = value * 2 * pi * freq
                 
--------------------------------------------------------------------------------
--- Build vector i of Yv = i
mk_i :: [([Int], [Double])] -> [Double] -> Array Int (Complex Double)    
mk_i (e:es) param = array (1, n) ([(j, z j) | j <- [1..n]])
          where z j | j <= m = ii_1 j (srcList (e:es)) 
                    | otherwise = ii_2 j (srcList (e:es))
                
                ii_1 _ [] = 0
                ii_1 j (e:es) | typeDev (fst e) == tIAC 
                                = isrc j (e:es) + ii_1 j es
                              | otherwise = ii_1 j es 
                
                ii_2 _ [] = 0
                ii_2 j (e:es) | typeDev (fst e) == tVAC 
                                = vsrc j (e:es) + ii_2 j es
                              | otherwise = ii_2 j es 
                
                vsrc j (v:vs) | length vs == (u - j) = ((snd v)!!0):+0      
                              | otherwise = 0 
                 
                isrc j (i:is) | j == nodez 1 (fst i) = (((snd i)!!0) :+ 0) 
                              | j == nodez 2 (fst i) = -1 * (((snd i)!!0) :+ 0)  
                              | otherwise = 0
                              
                m = ncol (e:es)
                n = nrow (e:es)
                u = m + length (srcList (e:es))
                          
--- Utilities
typeDev [] = 0
typeDev parm1 = parm1!!0

nodez _ [] = 0   
nodez n parm1 = parm1!!n
                
nodeList [] = []
nodeList (n:ns) = nub (union (tail $ fst n) (nodeList ns)) 

ncol netlist = length (nodeList netlist) - 1
nrow netlist = ncol netlist + length (srcList netlist) 
               - length (isrcList netlist)

srcList [] = []
srcList netlist = vsrcList netlist ++ isrcList netlist

admList [] = []
admList netlist = resList netlist ++ indList netlist ++ capList netlist 


resList :: [([Int], [Double])] -> [([Int], [Double])]
resList [] = []
resList (n:ns) | typeDev (fst n) == tRES = n : resList ns
               | otherwise = resList ns

capList :: [([Int], [Double])] -> [([Int], [Double])]
capList [] = []
capList (n:ns) | typeDev (fst n) == tCAP = n : capList ns
               | otherwise = capList ns

indList :: [([Int], [Double])] -> [([Int], [Double])]
indList [] = []
indList (n:ns) | typeDev (fst n) == tIND = n : indList ns
               | otherwise = indList ns


vsrcList :: [([Int], [Double])] -> [([Int], [Double])]
vsrcList [] = []
vsrcList (n:ns) | typeDev (fst n) == tVAC = n : vsrcList ns
                | otherwise = vsrcList ns

isrcList :: [([Int], [Double])] -> [([Int], [Double])]
isrcList [] = []
isrcList (n:ns) | typeDev (fst n) ==  tISRC = n : isrcList ns
                | otherwise = isrcList ns

-------------------------------------------------------------------------------
--- Solve v of Yv = i

pFREQ :: Int
pFREQ = 0

pPHASE :: Int
pPHASE = 1

acCalc :: [([Int], [Double])] -> [Double] -> Array Int (Complex Double)    
acCalc netlist param = psolve (mk_Y netlist param) (mk_i netlist param)