/* -*- C++ -*-
 
  This file is part of ViPEC
  Copyright (C) 1991-2001 Johan Rossouw (jrossouw@alcatel.altech.co.za)
 
  This program is free software; you can redistribute it and/or modify
  it under the terms of the GNU Library 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 Library General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 
*/

#include <TwoPorts.h>
#include <Setup.h>
#include <Matrix.h>
#include <Exception.h>

#include <math.h>
#include <iostream>

/* a few static constant variables used by the models */
#ifndef _WS_WIN_
static const TReal ZERO=1E-1000;
static const TReal BIGNUM=1E1000;
#else
static const TReal ZERO=1E-200;
static const TReal BIGNUM=1E200;
#endif

using namespace std;

static const TReal C = 2.997925E8; 	// Speed of light (m/s)
static const TReal muo =4*M_PI*1E-7;	//permeabilitty of vacuum (V.s.A^{-1}.m^{-1})
static const TReal KapaC = 58E6;	// conductivity of cooper  (Ohm^{-1}.m^{-1})

//---------------------------------------------------------------------------
void TwoPorts::getTLIN(TReal Freq, TReal Z, TReal E, TReal Fo, Matrix& Y)
{

  TReal Alpha = 0.00001;  //LOSS FACTOR

  TReal Beta = (M_PI*E/180) * (Freq/Fo);
  TComplex Gamma = TComplex(Alpha, Beta);


  TComplex Yi = cosh(Gamma) / (sinh(Gamma) * TComplex(Z,0));
  TComplex Yt = TComplex(1,0) / (sinh(Gamma) * TComplex(-Z,0));

  Y.resize(2);
  Y(0,0) = Yi;
  Y(0,1) = Yt;
  Y(1,0) = Yt;
  Y(1,1) = Yi;

}

//---------------------------------------------------------------------------
void TwoPorts::getTLSC(TReal Freq, TReal Z, TReal E, TReal Fo, TComplex &Y)
{
  TReal Alpha = 0.00001;  //LOSS FACTOR

  TReal Beta = (M_PI*E/180) * (Freq/Fo);
  TComplex Gamma = TComplex(Alpha, Beta);

  Y = cosh(Gamma) / (sinh(Gamma) * TComplex(Z,0));

}

//---------------------------------------------------------------------------
void TwoPorts::getTLOC(TReal Freq, TReal Z, TReal E, TReal Fo, TComplex &Y)
{
  TReal Alpha = 0.00001;  //LOSS FACTOR

  TReal Beta = (M_PI*E/180) * (Freq/Fo);
  TComplex Gamma = TComplex(Alpha, Beta);

  Y = sinh(Gamma) / (TComplex(Z,0) * cosh(Gamma));

}

//---------------------------------------------------------------------------
void TwoPorts::getCLIN(TReal Freq, TReal Ze, TReal Zo, TReal E, TReal Fo, Matrix& Y)
{
  TReal  d = cos(Freq/Fo);
  TReal Ye = 1/Ze;
  TReal Yo = 1/Zo;

  /* calculate beta - actually this is beta*length */
  TReal Beta = (M_PI*E/180) * (Freq/Fo);

  TComplex s,m;

  if (fabs(d)<ZERO)
    {
      if (d<0)
        s = TComplex(0,-BIGNUM);
      else
        s = TComplex(0,+BIGNUM);

      m = TComplex(0,-1);
    }
  else
    {
      s = TComplex(0,tan(Beta));
      m = sqrt(TComplex(1,0) - pow(s,2)) / s;
    }

  Y.resize(2);
  Y(0,0) = TComplex((Yo+Ye),0) / (TComplex(-2,0)*s);
  Y(0,1) = TComplex((Yo-Ye),0) / (TComplex(-2,0)*s);
  Y(1,0) = m*TComplex((Yo-Ye),0) / (TComplex(-2,0));
  Y(1,1) = m*TComplex((Yo+Ye),0) / (TComplex(-2,0));
}

//---------------------------------------------------------------------------

/* NOTE: - need to take effect of finite strip thickness into accout
 		 - losses still to be added */

void TwoPorts::getSCLIN(TReal Freq,
                        TReal W,  //Line width
                        TReal S,  //Line spacing
                        TReal L,  //Line length
                        TReal H,  //Substrate height
                        TReal Er, //Relative permittivity
                        Matrix& Y)
{
  /* k and kprime are functions of the physical dimensions */
  TReal ke = tanh(M_PI*W/(2*H))*tanh((M_PI/2)*(W+S)/H);
  TReal keP = sqrt(1 - pow(ke,2));
  TReal ko = tanh(M_PI*W/(2*H))/tanh((M_PI/2)*(W+S)/H);
  TReal koP = sqrt(1 - pow(ko,2));

  /* even mode impedance */
  TReal Zoe = 1/sqrt(Er);

  if (ke>=ZERO && ke<=0.7)
    {
      keP = log(2*((1+sqrt(keP))/(1-sqrt(keP))))/M_PI;
      Zoe = Zoe*30.0*M_PI*keP;
    }
  else if (ke>0.7 && ke<=1)
    {
      ke = log(2*(1+sqrt(ke))/(1-sqrt(ke)))/M_PI;
      Zoe = Zoe*30.0*M_PI/ke;
    }
  else
    {
      throw Exception::LineDimensionsOutOfRange();
    }

  /* odd mode impedance */
  TReal Zoo = 1/sqrt(Er);

  if (ko>=ZERO && ko<=0.7)
    {
      koP = log(2*((1+sqrt(koP))/(1-sqrt(koP))))/M_PI;
      Zoo = Zoo*30.0*M_PI*koP;
    }
  else if (ko>0.7 && ko<=1)
    {
      ko = log(2*(1+sqrt(ko))/(1-sqrt(ko)))/M_PI;
      Zoo = Zoo*30.0*M_PI/ko;
    }
  else
    {
      throw Exception::LineDimensionsOutOfRange();
    }

  /* Phase length at this freq computed by first getting beta=2*pi/wavelength (rad)
     with the freq in Hz!!  - electrical length in rad just by multiplying with L */
  TReal Beta = Freq*sqrt(Er)/(C);  //2*pi above and below line cancelled

  /* losses contained in alpha */
  //TReal Alpha = (27.3/2*M_PI)*Freq*tanRho*sqrt(Er)/C;

  //TComplex Gamma(Alpha, Beta);

  /* S = j*tan(beta*l) */
  TComplex s = TComplex(0,tan(Beta*L));
  //TComplex s = TComplex(0,tan(Gamma*L));
  TComplex m = sqrt(TComplex(1,0) - pow(s,2)) / s;

  TReal Ye = 1/Zoe;
  TReal Yo = 1/Zoo;

  Y.resize(2);
  Y(0,0) = TComplex((Yo+Ye),0) / (TComplex(-2,0)*s);
  Y(0,1) = TComplex((Yo-Ye),0) / (TComplex(-2,0)*s);
  Y(1,0) = m*TComplex((Yo-Ye),0) / (TComplex(-2,0));
  Y(1,1) = m*TComplex((Yo+Ye),0) / (TComplex(-2,0));

}

//---------------------------------------------------------------------------
void TwoPorts::getGYR(TReal R, Matrix& Y)
{
  Y.resize(2);
  Y(0,0) = TComplex(0,0);
  Y(0,1) = TComplex(1/R,0);
  Y(1,0) = TComplex(-1/R,0);
  Y(1,1) = TComplex(0,0);
}

//---------------------------------------------------------------------------
void TwoPorts::getTransformer(TReal f, TReal L1, TReal L2, TReal k, Matrix& Y)
{

  Y.resize(2);

  // small case, to restrict the users possibility of getting a singular matrix
  if  (k>=1)
    {
      k = 0.99999;
    }
  else if (k<0)
    {
      k = 0;
    }
  TReal L12 = sqrt(L1*L2);  // Maximal mutual  inductance
  TReal k2m1 = k*k-1;       //  k=M/L12, where M is the diagonal term of the impedance matrix
  Y(0,0) = TComplex(0, 1/(f*L1*k2m1));
  Y(0,1) = TComplex(0, -k/(f*L12*k2m1));
  Y(1,0) = Y(0,1);
  Y(1,1) = TComplex(0, 1/(f*L2*k2m1));
}

//---------------------------------------------------------------------------
void TwoPorts::getBIPB(TReal Freq, TReal B, TReal Fo, TReal Re, TReal C, Matrix& Y)
{
  // PARAMETER SEQUENCY B=beta F=frequency Re=resistance C=capac

  TComplex beta = TComplex(B,0) / TComplex(1,B*Freq/Fo);
  TComplex zr = (TComplex(1,0) + beta) * TComplex(Re,0);
  TComplex yc = TComplex(0,Freq*C);
  Y.resize(2);
  Y(0,0) = (TComplex(1,0) / zr) + yc;
  Y(1,0) = (beta / zr) - yc;
  Y(0,1) = TComplex(-1,0) * yc;
  Y(1,1) = yc;
}

//---------------------------------------------------------------------------
void TwoPorts::getSLIN(TReal Freq,
                       TReal W,  //Line width
                       TReal L,  //Line length
                       TReal T,  //Conductor thickness
                       TReal H,  //Substrate height
                       TReal Er, //Relative permittivity
                       TReal tanRho, //Loss
                       Matrix& Y)
{
  /* work with normalized Z for the rest of function */
  TReal Z = 1.0/sqrt(Er);	//normalized Z

  /* Phase length at this freq computed by first getting beta=2*pi/wavelength (rad)
     with the freq in Hz!!  - electrical length in rad just by multiplying with L */
  TReal Beta = Freq*sqrt(Er)/(C);  //2*pi above and below line cancelled out

  /* k and kprime are functions of W and height */
  TReal k = tanh(M_PI*W/(2*H));
  TReal kPrime = sqrt(1-pow(k,2));

  if (T<=ZERO)
    {   //case for zero thickness conductor
      if (k>=ZERO && k<=0.7)
        {
          kPrime = log(2*((1+sqrt(kPrime))/(1-sqrt(kPrime))))/M_PI;
          Z = Z*30.0*M_PI*kPrime;
        }
      else if (k>0.7 && k<=1)
        {
          k = log(2*(1+sqrt(k))/(1-sqrt(k)))/M_PI;
          Z = Z*30.0*M_PI/k;
        }
      else
        {
          throw Exception::LineDimensionsOutOfRange();
        }
    }//if T<=ZERO
  else
    {  //case for finite thickness conductor
      TReal x = T/H;
      TReal m = 2/(1+(2/3)*x/(1-x));
      TReal Wnorm = (x/M_PI*(1-x))*(1-0.5*log(pow(x/(2-x),2)+pow(0.0796*x/(W/H+1.1*x),m)));
      Wnorm = 1/(Wnorm + W/(H-T));
      Z = Z*30*log(1+(Wnorm*4/M_PI)*(Wnorm*8/M_PI+sqrt(pow(Wnorm*8/M_PI,2)+6.27)));
    }//else

  /* alpha = pi*tan(rho)/wavelength - note that freq is radians in ALL functions */
  //TReal Alpha = 0.5*Freq*tanRho*sqrt(Er)/C;

  /* another dielectric loss tangent is alpha = 27.3*sqrt(Er)*tanRho/(natural wavel)
     NOTE this one gives better results!! */
  TReal Alpha = (27.3/2*M_PI)*Freq*tanRho*sqrt(Er)/C;

  /* before calculating input admittance, multiply alpha, beta with length to get angle */
  Beta *= L;
  Alpha *= L;

  /* Now get gamma */
  TComplex Gamma = TComplex(Alpha, Beta);

  TComplex Yi = cosh(Gamma) / (sinh(Gamma) * TComplex(Z,0));
  TComplex Yt = TComplex(1,0) / (sinh(Gamma) * TComplex(-Z,0));

  Y.resize(2);
  Y(0,0) = Yi;
  Y(0,1) = Yt;
  Y(1,0) = Yt;
  Y(1,1) = Yi;

}//end of getSLIN()


//---------------------------------------------------------------------------
void TwoPorts::getMLIN(TReal Freq,
                       TReal W,  //Line width
                       TReal L,  //Line length
                       TReal T,  //Line thickness
                       TReal H,  //Substrate height
                       TReal Er, //Relative permittivity
                       TReal tanRho, //Loss
                       Matrix& Y)
{
  const TReal ZETA = 120*M_PI;	/* units of OHM */
  TReal Z=0;	/* variable used for line impedance */
  TReal Ere=0;	/* variable for the effective relabive permm */

  /* We take strip thickness into account... */
  if (T/H <=0.005)
    {		//small enough to assume T==0

      /* because the mode is not purely TEM but quasi-TEM the effective dielectric
      constant is lower than that given for the substrate and takes into account
      the fields external to the substrate */
      Ere = 0.5*(Er+1) + 0.5*(Er-1)*pow((1+10*H/W),-0.5);

      /* Z computed differently for different ratios of W/H */
      if (W/H <= 1)
        {
          Z = ZETA/(2*M_PI*sqrt(Ere)) * log( (8*H/W)+0.25*W/H );
        }
      else if (W/H > 1)
        {
          Z = (ZETA/sqrt(Ere)) / (W/H + 1.393 + 0.667*log(W/H+1.444));
        }
      else
        {
          throw Exception::LineDimensionsOutOfRange();
        }
    }//if T/H<=0.005
  else
    {
      Ere = 0.5*(Er+1) + 0.5*(Er-1)/sqrt(1+10*H/W) - (pow(W/H,-0.5)*T/H)*(Er-1)/4.6;

      /* ratio needed to calculate effective width */
      TReal dWonH;
      if (W/H <= 1/(2*M_PI))
        {
          dWonH = (1.25/M_PI)*(T/H)*(1 + log(4*M_PI*W/T));
        }
      else
        {
          dWonH = (1.25/M_PI)*(T/H)*(1 + log(2*H/T));
        }

      /* effective Width */
      TReal We = W + dWonH*H;

      /* and now we can get Z - different for ratios of W/H */
      if (W/H <= 1)
        {
          Z = (60/sqrt(Ere)) * log( (8*H/We)+0.25*We/H );
        }
      else if (W/H > 1)
        {
          Z = (376.7/sqrt(Ere)) / (We/H + 1.393 + 0.667*log(We/H+1.444));
        }
      else
        {
          throw Exception::LineDimensionsOutOfRange();
        }
    }//if-else T/H

  /* Phase length at this freq computed by first getting beta=2*pi/wavelength (rad)
     with the freq in Hz!!  - electrical length in rad just by multiplying with L */
  TReal Beta = Freq*sqrt(Ere)/(C);  //2*pi above and below line cancelled out

  /* alpha - note that freq is radians in ALL functions, therefore the 2*pi
     NOTE the mixture of effective and normal relative permm */

  // 27/02/2001 - Bug in calculation of Alpha, the commented line was the original
  //TMyReal Alpha = (27.3/2*M_PI)*Er/(Er-1)*((Ere-1)/sqrt(Ere))*Freq*tanRho*sqrt(Ere)/C;
  TReal Alpha = 0.5*Er/(Er-1)*((Ere-1)/sqrt(Ere))*Freq*tanRho/C;

  /* before calculating input admittance, multiply alpha, beta with length to get angle */
  Beta *= L;
  Alpha *= L;

  /* Now get gamma */
  TComplex Gamma = TComplex(Alpha, Beta);

  TComplex Yi = cosh(Gamma) / (sinh(Gamma) * TComplex(Z,0));
  TComplex Yt = TComplex(1,0) / (sinh(Gamma) * TComplex(-Z,0));

  Y.resize(2);
  Y(0,0) = Yi;
  Y(0,1) = Yt;
  Y(1,0) = Yt;
  Y(1,1) = Yi;

}

//-----------------------------------------------------------------------------------------
void TwoPorts::getCoaxLIN(TReal Freq,
                          TReal W,  //Outer Diameter of the inner conductor
                          TReal L,  //Line length
                          TReal /*T*/,  //Thickness of the inner conductor for further calculations (not used for the monent)
                          TReal H,  //Inner Diameter of the outer conductor (shield)
                          TReal Er, //Relative permittivity beween the shield and the inner conductor
                          TReal tanRho, //Loss
                          Matrix& Y)
{
  /* the theory behind this model is developped in the Book of Brian C.Wandell : Transmission Line Design Hadbook   */
  /* elements of this model are presented in efg, internal report plx2.doc  991202  */

  const TReal Zo = 377;	/* units of OHM */
  TReal freal = Freq/(2*M_PI);  /*  note that freq is in radians in all functions, therefore we have to consider 2*pi each time*/

  /* Caracteristic impedance (2.6) */
  TReal Z = Zo/(2*M_PI*sqrt(Er))*log(H/W);

  /* Computing the losses along TL */


  /* Conductive losses metal conductor */
  TReal alphaL = 1/(2*H*Z)*pow((freal*muo/(M_PI*KapaC)),0.5)*(1+H/W);      // the conductor is not paramagnetic mur=1

  /*  Dielectric losses */
  TReal alphaD = M_PI*freal*(sqrt(Er)*(tanRho)/C);

  /*Total damping = conductor and dielectric (radiation damping is assumed negligible) */
  TReal Alpha = alphaL+alphaD;

  /*   We have defined Ere, as an electric permittivity with E'+jE'',  */
  /* In our range of materials Er is not linked to frequency              */
  /* and therefore the TL will resonate at equally spaced harmonics */
  /* If tanRho is not always negligible, Ere is not real */

  //TComplex Ere = TComplex(Er,Er*tanRho);
  TReal Ere = Er;

  /* in the full TEM case */
  /*The propagation constant become */

  //TComplex Beta=2*M_PI/C*freal*pow(Ere,0.5);
  TReal Beta=2*M_PI*freal*(pow(Ere,0.5)/C);

  /* before calculating input admittance, multiply alpha, beta with length to get angle */
  // Beta = Beta*TComplex(L,0);
  Beta *= L;
  Alpha *= L;

  /* Now get gamma */
  //TComplex Gamma = TComplex(Alpha,0)+Beta;
  TComplex Gamma = TComplex(Alpha,Beta);

  TComplex Yi = cosh(Gamma) / (sinh(Gamma) * TComplex(Z,0));
  TComplex Yt = TComplex(1,0) / (sinh(Gamma) * TComplex(-Z,0));

  Y.resize(2);
  Y(0,0) = Yi;
  Y(0,1) = Yt;
  Y(1,0) = Yt;
  Y(1,1) = Yi;
}