--------------------------------------------------------------------------------
module DC where
import Circuit
import Array
import List
import Matrix.LU
--import Matrix.Pivot
import Maybe
---- Y node admittance matrix
---- v vector of unknown voltages
---- i equivalent current source vector
--- Build matrix Y of Yv = i
mk_Y :: [Device] -> Array (Int, Int) Double
mk_Y netlist = yy netlist
where yy (e:es) = 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 (mYList (e:es))
| i > m || j > m
= yy_2 i j (srcList (e:es))
| otherwise = 0
yy_2 _ _ [] = 0
yy_2 i j (e:es) | devtype e == VSRC
= vsrc i j (e:es) + yy_2 i j es
| devtype e == VCVS
= vcvs i j (e:es) + yy_2 i j es
| devtype e == CCCS
= cccs i j (e:es) + yy_2 i j es
| devtype e == CCVS0
= ccvs0 i j (e:es) + yy_2 i j es
| devtype e == CCVS1
= ccvs1 i j (e:es) + yy_2 i j es
| devtype e == VCCS
= vccs i j (e:es) + yy_2 i j es
| otherwise = yy_2 i j es
yy_1 _ _ [] = 0
yy_1 i j (e:es) | devtype e == ADM
= adm i j (e:es) + yy_1 i j es
| otherwise = yy_1 i j es
adm m n (e:es) | m == nodez 1 e &&
n == nodez 2 e = -1 * devvalue e
| n == nodez 1 e &&
m == nodez 2 e = -1 * devvalue e
| n == m && (n == nodez 1 e ||
m == nodez 2 e) = devvalue e
| otherwise = 0
vsrc i j (e:es) | i == m + (vseqz e netlist) &&
j == nodez 1 e = 1
| i == m + (vseqz e netlist) &&
j == (nodez 2 e) = -1
| j == m + (vseqz e netlist) &&
i == nodez 1 e = 1
| j == m + (vseqz e netlist) &&
i == nodez 2 e = -1
| otherwise = 0
vccs i j (e:es) | i == m + (vseqz e netlist) &&
j == nodez 3 e = 1
| i == m + (vseqz e netlist) &&
j == nodez 4 e = -1
| j == m + (vseqz e netlist) &&
i == nodez 1 e = 1
| j == m + (vseqz e netlist) &&
i == nodez 2 e = -1
| i == m + (vseqz e netlist) &&
j == m + (vseqz e netlist)
= -1 * devvalue e
| otherwise = 0
vcvs i j (e:es) | i == m + (vseqz e netlist) &&
j == nodez 1 e = -1
| i == m + (vseqz e netlist) &&
j == nodez 2 e = 1
| i == m + (vseqz e netlist) &&
j == nodez 3 e = 1 * devvalue e
| i == m + (vseqz e netlist) &&
j == nodez 4 e = -1 * devvalue e
| j == m + (vseqz e netlist) &&
i == nodez 1 e = -1
| j == m + (vseqz e netlist) &&
i == nodez 2 e = 1
| otherwise = 0
cccs i j (e:es) | i == m + (vseqz e netlist) &&
j == nodez 3 e = 1
| i == m + (vseqz e netlist) &&
j == nodez 4 e = -1
| j == m + (vseqz e netlist) &&
i == nodez 1 e = 1
| j == m + (vseqz e netlist) &&
i == nodez 2 e = -1
| j == m + (vseqz e netlist) &&
i == nodez 3 e = devvalue e
| j == m + (vseqz e netlist) &&
i == nodez 4 e = -1 * devvalue e
| otherwise = 0
ccvs0 i j (e:es) | j == m + (vseqz e netlist) &&
i == nodez 3 e = 1
| j == m + (vseqz e netlist) &&
i == nodez 4 e = -1
| i == m + (vseqz e netlist) &&
j == nodez 1 e = 1
| i == m + (vseqz e netlist) &&
j == nodez 2 e = -1
| i == m + (vseqz e netlist) &&
j == m + (vseqz e netlist)
= -1 * devvalue e
| otherwise = 0
ccvs1 i j (e:es) | j == m + (vseqz e netlist) &&
i == nodez 1 e = -1
| j == m + (vseqz e netlist) &&
i == nodez 2 e = 1
| i == m + (vseqz e netlist) &&
j == nodez 3 e = 1
| i == m + (vseqz e netlist) &&
j == nodez 4 e = -1
| otherwise = 0
m = ncol netlist
n = nrow netlist
--------------------------------------------------------------------------------
--- Build vector i of Yv = i
mk_i :: [Device] -> Array Int Double
mk_i netlist = ii netlist
where ii (e:es) = 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) | devtype e == ISRC
= isrc j (e:es) + ii_1 j es
| otherwise = ii_1 j es
ii_2 _ [] = 0
ii_2 j (e:es) | devtype e == VSRC
= vsrc j (e:es) + ii_2 j es
| otherwise = ii_2 j es
vsrc j (v:vs) | j == m + (vseqz v netlist) = devvalue v
| otherwise = 0
isrc j (i:is) | j == nodez 1 i = devvalue i
| j == nodez 2 i = -1 * devvalue i
| otherwise = 0
m = ncol netlist
n = nrow netlist
--- Utilities
vseqz n netlist = (fromJust (elemIndex n (srcList (netlist))) + 1)
nodez n e | length (devnode e) == 0 = 0
| otherwise = (devnode e)!!(n - 1)
nodeList [] = []
nodeList netlist = sort (nub (foldr1 (++) (map (node) netlist)))
where node dev = (devnode dev)
ncol netlist = length (nodeList netlist) - 1
nrow netlist = ncol netlist + length (srcList netlist)
- length (isrcList netlist)
srcList [] = []
srcList netlist = vsrcList netlist ++ vccsList netlist ++vcvsList netlist
++ cccsList netlist ++ ccvsList netlist ++ isrcList netlist
mYList [] = []
mYList netlist = admList netlist
admList :: [Device] -> [Device]
admList [] = []
admList (n:ns) | devtype n == ADM = n : admList ns
| otherwise = admList ns
vsrcList :: [Device] -> [Device]
vsrcList [] = []
vsrcList (n:ns) | devtype n == VSRC = n : vsrcList ns
| otherwise = vsrcList ns
vccsList :: [Device] -> [Device]
vccsList [] = []
vccsList (n:ns) | devtype n == VCCS = n : vccsList ns
| otherwise = vccsList ns
vcvsList :: [Device] -> [Device]
vcvsList [] = []
vcvsList (n:ns) | devtype n == VCVS = n : vcvsList ns
| otherwise = vcvsList ns
cccsList :: [Device] -> [Device]
cccsList [] = []
cccsList (n:ns) | devtype n == CCCS = n : cccsList ns
| otherwise = cccsList ns
ccvsList :: [Device] -> [Device]
ccvsList [] = []
ccvsList (n:ns) | devtype n == CCVS0 || devtype n == CCVS1
= n : ccvsList ns
| otherwise = ccvsList ns
isrcList :: [Device] -> [Device]
isrcList [] = []
isrcList (n:ns) | devtype n == ISRC = n : isrcList ns
| otherwise = isrcList ns
nodenames netlist = nodeList (netlist defaultESim)
srcnames netlist = map pp (srcList (netlist defaultESim))
where pp device = (devname device)
--- Solve v of Yv = i ---------------------------------------------------------
dcCalc :: [Device] -> Array Int Double
dcCalc netlist = psolve (mk_Y netlist) (mk_i netlist)
-------------------------------------------------------------------------------
dcInit netlist = map cc netlist
where cc dev = (Device (devtype dev)
(devname dev)
(node dev)
(devvalue dev))
node dev = map dd (devnode dev)
dd n = fromJust (elemIndex n (nodeList netlist))
dcPoint1 z0 netlist = (ESim p0 nrr trr zz t0)
where z' n e nn | n == 10000 = z
| e < 0 = z
| otherwise = z' (n + 1) zerr z
where z = elems (dcCalc (dcInit (netlist sim0)))
zerr = maximum (err nn z)
sim0 = ESim p0 (n + 1) trr nn t0
e = 1.0
zz = z' 0 e (nrData z0)
p0 = simInfo z0
t0 = trData z0
trr = trSteps z0
nrr = nrSteps z0
err (a:as) (b:bs) | as == [] && bs == [] = df a b : []
| otherwise = df a b : err as bs
df l m = (abs (l - m)) - ((abs (0.001 * l)) + 1e-10)
--- DC Operating Point
data OPData = OPData
{ opInfo :: SimInfo
, opvalue :: [Double]
} deriving (Show)
-- n0 = Initial node voltage guesses
-- t0 = Initial node time dependent node voltage
dcOP1 n0 t0 netlist = (OPData p0 dd)
where dd = nrData (dcPoint1 parm0 netlist)
p0 = (SimInfo OP n s [] 0)
parm0 = (ESim p0 0 0 nn0 tt0)
n = nodenames netlist
s = srcnames netlist
tt0 | length t0 == 0
= map (0*) [1..fromIntegral(nrow
(netlist defaultESim))]
| otherwise = t0
nn0 | length n0 == 0 = tt0
| otherwise = n0
--- transient analysis
data TrData = TrData
{ trInfo :: SimInfo
, trvalue :: [[Double]]
} deriving (Show)
dcTran1 tstep tend n0 t0 netlist = (TrData p0 (z' 0 zz))
where z' n j | n > steps = []
| otherwise = z : z' (n + 1) sim0
where z = nrData (dcPoint1 j netlist)
sim0 = (ESim p0 0 (n + 1) nz0 tz0)
nz0 = z --trData j
tz0 = z --trData j
steps = tend / tstep
zz = (ESim p0 0 0 nn0 tt0)
p0 = (SimInfo TRANS m s [tstep] 0)
m = nodenames netlist
s = srcnames netlist
tt0 | length t0 == 0
= opvalue $ dcOP1 [] [] netlist
| otherwise = t0
nn0 | length n0 == 0 = tt0
| otherwise = n0
|
