#!/usr/bin/python3

import numpy as np
import scipy.optimize as opt

np.seterr(under='ignore')
### STATIC  approx
def eps_re(u, e_r):
    '''
    Calculate effective dielectric constant
    using the expressions given in
    E. Hammerstad and O. Jensen, "Accurate Models for
    Microstrip Computer-Aided Design," 1980 IEEE MTT-S
    International Microwave symposium Digest, Washington,
    DC, USA, 1980, pp. 407-409.
    doi: 10.1109/MWSYM.1980.1124303

    e_re = eps_re(u, e_r):
    u <-  ratio W/h
    e_r <- substrate relative dielectric constant
    e_re -> effective relative dielectric constant
    of microstrip line
    '''
    a = 1 + np.log((u**4 + (u/52)**2)/(u**4+0.432))/49 +\
            np.log(1+(u/18.1)**3)/18.7
    b = 0.564*( (e_r-0.9)/(e_r+3) )**(0.053)

    e_re = 0.5*(e_r+1)+0.5*(e_r-1)*(1+10/u)**(-a*b)
    return e_re

def getZc(u, e_r):
    '''
    Calculate characteristic impedance Zc 
    using the expressions given in
    E. Hammerstad and O. Jensen, "Accurate Models for
    Microstrip Computer-Aided Design," 1980 IEEE MTT-S
    International Microwave symposium Digest, Washington,
    DC, USA, 1980, pp. 407-409.
    doi: 10.1109/MWSYM.1980.1124303
    (neglected thickness of copper)

    Zc = getZc(u, e_r)
    u <- ratio W/h
    e_r <- substrate relative dielectric constant
    Zc -> characteristic impedance of microstrip line
    '''
    F = 6 + (2*np.pi-6)*np.exp(-(30.666/u)**0.7528)
    e_re = eps_re(u, e_r)
    Zc = (120*np.pi)/(2*np.pi*np.sqrt(e_re))*np.log(\
            F/u + np.sqrt(1 + (2/u)**2))
    return Zc

### Dispersion correction
def eps_f(W,h,e_r,f):
    '''
    Calculate effective dieletric constant using formulas given by:
    M. Kirschning and R. H. Jansen, "Accurate Model for Effective
    Dielectric Constant of Microstrip with Validity up to Millimeter-Wave
    Frequencies," Electronics Letters, vol. 8, no. 6, pp. 272-273,
    Mar. 1982. 
    R. H. Jansen and M. Kirschning, "Arguments and an accurate Model
    for the Power-Current Formulation of Microstrip Characteristic
    Impedance," Archiv für Elektronik und Übertragungstechnik (AEÜ),
    vol. 37, pp. 108-112, 1983.
    http://qucs.sourceforge.net/tech/node75.html

    e_re = eps_f(W,h,e_r,f)
    W <- width of microstrip (mm)
    h <- height of substrate (mm)
    e_r <- substrate relative dielectric cosntatn
    f <- frequency (GHz)
    e_re -> effective dielectric constant of microstrip line
    '''
    fn = f*h
    u = W/h
    P1 = 0.27488 + (0.6315+0.525/(1+0.0157*fn)**20)*u - 0.065683*np.exp(\
            -8.7513*u)
    P2 = 0.33622*(1-np.exp(-0.03442*e_r))
    P3 = 0.0363*np.exp(-4.6*u)*(1-np.exp(-(fn/38.7)**4.97))
    P4 = 1+2.751*(1-np.exp(-(e_r/15.916)**8))
    Pf = P1*P2*((0.1844+P3*P4)*fn)**1.5763
    eps0 = eps_re(u,e_r)
    return e_r - (e_r - eps0)/(1+Pf)
    
def getZcf(W,h,e_r,f,R=False):
    '''
    Calculate characteristic impedance using formulas given by:
    M. Kirschning and R. H. Jansen, "Accurate Model for Effective
    Dielectric Constant of Microstrip with Validity up to Millimeter-Wave
    Frequencies," Electronics Letters, vol. 8, no. 6, pp. 272-273,
    Mar. 1982. 
    R. H. Jansen and M. Kirschning, "Arguments and an accurate Model
    for the Power-Current Formulation of Microstrip Characteristic
    Impedance," Archiv für Elektronik und Übertragungstechnik (AEÜ),
    vol. 37, pp. 108-112, 1983.
    http://qucs.sourceforge.net/tech/node75.html

    Zc = Zcf(W,h,e_r,f,R=False)
    W <- width of microstrip (mm)
    h <- height of substrate (mm)
    e_r <- substrate relative dielectric cosntatn
    f <- frequency (GHz)
    R=False <- flag telling whether to return R17 or not
    Zc -> characteristic impedance of microstrip line
    R17 -> parameter used to calculate impedance and coupled line
    return only if R=True
    '''
    fn = f*h
    u = W/h
    ef = eps_f(W,h,e_r,f)
    e0 = eps_re(u,e_r)

    R1 = 0.03891*e_r**1.4
    R2 = 0.267*u**7
    R3 = 4.766*np.exp(-3.228*u**0.641)
    R4 = 0.016 + (0.0514*e_r)**4.524
    R5 = (fn/28.843)**12.0
    R6 = 22.20*u**1.92

    R7 = 1.206-0.3144*np.exp(-R1)*(1-np.exp(-R2))
    R8 = 1+1.275*(1-np.exp(-0.004625*R3*e_r**1.674)*(fn/18.365)**2.745)
    R9 = 5.086*R4*R5*np.exp(-R6)*(e_r-1)**6/\
            ((0.3838+0.386*R4)*(1+1.2992*R5)*(1+10*(e_r-1)**6))

    R10 = 0.00044*e_r**2.136+0.0184
    R11 = (fn/19.47)**6/(1+0.0962*(fn/19.47)**6)
    R12 = 1/(1+0.00245*u**2)
    R13 = 0.9408*ef**R8-0.9603
    R14 = (0.9408-R9)*e0**R8-0.9603
    R15 = 0.707*R10*(fn/12.3)**1.097
    R16 = 1+0.0503*e_r**2*R11*(1-np.exp(-(u/15)**6))
    R17 = R7*(1-1.1241*R12*np.exp(-0.026*fn**1.15656-R15)/R16)
    
    Z0 = getZc(u,e_r)
    Zc = Z0*(R13/R14)**R17
    if R:
        return (Zc, R17)
    else:
        return Zc

def getWf(Zk,h,e_r,f):
    '''
    Calculate width of microstrip line when also
    taking in to account dispersion.

    W = getWf(Zk,h,e_r,f)
    Zk <- characteristic impedance (ohms)
    h <- substrate height
    e_r <- substrate relative di. costant
    f <- frequency (GHz)
    W -> widht of microstrip line (mm)
    '''
    u = opt.root_scalar(lambda X: Zk-getZcf(X,h,e_r,f),x0=h,\
            method='brentq', bracket=[0.01,100])
    return u.root

### coupled lines
def mcl_eps(W,h,s,e_r,f):
    '''
    Calculate effective odd and even permittivity of parallel
    coupled microstrip lines.
    M. Kirschning, R. H. Jansen, and N. H. L. Koster,
    "Coupled microstrip parallel-gap model for improved filter and coupler design,"
    Electronics Letters, vol. 19, no. 10, pp. 377–379, May 1983,
    doi: 10.1049/el:19830261.

    even, odd = mcl_eps(W,h,s,e_r,f):
    W <- width of microstrip (mm)
    h <- height of substrate (mm)
    s <- spacing between lines (mm)
    e_r <- substrate relative dielectric cosntatn
    f <- frequency (GHz)
    even -> even mode dielectric constant
    odd -> odd mode dielectric constant
    '''
    u = W/h
    g = s/h
    fn = f*h

    eps0 = eps_re(u,e_r)

    v = u * (20+g**2)/(10+g**2)+ g*np.exp(-g)
    ae = 1+np.log((v**4+(v/52)**2)/(v**4+0.432))/49+np.log(1+(v/18.1)**3)/18.7
    be = 0.564*((e_r-0.9)/(e_r+3))**0.053
    even0 = 0.5*(e_r+1)+0.5*(e_r-1)*(1+10/v)**(-ae*be)

    ao = 0.7287*(eps0-0.5*(e_r+1))*(1-np.exp(-0.179*u))
    bo = 0.747*e_r/(0.15+e_r)
    co = bo - (bo-0.207)*np.exp(-0.414*u)
    do = 0.593+0.694*np.exp(-0.562*u)
    odd0 = (0.5*(e_r+1)+ao-eps0)*np.exp(-co*g**do)+eps0

    P1 = 0.27488+(0.6315+0.525/(1+0.0157*fn)**20)*u-0.065683*np.exp(-8.7513*u)
    P2 = 0.33622*(1-np.exp(-0.03442*e_r))
    P3 = 0.0363*np.exp(-4.6*u)*(1-np.exp(-(fn/38.7)**4.97))
    P4 = 1+2.751*(1-np.exp(-(e_r/15.916)**8))
    P5 = 0.334*np.exp(-3.3*(e_r/15)**3)+0.746
    P6 = P5*np.exp(-(fn/18)**0.368)
    P7 = 1+4.069*P6*g**0.479*np.exp(-1.347*g**0.595-0.17*g**2.5)
    Fe = P1*P2*((P3*P4+0.1844*P7)*fn)**1.5763
    even = e_r - (e_r - even0)/(1+Fe)

    P8 = 0.7168*(1+1.076/(1+0.0576*(e_r-1)))
    P9 = P8-0.7913*(1-np.exp(-(fn/20)**1.424))*np.arctan(2.481*(e_r/8)**0.946)
    P10 = 0.242*(e_r - 1)**0.55
    P11 = 0.6366*(np.exp(-0.3401*fn)-1)*np.arctan(1.263*(u/3)**1.629)
    P12 = P9 + (1-P9)/(1+1.183*u**1.376)
    P13 = 1.695 * P10/(0.414+1.605*P10)
    P14 = 0.8928+0.1072*(1-np.exp(-0.42*(fn/20)**3.215))
    P15 = np.abs(1-0.8928*(1+P11)*np.exp(-P13*g**1.092)*P12/P14)
    Fo = P1*P2*((P3*P4+1.844)*fn*P15)**1.5763
    odd = e_r - (e_r - odd0)/(1+Fo)

    return even,odd,even0,odd0

def mclZcf(W,h,s,e_r,f):
    '''
    Calculate odd and even characteristic impedance of parallel
    coupled microstrip lines.
    M. Kirschning, R. H. Jansen, and N. H. L. Koster,
    "Coupled microstrip parallel-gap model for improved filter and coupler design,"
    Electronics Letters, vol. 19, no. 10, pp. 377–379, May 1983,
    doi: 10.1049/el:19830261.

    even, odd = mcl_eps(W,h,s,e_r,f):
    W <- width of microstrip (mm)
    h <- height of substrate (mm)
    s <- spacing between lines (mm)
    e_r <- substrate relative dielectric cosntatn
    f <- frequency (GHz)
    Ze -> even mode charachterstic impedance
    Zo -> odd mode characteristic impedance
    '''
    u = W/h
    g = s/h
    fn = f*h

    eps0 = eps_re(u,e_r)
    epsf = eps_f(W,h,e_r,f)
    even,odd,even0,odd0 = mcl_eps(W,h,s,e_r,f)
    Z0 = getZc(u,e_r)
    Z0f, Q0 = getZcf(W,h,e_r,f,R=True)

    Q1 = 0.8695*u**0.194
    Q2 = 1+0.7519*g+0.189*g**2.31
    Q3 = 0.1975+(16.6 + (8.4/g)**6)**-0.387 +np.log(g**10/(1+(g/3.4)**10))/241
    Q4 = 2*Q1/Q2/(np.exp(-g)*u**Q3+(2-np.exp(-g))*u**-Q3)

    Ze0 = np.sqrt(eps0/even0)*Z0/(1-Z0/377*np.sqrt(eps0)*Q4)

    Q5 = 1.794 + 1.14*np.log(1+0.638/(g+0.517*g**2.43))
    Q6 = 0.2305+np.log(g**10/(1+(g/5.8)**10))/281.3+np.log(1+0.598*g**1.154)/5.1
    Q7 = (10+190*g**2)/(1+82.3*g**3)
    #print(-6.5-0.95*np.log(g)-(g/0.15)**5)
    Q8 = np.exp(-6.5-0.95*np.log(g)-(g/0.15)**5)
    Q9 = np.log(Q7)*(Q8+1/16.5)
    Q10 = (Q2*Q4-Q5*np.exp(np.log(u)*Q6*u**-Q9))/Q2
    
    Zo0 = np.sqrt(eps0/odd0)*Z0/(1-Z0/377*np.sqrt(eps0)*Q10)

    # dispersion accounted
    Q11 = 0.893*(1-0.3/(1+0.7*(e_r-1)))
    Q12 = 2.121*(fn/20)**4.91*np.exp(-2.87*g)*g**0.902/(1+Q11*(fn/20)**4.91)
    Q13 = 1+0.038*(e_r/8)**5.1
    Q14 = 1+1.203*(e_r/15)**4/(1+(e_r/15)**4)
    Q15 = 1.887*np.exp(-1.5*g**0.84)*g**Q14/(1+0.41*(fn/15)**3*u**(2/Q13)/(0.125+u**(1.625/Q13)))
    Q16 = Q15 * (1+ 9/(1+0.403*(e_r-1)**2))
    Q17 = 0.394*(1-np.exp(-1.47*(u/7)**0.672))*(1-np.exp(-4.25*(fn/20)**1.87))
    Q18 = 0.61*(1-np.exp(-2.13*(u/8)**1.593))/(1+6.544*g**4.17)
    Q19 = 0.21*g**4/(1+0.18*g**4.9)/(1+0.1*u**2)/(1+(fn/24)**3)
    Q20 = Q19 * (0.09+1/(1+0.1*(e_r-1)**2.7))
    Q21 = np.abs(1-42.54*g**0.133*np.exp(-0.812*g)*u**2.5/(1+0.033*u**2.5))

    re = (fn/28.843)**12
    qe = 0.016+(0.0514*e_r*Q21)**4.524
    pe = 4.766*np.exp(-3.228*u**0.641)
    de = 5.086*qe*re*np.exp(-22.2*u**1.92)*(e_r-1)**6/(0.3838+0.386*qe)/(1+1.2992*re)/(1+10*(e_r-1)**6)

    Ce = 1+1.275*(1-np.exp(-0.004625*pe*e_r**1.674*(fn/18.365)**2.745))-Q12+Q16-Q17+Q18+Q20
    Ze = Ze0*( (0.9408*epsf**Ce-0.9603)/( (0.9408-de)*eps0**Ce-0.9603))**Q0

    Q29 = 15.16/(1+0.196*(e_r-1)**2)
    Q28 = 0.149*(e_r-1)**3/(94.5+0.038*(e_r-1)**3)
    Q27 = 0.4*g**0.84*(1+(2.5*(e_r-1)**1.5)/(5+(e_r-1)**1.5))
    Q26 = 30 - ( 22.2*((e_r-1)/13)**13)/(1+3*((e_r-1)/13)**12) - Q29
    Q25 = (0.3*fn**2)/(10+fn**2)*(1+(2.333*(e_r-1)**2)/(5+(e_r-1)**2))
    Q24 = 2.506*Q28*u**0.894/(3.575+u**0.894)*((1+1.3*u)*fn/99.25)**4.29
    Q23 = 1+0.005*fn*Q27/(1+0.812*(fn/15)**1.9)/(1+0.025*u**2)
    Q22 = 0.925*(fn/Q26)**1.536/(1+0.3*(fn/30)**1.536)

    Zo = Z0f+(Zo0*(odd/odd0)**Q22 - Z0f*Q23)/(1+Q24+(0.46*g)**2.2*Q25)
    return Ze,Zo

def get_mcl(Ze, Zo, h, e_r, f):
    '''
    Find width W and gap s of parallel coupled microstrip lines

    [W, s] = get_mcl(Ze, Zo, h, e_r, f):
    Ze <- even mode characteristic impedance (ohm)
    Zo <- odd mode ch. impedance (ohm)
    h <- height of substrate (mm)
    e_r <- relative permitivity of substrate
    f <- frequency (GHz)
    W -> strip width (mm)
    s -> gap (mm)
    '''
    Z = np.array([Ze,Zo])
    u = opt.root(lambda X: Z - np.array(mclZcf(X[0],h,X[1],e_r,f)), [h/2,h/5],method='lm')
    return u.x

def open_end(W, h, e_r, f):
    '''
    Calculate open end displacment dl of microstrip open end.
    Kirschning, M., R. H. Jansen, and N. H. L. Koster.
    "Accurate model for open end effect of microstrip lines."
    Electronics Letters 17.3 (1981): 123-125.

    dl = open_end(W,h,e_r,f)
    W <- width of microstrip line (mm)
    h <- height of substrate (mm)
    e_r <- relative dielectric constant of substrate
    f <- frequency (GHz)
    '''
    epsf = eps_f(W,h,e_r,f)
    u = W/h

    x1 = 0.434907*(epsf**0.81+0.26)/(epsf**0.81-0.189)*(u**0.8544+0.236)/(u**0.8544+0.87)
    x2 = 1 + u**0.371/(2.358*e_r+1)
    x3 = 1 + 0.5274*np.arctan(0.084*u**(1.9413/x2))/epsf**0.9236
    x4 = 1 + 0.0377*np.arctan(0.067*u**1.456)*(6-5*np.exp(0.036*(1-e_r)))
    x5 = 1 - 0.218*np.exp(-7.5*u)

    return x1*x3*x5/x4 *h