Qucs-core  0.0.19
eqndefined.cpp
Go to the documentation of this file.
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 };