Qucs-core
0.0.19
|
00001 /* 00002 * eqndefined.cpp - equation defined device class implementation 00003 * 00004 * Copyright (C) 2007 Stefan Jahn <stefan@lkcc.org> 00005 * 00006 * This is free software; you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation; either version 2, or (at your option) 00009 * any later version. 00010 * 00011 * This software is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this package; see the file COPYING. If not, write to 00018 * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor, 00019 * Boston, MA 02110-1301, USA. 00020 * 00021 * $Id$ 00022 * 00023 */ 00024 00025 #if HAVE_CONFIG_H 00026 # include <config.h> 00027 #endif 00028 00029 #include "component.h" 00030 #include "equation.h" 00031 #include "environment.h" 00032 #include "device.h" 00033 #include "eqndefined.h" 00034 00035 using namespace qucs; 00036 using namespace qucs::eqn; 00037 00038 // Constructor for the equation defined device. 00039 eqndefined::eqndefined () : circuit () { 00040 type = CIR_EQNDEFINED; 00041 setVariableSized (true); 00042 veqn = NULL; 00043 ieqn = NULL; 00044 qeqn = NULL; 00045 geqn = NULL; 00046 ceqn = NULL; 00047 _jstat = NULL; 00048 _jdyna = NULL; 00049 _charges = NULL; 00050 } 00051 00052 // Destructor deletes equation defined device object from memory. 00053 eqndefined::~eqndefined () { 00054 free (veqn); 00055 free (ieqn); 00056 free (geqn); 00057 free (qeqn); 00058 free (ceqn); 00059 free (_jstat); 00060 free (_jdyna); 00061 free (_charges); 00062 } 00063 00064 // Callback for initializing the DC analysis. 00065 void eqndefined::initDC (void) { 00066 allocMatrixMNA (); 00067 if (ieqn == NULL) initModel (); 00068 doHB = false; 00069 } 00070 00071 // Some help macros. 00072 #define A(a) ((assignment *) (a)) 00073 #define C(c) ((constant *) (c)) 00074 #define BP(n) real (getV (n * 2 + 0) - getV (n * 2 + 1)) 00075 00076 // Creates a variable name from the given arguments. 00077 char * eqndefined::createVariable (const char * base, int n, bool pfx) { 00078 const char * str = strchr (getName (), '.'); 00079 if (str != NULL) 00080 str = strrchr (str, '.') + 1; 00081 else 00082 str = getName (); 00083 char * txt = (char *) malloc (strlen (str) + strlen (base) + 3); 00084 if (pfx) 00085 sprintf (txt, "%s.%s%d", str, base, n); 00086 else 00087 sprintf (txt, "%s%d", base, n); 00088 return txt; 00089 } 00090 00091 // Creates also a variable name from the given arguments. 00092 char * eqndefined::createVariable (const char * base, int r, int c, bool pfx) { 00093 const char * str = strchr (getName (), '.'); 00094 if (str != NULL) 00095 str = strrchr (str, '.') + 1; 00096 else 00097 str = getName (); 00098 char * txt = (char *) malloc (strlen (str) + strlen (base) + 4); 00099 if (pfx) 00100 sprintf (txt, "%s.%s%d%d", str, base, r, c); 00101 else 00102 sprintf (txt, "%s%d%d", base, r, c); 00103 return txt; 00104 } 00105 00106 // Saves the given value into the equation result. 00107 void eqndefined::setResult (void * eqn, nr_double_t val) { 00108 A(eqn)->evaluate (); 00109 constant * c = A(eqn)->getResult (); 00110 c->d = val; 00111 } 00112 00113 // Returns the result of the equation. 00114 nr_double_t eqndefined::getResult (void * eqn) { 00115 A(eqn)->evaluate (); 00116 return A(eqn)->getResultDouble (); 00117 } 00118 00119 // Initializes the equation defined device. 00120 void eqndefined::initModel (void) { 00121 int i, j, k, branches = getSize () / 2; 00122 char * in, * qn, * vn, * gn, * cn, * vnold, * inold; 00123 eqn::node * ivalue, * qvalue, * diff; 00124 00125 // allocate space for equation pointers 00126 veqn = (void **) malloc (sizeof (assignment *) * branches); 00127 ieqn = (void **) malloc (sizeof (assignment *) * branches); 00128 geqn = (void **) malloc (sizeof (assignment *) * branches * branches); 00129 qeqn = (void **) malloc (sizeof (assignment *) * branches); 00130 ceqn = (void **) malloc (sizeof (assignment *) * branches * branches); 00131 00132 // allocate space for Jacobians and charges 00133 _jstat = (nr_double_t *) malloc (sizeof (nr_double_t) * branches * branches); 00134 _jdyna = (nr_double_t *) malloc (sizeof (nr_double_t) * branches * branches); 00135 _charges = (nr_double_t *) malloc (sizeof (nr_double_t) * branches); 00136 00137 // first create voltage variables 00138 for (i = 0; i < branches; i++) { 00139 vn = createVariable ("V", i + 1); 00140 if ((veqn[i] = getEnv()->getChecker()->findEquation (vn)) == NULL) { 00141 veqn[i] = getEnv()->getChecker()->addDouble ("#voltage", vn, 0); 00142 A(veqn[i])->evalType (); 00143 A(veqn[i])->skip = 1; 00144 } 00145 free (vn); 00146 } 00147 00148 // prepare current and charge equations 00149 for (i = 0; i < branches; i++) { 00150 00151 // fetch current and charge equations 00152 in = createVariable ("I", i + 1); 00153 ivalue = getEnv()->getChecker()->findEquation (in); 00154 if (!ivalue) { 00155 logprint (LOG_ERROR, "ERROR: current equation `%s' not found for " 00156 "EDD `%s'\n", in, getName ()); 00157 } 00158 qn = createVariable ("Q", i + 1); 00159 qvalue = getEnv()->getChecker()->findEquation (qn); 00160 if (!qvalue) { 00161 logprint (LOG_ERROR, "ERROR: charge equation `%s' not found for " 00162 "EDD `%s'\n", qn, getName ()); 00163 } 00164 free (in); 00165 free (qn); 00166 00167 // replace voltage and current references 00168 for (j = 0; j < branches; j++) { 00169 in = createVariable ("I", j + 1); 00170 inold = createVariable ("I", j + 1, false); 00171 vn = createVariable ("V", j + 1); 00172 vnold = createVariable ("V", j + 1, false); 00173 if (ivalue) { 00174 ivalue->replace (vnold, vn); 00175 ivalue->replace (inold, in); 00176 } 00177 if (qvalue) { 00178 qvalue->replace (vnold, vn); 00179 qvalue->replace (inold, in); 00180 } 00181 free (vnold); free (vn); 00182 free (inold); free (in); 00183 } 00184 00185 // setup and save equations for currents and charges 00186 ieqn[i] = ivalue; 00187 qeqn[i] = qvalue; 00188 } 00189 00190 // evaluate types of currents and charges 00191 for (i = 0; i < branches; i++) { 00192 if (ieqn[i]) { A(ieqn[i])->evalType (); A(ieqn[i])->skip = 1; } 00193 if (qeqn[i]) { A(qeqn[i])->evalType (); A(qeqn[i])->skip = 1; } 00194 } 00195 00196 // create derivatives of currents 00197 for (k = 0, i = 0; i < branches; i++) { 00198 ivalue = A(ieqn[i]); 00199 00200 // create "static" differentiations 00201 for (j = 0; j < branches; j++, k++) { 00202 vn = createVariable ("V", j + 1); 00203 00204 // create conductance equations 00205 if (ivalue) { 00206 gn = createVariable ("G", i + 1, j + 1); 00207 if ((geqn[k] = getEnv()->getChecker()->findEquation (gn)) == NULL) { 00208 diff = ivalue->differentiate (vn); 00209 getEnv()->getChecker()->addEquation (diff); 00210 diff->evalType (); 00211 diff->skip = 1; 00212 geqn[k] = diff; 00213 A(diff)->rename (gn); 00214 } 00215 free (gn); 00216 #if DEBUG 00217 logprint (LOG_STATUS, "DEBUG: %s\n", A(geqn[k])->toString ()); 00218 #endif 00219 } 00220 else geqn[k] = NULL; 00221 00222 free (vn); 00223 } 00224 } 00225 00226 // create derivatives of charges 00227 for (k = 0, i = 0; i < branches; i++) { 00228 qvalue = A(qeqn[i]); 00229 00230 // create "dynamic" differentiations 00231 for (j = 0; j < branches; j++, k++) { 00232 vn = createVariable ("V", j + 1); 00233 00234 // create capacitance equations 00235 if (qvalue) { 00236 cn = createVariable ("C", i + 1, j + 1); 00237 if ((ceqn[k] = getEnv()->getChecker()->findEquation (cn)) == NULL) { 00238 diff = qvalue->differentiate (vn); 00239 getEnv()->getChecker()->addEquation (diff); 00240 diff->evalType (); 00241 ceqn[k] = diff; 00242 A(diff)->rename (cn); 00243 00244 // apply dQ/dI * dI/dV => dQ/dV derivatives 00245 for (int l = 0; l < branches; l++) { 00246 in = createVariable ("I", l + 1); 00247 diff = qvalue->differentiate (in); 00248 A(diff)->mul (A(geqn[l * branches + j])); 00249 A(ceqn[k])->add (A(diff)); 00250 delete diff; 00251 free (in); 00252 } 00253 A(ceqn[k])->evalType (); 00254 A(ceqn[k])->skip = 1; 00255 } 00256 free (cn); 00257 #if DEBUG 00258 logprint (LOG_STATUS, "DEBUG: %s\n", A(ceqn[k])->toString ()); 00259 #endif 00260 } 00261 else ceqn[k] = NULL; 00262 00263 free (vn); 00264 } 00265 } 00266 } 00267 00268 // Update local variable equations. 00269 void eqndefined::updateLocals (void) { 00270 int i, branches = getSize () / 2; 00271 00272 // update voltages for equations 00273 for (i = 0; i < branches; i++) { 00274 setResult (veqn[i], BP (i)); 00275 } 00276 // get local subcircuit values 00277 getEnv()->passConstants (); 00278 getEnv()->equationSolver (); 00279 } 00280 00281 // Callback for DC analysis. 00282 void eqndefined::calcDC (void) { 00283 int i, j, k, branches = getSize () / 2; 00284 00285 // update local equations 00286 updateLocals (); 00287 00288 // calculate currents and put into right-hand side 00289 for (i = 0; i < branches; i++) { 00290 nr_double_t c = getResult (ieqn[i]); 00291 setI (i * 2 + 0, -c); 00292 setI (i * 2 + 1, +c); 00293 } 00294 00295 // calculate derivatives and put into Jacobian and right-hand side 00296 for (k = 0, i = 0; i < branches; i++) { 00297 nr_double_t gv = 0; 00298 // usual G (dI/dV) entries 00299 for (j = 0; j < branches; j++, k++) { 00300 nr_double_t g = getResult (geqn[k]); 00301 setY (i * 2 + 0, j * 2 + 0, +g); 00302 setY (i * 2 + 1, j * 2 + 1, +g); 00303 setY (i * 2 + 0, j * 2 + 1, -g); 00304 setY (i * 2 + 1, j * 2 + 0, -g); 00305 gv += g * BP (j); 00306 } 00307 // during HB 00308 if (doHB) { 00309 setGV (i * 2 + 0, +gv); 00310 setGV (i * 2 + 1, -gv); 00311 } 00312 // DC and transient analysis 00313 else { 00314 addI (i * 2 + 0, +gv); 00315 addI (i * 2 + 1, -gv); 00316 } 00317 } 00318 } 00319 00320 // Evaluate operating points. 00321 void eqndefined::evalOperatingPoints (void) { 00322 int i, j, k, branches = getSize () / 2; 00323 00324 // save values for charges, conductances and capacitances 00325 for (k = 0, i = 0; i < branches; i++) { 00326 nr_double_t q = getResult (qeqn[i]); 00327 _charges[i] = q; 00328 for (j = 0; j < branches; j++, k++) { 00329 nr_double_t g = getResult (geqn[k]); 00330 _jstat[k] = g; 00331 nr_double_t c = getResult (ceqn[k]); 00332 _jdyna[k] = c; 00333 } 00334 } 00335 } 00336 00337 // Saves operating points. 00338 void eqndefined::saveOperatingPoints (void) { 00339 00340 // update local equations 00341 updateLocals (); 00342 00343 // save values for charges, conductances and capacitances 00344 evalOperatingPoints (); 00345 } 00346 00347 // Callback for initializing the AC analysis. 00348 void eqndefined::initAC (void) { 00349 allocMatrixMNA (); 00350 doHB = false; 00351 } 00352 00353 // Callback for AC analysis. 00354 void eqndefined::calcAC (nr_double_t frequency) { 00355 setMatrixY (calcMatrixY (frequency)); 00356 } 00357 00358 // Computes Y-matrix for AC analysis. 00359 matrix eqndefined::calcMatrixY (nr_double_t frequency) { 00360 int i, j, k, branches = getSize () / 2; 00361 matrix y (2 * branches); 00362 nr_double_t o = 2 * pi * frequency; 00363 00364 // merge conductances and capacitances 00365 for (k = 0, i = 0; i < branches; i++) { 00366 for (j = 0; j < branches; j++, k++) { 00367 int r = i * 2; 00368 int c = j * 2; 00369 nr_complex_t val = nr_complex_t (_jstat[k], o * _jdyna[k]); 00370 y (r + 0, c + 0) = +val; 00371 y (r + 1, c + 1) = +val; 00372 y (r + 0, c + 1) = -val; 00373 y (r + 1, c + 0) = -val; 00374 } 00375 } 00376 00377 return y; 00378 } 00379 00380 // Callback for initializing the TR analysis. 00381 void eqndefined::initTR (void) { 00382 int branches = getSize () / 2; 00383 setStates (2 * branches); 00384 initDC (); 00385 } 00386 00387 // Callback for the TR analysis. 00388 void eqndefined::calcTR (nr_double_t) { 00389 int state, i, j, k, branches = getSize () / 2; 00390 00391 // run usual DC iteration, then save operating points 00392 calcDC (); 00393 00394 // calculate Q and C 00395 evalOperatingPoints (); 00396 00397 // charge integrations 00398 for (i = 0; i < branches; i++) { 00399 int r = i * 2; 00400 state = 2 * i; 00401 transientCapacitanceQ (state, r + 0, r + 1, _charges[i]); 00402 } 00403 00404 // charge: 2-node, voltage: 2-node 00405 for (k = 0, i = 0; i < branches; i++) { 00406 for (j = 0; j < branches; j++, k++) { 00407 int r = i * 2; 00408 int c = j * 2; 00409 nr_double_t v = BP (j); 00410 transientCapacitanceC (r + 0, r + 1, c + 0, c + 1, _jdyna[k], v); 00411 } 00412 } 00413 } 00414 00415 // Callback for initializing the S-parameter analysis. 00416 void eqndefined::initSP (void) { 00417 allocMatrixS (); 00418 doHB = false; 00419 } 00420 00421 // Callback for S-parameter analysis. 00422 void eqndefined::calcSP (nr_double_t frequency) { 00423 setMatrixS (ytos (calcMatrixY (frequency))); 00424 } 00425 00426 // Callback for initializing the HB analysis. 00427 void eqndefined::initHB (int) { 00428 allocMatrixHB (); 00429 if (ieqn == NULL) initModel (); 00430 doHB = true; 00431 } 00432 00433 // Callback for HB analysis. 00434 void eqndefined::calcHB (int) { 00435 int i, j, k, branches = getSize () / 2; 00436 00437 // G's (dI/dV) into Y-Matrix and I's into I-Vector 00438 calcDC (); 00439 00440 // calculate Q and C 00441 evalOperatingPoints (); 00442 00443 // setup additional charge right-hand side 00444 for (i = 0; i < branches; i++) { 00445 setQ (i * 2 + 0, -_charges[i]); 00446 setQ (i * 2 + 1, +_charges[i]); 00447 } 00448 00449 // fill in C's (dQ/dV) in QV-Matrix and CV into right hand side 00450 for (k = 0, i = 0; i < branches; i++) { 00451 nr_double_t cv = 0; 00452 for (j = 0; j < branches; j++, k++) { 00453 int r = i * 2; 00454 int c = j * 2; 00455 nr_double_t val = _jdyna[k]; 00456 setQV (r + 0, c + 0, +val); 00457 setQV (r + 1, c + 1, +val); 00458 setQV (r + 0, c + 1, -val); 00459 setQV (r + 1, c + 0, -val); 00460 cv += val * BP (j); 00461 } 00462 setCV (i * 2 + 0, +cv); 00463 setCV (i * 2 + 1, -cv); 00464 } 00465 } 00466 00467 // properties 00468 PROP_REQ [] = { 00469 { "I1", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00470 { "Q1", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00471 PROP_NO_PROP }; 00472 PROP_OPT [] = { 00473 { "I2", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00474 { "Q2", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00475 { "I3", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00476 { "Q3", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00477 { "I4", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00478 { "Q4", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00479 { "I5", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00480 { "Q5", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00481 { "I6", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00482 { "Q6", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00483 { "I7", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00484 { "Q7", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00485 { "I8", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00486 { "Q8", PROP_REAL, { 0, PROP_NO_STR }, PROP_NO_RANGE }, 00487 PROP_NO_PROP }; 00488 struct define_t eqndefined::cirdef = 00489 { "EDD", 00490 PROP_NODES, PROP_COMPONENT, PROP_NO_SUBSTRATE, PROP_NONLINEAR, PROP_DEF };