-------------------------------------------------------------------------------- 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