octave_packages/ocs-0.1.3/asm/asm_initialize_system.m

   1 ## Copyright (C) 2006,2007,2008  Carlo de Falco, Culpo Massimiliano
   2 ##
   3 ## This file is part of:
   4 ## OCS - A Circuit Simulator for Octave
   5 ##
   6 ## OCS is free software; you can redistribute it and/or modify
   7 ## it under the terms of the GNU General Public License as published by
   8 ## the Free Software Foundation.
   9 ##
  10 ## This program is distributed in the hope that it will be useful,
  11 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 ## GNU General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with this program (see the file LICENSE); if not,
  17 ## see <http://www.gnu.org/licenses/>.
  18 ##
  19 ## author: Carlo de Falco <cdf _AT_ users.sourceforge.net>
  20 ## author: culpo@math.uni-wuppertal.de
  21
  22 ## -*- texinfo -*-
  23 ##
  24 ## @deftypefn{Function File} {[@var{A},@var{B},@var{C}] =}@
  25 ## asm_initialize_system(@var{instruct},@var{x})
  26 ##
  27 ## Cycle through the circuit description structure @var{instruct}
  28 ## to build the system matrices @var{A}, @var{B}, @var{C} for
  29 ## the linear and time-invariant part of the system.
  30 ##
  31 ## @var{x} is the current value of the state variables.
  32 ##
  33 ## See the @cite{IFF file format specifications} for details about
  34 ## the output matrices.
  35 ##
  36 ## @seealso{asm_build_system,prs_iff}
  37 ##
  38 ## @end deftypefn
  39
  40 function  [A,B,C] = asm_initialize_system(instruct,x);
  41
  42   ## Check input
  43   if nargin != 2
  44     error("asm_initialize_system: wrong number of input parameters.");
  45   elseif !(isstruct(instruct) && isfield(instruct,"LCR") &&
  46            isfield(instruct,"NLC"))
  47     error("asm_initialize_system: first input is not a valid structure.");
  48   elseif !isvector(x)
  49     error("asm_initialize_system: second input is not a valid vector.");
  50   endif
  51
  52   ## Build linear part of the system
  53   n  = instruct.totextvar + instruct.totintvar; # Number of variables
  54
  55   ## Initialize to zero any state variable that is not in the input argument
  56   lx = length(x);
  57   if lx < n
  58     x(lx+1:n) = 0;
  59   endif
  60   ## FIXME: should a warning be passed if lx != n ?
  61
  62   A = sparse(n,n);
  63
  64   ## LCR section
  65   B = sparse(n,n);
  66   C = sparse(n,1);
  67
  68   nblocks = length(instruct.LCR);
  69
  70   for ibl = 1:nblocks
  71     for iel = 1:instruct.LCR(ibl).nrows
  72
  73       ## Evaluate element
  74       if instruct.LCR(ibl).nintvar(iel)
  75         intvars = instruct.totextvar + instruct.LCR(ibl).osintvar(iel) + [1:instruct.LCR(ibl).nintvar(iel)]';
  76       else
  77         intvars=[];
  78       endif
  79
  80       il   = instruct.LCR(ibl).vnmatrix(iel,:)';
  81       nzil = find(il!=0);
  82
  83       y       = zeros(size(il));
  84       y(nzil) = x(il(nzil));
  85       z       = x(intvars);
  86
  87       [a,b,c] = feval(instruct.LCR(ibl).func,
  88                       instruct.LCR(ibl).section,
  89                       instruct.LCR(ibl).pvmatrix(iel,:),
  90                       instruct.LCR(ibl).parnames,
  91                       y,z,0);
  92
  93       ## Assemble matrices
  94
  95       ## Global indexing
  96       vars = [il(nzil);intvars];
  97
  98       ## Local indexing
  99       lclvars = [nzil; instruct.LCR(ibl).nextvar + (1:length(intvars))' ];
 100
 101       ## Reshaping sparse stamps
 102       a = a(lclvars,lclvars);
 103       b = b(lclvars,lclvars);
 104       c = reshape(c(lclvars),[],1);
 105
 106       [na,ma,va] = find(a);
 107       [nb,mb,vb] = find(b);
 108       [nc,mc,vc] = find(c);
 109
 110       ## Stamping
 111       A += sparse(vars(na),vars(ma),va,n,n);
 112       B += sparse(vars(nb),vars(mb),vb,n,n);
 113       C += sparse(vars(nc),1,vc,n,1);
 114
 115     endfor
 116   endfor
 117
 118 endfunction