%load_ext autoreload
%autoreload 2
import skrf as rf
import skrf.mathFunctions as mf
import numpy as np
from numpy import real, log, log10, sum, absolute, pi, sqrt
import matplotlib.pyplot as plt
from matplotlib.ticker import AutoMinorLocator
from scipy.optimize import minimize

rf.stylely()

"""
Explanation: DefinedAEpTandZ0 media example
End of explanation
"""


# Load raw measurements
TL100 = rf.Network('CPWG100.s2p')
TL200 = rf.Network('CPWG200.s2p')
TL100_dc = TL100.extrapolate_to_dc(kind='linear')
TL200_dc = TL200.extrapolate_to_dc(kind='linear')

plt.figure()
plt.suptitle('Raw measurement')
TL100.plot_s_db()
TL200.plot_s_db()

plt.figure()
t0 = -2e-9
t1 = 4e-9
plt.suptitle('Time domain reflexion step response (DC extrapolation)')
ax = plt.subplot(1, 1, 1)
TL100_dc.s11.plot_z_time_step(pad=2000, window='hamming', z0=50, label='TL100', ax=ax, color='0.0')
TL200_dc.s11.plot_z_time_step(pad=2000, window='hamming', z0=50, label='TL200', ax=ax, color='0.2')
ax.set_xlim(t0, t1)
ax.xaxis.set_minor_locator(AutoMinorLocator(10))
ax.yaxis.set_minor_locator(AutoMinorLocator(5))
ax.patch.set_facecolor('1.0')
ax.grid(True, color='0.8', which='minor')
ax.grid(True, color='0.4', which='major')
plt.show()


"""
Explanation: Measurement of two CPWG lines with different lenghts
The measurement where performed the 21th March 2017 on a Anritsu MS46524B 20GHz Vector Network Analyser. The setup is a linear frequency sweep from 1MHz to 10GHz with 10'000 points. Output power is 0dBm, IF bandwidth is 1kHz and neither averaging nor smoothing are used.
CPWGxxx is a L long, W wide, with a G wide gap to top ground, T thick copper coplanar waveguide on ground on a H height substrate with top and bottom ground plane. A closely spaced via wall is placed on both side of the line and the top and bottom ground planes are connected by many vias.
| Name | L (mm) | W (mm) | G (mm) | H (mm) | T (um) | Substrate |
| :--- | ---: | ---: | ---: | ---: | ---: | :--- |
| MSL100 | 100 | 1.70 | 0.50 | 1.55 | 50 | FR-4 |
| MSL200 | 200 | 1.70 | 0.50 | 1.55 | 50 | FR-4 |
The milling of the artwork is performed mechanically with a lateral wall of 45°.
The relative permittivity of the dielectric was assumed to be approximatively 4.5 for design purpose.
End of explanation
"""


Z_conn = 53.2    # ohm, connector impedance
Z_line = 51.4    # ohm, line plateau impedance
d_conn = 0.05e-9 # s,   connector discontinuity delay
d_line = 0.45e-9 # s,   line plateau delay, without connectors

"""
Explanation: Impedance from the line and from the connector section may be estimated on the step response.
The line section is not flat, there is some variation in the impedance which may be induced by manufacturing tolerances and dielectric inhomogeneity.
Note that the delay on the reflexion plot are twice the effective section delays because the wave travel back and forth on the line.
Connector discontinuity is about 50 ps long. TL100 line plateau (flat impedance part) is about 450 ps long.
End of explanation
"""


#Make the missing reflect measurement
#Reflect only affects sign of the corrected
reflect = TL100.copy()
reflect.s[:,0,0] = 1
reflect.s[:,1,1] = 1
reflect.s[:,1,0] = 0
reflect.s[:,0,1] = 0

# Perform NISTMultilineTRL algorithm
cal = rf.NISTMultilineTRL([TL100, reflect, TL200], [1], [100e-3, 200e-3], er_est=3.0, refl_offset=[0])

plt.figure()
plt.title('Corrected lines')
cal.apply_cal(TL100).plot_s_db()
cal.apply_cal(TL200).plot_s_db()
plt.show()

"""
Explanation: Dielectric effective relative permittivity extraction by multiline method
End of explanation
"""


freq  = TL100.frequency
f     = TL100.frequency.f
f_ghz = TL100.frequency.f/1e9
L     = 0.1
A     = 0.0
f_A   = 1e9
ep_r0 = 2.0
tanD0 = 0.001
f_ep  = 1e9
x0 = [ep_r0, tanD0]

ep_r_mea = cal.er_eff.real
A_mea    = 20/log(10)*cal.gamma.real

def model(x, freq, ep_r_mea, A_mea, f_ep):
    ep_r, tanD = x[0], x[1]
    m = rf.media.DefinedAEpTandZ0(frequency=freq, ep_r=ep_r, tanD=tanD, Z0=50,
                              f_low=1e3, f_high=1e18, f_ep=f_ep, model='djordjevicsvensson')
    ep_r_mod = m.ep_r_f.real
    A_mod = m.alpha * log(10)/20
    return sum((ep_r_mod - ep_r_mea)**2)  + 0.001*sum((20/log(10)*A_mod - A_mea)**2)

res = minimize(model, x0, args=(TL100.frequency, ep_r_mea, A_mea, f_ep),
               bounds=[(2, 4), (0.001, 0.013)])
ep_r, tanD = res.x[0], res.x[1]

print('epr={:.3f}, tand={:.4f} at {:.1f} GHz.'.format(ep_r, tanD, f_ep * 1e-9))

m = rf.media.DefinedAEpTandZ0(frequency=freq, ep_r=ep_r, tanD=tanD, Z0=50,
                              f_low=1e3, f_high=1e18, f_ep=f_ep, model='djordjevicsvensson')

plt.figure()
plt.suptitle('Effective relative permittivity and attenuation')
plt.subplot(2,1,1)
plt.ylabel('$\epsilon_{r,eff}$')
plt.plot(f_ghz, ep_r_mea, label='measured')
plt.plot(f_ghz, m.ep_r_f.real, label='model')
plt.legend()

plt.subplot(2,1,2)
plt.xlabel('Frequency [GHz]')
plt.ylabel('A (dB/m)')
plt.plot(f_ghz, A_mea, label='measured')
plt.plot(f_ghz, 20/log(10)*m.alpha, label='model')
plt.legend()
plt.show()


"""
Explanation: Calibration results shows a very low residual noise floor. The error model is well fitted.
End of explanation
"""


# note: a half line is embedded in connector network
coefs = cal.coefs
r = mf.sqrt_phase_unwrap(coefs['forward reflection tracking'])
s1 = np.array([[coefs['forward directivity'],r],
        [r, coefs['forward source match']]]).transpose()

conn = TL100.copy()
conn.name = 'Connector'
conn.s = s1

# delay estimation,
phi_conn = (np.angle(conn.s[:500,1,0]))
z = np.polyfit(f[:500], phi_conn, 1)
p = np.poly1d(z)
delay = -z[0]/(2*np.pi)
print('Connector + half thru delay: {:.0f} ps'.format(delay * 1e12))
print('TDR readed half thru delay: {:.0f} ps'.format(d_line/2 * 1e12))
d_conn_p = delay - d_line/2
print('Connector delay: {:.0f} ps'.format(d_conn_p * 1e12))


# connector model with guessed loss
half = m.line(d_line/2, 's', z0=Z_line)
mc = rf.media.DefinedAEpTandZ0(m.frequency, ep_r=1, tanD=0.025, Z0=50,
                              f_low=1e3, f_high=1e18, f_ep=f_ep, model='djordjevicsvensson')
left = mc.line(d_conn_p, 's', z0=Z_conn)
right = left.flipped()
check = mc.thru() ** left ** half ** mc.thru()

plt.figure()
plt.suptitle('Connector + half thru comparison')
plt.subplot(2,1,1)
conn.plot_s_deg(1, 0, label='measured')
check.plot_s_deg(1, 0, label='model')
plt.ylabel('phase (rad)')
plt.legend()

plt.subplot(2,1,2)
conn.plot_s_db(1, 0, label='Measured')
check.plot_s_db(1, 0, label='Model')
plt.xlabel('Frequency (GHz)')
plt.ylabel('Insertion Loss (dB)')
plt.legend()
plt.show()

"""
Explanation: Relative permittivity $\epsilon_{e,eff}$ and attenuation $A$ shows a reasonnable agreement.
A better agreement could be achieved by implementing the Kirschning and Jansen miscrostripline dispersion model or using a linear correction.
Connectors effects estimation
End of explanation
"""


DUT = m.line(d_line, 's', Z_line)
DUT.name = 'model'

Check = m.thru() ** left ** DUT ** right ** m.thru()
Check.name = 'model with connectors'

plt.figure()
TL100.plot_s_db()
Check.plot_s_db(1,0, color='k')
Check.plot_s_db(0,0, color='k')
plt.show()

Check_dc = Check.extrapolate_to_dc(kind='linear')

plt.figure()
plt.suptitle('Time domain step-response')
ax = plt.subplot(1,1,1)
TL100_dc.s11.plot_z_time_step(pad=2000, window='hamming', label='Measured', ax=ax, color='k')
Check_dc.s11.plot_z_time_step(pad=2000, window='hamming', label='Model', ax=ax, color='b')
t0 = -2e-9
t1 = 4e-9
ax.set_xlim(t0, t1)
ax.xaxis.set_minor_locator(AutoMinorLocator(10))
ax.yaxis.set_minor_locator(AutoMinorLocator(5))
ax.patch.set_facecolor('1.0')
ax.grid(True, color='0.8', which='minor')
ax.grid(True, color='0.5', which='major')


"""
Explanation: Connector + thru plots shows a reasonable agreement between calibration results and model. There is a phase jump in the calibration results.
Final check
End of explanation
"""


