# coding: utf-8

import math, cmath

# der Wellenwiderstand einer verlustbehafteten Leitung, die Verluste näherungsweise nachgebildet
# durch den die Widerstandsverluste, Skineffekt s*f**0.5 und die dieelektrischen Verluste g*f
#
# zwave = sqrt((R' + j*2*pi*f*L')/(G' + j*2*pi*f*C'))
# mit L'    Induktivitätsbelag in H/m
# mit C'    Kapazitätsbelag in F/m
# mit R'    Widerstandsbelag in Ohm/m
# mit G'    dieelektrische Verluste in Siemens/m = 1/(Ohm *m)

# aus z0**2 = L'/C' und ceff**2 = 1/(L' * C') ergibt sich L' = z0/ceff und C' = 1/(ceff * z0)
# mit z0    ideller (verlustfreier) Wellenwiderstand = OffsetZ0 in der HP Nomenklatur
# mit ceff  effektive Lichtgeschwindigkeit, ceff = c0/sqrt(epsilonr*micronr) = c0 * vf, c0 = 2.997925e8 [m/s]
# ergibt sich:
#
# zwave = sqrt((R' + j*2pi*f*z0/ceff)/(G' + j*2*pi*f/(ceff*z0))) 

# durch Multiplikation mit ceff/ceff ergibt sich:
#
# zwave = sqrt((R' * ceff + j*2pi*f*z0)/(G' * ceff + j*2*pi*f/z0)) 
#
# mit
# r' = R' * ceff = R' * 0.2997925 * vf [GigaOhm/s]
# r' wird in der HP-Nomenclatur mit Offset loss angegeben und hat die Einheit GigaOhm/s
# will man den Skin-Effekt näherungsweise berücksichtigen ist auf r' noch der Term s*sqrt(f) aufzuaddieren, 
# mit s als eine Konstante, die sich aus Material und Formgebung ergibt, und mit der Einheit Ohm/sqrt(s).
# Der gesamte Ohmsche Anteil nähert sich jedoch der Wirklichkeit besser an, wenn man den quadratischen
# Mittelwert bildet: sqrt(r'**2 + (s*f**0.5)**2)
#
# und mit
# g * f = G' * ceff = G' * 2.997925e8 * vf [1/(Ohm*s)], g hat damit die Einheit 1/Ohm und ist die Geradensteigungung.
# Bei den SOL Standards wird weder s noch g verwendet, also s=0 und g=0. Es ergibt sich letztlich:
#
# zwave = sqrt((sqrt( r'**2.0 + (s * sqrt(f))**2.0 ) + j*2*pi*f*z0)/ (g*f + j*2*pi*f/z0))

def zwave(f, g, s, z0, r=0):
    return (complex(math.sqrt(r**2.0 + s**2.0 * f), 2*math.pi*f*z0)/complex(g*f, 2*math.pi*f/z0))**0.5
#end zwave



# der komplexe Ausbreitungskoeffizient 
# gamma = alpha + j* beta in e**(gamma*l) 
        = e**(gamma * ceff * delay) 
        = e**(gammadelay * delay)

def gammadelay(f, g, s, z0, r=0):
    return (complex(math.sqrt(r**2.0 + s**2.0 * f), 2*math.pi*f*z0)*complex(g*f, 2*math.pi*f/z0))**0.5
#end gammadelay



# Berechnung des transformierten s11 am TX-Messport (reflektierte Welle) des VNA
# za ist die Abschlussimpedanz der Leitung
# ZR ist die Systemimpedanz und in der Regel = z0
def s11(f, za, d, g, s, loss, z0):
    """
    zline = zwave(f, g, s, z0, loss)
    s111 = (za-zline)/(za+zline)
    s112 = s111 * math.e**(-2*gammadelay(f, g, s, z0, loss)*d)
    zin = zline *(1+s112)/(1-s112)
    return (zin-ZR)/(zin+ZR)
    Die vorstehende Berechnung liefert den selben Wert wie die nachfolgende, die in der HP-Literatur verwendet wird
    """
    gt = (za-ZR)/(za+ZR)
    zc = zwave(f, g, s, z0)
    g1 = (zc-ZR)/(zc+ZR)
    eg = math.e**(-2*gammadelay(f, g, s, z0, loss)*d)
    eg1 = 1.0-eg
    return (g1*(eg1-g1*gt)+eg*gt)/(1.0-g1*(eg*g1+gt*eg1))
#end s11

def s11open(f, d, loss, z0, C0, C1, C2, C3):
    c = C0*1e-15+C1*1e-27*f+C2*1e-36*f**2+C3*1e-45*f**3
    if c == 0.0:
        za = 1.0
    else:
        za = complex(0, -1/(2*math.pi*f*c)) 
    #end if
    return s11(f, za, d, 0.0, 0.0, loss, z0)
#end s11open

def s11short(f, d, loss, z0):
    return s11(f, 0.0, d, 0.0, 0.0, loss, z0)
#end s11short

# Parallelschaltung von R und parasitaerer Kapazitaet
def s11load(f, d, loss, z0, r, c):
    za = 1.0/complex(1.0/r, 2*math.pi*f *c)
    return s11(f, za, d, 0.0, 0,0, loss, z0)
#end s11load