Qucs-core  0.0.19
cpwline.cpp
Go to the documentation of this file.
00001 /*
00002  * cpwline.cpp - coplanar waveguide line class implementation
00003  *
00004  * Copyright (C) 2004, 2005 Vincent Habchi, F5RCS <10.50@free.fr>
00005  * Copyright (C) 2004, 2005, 2006, 2008 Stefan Jahn <stefan@lkcc.org>
00006  *
00007  * This is free software; you can redistribute it and/or modify
00008  * it under the terms of the GNU General Public License as published by
00009  * the Free Software Foundation; either version 2, or (at your option)
00010  * any later version.
00011  *
00012  * This software is distributed in the hope that it will be useful,
00013  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00014  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00015  * GNU General Public License for more details.
00016  *
00017  * You should have received a copy of the GNU General Public License
00018  * along with this package; see the file COPYING.  If not, write to
00019  * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
00020  * Boston, MA 02110-1301, USA.
00021  *
00022  * $Id$
00023  *
00024  */
00025 
00026 #if HAVE_CONFIG_H
00027 # include <config.h>
00028 #endif
00029 
00030 #include <limits>
00031 
00032 #include "component.h"
00033 #include "substrate.h"
00034 #include "cpwline.h"
00035 
00036 using namespace qucs;
00037 
00038 cpwline::cpwline () : circuit (2) {
00039   Zl = Er = 0;
00040   type = CIR_CPWLINE;
00041 }
00042 
00050 /* The function computes the complete elliptic integral of first kind
00051    K(k) using the arithmetic-geometric mean algorithm (AGM) found e.g.
00052    in Abramowitz and Stegun (17.6.1). 
00053    Note that the argument of the function here is the elliptic modulus k
00054    and not the parameter m=k^2 . */
00055 /* \todo move to common math */
00056 nr_double_t cpwline::ellipk (nr_double_t k) {
00057   if ((k < 0.0) || (k >= 1.0))
00058     // we use only the range from 0 <= k < 1
00059     return std::numeric_limits<nr_double_t>::quiet_NaN();
00060 
00061   nr_double_t a = 1.0;
00062   nr_double_t b = qucs::sqrt(1-k*k);
00063   nr_double_t c = k;
00064 
00065   while (c > std::numeric_limits<nr_double_t>::epsilon()) {
00066     nr_double_t tmp = (a + b) / 2.0;
00067     c = (a - b) / 2.0;
00068     b = qucs::sqrt(a * b);
00069     a = tmp;
00070   }
00071   return (pi_over_2 / a);
00072 }
00073 
00074 nr_double_t cpwline::KoverKp(nr_double_t k) {
00075   if ((k < 0.0) || (k >= 1.0))
00076     return std::numeric_limits<nr_double_t>::quiet_NaN();
00077 
00078   return (ellipk(k) / ellipk(qucs::sqrt(1-k*k)));
00079 }
00080 
00081 /* Approximation of K(k)/K'(k).
00082    First appeared in
00083    Hilberg, W., "From Approximations to Exact Relations for Characteristic
00084    Impedances," IEEE Trans. MTT, May 1969.
00085    More accurate expressions can be found in the above article and in
00086    Abbott, J. T., "Modeling the Capacitive Behavior of Coplanar Striplines
00087    and Coplanar Waveguides using Simple Functions", Rochester Institute of
00088    Technology, Rochester, New York, June 2011.
00089    The maximum relative error of the approximation implemented here is
00090    about 2 ppm, so good enough for any practical purpose.
00091  */
00092 nr_double_t cpwline::ellipa (nr_double_t k) {
00093   nr_double_t r, kp;
00094   if (k < sqrt1_2) {
00095     kp = qucs::sqrt (1 - k * k);
00096     r = pi / qucs::log (2 * (1 + qucs::sqrt (kp)) / (1 - qucs::sqrt (kp)));
00097   }
00098   else {
00099     r = qucs::log (2 * (1 + qucs::sqrt (k)) / (1 - qucs::sqrt (k))) / pi;
00100   }
00101   return r;
00102 }
00103 
00104 void cpwline::initSP (void) {
00105   // allocate S-parameter matrix
00106   allocMatrixS ();
00107   // pre-compute propagation factors
00108   initPropagation ();
00109 }
00110 
00111 void cpwline::initPropagation (void) {
00112   // get properties of substrate and coplanar line
00113   nr_double_t W =  getPropertyDouble ("W");
00114   nr_double_t s =  getPropertyDouble ("S");
00115   substrate * subst = getSubstrate ();
00116   nr_double_t er = subst->getPropertyDouble ("er");
00117   nr_double_t h  = subst->getPropertyDouble ("h");
00118   nr_double_t t  = subst->getPropertyDouble ("t");
00119   int backMetal  = !strcmp (getPropertyString ("Backside"), "Metal");
00120   int approx     = !strcmp (getPropertyString ("Approx"), "yes");
00121 
00122   tand = subst->getPropertyDouble ("tand");
00123   rho  = subst->getPropertyDouble ("rho");
00124   len  = getPropertyDouble ("L");
00125 
00126   // other local variables (quasi-static constants)
00127   nr_double_t k1, kk1, kpk1, k2, k3, q1, q2, q3 = 0, qz, er0 = 0;
00128 
00129   // compute the necessary quasi-static approx. (K1, K3, er(0) and Z(0))
00130   k1   = W / (W + s + s);
00131   kk1  = ellipk (k1);
00132   kpk1 = ellipk (qucs::sqrt (1 - k1 * k1));
00133   if (approx) {
00134     q1 = ellipa (k1);
00135   } else {
00136     q1 = kk1 / kpk1;
00137   }
00138 
00139   // backside is metal
00140   if (backMetal) {
00141     k3  = qucs::tanh ((pi / 4) * (W / h)) / qucs::tanh ((pi / 4) * (W + s + s) / h);
00142     if (approx) {
00143       q3 = ellipa (k3);
00144     } else {
00145       q3 = ellipk (k3) / ellipk (qucs::sqrt (1 - k3 * k3));
00146     }
00147     qz  = 1 / (q1 + q3);
00148     er0 = 1 + q3 * qz * (er - 1);
00149     zl_factor = Z0 / 2 * qz;
00150   }
00151   // backside is air
00152   else {
00153     k2  = qucs::sinh ((pi / 4) * (W / h)) / qucs::sinh ((pi / 4) * (W + s + s) / h);
00154     if (approx) {
00155       q2 = ellipa (k2);
00156     } else {
00157       q2 = ellipk (k2) / ellipk (qucs::sqrt (1 - k2 * k2));
00158     }
00159     er0 = 1 + (er - 1) / 2 * q2 / q1;
00160     zl_factor = Z0 / 4 / q1;
00161   }
00162 
00163   // adds effect of strip thickness
00164   if (t > 0) {
00165     nr_double_t d, ke, qe;
00166     d  = (t * 1.25 / pi) * (1 + qucs::log (4 * pi * W / t));
00167 
00168     // modifies k1 accordingly (k1 = ke)
00169     ke = k1 + (1 - k1 * k1) * d / 2 / s;
00170     if (approx) {
00171       qe = ellipa (ke);
00172     } else {
00173       qe = ellipk (ke) / ellipk (qucs::sqrt (1 - ke * ke));
00174     }
00175     // backside is metal
00176     if (backMetal) {
00177       qz  = 1 / (qe + q3);
00178       //er0 = 1 + q3 * qz * (er - 1);
00179       zl_factor = Z0 / 2 * qz;
00180     }
00181     // backside is air
00182     else {
00183       zl_factor = Z0 / 4 / qe;
00184     }
00185 
00186     // modifies er0 as well
00187     er0 = er0 - (0.7 * (er0 - 1) * t / s) / (q1 + (0.7 * t / s));
00188   }
00189 
00190   // pre-compute square roots
00191   sr_er = qucs::sqrt (er);
00192   sr_er0 = qucs::sqrt (er0);
00193 
00194   // cut-off frequency of the TE0 mode
00195   fte = (C0 / 4) / (h * qucs::sqrt (er - 1));
00196 
00197   // dispersion factor G
00198   nr_double_t p = qucs::log (W / h);
00199   nr_double_t u = 0.54 - (0.64 - 0.015 * p) * p;
00200   nr_double_t v = 0.43 - (0.86 - 0.54 * p) * p;
00201   G = qucs::exp (u * qucs::log (W / s) + v);
00202 
00203   // loss constant factors (computed only once for efficency sake)
00204   nr_double_t ac = 0;
00205   if (t > 0) {
00206     // equations by GHIONE
00207     nr_double_t n  = (1 - k1) * 8 * pi / (t * (1 + k1));
00208     nr_double_t a  = W / 2;
00209     nr_double_t b  = a + s;
00210     ac = (pi + qucs::log (n * a)) / a + (pi + qucs::log (n * b)) / b;
00211   }
00212   ac_factor  = ac / (4 * Z0 * kk1 * kpk1 * (1 - k1 * k1));
00213   ac_factor *= qucs::sqrt (pi * MU0 * rho); // Rs factor
00214   ad_factor  = (er / (er - 1)) * tand * pi / C0;
00215 
00216   // propagation constant (partial, final value computed in calcAB() )
00217   bt_factor  = 2 * pi / C0;
00218 }
00219 
00220 void cpwline::calcAB (nr_double_t f, nr_double_t& zl, nr_double_t& al,
00221                       nr_double_t& bt) {
00222   nr_double_t sr_er_f = sr_er0;
00223   nr_double_t ac = ac_factor;
00224   nr_double_t ad = ad_factor;
00225 
00226   // by initializing as much as possible outside this function, the
00227   // overhead is minimal
00228 
00229   // add the dispersive effects to er0
00230   sr_er_f += (sr_er - sr_er0) / (1 + G * qucs::pow (f / fte, -1.8));
00231 
00232   // computes impedance
00233   zl /= sr_er_f;
00234 
00235   // for now, the loss are limited to strip losses (no radiation
00236   // losses yet) losses in neper/length
00237   ad *= f * (sr_er_f - 1 / sr_er_f);
00238   ac *= qucs::sqrt (f) * sr_er0;
00239 
00240   al  = ac + ad;
00241   bt *= sr_er_f * f;
00242 
00243   Er = sqr (sr_er_f);
00244   Zl = zl;
00245 }
00246 
00247 void cpwline::saveCharacteristics (nr_double_t) {
00248   setCharacteristic ("Zl", Zl);
00249   setCharacteristic ("Er", Er);
00250 }
00251 
00252 void cpwline::calcSP (nr_double_t frequency) {
00253 
00254   nr_double_t zl = zl_factor;
00255   nr_double_t beta = bt_factor;
00256   nr_double_t alpha;
00257 
00258   calcAB (frequency, zl, alpha, beta);
00259 
00260   // calculate and set S-parameters
00261   nr_double_t z = zl / z0;
00262   nr_double_t y = 1 / z;
00263   nr_complex_t g = nr_complex_t (alpha, beta);
00264   nr_complex_t n = 2.0 * cosh (g * len) + (z + y) * sinh (g * len);
00265   nr_complex_t s11 = (z - y) * sinh (g * len) / n;
00266   nr_complex_t s21 = 2.0 / n;
00267 
00268   setS (NODE_1, NODE_1, s11); setS (NODE_2, NODE_2, s11);
00269   setS (NODE_1, NODE_2, s21); setS (NODE_2, NODE_1, s21);
00270 }
00271 
00272 /* FIXME : following function is unused? */
00273 /* The function calculates the quasi-static impedance of a coplanar
00274    waveguide line and the value of the effective dielectric constant
00275    for the given coplanar line and substrate properties. */
00276 void cpwline::analyseQuasiStatic (nr_double_t W, nr_double_t s, nr_double_t h,
00277                                   nr_double_t t, nr_double_t er, int backMetal,
00278                                   nr_double_t& ZlEff, nr_double_t& ErEff) {
00279 
00280   // local variables (quasi-static constants)
00281   nr_double_t k1, k2, k3, q1, q2, q3 = 0, qz;
00282 
00283   ErEff = er;
00284   ZlEff = 0;
00285 
00286   // compute the necessary quasi-static approx. (K1, K3, er(0) and Z(0))
00287   k1 = W / (W + s + s);
00288   q1 = ellipk (k1) / ellipk (qucs::sqrt (1 - k1 * k1));
00289 
00290   // backside is metal
00291   if (backMetal) {
00292     k3  = qucs::tanh ((pi / 4) * (W / h)) / qucs::tanh ((pi / 4) * (W + s + s) / h);
00293     q3 = ellipk (k3) / ellipk (qucs::sqrt (1 - k3 * k3));
00294     qz  = 1 / (q1 + q3);
00295     ErEff = 1 + q3 * qz * (er - 1);
00296     ZlEff = Z0 / 2 * qz;
00297   }
00298   // backside is air
00299   else {
00300     k2  = qucs::sinh ((pi / 4) * (W / h)) / qucs::sinh ((pi / 4) * (W + s + s) / h);
00301     q2 = ellipk (k2) / ellipk (qucs::sqrt (1 - k2 * k2));
00302     ErEff = 1 + (er - 1) / 2 * q2 / q1;
00303     ZlEff = Z0 / 4 / q1;
00304   }
00305 
00306   // adds effect of strip thickness
00307   if (t > 0) {
00308     nr_double_t d, ke, qe;
00309     d  = (t * 1.25 / pi) * (1 + qucs::log (4 * pi * W / t));
00310 
00311     // modifies k1 accordingly (k1 = ke)
00312     ke = k1 + (1 - k1 * k1) * d / 2 / s;
00313     qe = ellipk (ke) / ellipk (qucs::sqrt (1 - ke * ke));
00314 
00315     // backside is metal
00316     if (backMetal) {
00317       qz  = 1 / (qe + q3);
00318       //ErEff = 1 + q3 * qz * (er - 1);
00319       ZlEff = Z0 / 2 * qz;
00320     }
00321     // backside is air
00322     else {
00323       ZlEff = Z0 / 4 / qe;
00324     }
00325 
00326     // modifies ErEff as well
00327     ErEff = ErEff - (0.7 * (ErEff - 1) * t / s) / (q1 + (0.7 * t / s));
00328   }
00329   ErEff = qucs::sqrt (ErEff);
00330   ZlEff /= ErEff;
00331 }
00332 
00333 /* FIXME : following function is unused? */
00334 /* This function calculates the frequency dependent value of the
00335    effective dielectric constant and the coplanar line impedance for
00336    the given frequency. */
00337 void cpwline::analyseDispersion (nr_double_t W, nr_double_t s, nr_double_t h,
00338                                  nr_double_t er, nr_double_t ZlEff,
00339                                  nr_double_t ErEff, nr_double_t frequency,
00340                                  nr_double_t& ZlEffFreq,
00341                                  nr_double_t& ErEffFreq) {
00342   // local variables
00343   nr_double_t fte, G;
00344 
00345   ErEffFreq = ErEff;
00346   ZlEffFreq = ZlEff * ErEff;
00347 
00348   // cut-off frequency of the TE0 mode
00349   fte = (C0 / 4) / (h * qucs::sqrt (er - 1));
00350 
00351   // dispersion factor G
00352   nr_double_t p = qucs::log (W / h);
00353   nr_double_t u = 0.54 - (0.64 - 0.015 * p) * p;
00354   nr_double_t v = 0.43 - (0.86 - 0.54 * p) * p;
00355   G = qucs::exp (u * qucs::log (W / s) + v);
00356 
00357   // add the dispersive effects to er0
00358   ErEffFreq += (qucs::sqrt (er) - ErEff) / (1 + G * qucs::pow (frequency / fte, -1.8));
00359 
00360   // computes impedance
00361   ZlEffFreq /= ErEffFreq;
00362 }
00363 
00364 void cpwline::calcNoiseSP (nr_double_t) {
00365   // calculate noise using Bosma's theorem
00366   nr_double_t T = getPropertyDouble ("Temp");
00367   matrix s = getMatrixS ();
00368   matrix e = eye (getSize ());
00369   setMatrixN (celsius2kelvin (T) / T0 * (e - s * transpose (conj (s))));
00370 }
00371 
00372 void cpwline::initDC (void) {
00373   // a DC short (voltage source V = 0 volts)
00374   setVoltageSources (1);
00375   setInternalVoltageSource (1);
00376   allocMatrixMNA ();
00377   clearY ();
00378   voltageSource (VSRC_1, NODE_1, NODE_2);
00379 }
00380 
00381 void cpwline::initTR (void) {
00382   initDC ();
00383 }
00384 
00385 void cpwline::initAC (void) {
00386   setVoltageSources (0);
00387   allocMatrixMNA ();
00388   initPropagation ();
00389 }
00390 
00391 void cpwline::calcAC (nr_double_t frequency) {
00392 
00393   nr_double_t zl = zl_factor;
00394   nr_double_t beta = bt_factor;
00395   nr_double_t alpha;
00396 
00397   calcAB (frequency, zl, alpha, beta);
00398 
00399   // calculate and set Y-parameters
00400   nr_complex_t g = nr_complex_t (alpha, beta);
00401   nr_complex_t y11 = coth (g * len) / zl;
00402   nr_complex_t y21 = -cosech (g * len) / zl;
00403 
00404   setY (NODE_1, NODE_1, y11); setY (NODE_2, NODE_2, y11);
00405   setY (NODE_1, NODE_2, y21); setY (NODE_2, NODE_1, y21);
00406 }
00407 
00408 void cpwline::calcNoiseAC (nr_double_t) {
00409   // calculate noise using Bosma's theorem
00410   nr_double_t T = getPropertyDouble ("Temp");
00411   setMatrixN (4 * celsius2kelvin (T) / T0 * real (getMatrixY ()));
00412 }
00413 
00414 // properties
00415 PROP_REQ [] = {
00416   { "W", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_POS_RANGE },
00417   { "S", PROP_REAL, { 1e-3, PROP_NO_STR }, PROP_POS_RANGE },
00418   { "L", PROP_REAL, { 10e-3, PROP_NO_STR }, PROP_POS_RANGE },
00419   { "Subst", PROP_STR, { PROP_NO_VAL, "Subst1" }, PROP_NO_RANGE },
00420   PROP_NO_PROP };
00421 PROP_OPT [] = {
00422   { "Temp", PROP_REAL, { 26.85, PROP_NO_STR }, PROP_MIN_VAL (K) },
00423   { "Backside", PROP_STR, { PROP_NO_VAL, "Metal" },
00424     PROP_RNG_STR2 ("Metal", "Air") },
00425   { "Approx", PROP_STR, { PROP_NO_VAL, "no" }, PROP_RNG_YESNO },
00426   PROP_NO_PROP };
00427 struct define_t cpwline::cirdef =
00428   { "CLIN", 2, PROP_COMPONENT, PROP_NO_SUBSTRATE, PROP_LINEAR, PROP_DEF };