rectline.cpp

Go to the documentation of this file.

/*
 * rectline.cpp - rectangular waveguide class implementation
 *
 * Copyright (C) 2008 Bastien ROUCARIES <roucaries.bastien@gmail.com>
 * Copyright (C) 2008 Andrea Zonca <andrea.zonca@gmail.com>
 *
 * This is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2, or (at your option)
 * any later version.
 *
 * This software is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this package; see the file COPYING.  If not, write to
 * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
 * Boston, MA 02110-1301, USA.
 *
 * $Id$
 *
 */

#if HAVE_CONFIG_H
# include <config.h>
#endif

#include "component.h"
#include "rectline.h"

using namespace qucs;

rectline::rectline () : circuit (2) {
  alpha = beta = fc_low = fc_high = 0.0;
  zl = 0.0;
  type = CIR_RECTANGULAR;
}

void rectline::calcResistivity (const char * const Mat, nr_double_t T) {
  if (!strcmp (Mat, "Copper")) {
    if (T < 7) {
      rho = 2e-11;
    }
    else if (T < 15) {
      rho = 6.66667e-17 * qucs::pow (T, 5.0) - 3.88549e-15 * qucs::pow (T, 4.0) +
        9.82267e-14 * qucs::pow (T, 3.0) - 1.29684e-12 * qucs::pow (T, 2.0) +
        8.68341e-12 * T - 2.72120e-12;
    }
    else if (T < 45) {
      rho = 6.60731e-15 * qucs::pow (T, 3.0) - 1.14812e-13 * qucs::pow (T, 2.0) -
        1.11681e-12 * T + 4.23709e-11;
    }
    else if (T < 100) {
      rho = -6.53059e-15 * qucs::pow (T, 3.0) + 1.73783e-12 * qucs::pow (T, 2.0) -
        8.73888e-11 * T + 1.37016e-9;
    }
    else if (T < 350) {
      rho = 1.00018e-17 * qucs::pow (T, 3.0) - 8.72408e-15 * qucs::pow (T, 2.0) +
        7.06020e-11 * T - 3.51125e-9;
    }
    else {
      rho = 0.000000020628;
    }
    // in ADS iT_K is forced T_Ko Cu_300K:
    //rho = 1.7e-8;
  }
  else if (!strcmp (Mat, "StainlessSteel")) {
    rho = 7.4121e-17 * qucs::pow (T, 4.0) - 5.3504e-14 * qucs::pow (T, 3.0) +
      1.2902e-11 * qucs::pow (T, 2.0) - 2.9186e-10 * T +4.9320e-7;
  }
  else if (!strcmp (Mat, "Gold")) {
    if (T < 20) {
      rho = 0.00000000024;
    }
    else if (T < 65) {
      rho = 2e-12 * qucs::pow (T, 2.0) - 8e-11 * T + 1e-9;
    }
    else if (T < 80) {
      rho = 5e-13 * qucs::pow (T, 3.0) - 1e-10 * qucs::pow (T, 2.0) + 9e-9 * T - 2e-7;
    }
    else if (T < 300) {
      rho = 8e-11 * T - 1e-10;
    }
    else {
      rho = 2.4e-8;
    }
  }
}

void rectline::calcPropagation (nr_double_t frequency) {
  nr_double_t er   = getPropertyDouble ("er");
  nr_double_t mur  = getPropertyDouble ("mur");
  nr_double_t tand = getPropertyDouble ("tand");
  nr_double_t a    = getPropertyDouble ("a");
  nr_double_t b    = getPropertyDouble ("b");

  /* wave number */
  nr_double_t k0;
  nr_double_t kc;
  /* dielectric loss */
  nr_double_t ad, ac, rs;

  // check cutoff frequency
  if (frequency >= fc_high) {
    logprint (LOG_ERROR, "WARNING: Operating frequency (%g) outside TE10 "
              "band (%g <= TE10 <= %g) or outside non propagative mode "
              "<= %g\n", frequency, fc_low, fc_high, fc_low);
  }
  // calculate wave number
  k0 = (2.0 * pi * frequency * std::sqrt (er * mur)) / C0;
  kc = pi / a;

  // calculate losses only for propagative mode
  if (frequency >= fc_low) {
    // calculate beta
    beta = std::sqrt (sqr (k0) - sqr (kc));

    // dielectric
    ad = (sqr(k0) * tand) / (2.0 * beta);
    // resistive
    rs = std::sqrt (pi * frequency * mur * MU0 * rho);
    ac = rs * (2.0 * b * sqr (pi) + cubic (a) * sqr (k0)) /
      (cubic (a) * b * beta * k0 *  Z0);
    alpha = (ad + ac);

    // wave impedance
    zl = (k0 * Z0) / beta;

  } else {
    /* according to [2] eq 3.207 */
    beta = 0;
    alpha = -std::sqrt (- (sqr (k0) - sqr (kc)));
    // wave impedance
    zl = (k0 * Z0) / nr_complex_t (0, -alpha) ;
  }
}

void rectline::calcNoiseSP (nr_double_t) {
  nr_double_t l = getPropertyDouble ("L");
  if (l < 0) return;
  // calculate noise using Bosma's theorem
  nr_double_t T = getPropertyDouble ("Temp");
  matrix s = getMatrixS ();
  matrix e = eye (getSize ());
  setMatrixN (celsius2kelvin (T) / T0 * (e - s * transpose (conj (s))));
}

void rectline::initCheck (void) {
  nr_double_t a = getPropertyDouble ("a");
  nr_double_t b = getPropertyDouble ("b");
  nr_double_t epsr = getPropertyDouble ("er");
  nr_double_t mur = getPropertyDouble ("mur");
  // check validity
  if (a < b) {
    logprint (LOG_ERROR, "ERROR: a < b should be a >= b.\n");
  }
  nr_double_t c = std::sqrt (epsr * mur);
  fc_low =  C0 / (2.0 * a * c);
  /* min of second TE mode and first TM mode */
  fc_high = std::min (C0 / (a * c),  C0 / (2.0  * b * c));
  // calculation of resistivity
  rho  = getPropertyDouble ("rho");
  nr_double_t T = getPropertyDouble ("Temp");
  calcResistivity (getPropertyString ("Material"), celsius2kelvin (T));
}

void rectline::saveCharacteristics (nr_complex_t) {
  setCharacteristic ("Zl", real (zl));
}

void rectline::initSP (void) {
  // allocate S-parameter matrix
  allocMatrixS ();
  initCheck ();
}

void rectline::calcSP (nr_double_t frequency) {
  nr_double_t l = getPropertyDouble ("L");

  // calculate propagation constants
  calcPropagation (frequency);

  // calculate S-parameters
  nr_complex_t z = zl / z0;
  nr_complex_t y = 1.0 / z;
  nr_complex_t g = nr_complex_t (alpha, beta);
  nr_complex_t n = 2.0 * cosh (g * l) + (z + y) * sinh (g * l);
  nr_complex_t s11 = (z - y) * sinh (g * l) / n;
  nr_complex_t s21 = 2.0 / n;
  setS (NODE_1, NODE_1, s11); setS (NODE_2, NODE_2, s11);
  setS (NODE_1, NODE_2, s21); setS (NODE_2, NODE_1, s21);
}

/* ! Compute DC
     \note below cut off it is an open circuit
*/
void rectline::initDC (void) {
  allocMatrixMNA ();
  // open circuit
  clearY ();
}

void rectline::initAC (void) {
  setVoltageSources (0);
  allocMatrixMNA ();
  initCheck ();
}

void rectline::calcAC (nr_double_t frequency) {
  nr_double_t l = getPropertyDouble ("L");

  // calculate propagation constants
  calcPropagation (frequency);

  // calculate Y-parameters
  nr_complex_t g = nr_complex_t (alpha, beta);
  nr_complex_t y11 = coth (g * l) / zl;
  nr_complex_t y21 = -cosech (g * l) / zl;
  setY (NODE_1, NODE_1, y11); setY (NODE_2, NODE_2, y11);
  setY (NODE_1, NODE_2, y21); setY (NODE_2, NODE_1, y21);
}

void rectline::calcNoiseAC (nr_double_t) {
  nr_double_t l = getPropertyDouble ("L");
  if (l < 0) return;
  // calculate noise using Bosma's theorem
  nr_double_t T = getPropertyDouble ("Temp");
  setMatrixN (4.0 * celsius2kelvin (T) / T0 * real (getMatrixY ()));
}

// properties
PROP_REQ [] = {
  { "a", PROP_REAL, { 2.86e-2, PROP_NO_STR }, PROP_POS_RANGEX },
  { "b", PROP_REAL, { 1.016e-2, PROP_NO_STR }, PROP_POS_RANGEX },
  { "L", PROP_REAL, { 1500e-3, PROP_NO_STR }, PROP_NO_RANGE },
  { "er", PROP_REAL, { 1, PROP_NO_STR }, PROP_RNGII (1, 100) },
  { "mur", PROP_REAL, { 1, PROP_NO_STR }, PROP_RNGII (1, 100) },
  { "tand", PROP_REAL, { 4e-4, PROP_NO_STR }, PROP_POS_RANGE },
  { "rho", PROP_REAL, { 0.022e-6, PROP_NO_STR }, PROP_POS_RANGE },
  PROP_NO_PROP };
PROP_OPT [] = {
  { "Temp", PROP_REAL, { 26.85, PROP_NO_STR }, PROP_MIN_VAL (K) },
  { "Material", PROP_STR, { PROP_NO_VAL, "unspecified" },
    PROP_RNG_STR4 ("unspecified", "Copper", "StainlessSteel", "Gold") },
  PROP_NO_PROP };
struct define_t rectline::cirdef =
  { "RECTLINE", 2, PROP_COMPONENT, PROP_NO_SUBSTRATE, PROP_LINEAR, PROP_DEF };