share/contrib/state/tree.mac

   1 /*
   2 This subroutine computes a proper tree
   3 Copyright (C) 1999  Dan Stanger
   4
   5 This library is free software; you can redistribute it and/or modify it
   6 under the terms of the GNU Library General Public License as published
   7 by the Free Software Foundation; either version 2 of the License, or (at
   8 your option) any later version.
   9
  10 This library is distributed in the hope that it will be useful, but
  11 WITHOUT ANY WARRANTY; without even the implied warranty of
  12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 Library General Public License for more details.
  14
  15 You should have received a copy of the GNU Library General Public
  16 License along with this library; if not, write to the Free Software
  17 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  18
  19 Dan Stanger dan.stanger@ieee.org
  20
  21 */
  22 /* this function will return the a matrix
  23    due to no local arrays, this function will prefix all arrays used
  24    with the prefix pt, to avoid name clashes.  it will delete them at the end
  25    the reference node will be 0.
  26 */
  27 nameNode(l):=first(l)$
  28 fromNode(l):=third(l)$
  29 toNode(l):=fourth(l)$
  30 typeNode(l):=second(l)$
  31 valueNode(l):=fifth(l)$
  32
  33 countNodes(cl):=block([n:0],
  34    map(lambda([x],n:max(n,fromNode(x),toNode(x))),cl),n)$
  35
  36 printarray2(a,na,n,b):=block([ipw:integer_print_width],integer_print_width:4,
  37    apply(print, makelist(na[j],j,1,b)),
  38    for i:1 thru n do apply(print,makelist(a[i,j],j,1,b)),
  39    integer_print_width:ipw)$
  40
  41 swap(a,i,j)::=buildq([a,i,j],block([temp:a[i]],a[i]:a[j],a[j]:temp))$
  42
  43 properTree(filename):=block([cl, MB, MN, retvals, dbg:true],
  44 /* MB number of branches */
  45 /* MN number of nodes */
  46 if dbg then print("entry"),
  47     cl:readfile(filename),
  48 if dbg then print("before sort", cl),
  49     cl:sort(cl,
  50         lambda([x,y],
  51             getElementIndex(typeNode(x))<getElementIndex(typeNode(y)))),
  52 if dbg then print("after sort", cl),
  53     MB:length(cl),
  54     MN:countNodes(cl),
  55 if dbg then print("branches", MB, " nodes", MN),
  56     array(pt_A,MN,MB), array(pt_D, MN, MB),
  57     array([pt_NF,pt_NT,pt_V,pt_JT, pt_MT, pt_name, pt_aJT], MB),
  58 /*      DIMENSION A(10,10),D(10,10),NF(10),NT(10),V(10),JT(10),MT(10) */
  59 /*      READ 99, (NF(I),NT(I),V(I),JT(I),MT(I) I=1,MB)
  60 C       I IS BRANCH NUMBER, NF IS FROM NODE, NT IS TO NODE
  61 C       V IS ELEMENT VALUE, JT IS ELEMENT TYPE
  62 C       MT TELLS TREE OR LINK
  63 C       WHEN MT IS 0 BRANCH IS TREE, OTHER NUMBER BRANCH IS
  64 C       LINK
  65 C       MT(I) HAS ALL ZERO ELEMENTS INITIALLY */
  66 /* pt_name is the symbolic name of the circuit element, and if a element
  67    is not given a value, then it's value is its name.  maybe there should
  68    be a flag, to assign symbolic names even if there are numeric values,
  69    to be more compatible with spice.
  70    pt_aJT is the element type in the order of the pt_A array */
  71     block([i:0],
  72         for l in cl do (i:i+1,
  73             pt_NF[i]:fromNode(l), pt_NT[i]:toNode(l), pt_V[i]:valueNode(l),
  74             pt_JT[i]:typeNode(l), pt_MT[i]:0, pt_name[i]:nameNode(l),
  75             if pt_V[i] = false then pt_V[i]:pt_name[i])),
  76 /*      CONSTRUCT INCIDENCE MATRIX D */
  77 if dbg then print("constructing incidence matrix d"),
  78     fillarray(pt_D,[0]),
  79     for i:1 thru MB do (
  80         if pt_NF[i] > 0 then pt_d[pt_NF[i],i]:1,
  81         if pt_NT[i] > 0 then pt_d[pt_NT[i],i]:-1),
  82 if dbg then print("number of nodes", MN, "number of branches", MB),
  83     /* print the d matrix for debugging */
  84 if dbg then printarray2(pt_d,pt_name,MN,MB),
  85 /*      SELECT TREE BRANCHES AND LINKS
  86 C       KV IS TREE VOLTAGE SOURCE NUMBER
  87 C       KC IS TREE CAPACITOR NUMBER
  88 C       KR IS TREE RESISTOR NUMBER
  89 C       IR IS LINK RESISTOR NUMBER
  90 C       IL IS LINK INDUCTOR NUMBER
  91 C       IC IS LINK CURRENT SOURCE NUMBER
  92 C       IK IS TREE BRANCH NUMBER
  93 C       K IS LINK BRANCH NUMBER */
  94     block([KV:0, KC:0, KR:0, IR:0, IL:0, IC:0, IK:0, K:0],
  95         for i:1 thru MB do (
  96             if pt_MT[i] = 0 then (
  97                 if i < MB then block([continu,treebranch:false],
  98                     for j:(i+1) thru MB while (treebranch = false) do (
  99                         continu:false,
 100 if dbg then print("52", j, pt_MT[j]),
 101                         if pt_MT[j] = 0 then (
 102 if dbg then print("pt_MT[j]=0,i,pt_NT[i],j,pt_NT[j]", i, pt_NT[i], j, pt_NT[j]),
 103                             if pt_NT[i] = pt_NT[j] then ( /* 61 */
 104                                 if pt_NF[i] # pt_NF[j] then (/* 5 */
 105                                     pt_NT[j]:pt_NF[i], /* 43 */
 106 if dbg then print("1 j-MB",j-MB,52,102,102),
 107                                     /* IF(J-MB)52,102,102 * 220 */
 108                                     continu:true
 109                                 ) else (if dbg then print("should do 6"))
 110                             )
 111                             else (
 112 if dbg then print("33 i,pt_NT[i],j,pt_NF[j]",i,pt_NT[i],j,pt_NF[j]),
 113                                 if pt_NT[i] # pt_NF[j] then ( /* 33 */
 114 if dbg then print("2 j-MB",j-MB,52,102,102),
 115                                     /* IF(J-MB)52,102,102 * 220 */
 116                                     continu:true
 117                                 )
 118                                 else (
 119                                     if pt_NF[i] # pt_NT[j] then ( /* 63 */
 120                                         pt_NF[j]:pt_NF[i], /* 83 */
 121                                         /* IF(J-MB)52,102,102 * 220 */
 122 if dbg then print("3 j-MB",j-MB,52,102,102),
 123                                         continu:true
 124                                     ) else (if dbg then print("should do 6"))
 125                                 )
 126                             ),
 127                             if continu = true then (
 128 if dbg then print("goto 102"),
 129                                 if j >= MB then treebranch:true /* goto 102 */
 130                             )
 131                             else ( /* 6 branch j is a link */
 132                                 K:K+1,
 133 if dbg then print("k", k, " adding a link 6"),
 134                                 pt_MT[j]:MN + K,
 135                                 pt_aJT[pt_MT[j]]:pt_JT[j],
 136                                 for L:1 thru MN do pt_A[L,pt_MT[j]]:pt_D[L,j],
 137                                 if getElementIndex(pt_JT[j]) < 4 then IR:IR+1
 138                                 else (
 139                                     if getElementIndex(pt_JT[j]) > 4 then IC:IC+1
 140                                     else IL:IL+1
 141                                 ),
 142                                 /* IF(J-MB)52,102,103 */
 143 if dbg then print("4 j-MB",j-MB,52,102,103),
 144                                 if j = MB then treebranch:true,
 145                                 if j >= MB then return
 146                             )
 147                         )
 148                     ),
 149 if dbg then print("102 treebranch is ", treebranch),
 150                     if treebranch = true then ( /* 102 branch i is a tree branch */
 151                         IK:IK+1,
 152 if dbg then print("ik", ik,"i",i, " in 102 row swap"),
 153                         for M:1 thru MN do pt_A[M,IK]:pt_D[M,i],
 154                         pt_aJT[IK]:pt_JT[i],
 155                         swap(pt_name,i,IK),
 156                         if getElementIndex(pt_JT[i]) < 2 then KV:KV+1
 157                         else (
 158                             if getElementIndex(pt_JT[i]) > 2 then KR:KR+1
 159                             else KC:KC+1
 160                         )
 161                     )
 162                 )
 163             )
 164         ),
 165         retvals:[KC, KR, KV, IL, IR, IC, MN, MB]
 166     ),
 167 if dbg then print("103? MB MN",MB,MN),
 168 if dbg then printarray2(pt_a, pt_name, MN, MB),
 169     retvals
 170 )
 171 $