share/contrib/state/state.mac

   1 /*
   2 This subroutine computes state variable equations
   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 /* load in readfile and tree */
  22 /* get the filename from the user */
  23 /* this code works around problems caused by:
  24    matrix.scalar does not simplify properly
  25    1 dimensional matrixes do not simplify properly
  26    empty matrixes do not remember how they were created
  27 */
  28 f(m,n,p,q,r,s, MN, MB):=block([l,a_T, a_L,b_T,b_L,c_T,c_L],
  29     /* reorder the array a_pt and turn it into the a_T and a_L matrixes  */
  30     l:sortem(MN,MB),a_T:first(l),a_L:second(l),
  31     b_T:-transpose(a_L).invert(transpose(a_T)),
  32     b_L:ident(MN),
  33     c_T:ident(MN),
  34     c_L:-transpose(b_T),
  35     block([
  36         em:matrix(),
  37         c_LLC, c_LRC, c_LGC, c_LLR, c_LRR, c_LGR, c_LLG, c_LRG, c_LGG,
  38         b_TCL, b_TRL, b_TGL, b_TCR, b_TRR, b_TGR, b_TCG, b_TRG, b_TGG,
  39         x:transpose(append( /* state variable vector */
  40             makelist(concat('v,pt_name[i]),i,1,m),
  41             makelist(concat('i,pt_name[i]),i,m+n+p+1,m+n+p+q))),
  42         u:transpose(append(makelist(
  43                 if pt_V[i] = false then pt_name[i] else pt_V[i],
  44                 i,m+n+1,m+n+p),
  45             makelist(
  46                 if pt_V[i] = false then pt_name[i] else pt_V[i],
  47                 i,m+n+p+q+r+1,m+n+p+q+r+s))),
  48         r_L,g_T,se,sx,su,e1,e2,e3,e4],
  49         if length(x) = 1 then x:x[1],
  50         if length(u) = 1 then u:u[1],
  51         submat_m(c_LLC, c_L, 1 .. m,            1 .. q),
  52         submat_m(c_LRC, c_L, 1 .. m,            (q+1) .. (q+r)),
  53         submat_m(c_LGC, c_L, 1 .. m,            (q+r+1) .. (q+r+s)),
  54         submat_m(c_LLR, c_L, (m+1) .. (m+n),    1 .. q),
  55         submat_m(c_LRR, c_L, (m+1) .. (m+n),    (q+1) .. (q+r)),
  56         submat_m(c_LGR, c_L, (m+1) .. (m+n),    (q+r+1) .. (q+r+s)),
  57         submat_m(c_LLG, c_L, (m+n+1) .. (m+n+p),        1 .. q),
  58         submat_m(c_LRG, c_L, (m+n+1) .. (m+n+p),        (q+1) .. (q+r)),
  59         submat_m(c_LGG, c_L, (m+n+1) .. (m+n+p),        (q+r+1) .. (q+r+s)),
  60
  61         submat_m(b_TCL, b_T, 1 .. q,            1 .. m),
  62         submat_m(b_TRL, b_T, 1 .. q,            (m+1) .. (m+n)),
  63         submat_m(b_TGL, b_T, 1 .. q,            (m+n+1) .. (m+n+p)),
  64         submat_m(b_TCR, b_T, (q+1) .. (q+r),    1 .. m),
  65         submat_m(b_TRR, b_T, (q+1) .. (q+r),    (m+1) .. (m+n)),
  66         submat_m(b_TGR, b_T, (q+1) .. (q+r),    (m+n+1) .. (m+n+p)),
  67         submat_m(b_TCG, b_T, (q+r+1) .. (q+r+s),        1 .. m),
  68         submat_m(b_TRG, b_T, (q+r+1) .. (q+r+s),        (m+1) .. (m+n)),
  69         submat_m(b_TGG, b_T, (q+r+1) .. (q+r+s),        (m+n+1) .. (m+n+p)),
  70
  71         r_L:makelist(pt_name[i],i,m+n+p+q+1,m+n+p+q+r),
  72         if length(r_L) > 0 then r_L:apply(diag_matrix,r_L) else r_L:em,
  73         g_T:makelist(1/pt_name[i],i,m+1,m+n),
  74         if length(g_T) > 0 then g_T:apply(diag_matrix,g_T) else g_T:em,
  75         /* compute free response */
  76         e1:mat_unblocker(
  77             matrix2([zeromatrix2(mat_nrows_m(c_LRC),mat_ncols_m(b_TRL)),c_LRC],
  78                 [b_TRL,zeromatrix2(mat_nrows_m(b_TRL),mat_ncols_m(c_LRC))])),
  79         e2:block([e:matrix2([b_TRR,r_L], [g_T, c_LRR])],
  80             if matrixp(e) then invert(e) else 1/e),
  81         e3:mat_unblocker(
  82             matrix2( [b_TCR, zeromatrix2(mat_nrows_m(b_TCR),mat_ncols_m(c_LLR))],
  83                 [zeromatrix2(mat_nrows_m(c_LLR),mat_ncols_m(b_TCR)),c_LLR])),
  84         e4:mat_unblocker(
  85         matrix2([zeromatrix2(mat_nrows_m(c_LLC),mat_ncols_m(b_TCL)),c_LLC],
  86                 [b_TCL,zeromatrix2(mat_nrows_m(b_TCL),mat_ncols_m(c_LLC))])),
  87         sx:(e1.e2.e3),
  88         if e4 # em then sx:sx-e4, /* check no link inductors */
  89         sx:if matrixp(x) then sx.x else x*sx, /* if x is scalar use * mult */
  90         /* compute source response */
  91         e3:mat_unblocker(
  92             matrix2( [b_TGR, zeromatrix2(mat_nrows_m(b_TGR),mat_ncols_m(c_LGR))],
  93                 [zeromatrix2(mat_nrows_m(c_LGR),mat_ncols_m(b_TGR)),c_LGR])),
  94         e4:mat_unblocker(
  95         matrix2([zeromatrix2(mat_nrows_m(c_LGC),mat_ncols_m(b_TGL)),c_LGC],
  96                 [b_TGL,zeromatrix2(mat_nrows_m(b_TGL),mat_ncols_m(c_LGC))])),
  97         su:(e1.e2.e3),
  98         if e4 # em then su:su-e4,
  99         su:if matrixp(u) then su.u else u*su,
 100         se:sx+su,
 101         /* try to remove . if its there */
 102         if not matrixp(se) then subst("*",".",sx+su) else se
 103     )
 104 )$
 105 state():=
 106 block([filename:read("enter the filename"), se],
 107     /* compute the proper tree and use apply to destructure it */
 108     se:apply(f,propertree(filename)))
 109 $