/*
 * microstrip.cpp - microstrip class implementation
 * 
 * Copyright (C) 2001 Gopal Narayanan <gopal@astro.umass.edu>
 * Copyright (C) 2002 Claudio Girardi <claudio.girardi@ieee.org>
 * Copyright (C) 2005, 2006 Stefan Jahn <stefan@lkcc.org>
 * 
 * This program 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 of the License, or (at
 * your option) any later version.
 * 
 * This program 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.  
 *
 */


/* microstrip.c - Puts up window for microstrip and 
 * performs the associated calculations
 * Based on the original microstrip.c by Gopal Narayanan 
 */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <cmath>

#include "units.h"
#include "transline.h"
#include "microstrip.h"

microstrip::microstrip() : transline()
{
}

microstrip::~microstrip()
{
}

/*
 * Z0_homogeneous() - compute the impedance for a stripline in a
 * homogeneous medium, without cover effects
 */
double microstrip::Z0_homogeneous(double u)
{
  double f, Z0;
  f = 6.0 + (2.0 * pi - 6.0) * exp(-pow(30.666 / u, 0.7528));
  Z0 = (ZF0 / (2.0 * pi)) * log(f / u + sqrt(1.0 + 4.0 / (u * u)));
  return Z0;
}


/*
 * delta_Z0_cover() - compute the cover effect on impedance for a
 * stripline in a homogeneous medium
 */
double microstrip::delta_Z0_cover(double u, double h2h)
{
  double P, Q;
  double h2hp1;
  h2hp1 = 1.0 + h2h;
  P = 270.0 * (1.0 - tanh(1.192 + 0.706 * sqrt(h2hp1) - 1.389 / h2hp1));
  Q = 1.0109 - atanh((0.012 * u + 0.177 * u * u - 0.027 * u * u * u) / (h2hp1 * h2hp1));
  return (P * Q);
}


/*
 * filling_factor() - compute the filling factor for a microstrip
 * without cover and zero conductor thickness
 */
double microstrip::filling_factor(double u, double e_r)
{
  double a, b, q_inf;
  double u2, u3, u4;
  u2 = u * u;
  u3 = u2 * u;
  u4 = u3 * u;
  a = 1.0 + log((u4 + u2 / 2704) / (u4 + 0.432)) / 49.0 + log(1.0 + u3 / 5929.741) / 18.7;
  b = 0.564 * pow((e_r - 0.9) / (e_r + 3.0), 0.053);
  q_inf = pow(1.0 + 10.0 / u, -a * b);
  return q_inf;
}


/*
 * delta_q_cover() - compute the cover effect on filling factor
 */
double microstrip::delta_q_cover(double h2h)
{
  double q_c;
  q_c = tanh(1.043 + 0.121 * h2h - 1.164 / h2h);
  return q_c;
}


/*
 * delta_q_thickness() - compute the thickness effect on filling factor
 */
double microstrip::delta_q_thickness(double u, double t_h)
{
  double q_t;
  q_t = (2.0 * log(2.0) / pi) * (t_h / sqrt(u));
  return q_t;
}


/*
 * e_r_effective() - compute effective dielectric constant from
 * material e_r and filling factor
 */
double microstrip::e_r_effective(double e_r, double q)
{
  double e_r_eff;
  e_r_eff = 0.5 * (e_r + 1.0) + 0.5 * q * (e_r - 1.0);
  return e_r_eff;
}


/*
 * delta_u_thickness - compute the thickness effect on normalized width
 */
double microstrip::delta_u_thickness(double u, double t_h, double e_r)
{
  double delta_u;
  if (t_h > 0.0) {
    /* correction for thickness for a homogeneous microstrip */
    delta_u = (t_h / pi) * log(1.0 + (4.0 * e) * pow(tanh(sqrt(6.517 * u)), 2.0) / t_h);
    /* correction for strip on a substrate with relative permettivity e_r */
    delta_u = 0.5 * delta_u * (1.0 + 1.0 / cosh(sqrt(e_r - 1.0)));
  } else {
    delta_u = 0.0;
  }
  return delta_u;
}


/*
 * microstrip_Z0() - compute microstrip static impedance
 */
void microstrip::microstrip_Z0()
{
  double e_r, h2, h2h, u, t_h;
  double Z0_h_r, Z0;
  double delta_u_1, delta_u_r, q_inf, q_c, q_t, e_r_eff, e_r_eff_t, q;

  e_r = er;
  h2 = ht;
  h2h = h2 / h;
  u = w / h;
  t_h = t / h;

  /* compute normalized width correction for e_r = 1.0 */
  delta_u_1 = delta_u_thickness(u, t_h, 1.0);
  /* compute homogeneous stripline impedance */
  Z0_h_1 = Z0_homogeneous(u + delta_u_1);
  /* compute normalized width corection */
  delta_u_r = delta_u_thickness(u, t_h, e_r);
  u += delta_u_r;
  /* compute homogeneous stripline impedance */
  Z0_h_r = Z0_homogeneous(u);

  /* filling factor, with width corrected for thickness */
  q_inf = filling_factor(u, e_r);
  /* cover effect */
  q_c = delta_q_cover(h2h);
  /* thickness effect */
  q_t = delta_q_thickness(u, t_h);
  /* resultant filling factor */
  q = (q_inf - q_t) * q_c;

  /* e_r corrected for thickness and non homogeneous material */
  e_r_eff_t = e_r_effective(e_r, q);

  /* effective dielectric constant */
  e_r_eff = e_r_eff_t * pow(Z0_h_1 / Z0_h_r, 2.0);

  /* characteristic impedance, corrected for thickness, cover */
  /*   and non homogeneous material */
  Z0 = Z0_h_r / sqrt(e_r_eff_t);

  w_eff = u * h;
  er_eff_0 = e_r_eff;
  Z0_0 = Z0;
}


/*
 * e_r_dispersion() - computes the dispersion correction factor for
 * the effective permeability
 */
double microstrip::e_r_dispersion(double u, double e_r, double f_n)
{
  double P_1, P_2, P_3, P_4, P;

  P_1 = 0.27488 + u * (0.6315 + 0.525 / pow(1.0 + 0.0157 * f_n, 20.0)) - 0.065683 * exp(-8.7513 * u);
  P_2 = 0.33622 * (1.0 - exp(-0.03442 * e_r));
  P_3 = 0.0363 * exp(-4.6 * u) * (1.0 - exp(-pow(f_n / 38.7, 4.97)));
  P_4 = 1.0 + 2.751 * (1.0 - exp(-pow(e_r / 15.916, 8.0)));

  P = P_1 * P_2 * pow((P_3 * P_4 + 0.1844) * f_n, 1.5763);

  return P;
}


/*
 * Z0_dispersion() - computes the dispersion correction factor for the
 * characteristic impedance
 */
double microstrip::Z0_dispersion(double u, double e_r, double e_r_eff_0, double e_r_eff_f, double f_n)
{
  double R_1, R_2, R_3, R_4, R_5, R_6, R_7, R_8, R_9, R_10, R_11, R_12, R_13, R_14, R_15, R_16, R_17, D, tmpf;

  R_1 = 0.03891 * pow(e_r, 1.4);
  R_2 = 0.267 * pow(u, 7.0);
  R_3 = 4.766 * exp(-3.228 * pow(u, 0.641));
  R_4 = 0.016 + pow(0.0514 * e_r, 4.524);
  R_5 = pow(f_n / 28.843, 12.0);
  R_6 = 22.2 * pow(u, 1.92);
  R_7 = 1.206 - 0.3144 * exp(-R_1) * (1.0 - exp(-R_2));
  R_8 = 1.0 + 1.275 * (1.0 - exp(-0.004625 * R_3 * pow(e_r, 1.674) * pow(f_n / 18.365, 2.745)));
  tmpf = pow(e_r - 1.0, 6.0);
  R_9 = 5.086 * R_4 * (R_5 / (0.3838 + 0.386 * R_4)) * (exp(-R_6) / (1.0 + 1.2992 * R_5)) * (tmpf / (1.0 + 10.0 * tmpf));
  R_10 = 0.00044 * pow(e_r, 2.136) + 0.0184;
  tmpf = pow(f_n / 19.47, 6.0);
  R_11 = tmpf / (1.0 + 0.0962 * tmpf);
  R_12 = 1.0 / (1.0 + 0.00245 * u * u);
  R_13 = 0.9408 * pow(e_r_eff_f, R_8) - 0.9603;
  R_14 = (0.9408 - R_9) * pow(e_r_eff_0, R_8) - 0.9603;
  R_15 = 0.707 * R_10 * pow(f_n / 12.3, 1.097);
  R_16 = 1.0 + 0.0503 * e_r * e_r * R_11 * (1.0 - exp(-pow(u / 15.0, 6.0)));
  R_17 = R_7 * (1.0 - 1.1241 * (R_12 / R_16) * exp(-0.026 * pow(f_n, 1.15656) - R_15));

  D = pow(R_13 / R_14, R_17);

  return D;
}


/*
 * dispersion() - compute frequency dependent parameters of
 * microstrip
 */
void microstrip::dispersion()
{
  double e_r, e_r_eff_0;
  double u, f_n, P, e_r_eff_f, D, Z0_f;

  e_r = er;
  e_r_eff_0 = er_eff_0;
  u = w / h;

  /* normalized frequency [GHz * mm] */
  f_n = f * h / 1e06;

  P = e_r_dispersion(u, e_r, f_n);
  /* effective dielectric constant corrected for dispersion */
  e_r_eff_f = e_r - (e_r - e_r_eff_0) / (1.0 + P);

  D = Z0_dispersion(u, e_r, e_r_eff_0, e_r_eff_f, f_n);
  Z0_f = Z0_0 * D;

  er_eff = e_r_eff_f;
  Z0 = Z0_f;
}


/*
 * conductor_losses() - compute microstrip conductor losses per unit
 * length
 */
double microstrip::conductor_losses()
{
  double e_r_eff_0, delta;
  double K, R_s, Q_c, alpha_c;

  e_r_eff_0 = er_eff_0;
  delta = skindepth;

  if (f > 0.0) {
    /* current distribution factor */
    K = exp(-1.2 * pow(Z0_h_1 / ZF0, 0.7));
    /* skin resistance */
    R_s = 1.0 / (sigma * delta);
    
    /* correction for surface roughness */
    R_s *= 1.0 + ((2.0 / pi) * atan(1.40 * pow((rough / delta), 2.0)));
    /* strip inductive quality factor */
    Q_c = (pi * Z0_h_1 * w * f) / (R_s * C0 * K);
    alpha_c = (20.0 * pi / log(10.0)) * f * sqrt(e_r_eff_0) / (C0 * Q_c);
  } else {
    alpha_c = 0.0;
  }

  return alpha_c;
}


/*
 * dielectric_losses() - compute microstrip dielectric losses per unit
 * length
 */
double microstrip::dielectric_losses()
{
  double e_r, e_r_eff_0;
  double alpha_d;

  e_r = er;
  e_r_eff_0 = er_eff_0;

  alpha_d = (20.0 * pi / log(10.0)) * (f / C0) * (e_r / sqrt(e_r_eff_0)) * ((e_r_eff_0 - 1.0) / (e_r - 1.0)) * tand;

  return alpha_d;
}


/* 
 * attenuation() - compute attenuation of microstrip
 */
void microstrip::attenuation()
{
  skindepth = skin_depth();

  atten_cond = conductor_losses() * l;
  atten_dielectric = dielectric_losses() * l;
}


/*
 * mur_eff_ms() - returns effective magnetic permeability
 */
void microstrip::mur_eff_ms()
{
  double mureff;

  mureff = (2.0 * mur) / ((1.0 + mur) + ((1.0 - mur) * pow((1.0 + (10.0 * h / w)), -0.5)));

  mur_eff =  mureff;
}


/*
 * synth_width - calculate width given Z0 and e_r
 */
double microstrip::synth_width()
{
  double e_r, a, b;
  double w_h, w;


  e_r = er;


  a = ((Z0 / ZF0 / 2 / pi) * sqrt((e_r + 1) / 2.)) + ((e_r - 1) / (e_r + 1) * (0.23 + (0.11 / e_r)));
  b = ZF0 / 2 * pi / (Z0 * sqrt(e_r));

  if (a > 1.52) {
    w_h = 8 * exp(a) / (exp(2. * a) - 2);
  } else {
    w_h = (2. / pi) * (b - 1. - log((2 * b) - 1.) + ((e_r - 1) / (2 * e_r)) * (log(b - 1.) + 0.39 - 0.61 / e_r));
  }

  if (h > 0.0) {
    w = w_h * h;
    return w;
  } else {
    w = 0;
  }
  return w;
}


/*
 * line_angle() - calculate microstrip length in radians
 */
void microstrip::line_angle()
{
  double e_r_eff;
  double v, lambda_g;

  e_r_eff = er_eff;

  /* velocity */
  v = C0 / sqrt(e_r_eff * mur_eff);
  /* wavelength */
  lambda_g = v / f;
  /* electrical angles */
  ang_l = 2.0 * pi * l / lambda_g;	/* in radians */
}


void microstrip::calc()
{
  /* effective permeability */
  mur_eff_ms();
  /* static impedance */
  microstrip_Z0();
  /* calculate freq dependence of er and Z0 */
  dispersion();
  /* calculate electrical lengths */
  line_angle();
  /* calculate losses */
  attenuation();
}


/*
 * get_microstrip_sub () - get and assign microstrip substrate
 * parameters into microstrip structure
 */
void microstrip::get_microstrip_sub()
{
  er = getProperty ("Er");
  mur = getProperty ("Mur");
  h = getProperty ("H", UNIT_LENGTH, LENGTH_M);
  ht = getProperty ("H_t", UNIT_LENGTH, LENGTH_M);
  t = getProperty ("T", UNIT_LENGTH, LENGTH_M);
  sigma = getProperty ("Cond");
  tand = getProperty ("Tand");
  rough = getProperty ("Rough", UNIT_LENGTH, LENGTH_M);
}

/*
 * get_microstrip_comp() - get and assign microstrip component
 * parameters into microstrip structure
 */
void microstrip::get_microstrip_comp()
{
  f = getProperty ("Freq", UNIT_FREQ, FREQ_HZ);
}

/*
 * get_microstrip_elec() - get and assign microstrip electrical
 * parameters into microstrip structure
 */
void microstrip::get_microstrip_elec()
{
  Z0 = getProperty ("Z0", UNIT_RES, RES_OHM);
  ang_l = getProperty ("Ang_l", UNIT_ANG, ANG_RAD);
}


/*
 * get_microstrip_phys() - get and assign microstrip physical
 * parameters into microstrip structure
 */
void microstrip::get_microstrip_phys()
{
  w = getProperty ("W", UNIT_LENGTH, LENGTH_M);
  l = getProperty ("L", UNIT_LENGTH, LENGTH_M);
}


void microstrip::show_results()
{
  setProperty ("Z0", Z0, UNIT_RES, RES_OHM);
  setProperty ("Ang_l", ang_l, UNIT_ANG, ANG_RAD);

  setResult (0, er_eff, "");
  setResult (1, atten_cond, "dB");
  setResult (2, atten_dielectric, "dB");

  double val = convertProperty ("T", skindepth, UNIT_LENGTH, LENGTH_M);
  setResult (3, val, getUnit ("T"));
}

/*
 * analysis function
 */
void microstrip::analyze()
{
  /* Get and assign substrate parameters */
  get_microstrip_sub();

  /* Get and assign component parameters */
  get_microstrip_comp();

  /* Get and assign physical parameters */
  get_microstrip_phys();

  /* compute microstrip parameters */
  calc();

  /* print results in the subwindow */
  show_results();
}


/*
 * synthesis function
 */
int microstrip::synthesize()
{
  double Z0_dest, Z0_current, Z0_result, increment, slope, error;
  int iteration;
  const int maxiter = 100;

  /* Get and assign substrate parameters */
  get_microstrip_sub();

  /* Get and assign component parameters */
  get_microstrip_comp();

  /* Get and assign electrical parameters */
  get_microstrip_elec();

  /* Get and assign physical parameters */
  /* at present it is required only for getting strips length */
  get_microstrip_phys();


  /* calculate width and use for initial value in Newton's method */
  w = synth_width();

  /* required value of Z0 */
  Z0_dest = Z0;

  /* Newton's method */
  iteration = 0;

  /* compute microstrip parameters */
  calc();
  Z0_current = Z0;

  error = std::abs(Z0_dest - Z0_current);

  while (error > MAX_ERROR) {
    iteration++;
    increment = (w / 100.0);
    w += increment;
    /* compute microstrip parameters */
    calc();
    Z0_result = Z0;
    /* f(w(n)) = Z0 - Z0(w(n)) */
    /* f'(w(n)) = -f'(Z0(w(n))) */
    /* f'(Z0(w(n))) = (Z0(w(n)) - Z0(w(n+delw))/delw */
    /* w(n+1) = w(n) - f(w(n))/f'(w(n)) */
    slope = (Z0_result - Z0_current) / increment;
    /* printf("%g\n",slope); */
    w += (Z0_dest - Z0_current) / slope - increment;
    /*      printf("ms->w = %g\n", ms->w); */
    /* find new error */
    /* compute microstrip parameters */
    calc();
    Z0_current = Z0;
    error = std::abs(Z0_dest - Z0_current);
    /*      printf("Iteration = %d\n",iteration);
	    printf("w = %g\t Z0 = %g\n",ms->w, Z0_current); */
    if (iteration > maxiter)
      break;
  }

  setProperty ("W", w, UNIT_LENGTH, LENGTH_M);
  /* calculate physical length */
  ang_l = getProperty ("Ang_l", UNIT_ANG, ANG_RAD);
  l = C0 / f / sqrt(er_eff * mur_eff) * ang_l / 2.0 / pi;    /* in m */
  setProperty ("L", l, UNIT_LENGTH, LENGTH_M);

  /* compute microstrip parameters */
  calc();

  /* print results in the subwindow */
  show_results();

  if (iteration > maxiter)
    return -1;
  else
    return 0;
}