Branch data Line data Source code
1 : : /*
2 : : * eqndefined.cpp - equation defined device class implementation
3 : : *
4 : : * Copyright (C) 2007 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 "component.h"
30 : : #include "equation.h"
31 : : #include "environment.h"
32 : : #include "device.h"
33 : : #include "eqndefined.h"
34 : :
35 : : using namespace qucs;
36 : : using namespace qucs::eqn;
37 : :
38 : : // Constructor for the equation defined device.
39 : 0 : eqndefined::eqndefined () : circuit () {
40 : 0 : type = CIR_EQNDEFINED;
41 : 0 : setVariableSized (true);
42 : 0 : veqn = NULL;
43 : 0 : ieqn = NULL;
44 : 0 : qeqn = NULL;
45 : 0 : geqn = NULL;
46 : 0 : ceqn = NULL;
47 : 0 : _jstat = NULL;
48 : 0 : _jdyna = NULL;
49 : 0 : _charges = NULL;
50 : 0 : }
51 : :
52 : : // Destructor deletes equation defined device object from memory.
53 : 0 : eqndefined::~eqndefined () {
54 [ # # ][ # # ]: 0 : if (veqn) free (veqn);
55 [ # # ][ # # ]: 0 : if (ieqn) free (ieqn);
56 [ # # ][ # # ]: 0 : if (geqn) free (geqn);
57 [ # # ][ # # ]: 0 : if (qeqn) free (qeqn);
58 [ # # ][ # # ]: 0 : if (ceqn) free (ceqn);
59 [ # # ][ # # ]: 0 : if (_jstat) free (_jstat);
60 [ # # ][ # # ]: 0 : if (_jdyna) free (_jdyna);
61 [ # # ][ # # ]: 0 : if (_charges) free (_charges);
62 [ # # ][ # # ]: 0 : }
63 : :
64 : : // Callback for initializing the DC analysis.
65 : 0 : void eqndefined::initDC (void) {
66 : 0 : allocMatrixMNA ();
67 [ # # ]: 0 : if (ieqn == NULL) initModel ();
68 : 0 : doHB = false;
69 : 0 : }
70 : :
71 : : // Some help macros.
72 : : #define A(a) ((assignment *) (a))
73 : : #define C(c) ((constant *) (c))
74 : : #define BP(n) real (getV (n * 2 + 0) - getV (n * 2 + 1))
75 : :
76 : : // Creates a variable name from the given arguments.
77 : 0 : char * eqndefined::createVariable (const char * base, int n, bool pfx) {
78 : 0 : char * str = strchr (getName (), '.');
79 [ # # ]: 0 : if (str != NULL)
80 : 0 : str = strrchr (str, '.') + 1;
81 : : else
82 : 0 : str = getName ();
83 : 0 : char * txt = (char *) malloc (strlen (str) + strlen (base) + 3);
84 [ # # ]: 0 : if (pfx)
85 : 0 : sprintf (txt, "%s.%s%d", str, base, n);
86 : : else
87 : 0 : sprintf (txt, "%s%d", base, n);
88 : 0 : return txt;
89 : : }
90 : :
91 : : // Creates also a variable name from the given arguments.
92 : 0 : char * eqndefined::createVariable (const char * base, int r, int c, bool pfx) {
93 : 0 : char * str = strchr (getName (), '.');
94 [ # # ]: 0 : if (str != NULL)
95 : 0 : str = strrchr (str, '.') + 1;
96 : : else
97 : 0 : str = getName ();
98 : 0 : char * txt = (char *) malloc (strlen (str) + strlen (base) + 4);
99 [ # # ]: 0 : if (pfx)
100 : 0 : sprintf (txt, "%s.%s%d%d", str, base, r, c);
101 : : else
102 : 0 : sprintf (txt, "%s%d%d", base, r, c);
103 : 0 : return txt;
104 : : }
105 : :
106 : : // Saves the given value into the equation result.
107 : 0 : void eqndefined::setResult (void * eqn, nr_double_t val) {
108 : 0 : A(eqn)->evaluate ();
109 : 0 : constant * c = A(eqn)->getResult ();
110 : 0 : c->d = val;
111 : 0 : }
112 : :
113 : : // Returns the result of the equation.
114 : 0 : nr_double_t eqndefined::getResult (void * eqn) {
115 : 0 : A(eqn)->evaluate ();
116 : 0 : return A(eqn)->getResultDouble ();
117 : : }
118 : :
119 : : // Initializes the equation defined device.
120 : 0 : void eqndefined::initModel (void) {
121 : 0 : int i, j, k, branches = getSize () / 2;
122 : : char * in, * qn, * vn, * gn, * cn, * vnold, * inold;
123 : : eqn::node * ivalue, * qvalue, * diff;
124 : :
125 : : // allocate space for equation pointers
126 : 0 : veqn = (void **) malloc (sizeof (assignment *) * branches);
127 : 0 : ieqn = (void **) malloc (sizeof (assignment *) * branches);
128 : 0 : geqn = (void **) malloc (sizeof (assignment *) * branches * branches);
129 : 0 : qeqn = (void **) malloc (sizeof (assignment *) * branches);
130 : 0 : ceqn = (void **) malloc (sizeof (assignment *) * branches * branches);
131 : :
132 : : // allocate space for Jacobians and charges
133 : 0 : _jstat = (nr_double_t *) malloc (sizeof (nr_double_t) * branches * branches);
134 : 0 : _jdyna = (nr_double_t *) malloc (sizeof (nr_double_t) * branches * branches);
135 : 0 : _charges = (nr_double_t *) malloc (sizeof (nr_double_t) * branches);
136 : :
137 : : // first create voltage variables
138 [ # # ]: 0 : for (i = 0; i < branches; i++) {
139 : 0 : vn = createVariable ("V", i + 1);
140 [ # # ]: 0 : if ((veqn[i] = getEnv()->getChecker()->findEquation (vn)) == NULL) {
141 : 0 : veqn[i] = getEnv()->getChecker()->addDouble ("#voltage", vn, 0);
142 : 0 : A(veqn[i])->evalType ();
143 : 0 : A(veqn[i])->skip = 1;
144 : : }
145 : 0 : free (vn);
146 : : }
147 : :
148 : : // prepare current and charge equations
149 [ # # ]: 0 : for (i = 0; i < branches; i++) {
150 : :
151 : : // fetch current and charge equations
152 : 0 : in = createVariable ("I", i + 1);
153 : 0 : ivalue = getEnv()->getChecker()->findEquation (in);
154 [ # # ]: 0 : if (!ivalue) {
155 : : logprint (LOG_ERROR, "ERROR: current equation `%s' not found for "
156 : 0 : "EDD `%s'\n", in, getName ());
157 : : }
158 : 0 : qn = createVariable ("Q", i + 1);
159 : 0 : qvalue = getEnv()->getChecker()->findEquation (qn);
160 [ # # ]: 0 : if (!qvalue) {
161 : : logprint (LOG_ERROR, "ERROR: charge equation `%s' not found for "
162 : 0 : "EDD `%s'\n", qn, getName ());
163 : : }
164 : 0 : free (in);
165 : 0 : free (qn);
166 : :
167 : : // replace voltage and current references
168 [ # # ]: 0 : for (j = 0; j < branches; j++) {
169 : 0 : in = createVariable ("I", j + 1);
170 : 0 : inold = createVariable ("I", j + 1, false);
171 : 0 : vn = createVariable ("V", j + 1);
172 : 0 : vnold = createVariable ("V", j + 1, false);
173 [ # # ]: 0 : if (ivalue) {
174 : 0 : ivalue->replace (vnold, vn);
175 : 0 : ivalue->replace (inold, in);
176 : : }
177 [ # # ]: 0 : if (qvalue) {
178 : 0 : qvalue->replace (vnold, vn);
179 : 0 : qvalue->replace (inold, in);
180 : : }
181 : 0 : free (vnold); free (vn);
182 : 0 : free (inold); free (in);
183 : : }
184 : :
185 : : // setup and save equations for currents and charges
186 : 0 : ieqn[i] = ivalue;
187 : 0 : qeqn[i] = qvalue;
188 : : }
189 : :
190 : : // evaluate types of currents and charges
191 [ # # ]: 0 : for (i = 0; i < branches; i++) {
192 [ # # ]: 0 : if (ieqn[i]) { A(ieqn[i])->evalType (); A(ieqn[i])->skip = 1; }
193 [ # # ]: 0 : if (qeqn[i]) { A(qeqn[i])->evalType (); A(qeqn[i])->skip = 1; }
194 : : }
195 : :
196 : : // create derivatives of currents
197 [ # # ]: 0 : for (k = 0, i = 0; i < branches; i++) {
198 : 0 : ivalue = A(ieqn[i]);
199 : :
200 : : // create "static" differentiations
201 [ # # ]: 0 : for (j = 0; j < branches; j++, k++) {
202 : 0 : vn = createVariable ("V", j + 1);
203 : :
204 : : // create conductance equations
205 [ # # ]: 0 : if (ivalue) {
206 : 0 : gn = createVariable ("G", i + 1, j + 1);
207 [ # # ]: 0 : if ((geqn[k] = getEnv()->getChecker()->findEquation (gn)) == NULL) {
208 : 0 : diff = ivalue->differentiate (vn);
209 : 0 : getEnv()->getChecker()->addEquation (diff);
210 : 0 : diff->evalType ();
211 : 0 : diff->skip = 1;
212 : 0 : geqn[k] = diff;
213 : 0 : A(diff)->rename (gn);
214 : : }
215 : 0 : free (gn);
216 : : #if DEBUG
217 : 0 : logprint (LOG_STATUS, "DEBUG: %s\n", A(geqn[k])->toString ());
218 : : #endif
219 : : }
220 : 0 : else geqn[k] = NULL;
221 : :
222 : 0 : free (vn);
223 : : }
224 : : }
225 : :
226 : : // create derivatives of charges
227 [ # # ]: 0 : for (k = 0, i = 0; i < branches; i++) {
228 : 0 : qvalue = A(qeqn[i]);
229 : :
230 : : // create "dynamic" differentiations
231 [ # # ]: 0 : for (j = 0; j < branches; j++, k++) {
232 : 0 : vn = createVariable ("V", j + 1);
233 : :
234 : : // create capacitance equations
235 [ # # ]: 0 : if (qvalue) {
236 : 0 : cn = createVariable ("C", i + 1, j + 1);
237 [ # # ]: 0 : if ((ceqn[k] = getEnv()->getChecker()->findEquation (cn)) == NULL) {
238 : 0 : diff = qvalue->differentiate (vn);
239 : 0 : getEnv()->getChecker()->addEquation (diff);
240 : 0 : diff->evalType ();
241 : 0 : ceqn[k] = diff;
242 : 0 : A(diff)->rename (cn);
243 : :
244 : : // apply dQ/dI * dI/dV => dQ/dV derivatives
245 [ # # ]: 0 : for (int l = 0; l < branches; l++) {
246 : 0 : in = createVariable ("I", l + 1);
247 : 0 : diff = qvalue->differentiate (in);
248 : 0 : A(diff)->mul (A(geqn[l * branches + j]));
249 : 0 : A(ceqn[k])->add (A(diff));
250 [ # # ]: 0 : delete diff;
251 : 0 : free (in);
252 : : }
253 : 0 : A(ceqn[k])->evalType ();
254 : 0 : A(ceqn[k])->skip = 1;
255 : : }
256 : 0 : free (cn);
257 : : #if DEBUG
258 : 0 : logprint (LOG_STATUS, "DEBUG: %s\n", A(ceqn[k])->toString ());
259 : : #endif
260 : : }
261 : 0 : else ceqn[k] = NULL;
262 : :
263 : 0 : free (vn);
264 : : }
265 : : }
266 : 0 : }
267 : :
268 : : // Update local variable equations.
269 : 0 : void eqndefined::updateLocals (void) {
270 : 0 : int i, branches = getSize () / 2;
271 : :
272 : : // update voltages for equations
273 [ # # ]: 0 : for (i = 0; i < branches; i++) {
274 [ # # ][ # # ]: 0 : setResult (veqn[i], BP (i));
[ # # ]
275 : : }
276 : : // get local subcircuit values
277 : 0 : getEnv()->passConstants ();
278 : 0 : getEnv()->equationSolver ();
279 : 0 : }
280 : :
281 : : // Callback for DC analysis.
282 : 0 : void eqndefined::calcDC (void) {
283 : 0 : int i, j, k, branches = getSize () / 2;
284 : :
285 : : // update local equations
286 : 0 : updateLocals ();
287 : :
288 : : // calculate currents and put into right-hand side
289 [ # # ]: 0 : for (i = 0; i < branches; i++) {
290 : 0 : nr_double_t c = getResult (ieqn[i]);
291 [ # # ]: 0 : setI (i * 2 + 0, -c);
292 [ # # ]: 0 : setI (i * 2 + 1, +c);
293 : : }
294 : :
295 : : // calculate derivatives and put into Jacobian and right-hand side
296 [ # # ]: 0 : for (k = 0, i = 0; i < branches; i++) {
297 : 0 : nr_double_t gv = 0;
298 : : // usual G (dI/dV) entries
299 [ # # ]: 0 : for (j = 0; j < branches; j++, k++) {
300 : 0 : nr_double_t g = getResult (geqn[k]);
301 [ # # ]: 0 : setY (i * 2 + 0, j * 2 + 0, +g);
302 [ # # ]: 0 : setY (i * 2 + 1, j * 2 + 1, +g);
303 [ # # ]: 0 : setY (i * 2 + 0, j * 2 + 1, -g);
304 [ # # ]: 0 : setY (i * 2 + 1, j * 2 + 0, -g);
305 [ # # ][ # # ]: 0 : gv += g * BP (j);
306 : : }
307 : : // during HB
308 [ # # ]: 0 : if (doHB) {
309 [ # # ]: 0 : setGV (i * 2 + 0, +gv);
310 [ # # ]: 0 : setGV (i * 2 + 1, -gv);
311 : : }
312 : : // DC and transient analysis
313 : : else {
314 : 0 : addI (i * 2 + 0, +gv);
315 : 0 : addI (i * 2 + 1, -gv);
316 : : }
317 : : }
318 : 0 : }
319 : :
320 : : // Evaluate operating points.
321 : 0 : void eqndefined::evalOperatingPoints (void) {
322 : 0 : int i, j, k, branches = getSize () / 2;
323 : :
324 : : // save values for charges, conductances and capacitances
325 [ # # ]: 0 : for (k = 0, i = 0; i < branches; i++) {
326 : 0 : nr_double_t q = getResult (qeqn[i]);
327 : 0 : _charges[i] = q;
328 [ # # ]: 0 : for (j = 0; j < branches; j++, k++) {
329 : 0 : nr_double_t g = getResult (geqn[k]);
330 : 0 : _jstat[k] = g;
331 : 0 : nr_double_t c = getResult (ceqn[k]);
332 : 0 : _jdyna[k] = c;
333 : : }
334 : : }
335 : 0 : }
336 : :
337 : : // Saves operating points.
338 : 0 : void eqndefined::saveOperatingPoints (void) {
339 : :
340 : : // update local equations
341 : 0 : updateLocals ();
342 : :
343 : : // save values for charges, conductances and capacitances
344 : 0 : evalOperatingPoints ();
345 : 0 : }
346 : :
347 : : // Callback for initializing the AC analysis.
348 : 0 : void eqndefined::initAC (void) {
349 : 0 : allocMatrixMNA ();
350 : 0 : doHB = false;
351 : 0 : }
352 : :
353 : : // Callback for AC analysis.
354 : 0 : void eqndefined::calcAC (nr_double_t frequency) {
355 [ # # ]: 0 : setMatrixY (calcMatrixY (frequency));
356 : 0 : }
357 : :
358 : : // Computes Y-matrix for AC analysis.
359 : 0 : matrix eqndefined::calcMatrixY (nr_double_t frequency) {
360 : 0 : int i, j, k, branches = getSize () / 2;
361 : 0 : matrix y (2 * branches);
362 : 0 : nr_double_t o = 2 * M_PI * frequency;
363 : :
364 : : // merge conductances and capacitances
365 [ # # ]: 0 : for (k = 0, i = 0; i < branches; i++) {
366 [ # # ]: 0 : for (j = 0; j < branches; j++, k++) {
367 : 0 : int r = i * 2;
368 : 0 : int c = j * 2;
369 : 0 : nr_complex_t val = nr_complex_t (_jstat[k], o * _jdyna[k]);
370 : 0 : y (r + 0, c + 0) = +val;
371 : 0 : y (r + 1, c + 1) = +val;
372 : 0 : y (r + 0, c + 1) = -val;
373 : 0 : y (r + 1, c + 0) = -val;
374 : : }
375 : : }
376 : :
377 : 0 : return y;
378 : : }
379 : :
380 : : // Callback for initializing the TR analysis.
381 : 0 : void eqndefined::initTR (void) {
382 : 0 : int branches = getSize () / 2;
383 : 0 : setStates (2 * branches);
384 : 0 : initDC ();
385 : 0 : }
386 : :
387 : : // Callback for the TR analysis.
388 : 0 : void eqndefined::calcTR (nr_double_t) {
389 : 0 : int state, i, j, k, branches = getSize () / 2;
390 : :
391 : : // run usual DC iteration, then save operating points
392 : 0 : calcDC ();
393 : :
394 : : // calculate Q and C
395 : 0 : evalOperatingPoints ();
396 : :
397 : : // charge integrations
398 [ # # ]: 0 : for (i = 0; i < branches; i++) {
399 : 0 : int r = i * 2;
400 : 0 : state = 2 * i;
401 : 0 : transientCapacitanceQ (state, r + 0, r + 1, _charges[i]);
402 : : }
403 : :
404 : : // charge: 2-node, voltage: 2-node
405 [ # # ]: 0 : for (k = 0, i = 0; i < branches; i++) {
406 [ # # ]: 0 : for (j = 0; j < branches; j++, k++) {
407 : 0 : int r = i * 2;
408 : 0 : int c = j * 2;
409 [ # # ][ # # ]: 0 : nr_double_t v = BP (j);
410 : 0 : transientCapacitanceC (r + 0, r + 1, c + 0, c + 1, _jdyna[k], v);
411 : : }
412 : : }
413 : 0 : }
414 : :
415 : : // Callback for initializing the S-parameter analysis.
416 : 0 : void eqndefined::initSP (void) {
417 : 0 : allocMatrixS ();
418 : 0 : doHB = false;
419 : 0 : }
420 : :
421 : : // Callback for S-parameter analysis.
422 : 0 : void eqndefined::calcSP (nr_double_t frequency) {
423 [ # # ][ # # ]: 0 : setMatrixS (ytos (calcMatrixY (frequency)));
[ # # ][ # # ]
[ # # ]
424 : 0 : }
425 : :
426 : : // Callback for initializing the HB analysis.
427 : 0 : void eqndefined::initHB (int) {
428 : 0 : allocMatrixHB ();
429 [ # # ]: 0 : if (ieqn == NULL) initModel ();
430 : 0 : doHB = true;
431 : 0 : }
432 : :
433 : : // Callback for HB analysis.
434 : 0 : void eqndefined::calcHB (int) {
435 : 0 : int i, j, k, branches = getSize () / 2;
436 : :
437 : : // G's (dI/dV) into Y-Matrix and I's into I-Vector
438 : 0 : calcDC ();
439 : :
440 : : // calculate Q and C
441 : 0 : evalOperatingPoints ();
442 : :
443 : : // setup additional charge right-hand side
444 [ # # ]: 0 : for (i = 0; i < branches; i++) {
445 [ # # ]: 0 : setQ (i * 2 + 0, -_charges[i]);
446 [ # # ]: 0 : setQ (i * 2 + 1, +_charges[i]);
447 : : }
448 : :
449 : : // fill in C's (dQ/dV) in QV-Matrix and CV into right hand side
450 [ # # ]: 0 : for (k = 0, i = 0; i < branches; i++) {
451 : 0 : nr_double_t cv = 0;
452 [ # # ]: 0 : for (j = 0; j < branches; j++, k++) {
453 : 0 : int r = i * 2;
454 : 0 : int c = j * 2;
455 : 0 : nr_double_t val = _jdyna[k];
456 [ # # ]: 0 : setQV (r + 0, c + 0, +val);
457 [ # # ]: 0 : setQV (r + 1, c + 1, +val);
458 [ # # ]: 0 : setQV (r + 0, c + 1, -val);
459 [ # # ]: 0 : setQV (r + 1, c + 0, -val);
460 [ # # ][ # # ]: 0 : cv += val * BP (j);
461 : : }
462 [ # # ]: 0 : setCV (i * 2 + 0, +cv);
463 [ # # ]: 0 : setCV (i * 2 + 1, -cv);
464 : : }
465 : 0 : }
466 : :
467 : : // properties
468 : : PROP_REQ [] = {
469 : : { "I1", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
470 : : { "Q1", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
471 : : PROP_NO_PROP };
472 : : PROP_OPT [] = {
473 : : { "I2", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
474 : : { "Q2", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
475 : : { "I3", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
476 : : { "Q3", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
477 : : { "I4", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
478 : : { "Q4", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
479 : : { "I5", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
480 : : { "Q5", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
481 : : { "I6", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
482 : : { "Q6", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
483 : : { "I7", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
484 : : { "Q7", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
485 : : { "I8", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
486 : : { "Q8", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE },
487 : : PROP_NO_PROP };
488 : : struct define_t eqndefined::cirdef =
489 : : { "EDD",
490 : : PROP_NODES, PROP_COMPONENT, PROP_NO_SUBSTRATE, PROP_NONLINEAR, PROP_DEF };
|
