Branch data Line data Source code
1 : : /*
2 : : * hbsolver.cpp - harmonic balance solver class implementation
3 : : *
4 : : * Copyright (C) 2005, 2006, 2007, 2008 Stefan Jahn <stefan@lkcc.org>
5 : : *
6 : : * This 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; either version 2, or (at your option)
9 : : * any later version.
10 : : *
11 : : * This software is distributed in the hope that it will be useful,
12 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : : * GNU General Public License for more details.
15 : : *
16 : : * You should have received a copy of the GNU General Public License
17 : : * along with this package; see the file COPYING. If not, write to
18 : : * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
19 : : * Boston, MA 02110-1301, USA.
20 : : *
21 : : * $Id$
22 : : *
23 : : */
24 : :
25 : : #if HAVE_CONFIG_H
26 : : # include <config.h>
27 : : #endif
28 : :
29 : : #include <stdio.h>
30 : :
31 : : #include "object.h"
32 : : #include "logging.h"
33 : : #include "complex.h"
34 : : #include "vector.h"
35 : : #include "node.h"
36 : : #include "circuit.h"
37 : : #include "component_id.h"
38 : : #include "net.h"
39 : : #include "netdefs.h"
40 : : #include "strlist.h"
41 : : #include "ptrlist.h"
42 : : #include "tvector.h"
43 : : #include "tmatrix.h"
44 : : #include "eqnsys.h"
45 : : #include "analysis.h"
46 : : #include "dataset.h"
47 : : #include "fourier.h"
48 : : #include "hbsolver.h"
49 : :
50 : : #define HB_DEBUG 0
51 : :
52 : : namespace qucs {
53 : :
54 : : using namespace fourier;
55 : :
56 : : // Constructor creates an unnamed instance of the hbsolver class.
57 [ + - ][ + - ]: 1 : hbsolver::hbsolver () : analysis () {
[ + - ][ # # ]
[ # # ][ # # ]
58 : 1 : type = ANALYSIS_HBALANCE;
59 : 1 : frequency = 0;
60 : 1 : nlnodes = lnnodes = banodes = nanodes = NULL;
61 : 1 : Y = Z = A = NULL;
62 : 1 : NA = YV = JQ = JG = JF = NULL;
63 : 1 : OM = IR = QR = RH = IG = FQ = VS = VP = FV = IL = IN = IC = IS = NULL;
64 : 1 : vs = x = NULL;
65 : 1 : runs = 0;
66 : 1 : ndfreqs = NULL;
67 : 1 : }
68 : :
69 : : // Constructor creates a named instance of the hbsolver class.
70 [ # # ][ # # ]: 0 : hbsolver::hbsolver (char * n) : analysis (n) {
[ # # ][ # # ]
[ # # ][ # # ]
71 : 0 : type = ANALYSIS_HBALANCE;
72 : 0 : frequency = 0;
73 : 0 : nlnodes = lnnodes = banodes = nanodes = NULL;
74 : 0 : Y = Z = A = NULL;
75 : 0 : NA = YV = JQ = JG = JF = NULL;
76 : 0 : OM = IR = QR = RH = IG = FQ = VS = VP = FV = IL = IN = IC = IS = NULL;
77 : 0 : vs = x = NULL;
78 : 0 : runs = 0;
79 : 0 : ndfreqs = NULL;
80 : 0 : }
81 : :
82 : : // Destructor deletes the hbsolver class object.
83 : 1 : hbsolver::~hbsolver () {
84 : : // delete nodelists
85 [ + - ][ + - ]: 1 : if (nlnodes) delete nlnodes;
[ + - ][ # # ]
[ # # ][ # # ]
86 [ + - ][ + - ]: 1 : if (lnnodes) delete lnnodes;
[ + - ][ # # ]
[ # # ][ # # ]
87 [ + - ][ + - ]: 1 : if (banodes) delete banodes;
[ + - ][ # # ]
[ # # ][ # # ]
88 [ + - ][ + - ]: 1 : if (nanodes) delete nanodes;
[ + - ][ # # ]
[ # # ][ # # ]
89 : :
90 : : // delete temporary matrices
91 [ - + ][ # # ]: 1 : if (A) delete A;
[ # # ][ # # ]
92 [ - + ][ # # ]: 1 : if (Z) delete Z;
[ # # ][ # # ]
93 [ - + ][ # # ]: 1 : if (Y) delete Y;
[ # # ][ # # ]
94 : :
95 : : // delete matrices
96 [ + - ][ + - ]: 1 : if (NA) delete NA;
[ # # ][ # # ]
97 [ + - ][ + - ]: 1 : if (YV) delete YV;
[ # # ][ # # ]
98 [ + - ][ + - ]: 1 : if (JQ) delete JQ;
[ # # ][ # # ]
99 [ + - ][ + - ]: 1 : if (JG) delete JG;
[ # # ][ # # ]
100 [ + - ][ + - ]: 1 : if (JF) delete JF;
[ # # ][ # # ]
101 : :
102 : : // delete vectors
103 [ + - ][ + - ]: 1 : if (IC) delete IC;
[ # # ][ # # ]
104 [ + - ][ + - ]: 1 : if (IS) delete IS;
[ # # ][ # # ]
105 [ + - ][ + - ]: 1 : if (FV) delete FV;
[ # # ][ # # ]
106 [ + - ][ + - ]: 1 : if (IL) delete IL;
[ # # ][ # # ]
107 [ + - ][ + - ]: 1 : if (IN) delete IN;
[ # # ][ # # ]
108 [ + - ][ + - ]: 1 : if (IG) delete IG;
[ # # ][ # # ]
109 [ + - ][ + - ]: 1 : if (FQ) delete FQ;
[ # # ][ # # ]
110 [ + - ][ + - ]: 1 : if (VS) delete VS;
[ # # ][ # # ]
111 [ + - ][ + - ]: 1 : if (VP) delete VP;
[ # # ][ # # ]
112 [ + - ][ + - ]: 1 : if (vs) delete vs;
[ # # ][ # # ]
113 [ + - ][ + - ]: 1 : if (OM) delete OM;
[ # # ][ # # ]
114 [ + - ][ + - ]: 1 : if (IR) delete IR;
[ # # ][ # # ]
115 [ + - ][ + - ]: 1 : if (QR) delete QR;
[ # # ][ # # ]
116 [ + - ][ + - ]: 1 : if (RH) delete RH;
[ # # ][ # # ]
117 : :
118 [ + - ][ + - ]: 1 : if (x) delete x;
[ # # ][ # # ]
119 [ + - ][ + - ]: 1 : if (ndfreqs) delete[] ndfreqs;
[ # # ][ # # ]
120 [ - + ][ # # ]: 2 : }
121 : :
122 : : /* The copy constructor creates a new instance of the hbsolver class
123 : : based on the given hbsolver object. */
124 [ # # ][ # # ]: 0 : hbsolver::hbsolver (hbsolver & o) : analysis (o) {
[ # # ][ # # ]
[ # # ][ # # ]
125 : 0 : frequency = o.frequency;
126 [ # # ][ # # ]: 0 : negfreqs = o.negfreqs;
127 [ # # ][ # # ]: 0 : posfreqs = o.posfreqs;
128 : 0 : nlnodes = o.nlnodes;
129 : 0 : lnnodes = o.lnnodes;
130 : 0 : banodes = o.banodes;
131 : 0 : nanodes = o.nanodes;
132 : 0 : Y = Z = A = NULL;
133 : 0 : NA = YV = JQ = JG = JF = NULL;
134 : 0 : OM = IR = QR = RH = IG = FQ = VS = VP = FV = IL = IN = IC = IS = NULL;
135 : 0 : vs = x = NULL;
136 : 0 : runs = o.runs;
137 : 0 : ndfreqs = NULL;
138 : 0 : }
139 : :
140 : : #define VS_(r) (*VS) (r)
141 : : #define OM_(r) (*OM) (r)
142 : :
143 : : /* This is the HB netlist solver. It prepares the circuit list and
144 : : solves it then. */
145 : 1 : int hbsolver::solve (void) {
146 : :
147 : 1 : int iterations = 0, done = 0;
148 : 1 : int MaxIterations = getPropertyInteger ("MaxIter");
149 : :
150 : : // collect different parts of the circuit
151 : 1 : splitCircuits ();
152 : :
153 : : // create frequency array
154 : 1 : collectFrequencies ();
155 : :
156 : : // find interconnects between the linear and non-linear subcircuit
157 : 1 : getNodeLists ();
158 : :
159 : : // prepares the linear part --> 0 = IC + [YV] * VS
160 : 1 : prepareLinear ();
161 : :
162 : 1 : runs++;
163 : : logprint (LOG_STATUS, "NOTIFY: %s: solving for %d frequencies\n",
164 : 1 : getName (), lnfreqs);
165 : :
166 [ + - ]: 1 : if (nbanodes > 0) {
167 : :
168 : : // start balancing
169 : : logprint (LOG_STATUS, "NOTIFY: %s: balancing at %d nodes\n", getName (),
170 : 1 : nbanodes);
171 : :
172 : : // prepares the non-linear part
173 : 1 : prepareNonLinear ();
174 : :
175 : : #if HB_DEBUG
176 : : fprintf (stderr, "YV -- transY in f:\n"); YV->print ();
177 : : fprintf (stderr, "IC -- constant current in f:\n"); IC->print ();
178 : : #endif
179 : :
180 : : // start iteration
181 [ + - ][ + - ]: 10 : do {
[ + - ]
182 : 11 : iterations++;
183 : :
184 : : #if HB_DEBUG
185 : : fprintf (stderr, "\n -- iteration step: %d\n", iterations);
186 : : fprintf (stderr, "vs -- voltage in t:\n"); vs->print ();
187 : : #endif
188 : :
189 : : // evaluate component functionality and fill matrices and vectors
190 : 11 : loadMatrices ();
191 : :
192 : : #if HB_DEBUG
193 : : fprintf (stderr, "FQ -- charge in t:\n"); FQ->print ();
194 : : fprintf (stderr, "IG -- current in t:\n"); IG->print ();
195 : : #endif
196 : :
197 : : // currents into frequency domain
198 : 11 : VectorFFT (IG);
199 : :
200 : : // charges into frequency domain
201 : 11 : VectorFFT (FQ);
202 : :
203 : : // right hand side currents and charges into the frequency domain
204 : 11 : VectorFFT (IR);
205 : 11 : VectorFFT (QR);
206 : :
207 : : #if HB_DEBUG
208 : : fprintf (stderr, "VS -- voltage in f:\n"); VS->print ();
209 : : fprintf (stderr, "FQ -- charge in f:\n"); FQ->print ();
210 : : fprintf (stderr, "IG -- current in f:\n"); IG->print ();
211 : : fprintf (stderr, "IR -- corrected Jacobi current in f:\n"); IR->print ();
212 : : #endif
213 : :
214 : : // solve HB equation --> FV = IC + [YV] * VS + j[O] * FQ + IG
215 : 11 : solveHB ();
216 : :
217 : : #if HB_DEBUG
218 : : fprintf (stderr, "FV -- error vector F(V) in f:\n"); FV->print ();
219 : : fprintf (stderr, "IL -- linear currents in f:\n"); IL->print ();
220 : : fprintf (stderr, "IN -- non-linear currents in f:\n"); IN->print ();
221 : : fprintf (stderr, "RH -- right-hand side currents in f:\n"); RH->print ();
222 : : #endif
223 : :
224 : : // termination criteria met
225 [ + + ][ + + ]: 11 : if (iterations > 1 && checkBalance ()) {
[ + + ]
226 : 1 : done = 1;
227 : 1 : break;
228 : : }
229 : :
230 : : #if HB_DEBUG
231 : : fprintf (stderr, "JG -- G-Jacobian in t:\n"); JG->print ();
232 : : fprintf (stderr, "JQ -- C-Jacobian in t:\n"); JQ->print ();
233 : : #endif
234 : :
235 : : // G-Jacobian into frequency domain
236 : 10 : MatrixFFT (JG);
237 : :
238 : : // C-Jacobian into frequency domain
239 : 10 : MatrixFFT (JQ);
240 : :
241 : : #if HB_DEBUG
242 : : fprintf (stderr, "JQ -- dQ/dV C-Jacobian in f:\n"); JQ->print ();
243 : : fprintf (stderr, "JG -- dI/dV G-Jacobian in f:\n"); JG->print ();
244 : : #endif
245 : :
246 : : // calculate Jacobian --> JF = [YV] + j[O] * JQ + JG
247 : 10 : calcJacobian ();
248 : :
249 : : #if HB_DEBUG
250 : : fprintf (stderr, "JF -- full Jacobian in f:\n"); JF->print ();
251 : : #endif
252 : :
253 : : // solve equation system --> JF * VS(n+1) = JF * VS(n) - FV
254 : 10 : solveVoltages ();
255 : :
256 : : #if HB_DEBUG
257 : : fprintf (stderr, "VS -- next voltage in f:\n"); VS->print ();
258 : : #endif
259 : :
260 : : // inverse FFT of frequency domain voltage vector VS(n+1)
261 : 10 : VectorIFFT (vs);
262 : : }
263 : : // check termination criteria (balanced frequency domain currents)
264 : : while (!done && iterations < MaxIterations);
265 : :
266 [ - + ]: 1 : if (iterations >= MaxIterations) {
267 [ # # ]: 0 : qucs::exception * e = new qucs::exception (EXCEPTION_NO_CONVERGENCE);
268 : : e->setText ("no convergence in %s analysis after %d iterations",
269 : 0 : getName (), iterations);
270 : 0 : throw_exception (e);
271 : : logprint (LOG_ERROR, "%s: no convergence after %d iterations\n",
272 : 0 : getName (), iterations);
273 : : }
274 : : else {
275 : : logprint (LOG_STATUS, "%s: convergence reached after %d iterations\n",
276 : 1 : getName (), iterations);
277 : : }
278 : : }
279 : : else {
280 : : // no balancing necessary
281 : 0 : logprint (LOG_STATUS, "NOTIFY: %s: no balancing necessary\n", getName ());
282 : : }
283 : :
284 : : // print exception stack
285 : 1 : estack.print ();
286 : :
287 : : // apply AC analysis to the complete network in order to obtain the
288 : : // final results
289 : 1 : finalSolution ();
290 : :
291 : : // save results into dataset
292 : 1 : saveResults ();
293 : 1 : return 0;
294 : : }
295 : :
296 : : /* Goes through the list of circuit objects and runs its calcHB()
297 : : function. */
298 : 0 : void hbsolver::calc (hbsolver * self) {
299 : 0 : circuit * root = self->getNet()->getRoot ();
300 [ # # ]: 0 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) {
301 : 0 : c->calcHB (self->frequency);
302 : : }
303 : 0 : }
304 : :
305 : : /* Goes through the list of circuit objects and runs its initHB()
306 : : function. */
307 : 0 : void hbsolver::initHB (void) {
308 : 0 : circuit * root = subnet->getRoot ();
309 [ # # ]: 0 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) {
310 : 0 : c->initHB ();
311 : : }
312 : 0 : }
313 : :
314 : : /* Goes through the list of circuit objects and runs its initDC()
315 : : function. */
316 : 0 : void hbsolver::initDC (void) {
317 : 0 : circuit * root = subnet->getRoot ();
318 [ # # ]: 0 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) {
319 : 0 : c->initDC ();
320 : : }
321 : 0 : }
322 : :
323 : : // Returns true if circuit is a HB source.
324 : 4 : bool hbsolver::isExcitation (circuit * c) {
325 : 4 : int type = c->getType ();
326 [ + - ][ + + ]: 4 : if (type == CIR_PAC || type == CIR_VAC || type == CIR_VDC)
[ - + ]
327 : 1 : return true;
328 : 4 : return false;
329 : : }
330 : :
331 : : // Expands the frequency array using the given frequency and the order.
332 : 1 : void hbsolver::expandFrequencies (nr_double_t f, int n) {
333 [ + - ]: 1 : tvector<nr_double_t> nfreqs = negfreqs;
334 [ + - ]: 1 : tvector<nr_double_t> pfreqs = posfreqs;
335 : 1 : int i, k, len = nfreqs.getSize ();
336 : 1 : negfreqs.clear ();
337 : 1 : posfreqs.clear ();
338 [ - + ]: 1 : if (len > 0) {
339 : : // frequency expansion for full frequency sets
340 [ # # ]: 0 : for (i = 0; i <= n + 1; i++) {
341 [ # # ]: 0 : for (k = 0; k < len; k++) {
342 [ # # ]: 0 : negfreqs.add (i * f + nfreqs.get (k));
343 : : }
344 : : }
345 [ # # ]: 0 : for (i = -n; i < 0; i++) {
346 [ # # ]: 0 : for (k = 0; k < len; k++) {
347 [ # # ]: 0 : negfreqs.add (i * f + nfreqs.get (k));
348 : : }
349 : : }
350 [ # # ]: 0 : for (i = 0; i <= 2 * n + 1; i++) {
351 [ # # ]: 0 : for (k = 0; k < len; k++) {
352 [ # # ]: 0 : posfreqs.add (i * f + pfreqs.get (k));
353 : : }
354 : : }
355 : : }
356 : : else {
357 : : // first frequency
358 [ + + ]: 10 : for (i = 0; i <= n + 1; i++) negfreqs.add (i * f);
359 [ + + ]: 8 : for (i = -n; i < 0; i++) negfreqs.add (i * f);
360 [ + + ]: 17 : for (i = 0; i <= 2 * n + 1; i++) posfreqs.add (i * f);
361 : 1 : }
362 : 1 : }
363 : :
364 : : // Calculates an order fulfilling the "power of two" requirement.
365 : 1 : int hbsolver::calcOrder (int n) {
366 : 1 : int o, order = n * 2; // current order + DC + negative freqencies
367 [ + + ]: 5 : for (o = 1; o < order; o <<= 1) ; // a power of 2
368 : 1 : return o / 2 - 1;
369 : : }
370 : :
371 : : /* The function computes the harmonic frequencies excited in the
372 : : circuit list depending on the maximum number of harmonics per
373 : : exitation and saves its results into the 'negfreqs' vector. */
374 : 1 : void hbsolver::collectFrequencies (void) {
375 : :
376 : : // initialization
377 : 1 : negfreqs.clear ();
378 : 1 : posfreqs.clear ();
379 : 1 : rfreqs.clear ();
380 : 1 : dfreqs.clear ();
381 [ - + ][ # # ]: 1 : if (ndfreqs) delete[] ndfreqs;
382 : :
383 : : // obtain order
384 : 1 : int i, n = calcOrder (getPropertyInteger ("n"));
385 : :
386 : : // expand frequencies for each exitation
387 : : nr_double_t f;
388 [ + + ]: 2 : for (auto * c : excitations) {
389 [ + - ]: 1 : if (c->getType () != CIR_VDC) { // no extra DC sources
390 [ + - ][ + - ]: 1 : if ((f = c->getPropertyDouble ("f")) != 0.0) {
391 [ + - ]: 1 : if (!dfreqs.contains (f)) { // no double frequencies
392 : 1 : dfreqs.add (f);
393 [ + - ]: 1 : expandFrequencies (f, n);
394 : : }
395 : : }
396 : : }
397 : : }
398 : :
399 : : // no excitations
400 [ - + ]: 1 : if (negfreqs.getSize () == 0) {
401 : : // use specified frequency
402 : 0 : f = getPropertyDouble ("f");
403 : 0 : dfreqs.add (f);
404 : 0 : expandFrequencies (f, n);
405 : : }
406 : :
407 : : // build frequency dimension lengths
408 : 1 : ndfreqs = new int[dfreqs.getSize ()];
409 [ + + ]: 2 : for (i = 0; i < dfreqs.getSize (); i++) {
410 : 1 : ndfreqs[i] = (n + 1) * 2;
411 : : }
412 : :
413 : : #if HB_DEBUG
414 : : fprintf (stderr, "%d frequencies: [ ", negfreqs.getSize ());
415 : : for (i = 0; i < negfreqs.getSize (); i++) {
416 : : fprintf (stderr, "%g ", (double) real (negfreqs.get (i)));
417 : : }
418 : : fprintf (stderr, "]\n");
419 : : #endif /* HB_DEBUG */
420 : :
421 : : // build list of positive frequencies including DC
422 [ + + ]: 17 : for (n = 0; n < negfreqs.getSize (); n++) {
423 [ + + ]: 16 : if ((f = negfreqs (n)) < 0.0) continue;
424 : 9 : rfreqs.add (f);
425 : : }
426 : 1 : lnfreqs = rfreqs.getSize ();
427 : 1 : nlfreqs = negfreqs.getSize ();
428 : :
429 : : // pre-calculate the j[O] vector
430 [ + - ]: 1 : OM = new tvector<nr_complex_t> (nlfreqs);
431 [ + + ]: 17 : for (n = i = 0; n < nlfreqs; n++, i++)
432 : 16 : OM_(n) = nr_complex_t (0, 2 * M_PI * negfreqs (i));
433 : 1 : }
434 : :
435 : : // Split netlist into excitation, linear and non-linear part.
436 : 1 : void hbsolver::splitCircuits (void) {
437 : 1 : circuit * root = subnet->getRoot ();
438 [ + + ]: 6 : for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) {
439 [ + + ]: 5 : if (c->isNonLinear ()) {
440 : : // non-linear part
441 [ + - ]: 1 : nolcircuits.push_front(c);
442 : : }
443 [ + + ]: 4 : else if (isExcitation (c)) {
444 : : // get sinusoidal sources
445 [ + - ]: 1 : excitations.push_front(c);
446 : : }
447 [ + + ]: 3 : else if (c->getType () != CIR_GROUND) {
448 : : // linear part
449 [ + - ]: 2 : lincircuits.push_front (c);
450 : : }
451 : : }
452 : 1 : }
453 : :
454 : : // Creates a nodelist for the given circuit list.
455 : 3 : strlist * hbsolver::circuitNodes (ptrlist<circuit> circuits) {
456 [ + - ]: 3 : strlist * nodes = new strlist ();
457 [ + + ]: 7 : for (auto * c : circuits) {
458 [ + + ]: 12 : for (int i = 0; i < c->getSize (); i++) {
459 [ + - ]: 8 : char * n = c->getNode(i)->getName ();
460 [ + + ]: 8 : if (strcmp (n, "gnd")) {
461 [ + - ][ + + ]: 6 : if (!nodes->contains (n)) nodes->add (n);
[ + - ]
462 : : }
463 : : }
464 : : }
465 : 3 : return nodes;
466 : : }
467 : :
468 : : // Obtain node lists for linear and non-linear part.
469 : 1 : void hbsolver::getNodeLists (void) {
470 : : // non-linear nodes
471 [ + - ][ + - ]: 1 : nlnodes = circuitNodes (nolcircuits);
472 : : // linear nodes
473 [ + - ][ + - ]: 1 : lnnodes = circuitNodes (lincircuits);
474 : : // excitation nodes
475 [ + - ][ + - ]: 1 : exnodes = circuitNodes (excitations);
476 : :
477 : : #if DEBUG && 0
478 : : fprintf (stderr, "nonlinear nodes: [ %s ]\n", nlnodes->toString ());
479 : : fprintf (stderr, " linear nodes: [ %s ]\n", lnnodes->toString ());
480 : : #endif
481 : :
482 : : // organization of the nodes for the MNA:
483 : : // --------------------------------------
484 : : // 1.) balanced nodes: all connected to at least one non-linear device
485 : : // 2.a) the excitation nodes
486 : : // 2.b) all linear nodes not contained in 1. and 2.a
487 : : // 2.c) additional gyrators of linear nodes (built-in voltage sources)
488 : : // please note: excitation nodes also in 2.b; 1. and 2.a are 'ports'
489 : :
490 [ + - ][ + - ]: 1 : nanodes = new strlist (*nlnodes); // list 1.
491 [ + - ]: 1 : strlistiterator it;
492 : : // add excitation nodes; list 2.a
493 [ + - ][ + - ]: 2 : for (it = strlistiterator (exnodes); *it; ++it)
[ + - ][ + - ]
[ + + ]
494 [ + - ][ + - ]: 1 : nanodes->append (*it);
495 : : // add linear nodes; list 2.b
496 [ + - ][ + - ]: 4 : for (it = strlistiterator (lnnodes); *it; ++it) {
[ + - ][ + - ]
[ + + ]
497 [ + - ][ + - ]: 3 : if (!nanodes->contains (*it))
[ + + ]
498 [ + - ][ + - ]: 1 : nanodes->append (*it);
499 : : }
500 : :
501 [ + - ][ + - ]: 1 : banodes = new strlist (*nlnodes);
[ + - ]
502 : : #if DEBUG && 0
503 : : fprintf (stderr, " balanced nodes: [ %s ]\n", banodes->toString ());
504 : : fprintf (stderr, " exnodes nodes: [ %s ]\n", exnodes->toString ());
505 : : fprintf (stderr, " nanodes nodes: [ %s ]\n", nanodes->toString ());
506 : : #endif
507 : 1 : }
508 : :
509 : : /* The function applies unique voltage source identifiers to each
510 : : built in internal voltage source in the list of circuits. */
511 : 1 : int hbsolver::assignVoltageSources (ptrlist<circuit> circuits) {
512 : 1 : int sources = 0;
513 [ + + ]: 3 : for (auto *c: circuits) {
514 [ + - ][ + - ]: 2 : if (c->getVoltageSources () > 0) {
515 : 2 : c->setVoltageSource (sources);
516 [ + - ]: 2 : sources += c->getVoltageSources ();
517 : : }
518 : : }
519 : 1 : return sources;
520 : : }
521 : :
522 : : /* The function assigns unique node identifiers to each circuit node. */
523 : 3 : int hbsolver::assignNodes (ptrlist<circuit> circuits, strlist * nodes,
524 : : int offset) {
525 : : // through all collected nodes
526 [ + + ]: 12 : for (int nr = 0; nr < nodes->length (); nr++) {
527 : 9 : char * nn = nodes->get (nr);
528 : : // through all circuits
529 [ + + ]: 21 : for (auto *c : circuits) {
530 : : // assign current identifier if any of the circuit nodes equals
531 : : // the current one
532 [ + + ]: 36 : for (int i = 0; i < c->getSize (); i++) {
533 [ + - ]: 24 : node * n = c->getNode (i);
534 [ + + ]: 24 : if (!strcmp (n->getName (), nn)) n->setNode (offset + nr + 1);
535 : : }
536 : : }
537 : : }
538 : 3 : return nodes->length ();
539 : : }
540 : :
541 : : // Prepares the linear operations.
542 : 1 : void hbsolver::prepareLinear (void) {
543 [ + + ]: 3 : for (auto *lc : lincircuits)
544 [ + - ]: 2 : lc->initHB ();
545 [ + - ]: 1 : nlnvsrcs = assignVoltageSources (lincircuits);
546 : 1 : nnlvsrcs = excitations.size ();
547 : 1 : nnanodes = nanodes->length ();
548 : 1 : nexnodes = exnodes->length ();
549 : 1 : nbanodes = banodes->length ();
550 [ + - ]: 1 : assignNodes (lincircuits, nanodes);
551 [ + - ]: 1 : assignNodes (excitations, nanodes);
552 : 1 : createMatrixLinearA ();
553 : 1 : createMatrixLinearY ();
554 : 1 : calcConstantCurrent ();
555 : 1 : }
556 : :
557 : : /* The function creates the complex linear network MNA matrix. It
558 : : contains the MNA entries for all linear components for each
559 : : requested frequency. */
560 : 1 : void hbsolver::createMatrixLinearA (void) {
561 : 1 : int M = nlnvsrcs;
562 : 1 : int N = nnanodes;
563 : 1 : int f = 0;
564 : : nr_double_t freq;
565 : :
566 : : // create new MNA matrix
567 [ + - ]: 1 : A = new tmatrix<nr_complex_t> ((N + M) * lnfreqs);
568 : :
569 : : // through each frequency
570 [ + + ]: 10 : for (int i = 0; i < rfreqs.getSize (); i++) {
571 : 9 : freq = rfreqs (i);
572 : : // calculate components' MNA matrix for the given frequency
573 [ + + ]: 27 : for (auto *lc : lincircuits)
574 [ + - ]: 18 : lc->calcHB (freq);
575 : : // fill in all matrix entries for the given frequency
576 : 9 : fillMatrixLinearA (A, f++);
577 : : }
578 : :
579 : : // save a copy of the original MNA matrix
580 [ + - ]: 1 : NA = new tmatrix<nr_complex_t> (*A);
581 : 1 : }
582 : :
583 : : // some definitions for the linear matrix filler
584 : : #undef A_
585 : : #undef B_
586 : : #define A_(r,c) (*A) ((r)*lnfreqs+f,(c)*lnfreqs+f)
587 : : #define G_(r,c) A_(r,c)
588 : : #define B_(r,c) A_(r,c+N)
589 : : #define C_(r,c) A_(r+N,c)
590 : : #define D_(r,c) A_(r+N,c+N)
591 : :
592 : : /* This function fills in the MNA matrix entries into the A matrix for
593 : : a given frequency index. */
594 : 9 : void hbsolver::fillMatrixLinearA (tmatrix<nr_complex_t> * A, int f) {
595 : 9 : int N = nnanodes;
596 : :
597 : : // through each linear circuit
598 [ + + ]: 27 : for (auto *cir : lincircuits) {
599 : 18 : int s = cir->getSize ();
600 : : int nr, nc, r, c, v;
601 : :
602 : : // apply G-matrix entries
603 [ + + ]: 54 : for (r = 0; r < s; r++) {
604 [ + - ][ - + ]: 36 : if ((nr = cir->getNode(r)->getNode () - 1) < 0) continue;
605 [ + + ]: 108 : for (c = 0; c < s; c++) {
606 [ + - ][ - + ]: 72 : if ((nc = cir->getNode(c)->getNode () - 1) < 0) continue;
607 [ + - ][ + - ]: 72 : G_(nr, nc) += cir->getY (r, c);
608 : : }
609 : : }
610 : :
611 : : // augmented part -- built in voltage sources
612 [ + - ][ + - ]: 18 : if ((v = cir->getVoltageSources ()) > 0) {
613 : :
614 : : // apply B-matrix entries
615 [ + + ]: 54 : for (r = 0; r < s; r++) {
616 [ + - ][ - + ]: 36 : if ((nr = cir->getNode(r)->getNode () - 1) < 0) continue;
617 [ + + ]: 72 : for (c = 0; c < v; c++) {
618 : 36 : nc = cir->getVoltageSource () + c;
619 [ + - ][ + - ]: 36 : B_(nr, nc) += cir->getB (r, nc);
620 : : }
621 : : }
622 : :
623 : : // apply C-matrix entries
624 [ + + ]: 36 : for (r = 0; r < v; r++) {
625 : 18 : nr = cir->getVoltageSource () + r;
626 [ + + ]: 54 : for (c = 0; c < s; c++) {
627 [ + - ][ - + ]: 36 : if ((nc = cir->getNode(c)->getNode () - 1) < 0) continue;
628 [ + - ][ + - ]: 36 : C_(nr, nc) += cir->getC (nr, c);
629 : : }
630 : : }
631 : :
632 : : // apply D-matrix entries
633 [ + + ]: 36 : for (r = 0; r < v; r++) {
634 : 18 : nr = cir->getVoltageSource () + r;
635 [ + + ]: 36 : for (c = 0; c < v; c++) {
636 : 18 : nc = cir->getVoltageSource () + c;
637 [ + - ][ + - ]: 18 : D_(nr, nc) += cir->getD (nr, nc);
638 : : }
639 : : }
640 : : }
641 : : }
642 : 9 : }
643 : :
644 : : // The function inverts the given matrix A into the matrix H.
645 : 1 : void hbsolver::invertMatrix (tmatrix<nr_complex_t> * A,
646 : : tmatrix<nr_complex_t> * H) {
647 : 1 : eqnsys<nr_complex_t> eqns;
648 : 1 : int N = A->getCols ();
649 [ + - ][ + - ]: 1 : tvector<nr_complex_t> * x = new tvector<nr_complex_t> (N);
650 [ + - ][ + - ]: 1 : tvector<nr_complex_t> * z = new tvector<nr_complex_t> (N);
651 : :
652 : : try_running () {
653 : : // create LU decomposition of the A matrix
654 : 1 : eqns.setAlgo (ALGO_LU_FACTORIZATION_CROUT);
655 [ + - ]: 1 : eqns.passEquationSys (A, x, z);
656 [ + - ]: 1 : eqns.solve ();
657 : : }
658 : : // appropriate exception handling
659 [ + - ][ - + ]: 1 : catch_exception () {
[ # # ]
660 : : case EXCEPTION_PIVOT:
661 : : default:
662 [ # # ]: 0 : logprint (LOG_ERROR, "WARNING: %s: during TI inversion\n", getName ());
663 [ # # ]: 0 : estack.print ();
664 : : }
665 : :
666 : : // use the LU decomposition to obtain the inverse H
667 : 1 : eqns.setAlgo (ALGO_LU_SUBSTITUTION_CROUT);
668 [ + + ]: 19 : for (int c = 0; c < N; c++) {
669 : 18 : z->set (0.0);
670 [ + - ]: 18 : z->set (c, 1.0);
671 [ + - ]: 18 : eqns.passEquationSys (A, x, z);
672 [ + - ]: 18 : eqns.solve ();
673 [ + - ][ + - ]: 342 : for (int r = 0; r < N; r++) H->set (r, c, x->get (r));
[ + + ]
674 : : }
675 [ + - ]: 1 : delete x;
676 [ + - ]: 1 : delete z;
677 : 1 : }
678 : :
679 : : // Some defines for matrix element access.
680 : : #define V_(r) (*V) (r)
681 : : #define I_(r) (*I) (r)
682 : : #undef A_
683 : : #define A_(r,c) (*A) (r,c)
684 : :
685 : : #define Z_(r,c) (*Z) (r,c)
686 : : #define Y_(r,c) (*Y) (r,c)
687 : :
688 : : #define ZVU_(r,c) Z_(r,c)
689 : : #define ZVL_(r,c) Z_((r)*lnfreqs+f+sn,c)
690 : : #define ZCU_(r,c) Z_(r,(c)*lnfreqs+f+sn)
691 : : #define ZCL_(r,c) Z_((r)*lnfreqs+f+sn,(c)*lnfreqs+f+sn)
692 : :
693 : : #define YV_(r,c) (*YV) (r,c)
694 : : #define NA_(r,c) (*NA) (r,c)
695 : : #define JF_(r,c) (*JF) (r,c)
696 : :
697 : : /* The following function performs the following steps:
698 : : 1. form the MNA matrix A including all nodes (linear, non-linear and
699 : : excitations)
700 : : 2. compute the variable transimpedance matrix entries for the nodes
701 : : to be balanced
702 : : 3. compute the constant transimpedance matrix entries for the constant
703 : : current vector caused by the excitations
704 : : 4. invert this overall transimpedance matrix
705 : : 5. extract the variable transadmittance matrix entries
706 : : */
707 : 1 : void hbsolver::createMatrixLinearY (void) {
708 : 1 : int M = nlnvsrcs;
709 : 1 : int N = nnanodes;
710 : : int c, r, f;
711 : :
712 : : // size of overall MNA matrix
713 : 1 : int sa = (N + M) * lnfreqs;
714 : 1 : int sv = nbanodes;
715 : 1 : int se = nnlvsrcs;
716 : 1 : int sy = sv + se;
717 : :
718 : : // allocate new transimpedance matrix
719 [ + - ][ + - ]: 1 : Z = new tmatrix<nr_complex_t> (sy * lnfreqs);
720 : :
721 : : // prepare equation system
722 : 1 : eqnsys<nr_complex_t> eqns;
723 : : tvector<nr_complex_t> * V;
724 : : tvector<nr_complex_t> * I;
725 : :
726 : : // 1. create variable transimpedance matrix entries relating
727 : : // voltages at the balanced nodes to the currents through these
728 : : // nodes into the non-linear part
729 : 1 : int sn = sv * lnfreqs;
730 [ + - ][ + - ]: 1 : V = new tvector<nr_complex_t> (sa);
731 [ + - ][ + - ]: 1 : I = new tvector<nr_complex_t> (sa);
732 : :
733 : : // connect a 100 Ohm resistor (to ground) to balanced node in the MNA matrix
734 [ + - ][ + + ]: 10 : for (c = 0; c < sv * lnfreqs; c++) A_(c, c) += 0.01;
735 : :
736 : : // connect a 100 Ohm resistor (in parallel) to each excitation
737 [ + + ]: 2 : for (auto *vs : excitations) {
738 : : // get positive and negative node
739 [ + - ]: 1 : int pnode = vs->getNode(NODE_1)->getNode ();
740 [ + - ]: 1 : int nnode = vs->getNode(NODE_2)->getNode ();
741 [ + + ]: 10 : for (f = 0; f < lnfreqs; f++) { // for each frequency
742 : 9 : int pn = (pnode - 1) * lnfreqs + f;
743 : 9 : int nn = (nnode - 1) * lnfreqs + f;
744 [ + - ][ + - ]: 9 : if (pnode) A_(pn, pn) += 0.01;
745 [ - + ][ # # ]: 9 : if (nnode) A_(nn, nn) += 0.01;
746 [ + - ][ - + ]: 9 : if (pnode && nnode) {
747 [ # # ]: 0 : A_(pn, nn) -= 0.01;
748 [ # # ]: 0 : A_(nn, pn) -= 0.01;
749 : : }
750 : : }
751 : : }
752 : :
753 : : // LU decompose the MNA matrix
754 : : try_running () {
755 : 1 : eqns.setAlgo (ALGO_LU_FACTORIZATION_CROUT);
756 [ + - ]: 1 : eqns.passEquationSys (A, V, I);
757 [ + - ]: 1 : eqns.solve ();
758 : : }
759 : : // appropriate exception handling
760 [ + - ][ - + ]: 1 : catch_exception () {
[ # # ]
761 : : case EXCEPTION_PIVOT:
762 : : default:
763 [ # # ]: 0 : logprint (LOG_ERROR, "WARNING: %s: during A factorization\n", getName ());
764 [ # # ]: 0 : estack.print ();
765 : : }
766 : :
767 : : // aquire variable transimpedance matrix entries
768 : 1 : eqns.setAlgo (ALGO_LU_SUBSTITUTION_CROUT);
769 [ + + ]: 10 : for (c = 0; c < sn; c++) {
770 : 9 : I->set (0.0);
771 [ + - ]: 9 : I_(c) = 1.0;
772 [ + - ]: 9 : eqns.passEquationSys (A, V, I);
773 [ + - ]: 9 : eqns.solve ();
774 : : // ZV | ..
775 : : // ---+---
776 : : // .. | ..
777 [ + - ][ + - ]: 90 : for (r = 0; r < sn; r++) ZVU_(r, c) = V_(r);
[ + + ]
778 : : // .. | ..
779 : : // ---+---
780 : : // ZV | ..
781 : 9 : r = 0;
782 [ + + ]: 18 : for (auto *e : excitations) {
783 : : // lower part entries
784 [ + + ]: 90 : for (f = 0; f < lnfreqs; f++) {
785 [ + - ][ + - ]: 81 : ZVL_(r, c) = excitationZ (V, e, f);
786 : : }
787 : : }
788 : : }
789 : :
790 : : // create constant transimpedance matrix entries relating the
791 : : // source voltages to the interconnection currents
792 : 1 : int vsrc = 0;
793 : : // aquire constant transadmittance matrix entries
794 [ + + ]: 2 : for (auto it = excitations.begin(); it != excitations.end(); ++it, vsrc++) {
795 : 1 : circuit * vs = *it;
796 : : // get positive and negative node
797 [ + - ]: 1 : int pnode = vs->getNode(NODE_1)->getNode ();
798 [ + - ]: 1 : int nnode = vs->getNode(NODE_2)->getNode ();
799 [ + + ]: 10 : for (f = 0; f < lnfreqs; f++) { // for each frequency
800 : 9 : int pn = (pnode - 1) * lnfreqs + f;
801 : 9 : int nn = (nnode - 1) * lnfreqs + f;
802 : 9 : I->set (0.0);
803 [ + - ][ + - ]: 9 : if (pnode) I_(pn) = +1.0;
804 [ - + ][ # # ]: 9 : if (nnode) I_(nn) = -1.0;
805 [ + - ]: 9 : eqns.passEquationSys (A, V, I);
806 [ + - ]: 9 : eqns.solve ();
807 : : // .. | ZC
808 : : // ---+---
809 : : // .. | ..
810 [ + + ]: 90 : for (r = 0; r < sn; r++) {
811 : : // upper part of the entries
812 [ + - ][ + - ]: 81 : ZCU_(r, vsrc) = V_(r);
813 : : }
814 : : // .. | ..
815 : : // ---+---
816 : : // .. | ZC
817 : 9 : r = 0;
818 [ + + ]: 18 : for (auto ite = excitations.begin(); ite != excitations.end(); ++ite, r++) {
819 : : // lower part entries
820 [ + - ][ + - ]: 9 : ZCL_(r, vsrc) = excitationZ (V, *ite, f);
821 : : }
822 : : }
823 : : }
824 [ + - ]: 1 : delete I;
825 [ + - ]: 1 : delete V;
826 : :
827 : : // allocate new transadmittance matrix
828 [ + - ][ + - ]: 1 : Y = new tmatrix<nr_complex_t> (sy * lnfreqs);
829 : :
830 : : // invert the Z matrix to a Y matrix
831 [ + - ]: 1 : invertMatrix (Z, Y);
832 : :
833 : : // substract the 100 Ohm resistor
834 [ + - ][ + + ]: 19 : for (c = 0; c < sy * lnfreqs; c++) Y_(c, c) -= 0.01;
835 : :
836 : : // extract the variable transadmittance matrix
837 [ + - ][ + - ]: 1 : YV = new tmatrix<nr_complex_t> (sv * nlfreqs);
838 : :
839 : : // variable transadmittance matrix must be continued conjugately
840 [ + - ][ + - ]: 1 : *YV = expandMatrix (*Y, sv);
[ + - ]
841 : :
842 : : // delete overall temporary MNA matrix
843 [ + - ]: 1 : delete A; A = NULL;
844 : : // delete transimpedance matrix
845 [ + - ]: 1 : delete Z; Z = NULL;
846 : 1 : }
847 : :
848 : : /* Little helper function obtaining a transimpedance value for the
849 : : given voltage source (excitation) and for a given frequency
850 : : index. */
851 : 90 : nr_complex_t hbsolver::excitationZ (tvector<nr_complex_t> * V, circuit * vs,
852 : : int f) {
853 : : // get positive and negative node
854 [ + - ]: 90 : int pnode = vs->getNode(NODE_1)->getNode ();
855 [ + - ]: 90 : int nnode = vs->getNode(NODE_2)->getNode ();
856 : 90 : nr_complex_t z = 0.0;
857 [ + - ][ + - ]: 90 : if (pnode) z += V_((pnode - 1) * lnfreqs + f);
858 [ - + ][ # # ]: 90 : if (nnode) z -= V_((nnode - 1) * lnfreqs + f);
859 : 90 : return z;
860 : : }
861 : :
862 : : /* This function computes the constant current vectors using the
863 : : voltage of the excitations and the transadmittance matrix
864 : : entries. */
865 : 1 : void hbsolver::calcConstantCurrent (void) {
866 : 1 : int se = nnlvsrcs * lnfreqs;
867 : 1 : int sn = nbanodes * lnfreqs;
868 : 1 : int r, c, vsrc = 0;
869 : :
870 : : // collect excitation voltages
871 [ + - ]: 1 : tvector<nr_complex_t> VC (se);
872 [ + + ]: 2 : for (auto it = excitations.begin(); it != excitations.end(); ++it, vsrc++) {
873 : 1 : circuit * vs = *it;
874 [ + - ]: 1 : vs->initHB ();
875 : 1 : vs->setVoltageSource (0);
876 [ + + ]: 10 : for (int f = 0; f < rfreqs.getSize (); f++) { // for each frequency
877 [ + - ]: 9 : nr_double_t freq = rfreqs (f);
878 [ + - ]: 9 : vs->calcHB (freq);
879 [ + - ][ + - ]: 9 : VC (vsrc * lnfreqs + f) = vs->getE (VSRC_1);
880 : : }
881 : : }
882 : :
883 : : // compute constant current vector for balanced nodes
884 [ + - ][ + - ]: 1 : IC = new tvector<nr_complex_t> (sn);
885 : : // .. | YC * VC
886 : : // ---+---
887 : : // .. | ..
888 [ + + ]: 10 : for (r = 0; r < sn; r++) {
889 : 9 : nr_complex_t i = 0.0;
890 [ + + ]: 90 : for (c = 0; c < se; c++) {
891 [ + - ][ + - ]: 81 : i += Y_(r, c + sn) * VC (c);
[ + - ]
892 : : }
893 : 9 : int f = r % lnfreqs;
894 [ + + ][ + + ]: 9 : if (f != 0 && f != lnfreqs - 1) i /= 2;
895 [ + - ]: 9 : IC->set (r, i);
896 : : }
897 : : // expand the constant current conjugate
898 [ + - ][ + - ]: 1 : *IC = expandVector (*IC, nbanodes);
[ + - ]
899 : :
900 : : // compute constant current vector for sources itself
901 [ + - ][ + - ]: 1 : IS = new tvector<nr_complex_t> (se);
902 : : // .. | ..
903 : : // ---+---
904 : : // .. | YC * VC
905 [ + + ]: 10 : for (r = 0; r < se; r++) {
906 : 9 : nr_complex_t i = 0.0;
907 [ + + ]: 90 : for (c = 0; c < se; c++) {
908 [ + - ][ + - ]: 81 : i += Y_(r + sn, c + sn) * VC (c);
[ + - ]
909 : : }
910 [ + - ]: 9 : IS->set (r, i);
911 : : }
912 : :
913 : : // delete overall transadmittance matrix
914 [ + - ]: 1 : delete Y; Y = NULL;
915 : 1 : }
916 : :
917 : : /* Checks whether currents through the interconnects of the linear and
918 : : non-linear subcircuit (in the frequency domain) are equal. */
919 : 10 : int hbsolver::checkBalance (void) {
920 : 10 : nr_double_t iabstol = getPropertyDouble ("iabstol");
921 : 10 : nr_double_t vabstol = getPropertyDouble ("vabstol");
922 : 10 : nr_double_t reltol = getPropertyDouble ("reltol");
923 : 10 : int n, len = FV->getSize ();
924 [ + + ]: 35 : for (n = 0; n < len; n++) {
925 : : // check iteration voltages
926 [ + - ][ + - ]: 25 : nr_double_t v_abs = abs (VS->get (n) - VP->get (n));
927 [ + - ]: 25 : nr_double_t v_rel = abs (VS->get (n));
928 [ + + ]: 25 : if (v_abs >= vabstol + reltol * v_rel) return 0;
929 : : // check balanced currents
930 [ + - ]: 18 : nr_complex_t il = IL->get (n);
931 [ + - ]: 18 : nr_complex_t in = IN->get (n);
932 [ + - ][ + - ]: 18 : if (il != in) {
933 : 18 : nr_double_t i_abs = abs (il + in);
934 [ + - ]: 18 : nr_double_t i_rel = abs ((il + in) / (il - in));
935 [ + + ][ + + ]: 18 : if (i_abs >= iabstol && 2.0 * i_rel >= reltol) return 0;
936 : : }
937 : : }
938 : 10 : return 1;
939 : : }
940 : :
941 : : // some definitions for the non-linear matrix filler
942 : : #undef G_
943 : : #undef C_
944 : : #define G_(r,c) (*jg) ((r)*nlfreqs+f,(c)*nlfreqs+f)
945 : : #define C_(r,c) (*jq) ((r)*nlfreqs+f,(c)*nlfreqs+f)
946 : : #undef FI_
947 : : #undef FQ_
948 : : #define FI_(r) (*ig) ((r)*nlfreqs+f)
949 : : #define FQ_(r) (*fq) ((r)*nlfreqs+f)
950 : : #define IR_(r) (*ir) ((r)*nlfreqs+f)
951 : : #define QR_(r) (*qr) ((r)*nlfreqs+f)
952 : :
953 : : /* This function fills in the matrix and vector entries for the
954 : : non-linear HB equations for a given frequency index. */
955 : 176 : void hbsolver::fillMatrixNonLinear (tmatrix<nr_complex_t> * jg,
956 : : tmatrix<nr_complex_t> * jq,
957 : : tvector<nr_complex_t> * ig,
958 : : tvector<nr_complex_t> * fq,
959 : : tvector<nr_complex_t> * ir,
960 : : tvector<nr_complex_t> * qr,
961 : : int f) {
962 : : // through each linear circuit
963 [ + + ]: 352 : for (auto *cir: nolcircuits) {
964 : 176 : int s = cir->getSize ();
965 : : int nr, nc, r, c;
966 : :
967 [ + + ]: 528 : for (r = 0; r < s; r++) {
968 [ + - ][ + + ]: 352 : if ((nr = cir->getNode(r)->getNode () - 1) < 0) continue;
969 : : // apply G- and C-matrix entries
970 [ + + ]: 528 : for (c = 0; c < s; c++) {
971 [ + - ][ + + ]: 352 : if ((nc = cir->getNode(c)->getNode () - 1) < 0) continue;
972 [ + - ][ + - ]: 176 : G_(nr, nc) += cir->getY (r, c);
973 [ + - ][ + - ]: 176 : C_(nr, nc) += cir->getQV (r, c);
974 : : }
975 : : // apply I- and Q-vector entries
976 [ + - ][ + - ]: 176 : FI_(nr) -= cir->getI (r);
977 [ + - ][ + - ]: 176 : FQ_(nr) -= cir->getQ (r);
978 : : // ThinkME: positive or negative?
979 [ + - ][ + - ]: 176 : IR_(nr) += cir->getGV (r) + cir->getI (r);
[ + - ]
980 [ + - ][ + - ]: 176 : QR_(nr) += cir->getCV (r) + cir->getQ (r);
[ + - ]
981 : : }
982 : : }
983 : 176 : }
984 : :
985 : : /* The function initializes the non-linear part of the HB. */
986 : 1 : void hbsolver::prepareNonLinear (void) {
987 : 1 : int N = nbanodes;
988 : :
989 : : // allocate matrices and vectors
990 [ + - ]: 1 : if (FQ == NULL) {
991 [ + - ]: 1 : FQ = new tvector<nr_complex_t> (N * nlfreqs);
992 : : }
993 [ + - ]: 1 : if (IG == NULL) {
994 [ + - ]: 1 : IG = new tvector<nr_complex_t> (N * nlfreqs);
995 : : }
996 [ + - ]: 1 : if (IR == NULL) {
997 [ + - ]: 1 : IR = new tvector<nr_complex_t> (N * nlfreqs);
998 : : }
999 [ + - ]: 1 : if (QR == NULL) {
1000 [ + - ]: 1 : QR = new tvector<nr_complex_t> (N * nlfreqs);
1001 : : }
1002 [ + - ]: 1 : if (JG == NULL) {
1003 [ + - ]: 1 : JG = new tmatrix<nr_complex_t> (N * nlfreqs);
1004 : : }
1005 [ + - ]: 1 : if (JQ == NULL) {
1006 [ + - ]: 1 : JQ = new tmatrix<nr_complex_t> (N * nlfreqs);
1007 : : }
1008 [ + - ]: 1 : if (JF == NULL) {
1009 [ + - ]: 1 : JF = new tmatrix<nr_complex_t> (N * nlfreqs);
1010 : : }
1011 : :
1012 : : // voltage vector in frequency and time domain
1013 [ + - ]: 1 : if (VS == NULL) {
1014 [ + - ]: 1 : VS = new tvector<nr_complex_t> (N * nlfreqs);
1015 : : }
1016 [ + - ]: 1 : if (vs == NULL) {
1017 [ + - ]: 1 : vs = new tvector<nr_complex_t> (N * nlfreqs);
1018 : : }
1019 [ + - ]: 1 : if (VP == NULL) {
1020 [ + - ]: 1 : VP = new tvector<nr_complex_t> (N * nlfreqs);
1021 : : }
1022 : :
1023 : : // error vector
1024 [ + - ]: 1 : if (FV == NULL) {
1025 [ + - ]: 1 : FV = new tvector<nr_complex_t> (N * nlfreqs);
1026 : : }
1027 : : // right hand side vector
1028 [ + - ]: 1 : if (RH == NULL) {
1029 [ + - ]: 1 : RH = new tvector<nr_complex_t> (N * nlfreqs);
1030 : : }
1031 : :
1032 : : // linear and non-linear current vector
1033 [ + - ]: 1 : if (IL == NULL) {
1034 [ + - ]: 1 : IL = new tvector<nr_complex_t> (N * nlfreqs);
1035 : : }
1036 [ + - ]: 1 : if (IN == NULL) {
1037 [ + - ]: 1 : IN = new tvector<nr_complex_t> (N * nlfreqs);
1038 : : }
1039 : :
1040 : : // assign nodes
1041 [ + - ]: 1 : assignNodes (nolcircuits, nanodes);
1042 : :
1043 : : // initialize circuits
1044 [ + + ]: 2 : for (auto *cir : nolcircuits) {
1045 [ + - ]: 1 : cir->initHB (nlfreqs);
1046 : : }
1047 : 1 : }
1048 : :
1049 : : /* Saves the node voltages of the given circuit and for the given
1050 : : frequency entry into the circuit voltage vector. */
1051 : 176 : void hbsolver::saveNodeVoltages (circuit * cir, int f) {
1052 : 176 : int r, nr, s = cir->getSize ();
1053 [ + + ]: 528 : for (r = 0; r < s; r++) {
1054 [ + + ]: 352 : if ((nr = cir->getNode(r)->getNode () - 1) < 0) continue;
1055 : : // apply V-vector entries
1056 [ + - ]: 176 : cir->setV (r, real (vs->get (nr * nlfreqs + f)));
1057 : : }
1058 : 176 : }
1059 : :
1060 : : /* The function saves voltages into non-linear circuits, runs each
1061 : : non-linear components' HB calculator for each frequency and applies
1062 : : the matrix and vector entries appropriately. */
1063 : 11 : void hbsolver::loadMatrices (void) {
1064 : : // clear matrices and vectors before
1065 : 11 : IG->set (0.0);
1066 : 11 : FQ->set (0.0);
1067 : 11 : IR->set (0.0);
1068 : 11 : QR->set (0.0);
1069 : 11 : JG->set (0.0);
1070 : 11 : JQ->set (0.0);
1071 : : // through each frequency
1072 [ + + ]: 187 : for (int f = 0; f < nlfreqs; f++) {
1073 : : // calculate components' HB matrices and vector for the given frequency
1074 [ + + ]: 352 : for (auto *cir : nolcircuits) {
1075 [ + - ]: 176 : saveNodeVoltages (cir, f); // node voltages
1076 [ + - ]: 176 : cir->calcHB (f); // HB calculator
1077 : : }
1078 : : // fill in all matrix entries for the given frequency
1079 : 176 : fillMatrixNonLinear (JG, JQ, IG, FQ, IR, QR, f);
1080 : : }
1081 : 11 : }
1082 : :
1083 : : /* The following function transforms a vector using a Fast Fourier
1084 : : Transformation from the time domain to the frequency domain. */
1085 : 74 : void hbsolver::VectorFFT (tvector<nr_complex_t> * V, int isign) {
1086 : : int i, k, r;
1087 : 74 : int n = nlfreqs;
1088 : 74 : int nd = dfreqs.getSize ();
1089 : 74 : int nodes = V->getSize () / n;
1090 : 74 : nr_double_t * d = (nr_double_t *) V->getData ();
1091 : :
1092 [ + - ]: 74 : if (nd == 1) {
1093 : : // for each node a single 1d-FFT
1094 [ + + ]: 148 : for (k = i = 0; i < nodes; i++, k += 2 * n) {
1095 : 74 : nr_double_t * dst = &d[k];
1096 : 74 : _fft_1d (dst, n, isign);
1097 [ + + ][ + + ]: 2122 : if (isign > 0) for (r = 0; r < 2 * n; r++) *dst++ /= n;
1098 : : }
1099 : : }
1100 : : else {
1101 : : // for each node a single nd-FFT
1102 [ # # ]: 0 : for (k = i = 0; i < nodes; i++, k += 2 * n) {
1103 : 0 : nr_double_t * dst = &d[k];
1104 : 0 : _fft_nd (dst, ndfreqs, nd, isign);
1105 [ # # ][ # # ]: 0 : if (isign > 0) for (r = 0; r < 2 * n; r++) *dst++ /= ndfreqs[0];
1106 : : }
1107 : : }
1108 : 74 : }
1109 : :
1110 : : /* The following function transforms a vector using an Inverse Fast
1111 : : Fourier Transformation from the frequency domain to the domain
1112 : : time. */
1113 : 10 : void hbsolver::VectorIFFT (tvector<nr_complex_t> * V, int isign) {
1114 : 10 : VectorFFT (V, -isign);
1115 : 10 : }
1116 : :
1117 : : /* The following function transforms a matrix using a Fast Fourier
1118 : : Transformation from the time domain to the frequency domain. */
1119 : 20 : void hbsolver::MatrixFFT (tmatrix<nr_complex_t> * M) {
1120 : :
1121 : : #if THE_SLO_ALGO
1122 : : // each column FFT
1123 : : for (int c = 0; c < M->getCols (); c++) {
1124 : : tvector<nr_complex_t> V = M->getCol (c);
1125 : : VectorFFT (&V);
1126 : : M->setCol (c, V);
1127 : : }
1128 : :
1129 : : // each row IFFT
1130 : : for (int r = 0; r < M->getRows (); r++) {
1131 : : tvector<nr_complex_t> V = M->getRow (r);
1132 : : VectorIFFT (&V);
1133 : : M->setRow (r, V);
1134 : : }
1135 : : #else
1136 : : int c, r, nc, nr;
1137 : : // for each non-linear node block
1138 [ + + ]: 40 : for (nc = c = 0; c < nbanodes; c++, nc += nlfreqs) {
1139 [ + + ]: 40 : for (nr = r = 0; r < nbanodes; r++, nr += nlfreqs) {
1140 [ + - ]: 20 : tvector<nr_complex_t> V (nlfreqs);
1141 : : int fr, fc, fi;
1142 : : // transform the sub-diagonal only
1143 [ + - ][ + - ]: 340 : for (fc = 0; fc < nlfreqs; fc++) V (fc) = M->get (nr + fc, nc + fc);
[ + + ]
1144 [ + - ]: 20 : VectorFFT (&V);
1145 : : // fill in resulting sub-matrix for the node
1146 [ + + ]: 340 : for (fc = 0; fc < nlfreqs; fc++) {
1147 [ + + ]: 5440 : for (fi = nlfreqs - 1 - fc, fr = 0; fr < nlfreqs; fr++) {
1148 [ + + ]: 5120 : if (++fi >= nlfreqs) fi = 0;
1149 [ + - ][ + - ]: 5120 : M->set (nr + fr, nc + fc, V (fi));
1150 : : }
1151 : : }
1152 : 20 : }
1153 : : }
1154 : : #endif
1155 : 20 : }
1156 : :
1157 : : /* This function solves the actual HB equation in the frequency domain.
1158 : : F(V) = IC + [YV] * VS + j[O] * FQ + IG -> 0
1159 : : Care must be taken with indexing here: In the frequency domain only
1160 : : real positive frequencies are used and computed, but in the time
1161 : : domain we have more values because of the periodic continuation in
1162 : : the frequency domain.
1163 : : RHS = j[O] * CV + GV - (IC + j[O] * FQ + IG)
1164 : : Also the right hand side of the equation system for the new voltage
1165 : : vector is computed here. */
1166 : 11 : void hbsolver::solveHB (void) {
1167 : : // for each non-linear node
1168 [ + + ]: 22 : for (int r = 0; r < nbanodes * nlfreqs; ) {
1169 : : // for each frequency
1170 [ + + ]: 187 : for (int f = 0; f < nlfreqs; f++, r++) {
1171 : 176 : nr_complex_t il = 0.0, in = 0.0, ir = 0.0;
1172 : : // constant current vector due to sources
1173 [ + - ]: 176 : il += IC->get (r);
1174 : : // part 1 of right hand side vector
1175 : 176 : ir -= il;
1176 : : // transadmittance matrix multiplied by voltage vector
1177 [ + + ]: 2992 : for (int c = 0; c < nbanodes * nlfreqs; c++) {
1178 [ + - ][ + - ]: 2816 : il += YV_(r, c) * VS_(c);
[ + - ]
1179 : : }
1180 : : // charge vector
1181 [ + - ][ + - ]: 176 : in += OM_(f) * FQ->get (r);
[ + - ]
1182 : : // current vector
1183 [ + - ]: 176 : in += IG->get (r);
1184 : : // part 2, 3 and 4 of right hand side vector
1185 [ + - ]: 176 : ir += IR->get (r);
1186 [ + - ][ + - ]: 176 : ir += OM_(f) * QR->get (r);
[ + - ]
1187 : : // put values into result vectors
1188 [ + - ]: 176 : RH->set (r, ir);
1189 [ + - ]: 176 : FV->set (r, il + in);
1190 [ + - ]: 176 : IL->set (r, il);
1191 [ + - ]: 176 : IN->set (r, in);
1192 : : }
1193 : : }
1194 : 11 : }
1195 : :
1196 : : /* The function calculates the full Jacobian JF = [YV] + j[O] * JQ + JG */
1197 : 10 : void hbsolver::calcJacobian (void) {
1198 : : int c, r, fc, fr, rt, ct;
1199 : : /* add admittances of capacitance matrix JQ and non-linear
1200 : : admittances matrix JG into complete Jacobian JF */
1201 [ + + ]: 20 : for (c = 0; c < nbanodes; c++) {
1202 : 10 : ct = c * nlfreqs;
1203 [ + + ]: 170 : for (fc = 0; fc < nlfreqs; fc++, ct++) {
1204 [ + + ]: 320 : for (r = 0; r < nbanodes; r++) {
1205 : 160 : rt = r * nlfreqs;
1206 [ + + ]: 2720 : for (fr = 0; fr < nlfreqs; fr++, rt++) {
1207 [ + - ][ + - ]: 2560 : JF_(rt, ct) = JG->get (rt, ct) + JQ->get (rt, ct) * OM_(fr);
1208 : : }
1209 : : }
1210 : : }
1211 : : }
1212 [ + - ]: 10 : *JF += *YV; // add linear admittance matrix
1213 : 10 : }
1214 : :
1215 : : /* The function expands the given vector in the frequency domain to
1216 : : make it a real valued signal in the time domain. */
1217 : 1 : tvector<nr_complex_t> hbsolver::expandVector (tvector<nr_complex_t> V,
1218 : : int nodes) {
1219 : 1 : tvector<nr_complex_t> res (nodes * nlfreqs);
1220 : : int r, ff, rf, rt;
1221 [ + + ]: 2 : for (r = 0; r < nodes; r++) {
1222 : 1 : rt = r * nlfreqs;
1223 : 1 : rf = r * lnfreqs;
1224 : : // copy first part of vector
1225 [ + + ]: 10 : for (ff = 0; ff < lnfreqs; ff++, rf++, rt++) {
1226 [ + - ][ + - ]: 9 : res (rt) = V (rf);
1227 : : }
1228 : : // continue vector conjugated
1229 [ + + ]: 8 : for (rf -= 2; ff < nlfreqs; ff++, rf--, rt++) {
1230 [ + - ][ + - ]: 7 : res (rt) = conj (V (rf));
1231 : : }
1232 : : }
1233 : 1 : return res;
1234 : : }
1235 : :
1236 : : /* The function expands the given matrix in the frequency domain to
1237 : : make it a real valued signal in the time domain. */
1238 : 1 : tmatrix<nr_complex_t> hbsolver::expandMatrix (tmatrix<nr_complex_t> M,
1239 : : int nodes) {
1240 : 1 : tmatrix<nr_complex_t> res (nodes * nlfreqs);
1241 : : int r, c, rf, rt, cf, ct, ff;
1242 [ + + ]: 2 : for (r = 0; r < nodes; r++) {
1243 [ + + ]: 2 : for (c = 0; c < nodes; c++) {
1244 : 1 : rf = r * lnfreqs;
1245 : 1 : rt = r * nlfreqs;
1246 : 1 : cf = c * lnfreqs;
1247 : 1 : ct = c * nlfreqs;
1248 : : // copy first part of diagonal
1249 [ + + ]: 10 : for (ff = 0; ff < lnfreqs; ff++, cf++, ct++, rf++, rt++) {
1250 [ + - ][ + - ]: 9 : res (rt, ct) = M (rf, cf);
1251 : : }
1252 : : // continue diagonal conjugated
1253 [ + + ]: 8 : for (cf -= 2, rf -= 2; ff < nlfreqs; ff++, cf--, ct++, rf--, rt++) {
1254 [ + - ][ + - ]: 7 : res (rt, ct) = conj (M (rf, cf));
1255 : : }
1256 : : }
1257 : : }
1258 : 1 : return res;
1259 : : }
1260 : :
1261 : : /* This function solves the equation system
1262 : : JF * VS(n+1) = JF * VS(n) - FV
1263 : : in order to obtains a new voltage vector in the frequency domain. */
1264 : 10 : void hbsolver::solveVoltages (void) {
1265 : : // save previous iteration voltage
1266 [ + - ]: 10 : *VP = *VS;
1267 : :
1268 : : // setup equation system
1269 : 10 : eqnsys<nr_complex_t> eqns;
1270 : : try_running () {
1271 : : // use LU decomposition for solving
1272 : 10 : eqns.setAlgo (ALGO_LU_DECOMPOSITION);
1273 [ + - ]: 10 : eqns.passEquationSys (JF, VS, RH);
1274 [ + - ]: 10 : eqns.solve ();
1275 : : }
1276 : : // appropriate exception handling
1277 [ + - ][ - + ]: 10 : catch_exception () {
[ # # ]
1278 : : case EXCEPTION_PIVOT:
1279 : : default:
1280 [ # # ]: 0 : logprint (LOG_ERROR, "WARNING: %s: during NR iteration\n", getName ());
1281 [ # # ]: 0 : estack.print ();
1282 : : }
1283 : :
1284 : : // save new voltages in time domain vector
1285 [ + - ]: 10 : *vs = *VS;
1286 : 10 : }
1287 : :
1288 : : /* The following function extends the existing linear MNA matrix to
1289 : : contain the additional rows and columns for the excitation voltage
1290 : : sources. */
1291 : 1 : tmatrix<nr_complex_t> hbsolver::extendMatrixLinear (tmatrix<nr_complex_t> M,
1292 : : int nodes) {
1293 : 1 : int no = M.getCols ();
1294 : 1 : tmatrix<nr_complex_t> res (no + nodes * lnfreqs);
1295 : : // copy the existing part
1296 [ + + ]: 46 : for (int r = 0; r < no; r++) {
1297 [ + + ]: 2070 : for (int c = 0; c < no; c++) {
1298 [ + - ][ + - ]: 2025 : res (r, c) = M (r, c);
1299 : : }
1300 : : }
1301 : 1 : return res;
1302 : : }
1303 : :
1304 : : /* The function fills in the missing MNA entries for the excitation
1305 : : voltage sources into the extended rows and columns as well as the
1306 : : actual voltage values into the right hand side vector. */
1307 : 1 : void hbsolver::fillMatrixLinearExtended (tmatrix<nr_complex_t> * Y,
1308 : : tvector<nr_complex_t> * I) {
1309 : : // through each excitation source
1310 : 1 : int of = lnfreqs * (nlnvsrcs + nnanodes);
1311 : 1 : int sc = of;
1312 : :
1313 [ + + ]: 2 : for (auto *vs : excitations) {
1314 : : // get positive and negative node
1315 [ + - ]: 1 : int pnode = vs->getNode(NODE_1)->getNode ();
1316 [ + - ]: 1 : int nnode = vs->getNode(NODE_2)->getNode ();
1317 [ + + ]: 10 : for (int f = 0; f < lnfreqs; f++, sc++) { // for each frequency
1318 : : // fill right hand side vector
1319 [ + - ]: 9 : nr_double_t freq = rfreqs (f);
1320 [ + - ]: 9 : vs->calcHB (freq);
1321 [ + - ][ + - ]: 9 : I_(sc) = vs->getE (VSRC_1);
1322 : : // fill MNA entries
1323 : 9 : int pn = (pnode - 1) * lnfreqs + f;
1324 : 9 : int nn = (nnode - 1) * lnfreqs + f;
1325 [ + - ]: 9 : if (pnode) {
1326 [ + - ]: 9 : Y_(pn, sc) = +1.0;
1327 [ + - ]: 9 : Y_(sc, pn) = +1.0;
1328 : : }
1329 [ - + ]: 9 : if (nnode) {
1330 [ # # ]: 0 : Y_(nn, sc) = -1.0;
1331 [ # # ]: 0 : Y_(sc, nn) = -1.0;
1332 : : }
1333 : : }
1334 : : }
1335 : 1 : }
1336 : :
1337 : : /* The function calculates and saves the final solution. */
1338 : 1 : void hbsolver::finalSolution (void) {
1339 : :
1340 : : // extend the linear MNA matrix
1341 [ + - ][ + - ]: 1 : *NA = extendMatrixLinear (*NA, nnlvsrcs);
1342 : :
1343 : 1 : int S = NA->getCols ();
1344 : 1 : int N = nnanodes * lnfreqs;
1345 : :
1346 : : // right hand side vector
1347 [ + - ]: 1 : tvector<nr_complex_t> * I = new tvector<nr_complex_t> (S);
1348 : : // temporary solution
1349 [ + - ]: 1 : tvector<nr_complex_t> * V = new tvector<nr_complex_t> (S);
1350 : : // final solution
1351 [ + - ]: 1 : x = new tvector<nr_complex_t> (N);
1352 : :
1353 : : // fill in missing MNA entries
1354 : 1 : fillMatrixLinearExtended (NA, I);
1355 : :
1356 : : // put currents through balanced nodes into right hand side
1357 [ + + ]: 2 : for (int n = 0; n < nbanodes; n++) {
1358 [ + + ]: 10 : for (int f = 0; f < lnfreqs; f++) {
1359 [ + - ]: 9 : nr_complex_t i = IL->get (n * nlfreqs + f);
1360 [ + + ][ + + ]: 9 : if (f != 0 && f != lnfreqs - 1) i *= 2;
1361 [ + - ]: 9 : I_(n * lnfreqs + f) = i;
1362 : : }
1363 : : }
1364 : :
1365 : : // use QR decomposition for the final solution
1366 : : try_running () {
1367 : 1 : eqnsys<nr_complex_t> eqns;
1368 : 1 : eqns.setAlgo (ALGO_LU_DECOMPOSITION);
1369 [ + - ]: 1 : eqns.passEquationSys (NA, V, I);
1370 [ + - ]: 1 : eqns.solve ();
1371 : : }
1372 : : // appropriate exception handling
1373 [ - + ]: 1 : catch_exception () {
1374 : : case EXCEPTION_PIVOT:
1375 : : default:
1376 : : logprint (LOG_ERROR, "WARNING: %s: during final AC analysis\n",
1377 : 0 : getName ());
1378 : 0 : estack.print ();
1379 : : }
1380 [ + + ]: 28 : for (int i = 0; i < N; i++) x->set (i, V_(i));
1381 : 1 : }
1382 : :
1383 : : // Saves simulation results.
1384 : 1 : void hbsolver::saveResults (void) {
1385 : : vector * f;
1386 : : // add current frequency to the dependency of the output dataset
1387 [ + - ]: 1 : if ((f = data->findDependency ("hbfrequency")) == NULL) {
1388 [ + - ]: 1 : f = new vector ("hbfrequency");
1389 : 1 : data->addDependency (f);
1390 : : }
1391 : : // save frequency vector
1392 [ + - ]: 1 : if (runs == 1) {
1393 [ + - ][ + + ]: 10 : for (int i = 0; i < lnfreqs; i++) f->add (rfreqs (i));
1394 : : }
1395 : : // save node voltage vectors
1396 : 1 : int nanode = 0;
1397 [ + - ][ + - ]: 4 : for (strlistiterator it (nanodes); *it; ++it, nanode++) {
[ + - ][ + + ]
1398 [ + - ]: 3 : int l = strlen (*it);
1399 : 3 : char * n = (char *) malloc (l + 4);
1400 [ + - ]: 3 : sprintf (n, "%s.Vb", *it);
1401 [ + + ]: 30 : for (int i = 0; i < lnfreqs; i++) {
1402 [ + - ][ + - ]: 27 : saveVariable (n, x->get (i + nanode * lnfreqs), f);
1403 : : }
1404 [ + - ]: 1 : }
1405 : 1 : }
1406 : :
1407 : : // properties
1408 : : PROP_REQ [] = {
1409 : : { "n", PROP_INT, { 1, PROP_NO_STR }, PROP_MIN_VAL (1) },
1410 : : PROP_NO_PROP };
1411 : : PROP_OPT [] = {
1412 : : { "f", PROP_REAL, { 1e9, PROP_NO_STR }, PROP_POS_RANGEX },
1413 : : { "iabstol", PROP_REAL, { 1e-12, PROP_NO_STR }, PROP_RNG_X01I },
1414 : : { "vabstol", PROP_REAL, { 1e-6, PROP_NO_STR }, PROP_RNG_X01I },
1415 : : { "reltol", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_RNG_X01I },
1416 : : { "MaxIter", PROP_INT, { 150, PROP_NO_STR }, PROP_RNGII (2, 10000) },
1417 : : PROP_NO_PROP };
1418 : : struct define_t hbsolver::anadef =
1419 : : { "HB", 0, PROP_ACTION, PROP_NO_SUBSTRATE, PROP_LINEAR, PROP_DEF };
1420 : :
1421 : : } // namespace qucs
1422 : :
1423 : : /* TODO list for HB Solver:
1424 : : - Take care about nodes with non-linear components only.
1425 : : - AC Power Sources (extra Z and open loop voltage).
1426 : : - Current sources.
1427 : : - Balancing of multi-dimensional non-linear networks.
1428 : : - Sources directly connected to non-linear components and no other
1429 : : linear component (insert short).
1430 : : - Bug: With capacitors at hand there is voltage convergence but no
1431 : : current balancing.
1432 : : - Output currents through voltage sources.
1433 : : */
|
