Qucs-core  0.0.19
matrix.cpp
Go to the documentation of this file.
00001 /*
00002  * matrix.cpp - matrix class implementation
00003  *
00004  * Copyright (C) 2003-2009 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  */
00082 #if HAVE_CONFIG_H
00083 # include <config.h>
00084 #endif
00085 
00086 #include <assert.h>
00087 #include <stdio.h>
00088 #include <cstdlib>
00089 #include <string.h>
00090 #include <cmath>
00091 
00092 #include "logging.h"
00093 #include "object.h"
00094 #include "complex.h"
00095 #include "vector.h"
00096 #include "matrix.h"
00097 
00098 namespace qucs {
00099 
00104 matrix::matrix () {
00105   rows = 0;
00106   cols = 0;
00107   data = NULL;
00108 }
00109 
00117 matrix::matrix (int s)  {
00118   rows = cols = s;
00119   data = (s > 0) ? new nr_complex_t[s * s] : NULL;
00120 }
00121 
00122 /* \brief Creates a matrix
00123 
00124    Constructor creates an unnamed instance of the matrix class with a
00125    certain number of rows and columns.
00126    \param[in] r number of rows
00127    \param[in] c number of column
00128    \todo Why not r and c constant
00129    \todo Assert r >= 0 and c >= 0
00130 */
00131 matrix::matrix (int r, int c)  {
00132   rows = r;
00133   cols = c;
00134   data = (r > 0 && c > 0) ? new nr_complex_t[r * c] : NULL;
00135 }
00136 
00137 /* \brief copy constructor
00138 
00139    The copy constructor creates a new instance based on the given
00140    matrix object.
00141    \todo Add assert tests
00142 */
00143 matrix::matrix (const matrix & m) {
00144   rows = m.rows;
00145   cols = m.cols;
00146   data = NULL;
00147 
00148   // copy matrix elements
00149   if (rows > 0 && cols > 0) {
00150     data = new nr_complex_t[rows * cols];
00151     memcpy (data, m.data, sizeof (nr_complex_t) * rows * cols);
00152   }
00153 }
00154 
00164 const matrix& matrix::operator=(const matrix & m) {
00165   if (&m != this) {
00166     rows = m.rows;
00167     cols = m.cols;
00168     if (data) {
00169       delete[] data;
00170       data = NULL;
00171     }
00172     if (rows > 0 && cols > 0) {
00173       data = new nr_complex_t[rows * cols];
00174       memcpy (data, m.data, sizeof (nr_complex_t) * rows * cols);
00175     }
00176   }
00177   return *this;
00178 }
00179 
00184 matrix::~matrix () {
00185   delete[] data;
00186 }
00187 
00194 nr_complex_t matrix::get (int r, int c) {
00195   return data[r * cols + c];
00196 }
00197 
00205 void matrix::set (int r, int c, nr_complex_t z) {
00206   data[r * cols + c] = z;
00207 }
00208 
00209 #ifdef DEBUG
00210 
00211 void matrix::print (void) {
00212   for (int r = 0; r < rows; r++) {
00213     for (int c = 0; c < cols; c++) {
00214       fprintf (stderr, "%+.2e,%+.2e ", (double) real (get (r, c)),
00215                (double) imag (get (r, c)));
00216     }
00217     fprintf (stderr, "\n");
00218   }
00219 }
00220 #endif /* DEBUG */
00221 
00228 matrix operator + (matrix a, matrix b) {
00229   assert (a.getRows () == b.getRows () && a.getCols () == b.getCols ());
00230 
00231   matrix res (a.getRows (), a.getCols ());
00232   for (int r = 0; r < a.getRows (); r++)
00233     for (int c = 0; c < a.getCols (); c++)
00234       res.set (r, c, a.get (r, c) + b.get (r, c));
00235   return res;
00236 }
00237 
00243 matrix matrix::operator += (matrix a) {
00244   assert (a.getRows () == rows && a.getCols () == cols);
00245 
00246   int r, c, i;
00247   for (i = 0, r = 0; r < a.getRows (); r++)
00248     for (c = 0; c < a.getCols (); c++, i++)
00249       data[i] += a.get (r, c);
00250   return *this;
00251 }
00252 
00259 matrix operator - (matrix a, matrix b) {
00260   assert (a.getRows () == b.getRows () && a.getCols () == b.getCols ());
00261 
00262   matrix res (a.getRows (), a.getCols ());
00263   for (int r = 0; r < a.getRows (); r++)
00264     for (int c = 0; c < a.getCols (); c++)
00265       res.set (r, c, a.get (r, c) - b.get (r, c));
00266   return res;
00267 }
00268 
00270 matrix matrix::operator - () {
00271   matrix res (getRows (), getCols ());
00272   int r, c, i;
00273   for (i = 0, r = 0; r < getRows (); r++)
00274     for (c = 0; c < getCols (); c++, i++)
00275       res.set (r, c, -data[i]);
00276   return res;
00277 }
00278 
00283 matrix matrix::operator -= (matrix a) {
00284   assert (a.getRows () == rows && a.getCols () == cols);
00285   int r, c, i;
00286   for (i = 0, r = 0; r < a.getRows (); r++)
00287     for (c = 0; c < a.getCols (); c++, i++)
00288       data[i] -= a.get (r, c);
00289   return *this;
00290 }
00291 
00298 matrix operator * (matrix a, nr_complex_t z) {
00299   matrix res (a.getRows (), a.getCols ());
00300   for (int r = 0; r < a.getRows (); r++)
00301     for (int c = 0; c < a.getCols (); c++)
00302       res.set (r, c, a.get (r, c) * z);
00303   return res;
00304 }
00305 
00313 matrix operator * (nr_complex_t z, matrix a) {
00314   return a * z;
00315 }
00316 
00323 matrix operator * (matrix a, nr_double_t d) {
00324   matrix res (a.getRows (), a.getCols ());
00325   for (int r = 0; r < a.getRows (); r++)
00326     for (int c = 0; c < a.getCols (); c++)
00327       res.set (r, c, a.get (r, c) * d);
00328   return res;
00329 }
00330 
00338 matrix operator * (nr_double_t d, matrix a) {
00339   return a * d;
00340 }
00341 
00348 matrix operator / (matrix a, nr_complex_t z) {
00349   matrix res (a.getRows (), a.getCols ());
00350   for (int r = 0; r < a.getRows (); r++)
00351     for (int c = 0; c < a.getCols (); c++)
00352       res.set (r, c, a.get (r, c) / z);
00353   return res;
00354 }
00355 
00362 matrix operator / (matrix a, nr_double_t d) {
00363   matrix res (a.getRows (), a.getCols ());
00364   for (int r = 0; r < a.getRows (); r++)
00365     for (int c = 0; c < a.getCols (); c++)
00366       res.set (r, c, a.get (r, c) / d);
00367   return res;
00368 }
00369 
00378 matrix operator * (matrix a, matrix b) {
00379   assert (a.getCols () == b.getRows ());
00380 
00381   int r, c, i, n = a.getCols ();
00382   nr_complex_t z;
00383   matrix res (a.getRows (), b.getCols ());
00384   for (r = 0; r < a.getRows (); r++) {
00385     for (c = 0; c < b.getCols (); c++) {
00386       for (i = 0, z = 0; i < n; i++) z += a.get (r, i) * b.get (i, c);
00387       res.set (r, c, z);
00388     }
00389   }
00390   return res;
00391 }
00392 
00399 matrix operator + (matrix a, nr_complex_t z) {
00400   matrix res (a.getRows (), a.getCols ());
00401   for (int r = 0; r < a.getRows (); r++)
00402     for (int c = 0; c < a.getCols (); c++)
00403       res.set (r, c, a.get (r, c) + z);
00404   return res;
00405 }
00406 
00414 matrix operator + (nr_complex_t z, matrix a) {
00415   return a + z;
00416 }
00417 
00424 matrix operator + (matrix a, nr_double_t d) {
00425   matrix res (a.getRows (), a.getCols ());
00426   for (int r = 0; r < a.getRows (); r++)
00427     for (int c = 0; c < a.getCols (); c++)
00428       res.set (r, c, a.get (r, c) + d);
00429   return res;
00430 }
00431 
00439 matrix operator + (nr_double_t d, matrix a) {
00440   return a + d;
00441 }
00442 
00450 matrix operator - (matrix a, nr_complex_t z) {
00451   return -z + a;
00452 }
00453 
00461 matrix operator - (nr_complex_t z, matrix a) {
00462   return -a + z;
00463 }
00464 
00472 matrix operator - (matrix a, nr_double_t d) {
00473   return -d + a;
00474 }
00475 
00483 matrix operator - (nr_double_t d, matrix a) {
00484   return -a + d;
00485 }
00486 
00492 matrix transpose (matrix a) {
00493   matrix res (a.getCols (), a.getRows ());
00494   for (int r = 0; r < a.getRows (); r++)
00495     for (int c = 0; c < a.getCols (); c++)
00496       res.set (c, r, a.get (r, c));
00497   return res;
00498 }
00499 
00505 matrix conj (matrix a) {
00506   matrix res (a.getRows (), a.getCols ());
00507   for (int r = 0; r < a.getRows (); r++)
00508     for (int c = 0; c < a.getCols (); c++)
00509       res.set (r, c, conj (a.get (r, c)));
00510   return res;
00511 }
00512 
00522 matrix adjoint (matrix a) {
00523   return transpose (conj (a));
00524 }
00525 
00531 matrix abs (matrix a) {
00532   matrix res (a.getRows (), a.getCols ());
00533   for (int r = 0; r < a.getRows (); r++)
00534     for (int c = 0; c < a.getCols (); c++)
00535       res.set (r, c, abs (a.get (r, c)));
00536   return res;
00537 }
00538 
00542 matrix dB (matrix a) {
00543   matrix res (a.getRows (), a.getCols ());
00544   for (int r = 0; r < a.getRows (); r++)
00545     for (int c = 0; c < a.getCols (); c++)
00546       res.set (r, c, dB (a.get (r, c)));
00547   return res;
00548 }
00549 
00555 matrix arg (matrix a) {
00556   matrix res (a.getRows (), a.getCols ());
00557   for (int r = 0; r < a.getRows (); r++)
00558     for (int c = 0; c < a.getCols (); c++)
00559       res.set (r, c, arg (a.get (r, c)));
00560   return res;
00561 }
00562 
00568 matrix real (matrix a) {
00569   matrix res (a.getRows (), a.getCols ());
00570   for (int r = 0; r < a.getRows (); r++)
00571     for (int c = 0; c < a.getCols (); c++)
00572       res.set (r, c, real (a.get (r, c)));
00573   return res;
00574 }
00575 
00581 matrix imag (matrix a) {
00582   matrix res (a.getRows (), a.getCols ());
00583   for (int r = 0; r < a.getRows (); r++)
00584     for (int c = 0; c < a.getCols (); c++)
00585       res.set (r, c, imag (a.get (r, c)));
00586   return res;
00587 }
00588 
00592 matrix sqr (matrix a) {
00593   return a * a;
00594 }
00595 
00603 matrix eye (int rs, int cs) {
00604   matrix res (rs, cs);
00605   for (int r = 0; r < res.getRows (); r++)
00606     for (int c = 0; c < res.getCols (); c++)
00607       if (r == c) res.set (r, c, 1);
00608   return res;
00609 }
00610 
00616 matrix eye (int s) {
00617   return eye (s, s);
00618 }
00619 
00624 matrix diagonal (qucs::vector diag) {
00625   int size = diag.getSize ();
00626   matrix res (size);
00627   for (int i = 0; i < size; i++) res (i, i) = diag (i);
00628   return res;
00629 }
00630 
00631 // Compute n-th power of the given matrix.
00632 matrix pow (matrix a, int n) {
00633   matrix res;
00634   if (n == 0) {
00635     res = eye (a.getRows (), a.getCols ());
00636   }
00637   else {
00638     res = a = n < 0 ? inverse (a) : a;
00639     for (int i = 1; i < std::abs (n); i++)
00640       res = res * a;
00641   }
00642   return res;
00643 }
00644 
00657 nr_complex_t cofactor (matrix a, int u, int v) {
00658   matrix res (a.getRows () - 1, a.getCols () - 1);
00659   int r, c, ra, ca;
00660   for (ra = r = 0; r < res.getRows (); r++, ra++) {
00661     if (ra == u) ra++;
00662     for (ca = c = 0; c < res.getCols (); c++, ca++) {
00663       if (ca == v) ca++;
00664       res.set (r, c, a.get (ra, ca));
00665     }
00666   }
00667   nr_complex_t z = detLaplace (res);
00668   return ((u + v) & 1) ? -z : z;
00669 }
00670 
00686 nr_complex_t detLaplace (matrix a) {
00687   assert (a.getRows () == a.getCols ());
00688   int s = a.getRows ();
00689   nr_complex_t res = 0;
00690   if (s > 1) {
00691     /* always use the first row for sub-determinant, but you should
00692        use the row or column with most zeros in it */
00693     int r = 0;
00694     for (int i = 0; i < s; i++) {
00695       res += a.get (r, i) * cofactor (a, r, i);
00696     }
00697     return res;
00698   }
00699   /* 1 by 1 matrix */
00700   else if (s == 1) {
00701     return a (0, 0);
00702   }
00703   /* 0 by 0 matrix */
00704   return 1;
00705 }
00706 
00717 nr_complex_t detGauss (matrix a) {
00718   assert (a.getRows () == a.getCols ());
00719   nr_double_t MaxPivot;
00720   nr_complex_t f, res;
00721   matrix b;
00722   int i, c, r, pivot, n = a.getCols ();
00723 
00724   // return special matrix cases
00725   if (n == 0) return 1;
00726   if (n == 1) return a (0, 0);
00727 
00728   // make copy of original matrix
00729   b = matrix (a);
00730 
00731   // triangulate the matrix
00732   for (res = 1, i = 0; i < n; i++) {
00733     // find maximum column value for pivoting
00734     for (MaxPivot = 0, pivot = r = i; r < n; r++) {
00735       if (abs (b.get (r, i)) > MaxPivot) {
00736         MaxPivot = abs (b.get (r, i));
00737         pivot = r;
00738       }
00739     }
00740     // exchange rows if necessary
00741     assert (MaxPivot != 0);
00742     if (i != pivot) { b.exchangeRows (i, pivot); res = -res; }
00743     // compute new rows and columns
00744     for (r = i + 1; r < n; r++) {
00745       f = b.get (r, i) / b.get (i, i);
00746       for (c = i + 1; c < n; c++) {
00747         b.set (r, c, b.get (r, c) - f * b.get (i, c));
00748       }
00749     }
00750   }
00751 
00752   // now compute determinant by multiplying diagonal
00753   for (i = 0; i < n; i++) res *= b.get (i, i);
00754   return res;
00755 }
00756 
00762 nr_complex_t det (matrix a) {
00763 #if 0
00764   return detLaplace (a);
00765 #else
00766   return detGauss (a);
00767 #endif
00768 }
00769 
00779 matrix inverseLaplace (matrix a) {
00780   matrix res (a.getRows (), a.getCols ());
00781   nr_complex_t d = detLaplace (a);
00782   assert (abs (d) != 0); // singular matrix
00783   for (int r = 0; r < a.getRows (); r++)
00784     for (int c = 0; c < a.getCols (); c++)
00785       res.set (r, c, cofactor (a, c, r) / d);
00786   return res;
00787 }
00788 
00798 matrix inverseGaussJordan (matrix a) {
00799   nr_double_t MaxPivot;
00800   nr_complex_t f;
00801   matrix b, e;
00802   int i, c, r, pivot, n = a.getCols ();
00803 
00804   // create temporary matrix and the result matrix
00805   b = matrix (a);
00806   e = eye (n);
00807 
00808   // create the eye matrix in 'b' and the result in 'e'
00809   for (i = 0; i < n; i++) {
00810     // find maximum column value for pivoting
00811     for (MaxPivot = 0, pivot = r = i; r < n; r++) {
00812       if (abs (b (r, i)) > MaxPivot) {
00813         MaxPivot = abs (b (r, i));
00814         pivot = r;
00815       }
00816     }
00817     // exchange rows if necessary
00818     assert (MaxPivot != 0); // singular matrix
00819     if (i != pivot) {
00820       b.exchangeRows (i, pivot);
00821       e.exchangeRows (i, pivot);
00822     }
00823 
00824     // compute current row
00825     for (f = b (i, i), c = 0; c < n; c++) {
00826       b (i, c) /= f;
00827       e (i, c) /= f;
00828     }
00829 
00830     // compute new rows and columns
00831     for (r = 0; r < n; r++) {
00832       if (r != i) {
00833         for (f = b (r, i), c = 0; c < n; c++) {
00834           b (r, c) -= f * b (i, c);
00835           e (r, c) -= f * e (i, c);
00836         }
00837       }
00838     }
00839   }
00840   return e;
00841 }
00842 
00847 matrix inverse (matrix a) {
00848 #if 0
00849   return inverseLaplace (a);
00850 #else
00851   return inverseGaussJordan (a);
00852 #endif
00853 }
00854 
00890 matrix stos (matrix s, qucs::vector zref, qucs::vector z0) {
00891   int d = s.getRows ();
00892   matrix e, r, a;
00893 
00894   assert (d == s.getCols () && d == z0.getSize () && d == zref.getSize ());
00895 
00896   e = eye (d);
00897   r = diagonal ((z0 - zref) / (z0 + zref));
00898   a = diagonal (sqrt (z0 / zref) / (z0 + zref));
00899   return inverse (a) * (s - r) * inverse (e - r * s) * a;
00900 }
00901 
00910 matrix stos (matrix s, nr_complex_t zref, nr_complex_t z0) {
00911   int d = s.getRows ();
00912   return stos (s, qucs::vector (d, zref), qucs::vector (d, z0));
00913 }
00914 
00923 matrix stos (matrix s, nr_double_t zref, nr_double_t z0) {
00924   return stos (s, nr_complex_t (zref, 0), nr_complex_t (z0, 0));
00925 }
00926 
00935 matrix stos (matrix s, qucs::vector zref, nr_complex_t z0) {
00936   return stos (s, zref, qucs::vector (zref.getSize (), z0));
00937 }
00938 
00947 matrix stos (matrix s, nr_complex_t zref, qucs::vector z0) {
00948   return stos (s, qucs::vector (z0.getSize (), zref), z0);
00949 }
00950 
00973 matrix stoz (matrix s, qucs::vector z0) {
00974   int d = s.getRows ();
00975   matrix e, zref, gref;
00976 
00977   assert (d == s.getCols () && d == z0.getSize ());
00978 
00979   e = eye (d);
00980   zref = diagonal (z0);
00981   gref = diagonal (sqrt (real (1 / z0)));
00982   return inverse (gref) * inverse (e - s) * (s * zref + zref) * gref;
00983 }
00984 
00992 matrix stoz (matrix s, nr_complex_t z0) {
00993   return stoz (s, qucs::vector (s.getRows (), z0));
00994 }
00995 
01018 matrix ztos (matrix z, qucs::vector z0) {
01019   int d = z.getRows ();
01020   matrix e, zref, gref;
01021 
01022   assert (d == z.getCols () && d == z0.getSize ());
01023 
01024   e = eye (d);
01025   zref = diagonal (z0);
01026   gref = diagonal (sqrt (real (1 / z0)));
01027   return gref * (z - zref) * inverse (z + zref) * inverse (gref);
01028 }
01029 
01037 matrix ztos (matrix z, nr_complex_t z0) {
01038   return ztos (z, qucs::vector (z.getRows (), z0));
01039 }
01040 
01050 matrix ztoy (matrix z) {
01051   assert (z.getRows () == z.getCols ());
01052   return inverse (z);
01053 }
01054 
01082 matrix stoy (matrix s, qucs::vector z0) {
01083   int d = s.getRows ();
01084   matrix e, zref, gref;
01085 
01086   assert (d == s.getCols () && d == z0.getSize ());
01087 
01088   e = eye (d);
01089   zref = diagonal (z0);
01090   gref = diagonal (sqrt (real (1 / z0)));
01091   return inverse (gref) * inverse (s * zref + zref) * (e - s) * gref;
01092 }
01093 
01101 matrix stoy (matrix s, nr_complex_t z0) {
01102   return stoy (s, qucs::vector (s.getRows (), z0));
01103 }
01104 
01133 matrix ytos (matrix y, qucs::vector z0) {
01134   int d = y.getRows ();
01135   matrix e, zref, gref;
01136 
01137   assert (d == y.getCols () && d == z0.getSize ());
01138 
01139   e = eye (d);
01140   zref = diagonal (z0);
01141   gref = diagonal (sqrt (real (1 / z0)));
01142   return gref * (e - zref * y) * inverse (e + zref * y) * inverse (gref);
01143 }
01151 matrix ytos (matrix y, nr_complex_t z0) {
01152   return ytos (y, qucs::vector (y.getRows (), z0));
01153 }
01181 matrix stoa (matrix s, nr_complex_t z1, nr_complex_t z2) {
01182   nr_complex_t d = s (0, 0) * s (1, 1) - s (0, 1) * s (1, 0);
01183   nr_complex_t n = 2.0 * s (1, 0) * sqrt (fabs (real (z1) * real (z2)));
01184   matrix a (2);
01185 
01186   assert (s.getRows () >= 2 && s.getCols () >= 2);
01187 
01188   a.set (0, 0, (conj (z1) + z1 * s (0, 0) -
01189                 conj (z1) * s (1, 1) - z1 * d) / n);
01190   a.set (0, 1, (conj (z1) * conj (z2) + z1 * conj (z2) * s (0, 0) +
01191                 conj (z1) * z2 * s (1, 1) + z1 * z2 * d) / n);
01192   a.set (1, 0, (1.0 - s (0, 0) - s (1, 1) + d) / n);
01193   a.set (1, 1, (conj (z2) - conj (z2) * s (0, 0) +
01194                 z2 * s (1, 1) - z2 * d) / n);
01195   return a;
01196 }
01197 
01198 
01222 matrix atos (matrix a, nr_complex_t z1, nr_complex_t z2) {
01223   nr_complex_t d = 2.0 * sqrt (fabs (real (z1) * real (z2)));
01224   nr_complex_t n = a (0, 0) * z2 + a (0, 1) +
01225     a (1, 0) * z1 * z2 + a (1, 1) * z1;
01226   matrix s (2);
01227 
01228   assert (a.getRows () >= 2 && a.getCols () >= 2);
01229 
01230   s.set (0, 0, (a (0, 0) * z2 + a (0, 1)
01231                 - a (1, 0) * conj (z1) * z2 - a (1, 1) * conj (z1)) / n);
01232   s.set (0, 1, (a (0, 0) * a (1, 1) -
01233                 a (0, 1) * a (1, 0)) * d / n);
01234   s.set (1, 0, d / n);
01235   s.set (1, 1, (a (1, 1) * z1 - a (0, 0) * conj (z2) +
01236                 a (0, 1) - a (1, 0) * z1 * conj (z2)) / n);
01237   return s;
01238 }
01239 
01269 matrix stoh (matrix s, nr_complex_t z1, nr_complex_t z2) {
01270   nr_complex_t n = s (0, 1) * s (1, 0);
01271   nr_complex_t d = (1.0 - s (0, 0)) * (1.0 + s (1, 1)) + n;
01272   matrix h (2);
01273 
01274   assert (s.getRows () >= 2 && s.getCols () >= 2);
01275 
01276   h.set (0, 0, ((1.0 + s (0, 0)) * (1.0 + s (1, 1)) - n) * z1 / d);
01277   h.set (0, 1, +2.0 * s (0, 1) / d);
01278   h.set (1, 0, -2.0 * s (1, 0) / d);
01279   h.set (1, 1, ((1.0 - s (0, 0)) * (1.0 - s (1, 1)) - n) / z2 / d);
01280   return h;
01281 }
01282 
01307 matrix htos (matrix h, nr_complex_t z1, nr_complex_t z2) {
01308   nr_complex_t n = h (0, 1) * h (1, 0);
01309   nr_complex_t d = (1.0 + h (0, 0) / z1) * (1.0 + z2 * h (1, 1)) - n;
01310   matrix s (2);
01311 
01312   assert (h.getRows () >= 2 && h.getCols () >= 2);
01313 
01314   s.set (0, 0, ((h (0, 0) / z1 - 1.0) * (1.0 + z2 * h (1, 1)) - n) / d);
01315   s.set (0, 1, +2.0 * h (0, 1) / d);
01316   s.set (1, 0, -2.0 * h (1, 0) / d);
01317   s.set (1, 1, ((1.0 + h (0, 0) / z1) * (1.0 - z2 * h (1, 1)) + n) / d);
01318   return s;
01319 }
01320 
01321 /*\brief Converts scattering parameters to second hybrid matrix.
01322   \bug{Programmed formulae are valid only for Z real}
01323   \bug{Not documented and references}
01324   \param[in] s Scattering matrix
01325   \param[in] z1 impedance at input 1
01326   \param[in] z2 impedance at input 2
01327   \return second hybrid matrix
01328   \note Assert 2 by 2 matrix
01329   \todo Why not s,z1,z2 const
01330 */
01331 matrix stog (matrix s, nr_complex_t z1, nr_complex_t z2) {
01332   nr_complex_t n = s (0, 1) * s (1, 0);
01333   nr_complex_t d = (1.0 + s (0, 0)) * (1.0 - s (1, 1)) + n;
01334   matrix g (2);
01335 
01336   assert (s.getRows () >= 2 && s.getCols () >= 2);
01337 
01338   g.set (0, 0, ((1.0 - s (0, 0)) * (1.0 - s (1, 1)) - n) / z1 / d);
01339   g.set (0, 1, -2.0 * s (0, 1) / d);
01340   g.set (1, 0, +2.0 * s (1, 0) / d);
01341   g.set (1, 1, ((1.0 + s (0, 0)) * (1.0 + s (1, 1)) - n) * z2 / d);
01342   return g;
01343 }
01344 
01345 /*\brief Converts second hybrid matrix to scattering parameters.
01346   \bug{Programmed formulae are valid only for Z real}
01347   \bug{Not documented and references}
01348   \param[in] g second hybrid matrix
01349   \param[in] z1 impedance at input 1
01350   \param[in] z2 impedance at input 2
01351   \return scattering matrix
01352   \note Assert 2 by 2 matrix
01353   \todo Why not g,z1,z2 const
01354 */
01355 matrix gtos (matrix g, nr_complex_t z1, nr_complex_t z2) {
01356   nr_complex_t n = g (0, 1) * g (1, 0);
01357   nr_complex_t d = (1.0 + g (0, 0) * z1) * (1.0 + g (1, 1) / z2) - n;
01358   matrix s (2);
01359 
01360   assert (g.getRows () >= 2 && g.getCols () >= 2);
01361 
01362   s.set (0, 0, ((1.0 - g (0, 0) * z1) * (1.0 + g (1, 1) / z2) + n) / d);
01363   s.set (0, 1, -2.0 * g (0, 1) / d);
01364   s.set (1, 0, +2.0 * g (1, 0) / d);
01365   s.set (1, 1, ((g (0, 0) * z1 + 1.0) * (g (1, 1) / z2 - 1.0) - n) / d);
01366   return s;
01367 }
01368 
01380 matrix ytoz (matrix y) {
01381   assert (y.getRows () == y.getCols ());
01382   return inverse (y);
01383 }
01384 
01404 matrix cytocs (matrix cy, matrix s) {
01405   matrix e = eye (s.getRows ());
01406 
01407   assert (cy.getRows () == cy.getCols () && s.getRows () == s.getCols () &&
01408           cy.getRows () == s.getRows ());
01409 
01410   return (e + s) * cy * adjoint (e + s) / 4;
01411 }
01412 
01430 matrix cstocy (matrix cs, matrix y) {
01431   matrix e = eye (y.getRows ());
01432 
01433   assert (cs.getRows () == cs.getCols () && y.getRows () == y.getCols () &&
01434           cs.getRows () == y.getRows ());
01435 
01436   return (e + y) * cs * adjoint (e + y);
01437 }
01438 
01457 matrix cztocs (matrix cz, matrix s) {
01458   matrix e = eye (s.getRows ());
01459 
01460   assert (cz.getRows () == cz.getCols () && s.getRows () == s.getCols () &&
01461           cz.getRows () == s.getRows ());
01462 
01463   return (e - s) * cz * adjoint (e - s) / 4;
01464 }
01465 
01483 matrix cstocz (matrix cs, matrix z) {
01484   assert (cs.getRows () == cs.getCols () && z.getRows () == z.getCols () &&
01485           cs.getRows () == z.getRows ());
01486   matrix e = eye (z.getRows ());
01487   return (e + z) * cs * adjoint (e + z);
01488 }
01489 
01507 matrix cztocy (matrix cz, matrix y) {
01508   assert (cz.getRows () == cz.getCols () && y.getRows () == y.getCols () &&
01509           cz.getRows () == y.getRows ());
01510 
01511   return y * cz * adjoint (y);
01512 }
01513 
01531 matrix cytocz (matrix cy, matrix z) {
01532   assert (cy.getRows () == cy.getCols () && z.getRows () == z.getCols () &&
01533           cy.getRows () == z.getRows ());
01534   return z * cy * adjoint (z);
01535 }
01536 
01543 void matrix::exchangeRows (int r1, int r2) {
01544   nr_complex_t * s = new nr_complex_t[cols];
01545   int len = sizeof (nr_complex_t) * cols;
01546 
01547   assert (r1 >= 0 && r2 >= 0 && r1 < rows && r2 < rows);
01548 
01549   memcpy (s, &data[r1 * cols], len);
01550   memcpy (&data[r1 * cols], &data[r2 * cols], len);
01551   memcpy (&data[r2 * cols], s, len);
01552   delete[] s;
01553 }
01554 
01561 void matrix::exchangeCols (int c1, int c2) {
01562   nr_complex_t s;
01563 
01564   assert (c1 >= 0 && c2 >= 0 && c1 < cols && c2 < cols);
01565 
01566   for (int r = 0; r < rows * cols; r += cols) {
01567     s = data[r + c1];
01568     data[r + c1] = data[r + c2];
01569     data[r + c2] = s;
01570   }
01571 }
01572 
01594 matrix twoport (matrix m, char in, char out) {
01595   assert (m.getRows () >= 2 && m.getCols () >= 2);
01596   nr_complex_t d;
01597   matrix res (2);
01598 
01599   switch (in) {
01600   case 'Y':
01601     switch (out) {
01602     case 'Y': // Y to Y
01603       res = m;
01604       break;
01605     case 'Z': // Y to Z
01606       d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0);
01607       res.set (0, 0, m (1, 1) / d);
01608       res.set (0, 1, -m (0, 1) / d);
01609       res.set (1, 0, -m (1, 0) / d);
01610       res.set (1, 1, m (0, 0) / d);
01611       break;
01612     case 'H': // Y to H
01613       d = m (0, 0);
01614       res.set (0, 0, 1.0 / d);
01615       res.set (0, 1, -m (0, 1) / d);
01616       res.set (1, 0, m (1, 0) / d);
01617       res.set (1, 1, m (1, 1) - m (0, 1) * m (1, 0) / d);
01618       break;
01619     case 'G': // Y to G
01620       d = m (1, 1);
01621       res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d);
01622       res.set (0, 1, m (0, 1) / d);
01623       res.set (1, 0, -m (1, 0) / d);
01624       res.set (1, 1, 1.0 / d);
01625       break;
01626     case 'A': // Y to A
01627       d = m (1, 0);
01628       res.set (0, 0, -m (1, 1) / d);
01629       res.set (0, 1, -1.0 / d);
01630       res.set (1, 0, m (0, 1) - m (1, 1) * m (0, 0) / d);
01631       res.set (1, 1, -m (0, 0) / d);
01632       break;
01633     case 'S': // Y to S
01634       res = ytos (m);
01635       break;
01636     }
01637     break;
01638   case 'Z':
01639     switch (out) {
01640     case 'Y': // Z to Y
01641       d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0);
01642       res.set (0, 0, m (1, 1) / d);
01643       res.set (0, 1, -m (0, 1) / d);
01644       res.set (1, 0, -m (1, 0) / d);
01645       res.set (1, 1, m (0, 0) / d);
01646       break;
01647     case 'Z': // Z to Z
01648       res = m;
01649       break;
01650     case 'H': // Z to H
01651       d = m (1, 1);
01652       res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d);
01653       res.set (0, 1, m (0, 1) / d);
01654       res.set (1, 0, -m (1, 0) / d);
01655       res.set (1, 1, 1.0 / d);
01656       break;
01657     case 'G': // Z to G
01658       d = m (0, 0);
01659       res.set (0, 0, 1.0 / d);
01660       res.set (0, 1, -m (0, 1) / d);
01661       res.set (1, 0, m (1, 0) / d);
01662       res.set (1, 1, m (1, 1) - m (0, 1) * m (1, 0) / d);
01663       break;
01664     case 'A': // Z to A
01665       d = m (1, 0);
01666       res.set (0, 0, m (0, 0) / d);
01667       res.set (0, 1, m (0, 0) * m (1, 1) / d - m (0, 1));
01668       res.set (1, 0, 1.0 / d);
01669       res.set (1, 1, m (1, 1) / d);
01670       break;
01671     case 'S': // Z to S
01672       res = ztos (m);
01673       break;
01674     }
01675     break;
01676   case 'H':
01677     switch (out) {
01678     case 'Y': // H to Y
01679       d = m (0, 0);
01680       res.set (0, 0, 1.0 / d);
01681       res.set (0, 1, -m (0, 1) / d);
01682       res.set (1, 0, m (1, 0) / d);
01683       res.set (1, 1, m (1, 1) - m (0, 1) * m.get(2, 1) / d);
01684       break;
01685     case 'Z': // H to Z
01686       d = m (1, 1);
01687       res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d);
01688       res.set (0, 1, m (0, 1) / d);
01689       res.set (1, 0, -m (1, 0) / d);
01690       res.set (1, 1, 1.0 / d);
01691       break;
01692     case 'H': // H to H
01693       res = m;
01694       break;
01695     case 'G': // H to G
01696       d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0);
01697       res.set (0, 0, m (1, 1) / d);
01698       res.set (0, 1, -m (0, 1) / d);
01699       res.set (1, 0, -m (1, 0) / d);
01700       res.set (1, 1, m (0, 0) / d);
01701       break;
01702     case 'A': // H to A
01703       d = m (1, 0);
01704       res.set (0, 0, m (0, 1) - m (0, 0) * m (1, 1) / d);
01705       res.set (0, 1, -m (0, 0) / d);
01706       res.set (1, 0, -m (1, 1) / d);
01707       res.set (1, 1, -1.0 / d);
01708       break;
01709     case 'S': // H to S
01710       res = htos (m);
01711       break;
01712     }
01713     break;
01714   case 'G':
01715     switch (out) {
01716     case 'Y': // G to Y
01717       d = m (1, 1);
01718       res.set (0, 0, m (0, 0) - m (0, 1) * m (1, 0) / d);
01719       res.set (0, 1, m (0, 1) / d);
01720       res.set (1, 0, -m (1, 0) / d);
01721       res.set (1, 1, 1.0 / d);
01722       break;
01723     case 'Z': // G to Z
01724       d = m (0, 0);
01725       res.set (0, 0, 1.0 / d);
01726       res.set (0, 1, -m (0, 1) / d);
01727       res.set (1, 0, m (1, 0) / d);
01728       res.set (1, 1, m (1, 1) - m (0, 1) * m (1, 0) / d);
01729       break;
01730     case 'H': // G to H
01731       d = m (0, 0) * m (1, 1) - m (0, 1) * m (1, 0);
01732       res.set (0, 0, m (1, 1) / d);
01733       res.set (0, 1, -m (0, 1) / d);
01734       res.set (1, 0, -m (1, 0) / d);
01735       res.set (1, 1, m (0, 0) / d);
01736       break;
01737     case 'G': // G to G
01738       res = m;
01739       break;
01740     case 'A': // G to A
01741       d = m (1, 0);
01742       res.set (0, 0, 1.0 / d);
01743       res.set (0, 1, m (1, 1) / d);
01744       res.set (1, 0, m (0, 0) / d);
01745       res.set (1, 1, m (0, 0) * m (1, 1) / d - m (0, 1));
01746       break;
01747     case 'S': // G to S
01748       res = gtos (m);
01749       break;
01750     }
01751     break;
01752   case 'A':
01753     switch (out) {
01754     case 'Y': // A to Y
01755       d = m (0, 1);
01756       res.set (0, 0, m (1, 1) / d);
01757       res.set (0, 1, m (1, 0) - m (0, 0) * m (1, 1) / d);
01758       res.set (1, 0, -1.0 / d);
01759       res.set (1, 1, m (0, 0) / d);
01760       break;
01761     case 'Z': // A to Z
01762       d = m (1, 0);
01763       res.set (0, 0, m (0, 0) / d);
01764       res.set (0, 1, m (0, 0) * m (1, 1) / d - m (0, 1));
01765       res.set (1, 0, 1.0 / d);
01766       res.set (1, 1, m (1, 1) / d);
01767       break;
01768     case 'H': // A to H
01769       d = m (1, 1);
01770       res.set (0, 0, m (0, 1) / d);
01771       res.set (0, 1, m (0, 0) - m (0, 1) * m (1, 0) / d);
01772       res.set (1, 0, -1.0 / d);
01773       res.set (1, 1, m (1, 0) / d);
01774       break;
01775     case 'G': // A to G
01776       d = m (0, 0);
01777       res.set (0, 0, m (1, 0) / d);
01778       res.set (0, 1, m (1, 0) * m (0, 1) / d - m (1, 1));
01779       res.set (1, 0, 1.0 / d);
01780       res.set (1, 1, m (0, 1) / d);
01781       break;
01782     case 'A': // A to A
01783       res = m;
01784       break;
01785     case 'S': // A to S
01786       res = atos (m);
01787       break;
01788     }
01789     break;
01790   case 'S':
01791     switch (out) {
01792     case 'S': // S to S
01793       res = m;
01794       break;
01795     case 'T': // S to T
01796       d = m (1, 0);
01797       res.set (0, 0, m (0, 1) - m (0, 0) * m (1, 1) / d);
01798       res.set (0, 1, m (0, 0) / d);
01799       res.set (1, 0, -m (1, 1) / d);
01800       res.set (0, 1, 1.0 / d);
01801       break;
01802     case 'A': // S to A
01803       res = stoa (m);
01804       break;
01805     case 'H': // S to H
01806       res = stoh (m);
01807       break;
01808     case 'G': // S to G
01809       res = stog (m);
01810       break;
01811     case 'Y': // S to Y
01812       res = stoy (m);
01813       break;
01814     case 'Z': // S to Z
01815       res = stoz (m);
01816       break;
01817     }
01818     break;
01819   case 'T':
01820     switch (out) {
01821     case 'S': // T to S
01822       d = m (1, 1);
01823       res.set (0, 0, m (0, 1) / d);
01824       res.set (0, 1, m (0, 0) - m (0, 1) * m (1, 0) / d);
01825       res.set (1, 0, 1.0 / d);
01826       res.set (0, 1, -m (1, 0) / d);
01827       break;
01828     case 'T': // T to T
01829       res = m;
01830       break;
01831     }
01832     break;
01833   }
01834   return res;
01835 }
01836 
01854 nr_double_t rollet (matrix m) {
01855   assert (m.getRows () >= 2 && m.getCols () >= 2);
01856   nr_double_t res;
01857   res = (1 - norm (m (0, 0)) - norm (m (1, 1)) + norm (det (m))) /
01858     2 / abs (m (0, 1) * m (1, 0));
01859   return res;
01860 }
01861 
01862 /* Computes stability measure B1 of the given S-parameter matrix. */
01863 nr_double_t b1 (matrix m) {
01864   assert (m.getRows () >= 2 && m.getCols () >= 2);
01865   nr_double_t res;
01866   res = 1 + norm (m (0, 0)) - norm (m (1, 1)) - norm (det (m));
01867   return res;
01868 }
01869 
01870 
01871 matrix rad2deg (matrix a) {
01872   matrix res (a.getRows (), a.getCols ());
01873   for (int r = 0; r < a.getRows (); r++)
01874     for (int c = 0; c < a.getCols (); c++)
01875       res.set (r, c, rad2deg (a.get (r, c)));
01876   return res;
01877 }
01878 
01879 matrix deg2rad (matrix a) {
01880   matrix res (a.getRows (), a.getCols ());
01881   for (int r = 0; r < a.getRows (); r++)
01882     for (int c = 0; c < a.getCols (); c++)
01883       res.set (r, c, deg2rad (a.get (r, c)));
01884   return res;
01885 }
01886 
01887 } // namespace qucs