Plan 9 from Bell Labs’s /usr/web/sources/contrib/fernan/escomma/DC.hs

Copyright © 2021 Plan 9 Foundation.
Distributed under the MIT License.
Download the Plan 9 distribution.


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









Bell Labs OSI certified Powered by Plan 9

(Return to Plan 9 Home Page)

Copyright © 2021 Plan 9 Foundation. All Rights Reserved.
Comments to [email protected].