Qucs-core  0.0.19
spsolver.cpp
Go to the documentation of this file.
00001 /*
00002  * spsolver.cpp - S-parameter solver class implementation
00003  *
00004  * Copyright (C) 2003-2008 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 <stdio.h>
00030 #include <stdlib.h>
00031 #include <string.h>
00032 
00033 #include "logging.h"
00034 #include "complex.h"
00035 #include "object.h"
00036 #include "node.h"
00037 #include "circuit.h"
00038 #include "strlist.h"
00039 #include "vector.h"
00040 #include "matvec.h"
00041 #include "dataset.h"
00042 #include "net.h"
00043 #include "analysis.h"
00044 #include "sweep.h"
00045 #include "nodelist.h"
00046 #include "netdefs.h"
00047 #include "characteristic.h"
00048 #include "spsolver.h"
00049 #include "constants.h"
00050 #include "components/component_id.h"
00051 #include "components/tee.h"
00052 #include "components/open.h"
00053 #include "components/itrafo.h"
00054 #include "components/cross.h"
00055 #include "components/ground.h"
00056 
00057 /* Evolved optimization flags. */
00058 #define USE_GROUNDS 1   // use extra grounds ?
00059 #define USE_CROSSES 1   // use additional cross connectors ?
00060 #define SORTED_LIST 1   // use sorted node list?
00061 
00062 #define TINYS (NR_TINY * 1.235) // 'tiny' value for singularities
00063 
00064 namespace qucs {
00065 
00066 // Constructor creates an unnamed instance of the spsolver class.
00067 spsolver::spsolver () : analysis () {
00068   type = ANALYSIS_SPARAMETER;
00069   swp = NULL;
00070   saveCVs = 0;
00071   noise = 0;
00072   nlist = NULL;
00073   tees = crosses = opens = grounds = 0;
00074   gnd = NULL;
00075 }
00076 
00077 // Constructor creates a named instance of the spsolver class.
00078 spsolver::spsolver (char * n) : analysis (n) {
00079   type = ANALYSIS_SPARAMETER;
00080   swp = NULL;
00081   saveCVs = 0;
00082   noise = 0;
00083   nlist = NULL;
00084   tees = crosses = opens = grounds = 0;
00085   gnd = NULL;
00086 }
00087 
00088 // Destructor deletes the spsolver class object.
00089 spsolver::~spsolver () {
00090   delete swp;
00091   delete nlist;
00092 }
00093 
00094 /* The copy constructor creates a new instance of the spsolver class
00095    based on the given spsolver object. */
00096 spsolver::spsolver (spsolver & n) : analysis (n) {
00097   tees = n.tees;
00098   crosses = n.crosses;
00099   opens = n.opens;
00100   grounds = n.grounds;
00101   noise = n.noise;
00102   saveCVs = n.saveCVs;
00103   swp = n.swp ? new sweep (*n.swp) : NULL;
00104   nlist = n.nlist ? new nodelist (*n.nlist) : NULL;
00105   gnd = n.gnd;
00106 }
00107 
00108 /* This function joins two nodes of a single circuit (interconnected
00109    nodes) and returns the resulting circuit. */
00110 circuit * spsolver::interconnectJoin (node * n1, node * n2) {
00111 
00112   circuit * s = n1->getCircuit ();
00113   circuit * result = new circuit (s->getSize () - 2);
00114   nr_complex_t p;
00115 
00116   // allocate S-parameter and noise corellation matrices
00117   result->initSP (); if (noise) result->initNoiseSP ();
00118 
00119   // interconnected port numbers
00120   int k = n1->getPort (), l = n2->getPort ();
00121 
00122   // denominator needs to be calculated only once
00123   nr_complex_t d = (1.0 - s->getS (k, l)) * (1.0 - s->getS (l, k)) -
00124     s->getS (k, k) * s->getS (l, l);
00125 
00126   // avoid singularity when two full reflective ports are interconnected
00127   nr_double_t tiny1 = (d == 0) ? 1.0 - TINYS : 1.0;
00128   nr_double_t tiny2 = tiny1 * tiny1;
00129   nr_double_t tiny3 = tiny1 * tiny2;
00130   d = (1.0 - s->getS (k, l) * tiny1) * (1.0 - s->getS (l, k) * tiny1) -
00131     s->getS (k, k) * s->getS (l, l) * tiny2;
00132 
00133   int j2; // column index for resulting matrix
00134   int i2; // row index for resulting matrix
00135   int j1; // column index for S matrix
00136   int i1; // row index for S matrix
00137 
00138   // handle single S block only
00139   i2 = j2 = 0;
00140   for (j1 = 0; j1 < s->getSize (); j1++) {
00141 
00142     // skip connected node
00143     if (j1 == k || j1 == l) continue;
00144 
00145     // assign node name of resulting circuit
00146     result->setNode (j2, s->getNode(j1)->getName ());
00147 
00148     // inside S only
00149     for (i1 = 0; i1 < s->getSize (); i1++) {
00150 
00151       // skip connected node
00152       if (i1 == k || i1 == l) continue;
00153 
00154       // compute S'ij
00155       p = s->getS (i1, j1);
00156       p +=
00157         (s->getS (k, j1) * s->getS (i1, l) * (1.0 - s->getS (l, k)) * tiny3 +
00158          s->getS (l, j1) * s->getS (i1, k) * (1.0 - s->getS (k, l)) * tiny3 +
00159          s->getS (k, j1) * s->getS (l, l) * s->getS (i1, k) * tiny3 +
00160          s->getS (l, j1) * s->getS (k, k) * s->getS (i1, l) * tiny3) / d;
00161       result->setS (i2++, j2, p);
00162     }
00163 
00164     // next column
00165     j2++; i2 = 0;
00166   }
00167   return result;
00168 }
00169 
00170 /* This function joins two nodes of two different circuits (connected
00171    nodes) and returns the resulting circuit. */
00172 circuit * spsolver::connectedJoin (node * n1, node * n2) {
00173 
00174   circuit * s = n1->getCircuit ();
00175   circuit * t = n2->getCircuit ();
00176   circuit * result = new circuit (s->getSize () + t->getSize () - 2);
00177   nr_complex_t p;
00178 
00179   // allocate S-parameter and noise corellation matrices
00180   result->initSP (); if (noise) result->initNoiseSP ();
00181 
00182   // connected port numbers
00183   int k = n1->getPort (), l = n2->getPort ();
00184 
00185   // denominator needs to be calculated only once
00186   nr_complex_t d = 1.0 - s->getS (k, k) * t->getS (l, l);
00187 
00188   // avoid singularity when two full reflective ports are connected
00189   nr_double_t tiny1 = (d == 0) ? 1.0 - TINYS : 1.0;
00190   nr_double_t tiny2 = tiny1 * tiny1;
00191   nr_double_t tiny3 = tiny1 * tiny2;
00192   d = 1.0 - s->getS (k, k) * t->getS (l, l) * tiny2;
00193 
00194   int j2; // column index for resulting matrix
00195   int i2; // row index for resulting matrix
00196   int j1; // column index for S matrix
00197   int i1; // row index for S matrix
00198 
00199   // handle S block
00200   i2 = j2 = 0;
00201   for (j1 = 0; j1 < s->getSize (); j1++) {
00202 
00203     // skip connected node
00204     if (j1 == k) continue;
00205 
00206     // assign node name of resulting circuit
00207     result->setNode (j2, s->getNode(j1)->getName());
00208 
00209     // inside S
00210     for (i1 = 0; i1 < s->getSize (); i1++) {
00211 
00212       // skip connected node
00213       if (i1 == k) continue;
00214 
00215       // compute S'ij
00216       p  = s->getS (i1, j1);
00217       p += s->getS (k, j1) * t->getS (l, l) * s->getS (i1, k) * tiny3 / d;
00218       result->setS (i2++, j2, p);
00219     }
00220 
00221     // across S and T
00222     for (i1 = 0; i1 < t->getSize (); i1++) {
00223 
00224       // skip connected node
00225       if (i1 == l) continue;
00226 
00227       // compute S'mj
00228       p = s->getS (k, j1) * t->getS (i1, l) * tiny2 / d;
00229       result->setS (i2++, j2, p);
00230     }
00231     // next column
00232     j2++; i2 = 0;
00233   }
00234 
00235   // handle T block
00236   for (j1 = 0; j1 < t->getSize (); j1++) {
00237 
00238     // skip connected node
00239     if (j1 == l) continue;
00240 
00241     // assign node name of resulting circuit
00242     result->setNode (j2, t->getNode(j1)->getName ());
00243 
00244     // across T and S
00245     for (i1 = 0; i1 < s->getSize (); i1++) {
00246 
00247       // skip connected node
00248       if (i1 == k) continue;
00249 
00250       // compute S'mj
00251       p = t->getS (l, j1) * s->getS (i1, k) * tiny2 / d;
00252       result->setS (i2++, j2, p);
00253     }
00254 
00255     // inside T
00256     for (i1 = 0; i1 < t->getSize (); i1++) {
00257 
00258       // skip connected node
00259       if (i1 == l) continue;
00260 
00261       // compute S'ij
00262       p  = t->getS (i1, j1);
00263       p += t->getS (l, j1) * s->getS (k, k) * t->getS (i1, l) * tiny3 / d;
00264       result->setS (i2++, j2, p);
00265     }
00266 
00267     // next column
00268     j2++; i2 = 0;
00269   }
00270 
00271   return result;
00272 }
00273 
00274 /* This function joins the two given nodes of a single circuit
00275    (interconnected nodes) and modifies the resulting circuit
00276    appropriately. */
00277 void spsolver::noiseInterconnect (circuit * result, node * n1, node * n2) {
00278 
00279   circuit * c = n1->getCircuit ();
00280   nr_complex_t p, k1, k2, k3, k4;
00281 
00282   // interconnected port numbers
00283   int k = n1->getPort (), l = n2->getPort ();
00284 
00285   // denominator needs to be calculated only once
00286   nr_complex_t t = (1.0 - c->getS (k, l)) * (1.0 - c->getS (l, k)) -
00287     c->getS (k, k) * c->getS (l, l);
00288 
00289   // avoid singularity when two full reflective ports are interconnected
00290   nr_double_t tiny1 = (t == 0) ? 1.0 - TINYS : 1.0;
00291   nr_double_t tiny2 = tiny1 * tiny1;
00292   t = (1.0 - c->getS (k, l) * tiny1) * (1.0 - c->getS (l, k) * tiny1) -
00293     c->getS (k, k) * c->getS (l, l) * tiny2;
00294 
00295   int j2; // column index for resulting matrix
00296   int i2; // row index for resulting matrix
00297   int j1; // column index for S matrix
00298   int i1; // row index for S matrix
00299 
00300   // handle single C block only
00301   i2 = j2 = 0;
00302   for (j1 = 0; j1 < c->getSize (); j1++) {
00303 
00304     // skip connected node
00305     if (j1 == k || j1 == l) continue;
00306 
00307     // inside C only
00308     for (i1 = 0; i1 < c->getSize (); i1++) {
00309 
00310       // skip connected node
00311       if (i1 == k || i1 == l) continue;
00312 
00313       k1 = (c->getS (i1, l) * (1.0 - c->getS (l, k)) +
00314             c->getS (l, l) * c->getS (i1, k)) * tiny2 / t;
00315       k2 = (c->getS (i1, k) * (1.0 - c->getS (k, l)) +
00316             c->getS (k, k) * c->getS (i1, l)) * tiny2 / t;
00317       k3 = (c->getS (j1, l) * (1.0 - c->getS (l, k)) +
00318             c->getS (l, l) * c->getS (j1, k)) * tiny2 / t;
00319       k4 = (c->getS (j1, k) * (1.0 - c->getS (k, l)) +
00320             c->getS (k, k) * c->getS (j1, l)) * tiny2 / t;
00321 
00322       p =
00323         c->getN (i1, j1) + c->getN (k, j1) * k1 + c->getN (l, j1) * k2 +
00324         conj (k3) * (c->getN (i1, k) + c->getN (k, k) * k1 +
00325                      c->getN (l, k) * k2) +
00326         conj (k4) * (c->getN (i1, l) + c->getN (k, l) * k1 +
00327                      c->getN (l, l) * k2);
00328       result->setN (i2, j2, p);
00329 
00330       if (i2 >= j2) break; // the other half need not be computed
00331       result->setN (j2, i2, conj (p));
00332       i2++;
00333     }
00334 
00335     // next column
00336     j2++; i2 = 0;
00337   }
00338 }
00339 
00340 
00341 /* The following function joins two nodes of two different circuits
00342    and saves the noise wave correlation matrix in the resulting
00343    circuit. */
00344 void spsolver::noiseConnect (circuit * result, node * n1, node * n2) {
00345   circuit * c = n1->getCircuit ();
00346   circuit * d = n2->getCircuit ();
00347   nr_complex_t p;
00348 
00349   // connected port numbers
00350   int k = n1->getPort (), l = n2->getPort ();
00351 
00352   // denominator needs to be calculated only once
00353   nr_complex_t t = 1.0 - c->getS (k, k) * d->getS (l, l);
00354 
00355   // avoid singularity when two full reflective ports are connected
00356   nr_double_t tiny1 = (t == 0) ? 1.0 - TINYS : 1.0;
00357   nr_double_t tiny2 = tiny1 * tiny1;
00358   nr_double_t tiny3 = tiny1 * tiny2;
00359   nr_double_t tiny4 = tiny1 * tiny3;
00360   t = 1.0 - c->getS (k, k) * d->getS (l, l) * tiny2;
00361 
00362   int j2; // column index for resulting matrix
00363   int i2; // row index for resulting matrix
00364   int j1; // column index for S matrix
00365   int i1; // row index for S matrix
00366 
00367   // handle C block
00368   i2 = j2 = 0;
00369   for (j1 = 0; j1 < c->getSize (); j1++) {
00370 
00371     // skip connected node
00372     if (j1 == k) continue;
00373 
00374     // inside C
00375     for (i1 = 0; i1 < c->getSize (); i1++) {
00376 
00377       // skip connected node
00378       if (i1 == k) continue;
00379 
00380       // compute C'ij
00381       p = c->getN (i1, j1) +
00382         c->getN (k, j1) * d->getS (l, l) * c->getS (i1, k) * tiny2 / t +
00383         c->getN (i1, k) * conj (d->getS (l, l) * c->getS (j1, k) * tiny2 / t) +
00384         (c->getN (k, k) * norm (d->getS (l, l)) + d->getN (l, l)) *
00385         c->getS (i1, k) * conj (c->getS (j1, k)) * tiny4 / norm (t);
00386 
00387       result->setN (i2, j2, p);
00388       if (i2 >= j2) break; // the other half need not be computed
00389       result->setN (j2, i2, conj (p));
00390       i2++;
00391     }
00392 
00393     /* The formulas "across C and D" are calculated elsewhere by the
00394        other half of the matrix (conjugate complex).  Therefore, they
00395        are missing here. */
00396 
00397     // next column
00398     j2++; i2 = 0;
00399   }
00400 
00401   // handle D block
00402   for (j1 = 0; j1 < d->getSize (); j1++) {
00403 
00404     // skip connected node
00405     if (j1 == l) continue;
00406 
00407     // across D and C
00408     for (i1 = 0; i1 < c->getSize (); i1++) {
00409 
00410       // skip connected node
00411       if (i1 == k) continue;
00412 
00413       // compute C'ij
00414       p = (c->getN (k, k) * d->getS (l, l) +
00415            d->getN (l, l) * conj (c->getS (k, k))) *
00416         c->getS (i1, k) * conj (d->getS (j1, l)) * tiny3 / norm (t) +
00417         d->getN (l, j1) * c->getS (i1, k) * tiny1 / t +
00418         c->getN (i1, k) * conj (d->getS (j1, l) * tiny1 / t);
00419       result->setN (i2, j2, p);
00420       result->setN (j2, i2, conj (p));
00421       i2++;
00422     }
00423 
00424     // inside D
00425     for (i1 = 0; i1 < d->getSize (); i1++) {
00426 
00427       // skip connected node
00428       if (i1 == l) continue;
00429 
00430       // compute C'ij
00431       p = d->getN (i1, j1) +
00432         (d->getN (l, l) * norm (c->getS (k, k)) + c->getN (k, k)) *
00433         d->getS (i1, l) * conj (d->getS (j1, l)) * tiny4 / norm (t) +
00434         d->getN (i1, l) * conj (c->getS (k, k) * d->getS (j1, l) * tiny2 / t) +
00435         d->getN (l, j1) * c->getS (k, k) * d->getS (i1, l) * tiny2 / t;
00436       result->setN (i2, j2, p);
00437       if (i2 >= j2) break; // the other half need not be computed
00438       result->setN (j2, i2, conj (p));
00439       i2++;
00440     }
00441 
00442     // next column
00443     j2++; i2 = 0;
00444   }
00445 }
00446 
00447 /* Goes through the list of circuit objects and runs its frequency
00448    dependent calcSP() function. */
00449 void spsolver::calc (nr_double_t freq) {
00450   circuit * root = subnet->getRoot ();
00451   for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) {
00452     c->calcSP (freq);
00453     if (noise) c->calcNoiseSP (freq);
00454   }
00455 }
00456 
00457 /* Go through each registered circuit object in the list and find the
00458    connection which results in a new subnetwork with the smallest
00459    number of s-parameters to calculate. */
00460 void spsolver::reduce (void) {
00461 
00462 #if SORTED_LIST
00463   node * n1, * n2;
00464   circuit * result, * cand1, * cand2;
00465 
00466   nlist->sortedNodes (&n1, &n2);
00467   cand1 = n1->getCircuit ();
00468   cand2 = n2->getCircuit ();
00469 #else /* !SORTED_LIST */
00470   node * n1, * n2, * cand;
00471   circuit * result, * c1, * c2, * cand1, * cand2;
00472   int ports;
00473   circuit * root = subnet->getRoot ();
00474 
00475   // initialize local variables
00476   result = c1 = c2 = cand1 = cand2 = NULL;
00477   n1 = n2 = cand = NULL;
00478   ports = 10000; // huge
00479 
00480   // go through the circuit list
00481   for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) {
00482 
00483     // skip signal ports
00484     if (c->getPort ()) continue;
00485 
00486     // and each node in the circuit
00487     for (int i = 0; i < c->getSize (); i++) {
00488 
00489       // find duplicate node
00490       if ((cand = subnet->findConnectedCircuitNode (c->getNode (i))) != NULL) {
00491 
00492         // save both candidates
00493         c1 = c; c2 = cand->getCircuit ();
00494         // connected
00495         if (c1 != c2) {
00496           if (c1->getSize () + c2->getSize () - 2 < ports) {
00497             ports = c1->getSize () + c2->getSize () - 2;
00498             cand1 = c1; cand2 = c2; n1 = c1->getNode (i); n2 = cand;
00499           }
00500         }
00501         // interconnect
00502         else {
00503           if (c1->getSize () - 2 < ports) {
00504             ports = c1->getSize () - 2;
00505             cand1 = c1; cand2 = c2; n1 = c1->getNode (i); n2 = cand;
00506           }
00507         }
00508       }
00509     }
00510   }
00511 #endif /* !SORTED_LIST */
00512 
00513   // found a connection ?
00514   if (cand1 != NULL && cand2 != NULL) {
00515     // connected
00516     if (cand1 != cand2) {
00517 #if DEBUG && 0
00518       logprint (LOG_STATUS, "DEBUG: connected node (%s): %s - %s\n",
00519                 n1->getName (), cand1->getName (), cand2->getName ());
00520 #endif /* DEBUG */
00521       result = connectedJoin (n1, n2);
00522       if (noise) noiseConnect (result, n1, n2);
00523       subnet->reducedCircuit (result);
00524 #if SORTED_LIST
00525       nlist->remove (cand1);
00526       nlist->remove (cand2);
00527       nlist->insert (result);
00528 #endif /* SORTED_LIST */
00529       subnet->removeCircuit (cand1);
00530       subnet->removeCircuit (cand2);
00531       subnet->insertCircuit (result);
00532       result->setOriginal (0);
00533     }
00534     // interconnect
00535     else {
00536 #if DEBUG && 0
00537       logprint (LOG_STATUS, "DEBUG: interconnected node (%s): %s\n",
00538                 n1->getName (), cand1->getName ());
00539 #endif
00540       result = interconnectJoin (n1, n2);
00541       if (noise) noiseInterconnect (result, n1, n2);
00542       subnet->reducedCircuit (result);
00543 #if SORTED_LIST
00544       nlist->remove (cand1);
00545       nlist->insert (result);
00546 #endif /* SORTED_LIST */
00547       subnet->removeCircuit (cand1);
00548       subnet->insertCircuit (result);
00549       result->setOriginal (0);
00550     }
00551   }
00552 }
00553 
00554 /* Goes through the list of circuit objects and runs initializing
00555    functions if necessary. */
00556 void spsolver::init (void) {
00557   circuit * root = subnet->getRoot ();
00558   for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) {
00559     if (c->isNonLinear ()) c->calcOperatingPoints ();
00560     c->initSP ();
00561     if (noise) c->initNoiseSP ();
00562   }
00563 }
00564 
00565 /* This is the netlist solver.  It prepares the circuit list for each
00566    requested frequency and solves it then. */
00567 int spsolver::solve (void) {
00568   nr_double_t freq;
00569   int ports;
00570   runs++;
00571 
00572   // fetch simulation properties
00573   saveCVs |= !strcmp (getPropertyString ("saveCVs"), "yes") ? SAVE_CVS : 0;
00574   saveCVs |= !strcmp (getPropertyString ("saveAll"), "yes") ? SAVE_ALL : 0;
00575 
00576   // run additional noise analysis ?
00577   noise = !strcmp (getPropertyString ("Noise"), "yes") ? 1 : 0;
00578 
00579   // create frequency sweep if necessary
00580   if (swp == NULL) {
00581     swp = createSweep ("frequency");
00582   }
00583 
00584   init ();
00585   insertConnections ();
00586 
00587 #if SORTED_LIST
00588 #if DEBUG
00589   logprint (LOG_STATUS, "NOTIFY: %s: creating sorted nodelist for "
00590             "SP analysis\n", getName ());
00591 #endif
00592   nlist = new nodelist (subnet);
00593   nlist->sort ();
00594 #endif /* SORTED_LIST */
00595 
00596 #if DEBUG
00597   logprint (LOG_STATUS, "NOTIFY: %s: solving SP netlist\n", getName ());
00598 #endif
00599 
00600   swp->reset ();
00601   for (int i = 0; i < swp->getSize (); i++) {
00602     freq = swp->next ();
00603     if (progress) logprogressbar (i, swp->getSize (), 40);
00604 
00605     ports = subnet->countNodes ();
00606     subnet->setReduced (0);
00607     calc (freq);
00608 
00609 #if DEBUG && 0
00610     logprint (LOG_STATUS, "NOTIFY: %s: solving netlist for f = %e\n",
00611               getName (), (double) freq);
00612 #endif
00613 
00614     while (ports > subnet->getPorts ()) {
00615       reduce ();
00616       ports -= 2;
00617     }
00618 
00619     saveResults (freq);
00620     subnet->getDroppedCircuits (nlist);
00621     subnet->deleteUnusedCircuits (nlist);
00622     if (saveCVs & SAVE_CVS) saveCharacteristics (freq);
00623   }
00624   if (progress) logprogressclear (40);
00625   dropConnections ();
00626 #if SORTED_LIST
00627   delete nlist; nlist = NULL;
00628 #endif
00629   return 0;
00630 }
00631 
00632 /* The function goes through the list of circuit objects and creates
00633    tee and cross circuits if necessary.  It looks for nodes in the
00634    circuit list connected to the given node. */
00635 void spsolver::insertConnectors (node * n) {
00636 
00637   int count = 0;
00638   node * nodes[4], * _node;
00639   const char * _name = n->getName ();
00640   circuit * root = subnet->getRoot ();
00641 
00642 #if USE_GROUNDS
00643   if (!strcmp (_name, "gnd")) return;
00644 #endif /* USE_GROUNDS */
00645 
00646   nodes[0] = n;
00647 
00648   // go through list of circuit objects
00649   for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) {
00650     // and each node in a circuit
00651     for (int i = 0; i < c->getSize (); i++) {
00652       _node = c->getNode (i);
00653       if (!strcmp (_node->getName (), _name)) {
00654         if (_node != n) {
00655 
00656           // found a connected node
00657           nodes[++count] = _node;
00658 #if USE_CROSSES
00659           if (count == 3) {
00660             // create an additional cross and assign its nodes
00661             insertCross (nodes, _name);
00662             count = 1;
00663           }
00664 #else /* !USE_CROSSES */
00665           if (count == 2) {
00666             // create an additional tee and assign its nodes
00667             insertTee (nodes, _name);
00668             count = 1;
00669           }
00670 #endif /* !USE_CROSSES */
00671         }
00672       }
00673     }
00674   }
00675 #if USE_CROSSES
00676   /* if using crosses there can be a tee left here */
00677   if (count == 2) {
00678     insertTee (nodes, _name);
00679   }
00680 #endif /* USE_CROSSES */
00681 }
00682 
00683 /* The following function creates a tee circuit with the given nodes
00684    and the node name.  The tee's node names are adjusted to be
00685    internal nodes. */
00686 void spsolver::insertTee (node ** nodes, const char * name) {
00687   circuit * result;
00688   // create a tee and assign its node names
00689   result = new tee ();
00690   subnet->insertedCircuit (result);
00691   result->setNode (0, name);
00692   subnet->insertedNode (result->getNode (1));
00693   subnet->insertedNode (result->getNode (2));
00694   // rename the nodes connected to the tee
00695   nodes[1]->setName (result->getNode(1)->getName ());
00696   nodes[2]->setName (result->getNode(2)->getName ());
00697   // complete the nodes of the tee
00698   result->getNode(1)->setCircuit (result);
00699   result->getNode(2)->setCircuit (result);
00700   result->getNode(1)->setPort (1);
00701   result->getNode(2)->setPort (2);
00702   // put the tee into the circuit list and initialize it
00703   subnet->insertCircuit (result);
00704   result->initSP (); if (noise) result->initNoiseSP ();
00705   // put the tee's first node into the node collection
00706   nodes[1] = result->getNode (0);
00707   tees++;
00708 }
00709 
00710 /* The following function creates a cross circuit with the given nodes
00711    and the node name.  The cross's node names are adjusted to be
00712    internal nodes. */
00713 void spsolver::insertCross (node ** nodes, const char * name) {
00714   circuit * result;
00715   // create a cross and assign its node names
00716   result = new cross ();
00717   subnet->insertedCircuit (result);
00718   result->setNode (0, name);
00719   subnet->insertedNode (result->getNode (1));
00720   subnet->insertedNode (result->getNode (2));
00721   subnet->insertedNode (result->getNode (3));
00722   // rename the nodes connected to the cross
00723   nodes[1]->setName (result->getNode(1)->getName ());
00724   nodes[2]->setName (result->getNode(2)->getName ());
00725   nodes[3]->setName (result->getNode(3)->getName ());
00726   // complete the nodes of the cross
00727   result->getNode(1)->setCircuit (result);
00728   result->getNode(2)->setCircuit (result);
00729   result->getNode(3)->setCircuit (result);
00730   result->getNode(1)->setPort (1);
00731   result->getNode(2)->setPort (2);
00732   result->getNode(3)->setPort (3);
00733   // put the cross into the circuit list and initialize it
00734   subnet->insertCircuit (result);
00735   result->initSP (); if (noise) result->initNoiseSP ();
00736   // put the cross's first node into the node collection
00737   nodes[1] = result->getNode (0);
00738   crosses++;
00739 }
00740 
00741 /* This function removes an inserted tee from the netlist and restores
00742    the original node names. */
00743 void spsolver::dropTee (circuit * c) {
00744   node * n;
00745   if (c->getType () == CIR_TEE) {
00746     const char * name = c->getNode(0)->getName ();
00747     n = subnet->findConnectedNode (c->getNode (1)); n->setName (name);
00748     n = subnet->findConnectedNode (c->getNode (2)); n->setName (name);
00749     c->setOriginal (0);
00750     subnet->removeCircuit (c);
00751   }
00752 }
00753 
00754 /* This function removes an inserted cross from the netlist and restores
00755    the original node names. */
00756 void spsolver::dropCross (circuit * c) {
00757   node * n;
00758   if (c->getType () == CIR_CROSS) {
00759     const char * name = c->getNode(0)->getName ();
00760     n = subnet->findConnectedNode (c->getNode (1)); n->setName (name);
00761     n = subnet->findConnectedNode (c->getNode (2)); n->setName (name);
00762     n = subnet->findConnectedNode (c->getNode (3)); n->setName (name);
00763     c->setOriginal (0);
00764     subnet->removeCircuit (c);
00765   }
00766 }
00767 
00768 /* The function adds an open to the circuit list if the given node is
00769    unconnected. */
00770 void spsolver::insertOpen (node * n) {
00771   if (strcmp (n->getName (), "gnd") &&
00772       subnet->findConnectedNode (n) == NULL) {
00773     circuit * result = new open ();
00774     subnet->insertedCircuit (result);
00775     result->setNode (0, n->getName ());
00776     subnet->insertCircuit (result);
00777     result->initSP (); if (noise) result->initNoiseSP ();
00778     opens++;
00779   }
00780 }
00781 
00782 // This function removes an inserted open from the netlist.
00783 void spsolver::dropOpen (circuit * c) {
00784   if (c->getType () == CIR_OPEN) {
00785     c->setOriginal (0);
00786     subnet->removeCircuit (c);
00787   }
00788 }
00789 
00790 /* The function adds a ground circuit to the circuit list if the given
00791    node is a ground connection. */
00792 void spsolver::insertGround (node * n) {
00793   if (!strcmp (n->getName (), "gnd") && !n->getCircuit()->getPort () &&
00794       n->getCircuit()->getType () != CIR_GROUND) {
00795     circuit * result = new ground ();
00796     subnet->insertedCircuit (result);
00797     subnet->insertedNode (result->getNode (0));
00798     result->getNode(0)->setCircuit (result);
00799     result->getNode(0)->setPort (0);
00800     n->setName (result->getNode(0)->getName ());
00801     subnet->insertCircuit (result);
00802     result->initSP (); if (noise) result->initNoiseSP ();
00803     grounds++;
00804   }
00805 }
00806 
00807 // This function removes an inserted ground from the netlist.
00808 void spsolver::dropGround (circuit * c) {
00809   if (c->getType () == CIR_GROUND) {
00810     node * n = subnet->findConnectedNode (c->getNode (0));
00811     n->setName ("gnd");
00812     c->setOriginal (0);
00813     subnet->removeCircuit (c);
00814   }
00815 }
00816 
00817 /* This function prepares the circuit list by adding Ts and opens to
00818    the circuit list.  With this adjustments the solver is able to
00819    solve the circuit. */
00820 void spsolver::insertConnections (void) {
00821 
00822   circuit * root, * c;
00823 #if DEBUG
00824   logprint (LOG_STATUS, "NOTIFY: %s: preparing circuit for analysis\n",
00825             getName ());
00826 #endif /* DEBUG */
00827 
00828 #if USE_GROUNDS
00829   // remove original ground circuit from netlist
00830   root = subnet->getRoot ();
00831   for (c = root; c != NULL; c = (circuit *) c->getNext ()) {
00832     if (c->getType () == CIR_GROUND) {
00833       gnd = c;
00834       subnet->removeCircuit (c, 0);
00835       break;
00836     }
00837   }
00838 #endif /* USE_GROUNDS */
00839 
00840   // insert opens, tee and crosses where necessary
00841   tees = crosses = opens = grounds = 0;
00842   root = subnet->getRoot ();
00843   for (c = root; c != NULL; c = (circuit *) c->getNext ()) {
00844     for (int i = 0; i < c->getSize (); i++) {
00845       insertConnectors (c->getNode (i));
00846       insertOpen (c->getNode (i));
00847     }
00848   }
00849 
00850   // insert S-parameter port transformers
00851   insertDifferentialPorts ();
00852 
00853 #if USE_GROUNDS
00854   // insert grounds where necessary
00855   root = subnet->getRoot ();
00856   for (c = root; c != NULL; c = (circuit *) c->getNext ()) {
00857     for (int i = 0; i < c->getSize (); i++) {
00858       insertGround (c->getNode (i));
00859     }
00860   }
00861 #endif /* USE_GROUNDS */
00862 
00863 #if DEBUG
00864   logprint (LOG_STATUS, "NOTIFY: %s: inserted %d tees, %d crosses, %d opens "
00865             "and %d grounds\n",
00866             getName (), tees, crosses, opens, grounds);
00867 #endif /* DEBUG */
00868 }
00869 
00870 /* The function is the counterpart of insertConnections().  It removes
00871    all additional circuits from the netlist which were necessary to
00872    run the analysis algorithm. */
00873 void spsolver::dropConnections (void) {
00874   circuit * next, * cand;
00875   int inserted;
00876 
00877   // drop all additional inserted circuits in correct order
00878   do {
00879     // find last inserted circuit
00880     inserted = -1;
00881     cand = NULL;
00882     for (circuit * c = subnet->getRoot (); c != NULL; c = next) {
00883       next = (circuit *) c->getNext ();
00884       if (c->getInserted () > inserted) {
00885         inserted = c->getInserted ();
00886         cand = c;
00887       }
00888     }
00889     // if found, then drop that circuit
00890     if (cand != NULL) {
00891       switch (cand->getType ()) {
00892       case CIR_OPEN:
00893         dropOpen (cand);
00894         break;
00895       case CIR_TEE:
00896         dropTee (cand);
00897         break;
00898       case CIR_CROSS:
00899         dropCross (cand);
00900         break;
00901       case CIR_GROUND:
00902         dropGround (cand);
00903         break;
00904       case CIR_ITRAFO:
00905         dropDifferentialPort (cand);
00906         break;
00907       }
00908     }
00909   } while (cand != NULL);
00910 
00911 #if USE_GROUNDS
00912   // attach the original ground to the netlist
00913   subnet->insertCircuit (gnd);
00914 #endif /* USE_GROUNDS */
00915 }
00916 
00917 /* This function inserts an ideal transformer before an AC power
00918    source in order to allow differential S parameter ports.  */
00919 void spsolver::insertDifferentialPorts (void) {
00920   circuit * root = subnet->getRoot ();
00921   for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) {
00922     if (c->getPort ()) {
00923 
00924       // create an ideal transformer and assign its node names
00925       circuit * result = new itrafo ();
00926       subnet->insertedCircuit (result);
00927       subnet->insertedNode (result->getNode (0));
00928       result->setNode (1, c->getNode(0)->getName ());
00929       result->setNode (2, c->getNode(1)->getName ());
00930 
00931       // rename the nodes connected to the trafo
00932       c->getNode(0)->setName (result->getNode(0)->getName ());
00933       c->getNode(1)->setName ("PacGround");
00934 
00935       // complete the nodes of the trafo
00936       result->getNode(0)->setCircuit (result);
00937       result->getNode(0)->setPort (0);
00938 
00939       // pass the port impedance to the ideal trafo
00940       result->addProperty ("Z", c->getPropertyDouble ("Z"));
00941 
00942       // put the trafo in the circuit list
00943       subnet->insertCircuit (result);
00944 
00945       // allocate S-parameter and noise correlation matrices
00946       result->initSP (); if (noise) result->initNoiseSP ();
00947     }
00948   }
00949 }
00950 
00951 /* This function removes an ideal transformer which was necessary to
00952    be placed in front of a s-parameter port in order to allow
00953    differential s-parameters.  It also restores the original node
00954    names. */
00955 void spsolver::dropDifferentialPort (circuit * c) {
00956   circuit * pac;
00957   node * n;
00958   if (c->getType () == CIR_ITRAFO) {
00959     n = subnet->findConnectedNode (c->getNode (0));
00960     pac = n->getCircuit ();
00961     pac->getNode(0)->setName (c->getNode(1)->getName ());
00962     pac->getNode(1)->setName (c->getNode(2)->getName ());
00963     c->setOriginal (0);
00964     subnet->removeCircuit (c);
00965   }
00966 }
00967 
00968 /* This function saves the results of a single solve() functionality
00969    (for the given frequency) into the output dataset. */
00970 void spsolver::saveResults (nr_double_t freq) {
00971 
00972   vector * f;
00973   node * sig_i, * sig_j;
00974   char * n;
00975   int res_i, res_j;
00976   circuit * root = subnet->getRoot ();
00977 
00978   // temporary noise matrices and input port impedance
00979   nr_complex_t noise_c[4], noise_s[4];
00980   nr_double_t z0 = circuit::z0;
00981 
00982   // add current frequency to the dependency of the output dataset
00983   if ((f = data->findDependency ("frequency")) == NULL) {
00984     f = new vector ("frequency");
00985     data->addDependency (f);
00986   }
00987   if (runs == 1) f->add (freq);
00988 
00989   // go through the list of remaining circuits
00990   for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) {
00991     // skip signals
00992     if (!c->getPort ()) {
00993       // handle each s-parameter
00994       for (int i = 0; i < c->getSize (); i++) {
00995         for (int j = 0; j < c->getSize (); j++) {
00996 
00997           // generate the appropriate variable name
00998           sig_i = subnet->findConnectedNode (c->getNode (i));
00999           sig_j = subnet->findConnectedNode (c->getNode (j));
01000           res_i = sig_i->getCircuit()->getPropertyInteger ("Num");
01001           res_j = sig_j->getCircuit()->getPropertyInteger ("Num");
01002           n = createSP (res_i, res_j);
01003 
01004           // add variable data item to dataset
01005           saveVariable (n, c->getS (i, j), f);
01006 
01007           // if noise analysis is requested
01008           if (noise) {
01009             int ro, co;
01010             int ni = getPropertyInteger ("NoiseIP");
01011             int no = getPropertyInteger ("NoiseOP");
01012             if ((res_i == ni || res_i == no) && (res_j == ni || res_j == no)) {
01013               if (ni == res_i) {
01014                 // assign input port impedance
01015                 z0 = sig_i->getCircuit()->getPropertyDouble ("Z");
01016               }
01017               ro = (res_i == ni) ? 0 : 1;
01018               co = (res_j == ni) ? 0 : 1;
01019               // save results in temporary data items
01020               noise_c[co + ro * 2] = c->getN (i, j);
01021               noise_s[co + ro * 2] = c->getS (i, j);
01022             }
01023           }
01024         }
01025       }
01026     }
01027   }
01028 
01029   // finally compute and save noise parameters
01030   if (noise) {
01031     saveNoiseResults (noise_s, noise_c, z0, f);
01032   }
01033 }
01034 
01035 /* This function takes the s-parameter matrix and noise wave
01036    correlation matrix and computes the noise parameters based upon
01037    these values.  Then it save the results into the dataset. */
01038 void spsolver::saveNoiseResults (nr_complex_t s[4], nr_complex_t c[4],
01039                                  nr_double_t z0, vector * f) {
01040   nr_complex_t c22 = c[3], c11 = c[0], c12 = c[1];
01041   nr_complex_t s11 = s[0], s21 = s[2];
01042   nr_complex_t n1, n2, F, Sopt, Fmin, Rn;
01043 
01044   // linear noise figure
01045   F    = real (1.0 + c22 / norm (s21));
01046   n1   =
01047     c11 * norm (s21) - 2.0 * real (c12 * s21 * conj (s11)) +
01048     c22 * norm (s11);
01049   n2   = 2.0 * (c22 * s11 - c12 * s21) / (c22 + n1);
01050 
01051   // optimal source reflection coefficient
01052   Sopt = 1.0 - norm (n2);
01053   if (real (Sopt) < 0.0)
01054     Sopt = (1.0 + sqrt (Sopt)) / n2;  // avoid a negative radicant
01055   else
01056     Sopt = (1.0 - sqrt (Sopt)) / n2;
01057 
01058   // minimum noise figure
01059   Fmin = real (1.0 + (c22 - n1 * norm (Sopt)) /
01060                norm (s21) / (1.0 + norm (Sopt)));
01061 
01062   // equivalent noise resistance
01063   Rn   = real ((c11 - 2.0 * real (c12 * conj ((1.0 + s11) / s21)) +
01064                 c22 * norm ((1.0 + s11) / s21)) / 4.0);
01065   Rn   = Rn * z0;
01066 
01067   // add variable data items to dataset
01068   saveVariable ("F", F, f);
01069   saveVariable ("Sopt", Sopt, f);
01070   saveVariable ("Fmin", Fmin, f);
01071   saveVariable ("Rn", Rn, f);
01072 }
01073 
01074 // Create an appropriate variable name.
01075 char * spsolver::createSP (int i, int j) {
01076   return matvec::createMatrixString ("S", i - 1, j - 1);
01077 }
01078 
01079 /* Create an appropriate variable name for characteristic values.  The
01080    caller is responsible to free() the returned string. */
01081 const char * spsolver::createCV (const std::string &c, const std::string &n) {
01082   return (c+"."+n).c_str();
01083 }
01084 
01085 /* Goes through the list of circuit objects and runs its
01086    saveCharacteristics() function.  Then puts these values into the
01087    dataset. */
01088 void spsolver::saveCharacteristics (nr_double_t freq) {
01089   circuit * root = subnet->getRoot ();
01090   const char * n;
01091   vector * f = data->findDependency ("frequency");
01092   for (circuit * c = root; c != NULL; c = (circuit *) c->getNext ()) {
01093     c->saveCharacteristics (freq);
01094     if (!c->getSubcircuit ().empty() && !(saveCVs & SAVE_ALL)) continue;
01095     c->calcCharacteristics (freq);
01096     for (auto ps: c->getCharacteristics ()) {
01097       characteristic &p = ps.second;
01098       n = createCV (c->getName (), p.getName ());
01099       saveVariable (n, p.getValue (), f);
01100     }
01101   }
01102 }
01103 
01104 // properties
01105 PROP_REQ [] = {
01106   { "Type", PROP_STR, { PROP_NO_VAL, "lin" }, PROP_RNG_TYP },
01107   PROP_NO_PROP };
01108 PROP_OPT [] = {
01109   { "Noise", PROP_STR, { PROP_NO_VAL, "no" }, PROP_RNG_YESNO },
01110   { "NoiseIP", PROP_INT, { 1, PROP_NO_STR }, PROP_RNGII (1, MAX_PORTS) },
01111   { "NoiseOP", PROP_INT, { 2, PROP_NO_STR }, PROP_RNGII (1, MAX_PORTS) },
01112   { "Start", PROP_REAL, { 1e9, PROP_NO_STR }, PROP_POS_RANGE },
01113   { "Stop", PROP_REAL, { 10e9, PROP_NO_STR }, PROP_POS_RANGE },
01114   { "Points", PROP_INT, { 10, PROP_NO_STR }, PROP_MIN_VAL (2) },
01115   { "Values", PROP_LIST, { 10, PROP_NO_STR }, PROP_POS_RANGE },
01116   { "saveCVs", PROP_STR, { PROP_NO_VAL, "no" }, PROP_RNG_YESNO },
01117   { "saveAll", PROP_STR, { PROP_NO_VAL, "no" }, PROP_RNG_YESNO },
01118   PROP_NO_PROP };
01119 struct define_t spsolver::anadef =
01120   { "SP", 0, PROP_ACTION, PROP_NO_SUBSTRATE, PROP_LINEAR, PROP_DEF };
01121 
01122 } // namespace qucs