<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
Step1:  Power Handling
For coaxial, the electric field is 
Step2:  This is maximum for $\rho=a$. Let us suppose that the peak voltage $V_p$ (or equivalently the maximum RF power $P_p$ the line can handle) is set by the breakdown voltage of the insulator. The electric field strength at which breakdown occurs depends on the respective geometries of the insulator and the electrodes with which the electric field is applied. In the case of the air, the dielectric strength is $E_p$=3.0 MV/m.
$$
E_{max} = E(\rho=a) = V_p / (a \ln(b/a))
$$
Thus,
$$
V_p = E_{max} a \ln(b/a)
$$
Step3:  One get a peak voltage of 104kV. In terms of power, for a 30 Ohm characteric impedance, this leads to
Step4:  A maximum power of 181 MW.
Now let's look how the power handling (ie the max power) evolves vs the geometrical parameters $a$ and $b$. As seen above, the maximum power is
Step6:  The second term is maximum at $b/a$ is around 1.65.
Step7:  This maximum corresponds to an impedance of 
Step8:  Probe equations
We measure the $S_{21}$ term, where port 1 energizes the rigid coaxial line and port 2 the probe. By definition
Step9:  The voltage in the rigid line is thus given from the [0-10V] signal voltage by
Step10:  As we want to measure voltage between 5 and 55 kV, using the previous expressions lead to
Step11:  The latter equation is satisfied for $b/a=3.6$, which gives an optimum characteristic impedance of 77 Ohm
Step12:  Optimum ?
Let's plot together the attenuation and the power handling vs characteristic impedance vs Z0
Step13:  Low Frequency Measurements
We load the $b_2/b_1$ measurement which is performed at low frequencies
Step14:  Let's plot the results versus the frequency
Step15:  Let's do the same into a log scale
Step16:  The plot look like pretty linear. Let's fit with a linear trends to see how it evolves with respect to the log10 of the frequency
Step17:  One can clearly see that the $b_2/b_1$ ratio increases linearly with respect to frequency decades (a decade is a factor of 10 difference between two numbers). More precisely, it increases by 20dB per decade.  
So, if we assumes that this gain $G=b_2/b_1$ can be extrapolated from the low frequencies $f_{low}$ to the working frequency $f$ (60 MHz), one should get 
Step18:  The gain at 60 MHz extrapolated from low frequency measurement is thus G=-52.99 dB $\pm$ 0.25 dB
Let's look to the welded probe results
CEA-side and DUT-side probes
Step19:  A little bit of theory
So, why does the gain evolves that way? Let's look to the probe configuration 
The transfer function is such
<ASSISTANT_TASK:>
Python Code:
a = 140e-3 /2 # inner conductor radius
b = 230e-3 /2 # inner conductor radius
def coax_electric_field(rho, V, a, b):
    
    Returns the electric field in a coaxial line.
    
    return 1/rho*V/log(b/a)
rho = linspace(a, b, 101)
V = 1 # V
plot(rho/a, coax_electric_field(rho, V, a, b)/V)
xlabel('rho/a')
ylabel('E_{rho} / V')
Explanation: Power Handling
For coaxial, the electric field is :
$$
E(\rho)=\frac{1}{\rho}\frac{V}{\ln\left(b/a\right)}
$$
where $a<\rho<b$.
This is maximum
End of explanation
Ep = 3.0e6 # V/m
Vp = Ep*a*log(b/a)
Vp
Explanation: This is maximum for $\rho=a$. Let us suppose that the peak voltage $V_p$ (or equivalently the maximum RF power $P_p$ the line can handle) is set by the breakdown voltage of the insulator. The electric field strength at which breakdown occurs depends on the respective geometries of the insulator and the electrodes with which the electric field is applied. In the case of the air, the dielectric strength is $E_p$=3.0 MV/m.
$$
E_{max} = E(\rho=a) = V_p / (a \ln(b/a))
$$
Thus,
$$
V_p = E_{max} a \ln(b/a)
$$
End of explanation
Z0 = 30
Pd = Vp**2 / (2*Z0)
print(Pd/1e6)
Explanation: One get a peak voltage of 104kV. In terms of power, for a 30 Ohm characteric impedance, this leads to:
$$
P_d = \frac{V^2}{2 Z_0}
$$
End of explanation
x = linspace(1, 10, 1001)
y = log(x)/x**2
plot(x, y, lw=2)
xlabel('b/a')
axvline(x=exp(1/2), color='r')
grid(True)
title('$\ln(b/a) / (b/a)^2$')
print('The maximum is for x=', x[argmax(y)])
Explanation: A maximum power of 181 MW.
Now let's look how the power handling (ie the max power) evolves vs the geometrical parameters $a$ and $b$. As seen above, the maximum power is: 
$$
P_{max} 
= \frac{V^2_{max}}{2 Z_0}
$$
with 
$$
    Z_0 = \frac{377}{2\pi}\sqrt{\frac{\mu_r}{\varepsilon_r}} \ln\left(\frac{b}{a}\right)
$$
and
$$
    V_{max} = E_{max} a \ln\left(\frac{b}{a}\right)
$$
Thus,
$$
P_{max} \approx \frac{E^2_{max} a^2}{120} \ln\left(\frac{b}{a}\right)
$$
Or, rewriting in terms of $b/a$ :
$$
P_{max} \approx \frac{E^2_{max} b^2}{120} \frac{ \ln\left(\frac{b}{a}\right) }{\left(b/a\right)^2}
$$
$$
\frac{\partial P_{max}}{\partial a}
= 
0
\rightarrow
\frac{\partial}{\partial a}
\left( \frac{\ln(b/a)}{(b/a)^2}\right)
=
0
$$
ie
$$
2\ln(b/a) - 1 = 0
$$
which is satisfied for 
$$
b/a=\exp(1/2)\approx 1.65
$$
End of explanation
def coax_char_impedance(a, b, eps_r=1, mu_r=1):
    
    Returns the characteristic impedance of a coaxial line
    
    return 377/(2*pi)*sqrt(mu_r/eps_r)*log(b/a)
Explanation: The second term is maximum at $b/a$ is around 1.65.
End of explanation
coax_char_impedance(1, 1.65)
Explanation: This maximum corresponds to an impedance of :
End of explanation
30/50 * 1/10
(sqrt(6/100) - sqrt(1/10))/sqrt(1/10)
20*log10(sqrt(50/30)*10**(-80/20))
Explanation: Probe equations
We measure the $S_{21}$ term, where port 1 energizes the rigid coaxial line and port 2 the probe. By definition:
$$
S_{21} = \frac{b_2}{a_1}
$$
Our measurement is the combination of transmission losses due to the impedance transformers and the rigid line up to the probe duct, i.e. $S_{21} = S_{21,loss} + S_{21,probe}$, with: 
$$
S_{21,probe} = \frac{b_2}{a_l}
$$
where $a_l$ is the power wave in the rigid line at the probe location. 
If we assume that the rigid line is perfectly matched, one gets:
$$
V_{l} = \sqrt{30} (a_l + b_l) = \sqrt{30} a_l \approx \sqrt{30} a_1 \rightarrow  a_1 = \frac{V_l}{\sqrt{30}}
$$
and if the probe is not energized from port 2:
$$
V_2 = \sqrt{50} b_2 \rightarrow b_2 = \frac{V_2}{\sqrt{50}}
$$
Thus,
$$
S_{21,probe} = \frac{\sqrt{30}}{\sqrt{50}} \frac{V_2}{V_l} 
$$
The voltage in the rigid line $V_l$ can be deduced from the voltage at the probe $V_2$ and the transmission loss of the setup $S_{21,probe}$:
$$
\boxed{
V_l = \frac{\sqrt{30}}{\sqrt{50}} \frac{V_2}{S_{21,probe}} 
}$$
End of explanation
S21 = -80 # Example of attenuation factor [dB]
V_line = 55e3 # Voltage in the transmission line [V]
V_probe = V_line * 10**(S21/20) # voltage at probe in [V]
print('Voltage at probe: {} V'.format(V_probe))
Z0 = 50 # Characteristic Impedance of the probe cable and amplifier [Ohm]
Power_probe = 10*np.log10(V_probe**2 / (2*Z0) * 1000) # power at the probe in [dBm]
print('Power at probe: {} dBm'.format(Power_probe)) 
Explanation: The voltage in the rigid line is thus given from the [0-10V] signal voltage by:
$$
V_{line} = \sqrt{2 Z_0 P_{line}}
$$
where the power in the line $P_{line}$ is given from the power in the probe $P_{probe}$ by:
$$
P_{line} [W] = P_{probe} [W] \times 10^{|G|/10}
$$
with $G$ the probe gain ($G=|S_{21,probe,dB}| - losses$) in dB. Finally, if the probe power is expressed in dBm, we have:
$$
\boxed{
V_{line} = \sqrt{2 \times 30 \times \frac{10^{(P_{probe}[dBm]+|G])/10}}{1000} }
}
$$
Signal Conditioning
The probe input power $P_{probe}$ is generally expressed in dBm and is rectified into a 0-10V voltage $V_s$ signal using a log amplifier. The output signal voltage is then proportional to the power input signal:
$$
V_{s} [V] = a \times P_{probe} [dBm] + b
$$
The calibration process which gives coefficients $a$ and $b$ is performed for different frequencies. Once calibrated, these coefficients are used to invert the signal voltage to retrieve the power in the prove in the above equation:
$$
P_{probe} [dBm] = \frac{V_s}{a} - \frac{b}{a}
$$
Probe Goal Attenuation
We recall here that we search for the ideal probe attenuation in order that the power ouput by the probe should be between [-4, 25] dBm and gives a [0,10] V signal, for an input voltage in the line between [5, 55] kV
Let's first illustrate this by an example, for a voltage in the line of 55 kV and -80 dB probe attenuation:
End of explanation
from scipy.optimize import minimize
def fun(u):
    return abs(1+1/u-log(u))
sol = minimize(fun, x0=3)
print(sol.x)
u = linspace(0.1, 10, 501)
test = 1+1/u-log(u)
plot(u, test, lw=2)
xlabel('b/a')
title('$1+1/(b/a)-ln(b/a)$')
#u_0 = u[argmin(abs(test))]
axhline(y=0, color='k')
grid()
axvline(x=sol.x, color='r')
Explanation: As we want to measure voltage between 5 and 55 kV, using the previous expressions lead to: 
$$
10 \log\left( 
\frac{10^{P_{probe,[dBm]}/10} 2 Z_0 }{ 1000  V_{line}^2 }
\right)
$$
Coaxial Line Attenuation
Let's fix the frequency and vary the dimensions to see how evolve the attenutation
$$
\frac{\partial \alpha}{\partial a} = 0 
\rightarrow 
\frac{R_s}{2} \sqrt{\frac{\varepsilon}{\mu}} 
\frac{\partial }{\partial a}
\left( 
\frac{1/a+1/b}{\ln(b/a)}
\right)
=0
$$
$$
\rightarrow 1 + a/b - \ln\left(\frac{b}{a}\right) = 0
$$
Let's see where this equation is satisfied :
End of explanation
a = 200e-3
coax_char_impedance(a, 3.6*a)
Explanation: The latter equation is satisfied for $b/a=3.6$, which gives an optimum characteristic impedance of 77 Ohm:
End of explanation
def coax_maximum_power(a, b, Emax=3.0e6):
    # deduces the maximum electric field from voltage breakdown value
    Vp = Emax*a*log(b/a)
    Z0 = coax_char_impedance(a, b)
    return Vp**2 / (2*Z0)
def coax_alpha(a, b, f, sigma=6e7):
    # sheet resistance [Ohm]
    Rs = sqrt(2*pi*f*mu_0/(2*sigma))
    # Resistance per meter [Ohm/m]
    R = Rs/(2*pi)*(1/a + 1/b)
    # Characteristic Impedance [Ohm]
    Z0 = coax_char_impedance(a, b)
    # Attenuation [1/m]
    return R/(2*Z0)
a = linspace(20e-3, 450e-3) # inner conductor radius
b = 500e-3
f = 60e6
Z0 = coax_char_impedance(a, b)
alpha = coax_alpha(a, b, f)
alpha_norm = alpha/min(alpha)
max_pow = coax_maximum_power(a, b)
max_pow_norm = max_pow/max(max_pow)
semilogx(Z0, alpha_norm, lw=2)
semilogx(Z0, max_pow_norm, '--', lw=2)
legend(('Attenuation', 'Power handling'), fontsize=14, loc='best')
grid(True, which='major')
grid(True, which='minor', axis='x')
xlabel('$Z_0$ [$\Omega$]', fontsize=14)
ylabel('Normalized Values', fontsize=14)
axis([10, 200, 0, 2])
ax=gca()
ax.annotate('Maximum at $30\Omega$', xy=(30,1), xytext=(15,0.60), 
            color='g', arrowprops=dict(facecolor='g', width=5),
            fontsize=14)
ax.annotate('Maximum at $77\Omega$', xy=(77,1), xytext=(50,1.5), 
            color='b', arrowprops=dict(facecolor='b', width=5),
            fontsize=14)
Explanation: Optimum ?
Let's plot together the attenuation and the power handling vs characteristic impedance vs Z0
End of explanation
freq, reB2overB1, imB2overB1 = loadtxt('ICRH_Probe-Calibration/SLIDING_BF.ASC', skiprows=14, delimiter=';', unpack=True)
B2overB1 = reB2overB1 + 1j*imB2overB1
B2overB1_dB = 20*log10(abs(B2overB1))
Explanation: Low Frequency Measurements
We load the $b_2/b_1$ measurement which is performed at low frequencies
End of explanation
plot(freq/1e6, B2overB1_dB)
xlabel('f [MHz]')
ylabel('$b_2/b_1$ [dB]')
grid()
xlim(min(freq/1e6), max(freq/1e6))
Explanation: Let's plot the results versus the frequency
End of explanation
semilogx(freq, B2overB1_dB)
xlabel('f [Hz]')
ylabel('$b_2/b_1$ [dB]')
grid()
xlim(min(freq), max(freq))
Explanation: Let's do the same into a log scale:
End of explanation
freq_log = log10(freq)
# fits the curve with a linear polynom y=a*x+b
a,b=polyfit(freq_log, B2overB1_dB, deg=1)
print(a,b)
semilogx(freq, B2overB1_dB, lw=2)
semilogx(freq, a*freq_log + b, '--r', lw=3)
grid()
grid(which='minor', axis='x')
xlim(min(freq), max(freq))
xlabel('f [Hz]')
ylabel('$b_2/b_1$ [dB]')
Explanation: The plot look like pretty linear. Let's fit with a linear trends to see how it evolves with respect to the log10 of the frequency:
End of explanation
f = 62.5e6
G = B2overB1_dB + 20*log10(f/freq)
G_mean = mean(G)
G_std = std(G)
plot(G)
axhline(y=G_mean, color='r')
grid()
ylabel('G [dB]')
xlabel('measurement points')
print(G_mean, G_std)
Explanation: One can clearly see that the $b_2/b_1$ ratio increases linearly with respect to frequency decades (a decade is a factor of 10 difference between two numbers). More precisely, it increases by 20dB per decade.  
So, if we assumes that this gain $G=b_2/b_1$ can be extrapolated from the low frequencies $f_{low}$ to the working frequency $f$ (60 MHz), one should get :
$$
G_{[dB]}(f) = G(f_{low}) + 20 \log10 \left( f / f_{low} \right) 
$$
Let's make an exemple to get the gain at 60 MHz from the low frequency measurements:
End of explanation
# Band sweep
CEA1_freq, CEA1_reB2overB1, CEA1_imB2overB1 = loadtxt('ICRH_Probe-Calibration/CEA_BF1.ASC', skiprows=14, delimiter=';', unpack=True)
CEA1_B2overB1_dB = 20*log10(abs(CEA1_reB2overB1 + 1j*CEA1_imB2overB1))
# Point sweep
CEA2_freq, CEA2_reB2overB1, CEA2_imB2overB1 = loadtxt('ICRH_Probe-Calibration/CEA_BF2.ASC', skiprows=14, delimiter=';', unpack=True)
CEA2_B2overB1_dB = 20*log10(abs(CEA2_reB2overB1 + 1j*CEA2_imB2overB1))
DUT1_freq, DUT1_reB2overB1, DUT1_imB2overB1 = loadtxt('ICRH_Probe-Calibration/DUT_BF1.ASC', skiprows=14, delimiter=';', unpack=True)
DUT1_B2overB1_dB = 20*log10(abs(DUT1_reB2overB1 + 1j*DUT1_imB2overB1))
DUT2_freq, DUT2_reB2overB1, DUT2_imB2overB1 = loadtxt('ICRH_Probe-Calibration/DUT_BF2.ASC', skiprows=14, delimiter=';', unpack=True)
DUT2_B2overB1_dB = 20*log10(abs(DUT2_reB2overB1 + 1j*DUT2_imB2overB1))
semilogx(CEA1_freq, CEA1_B2overB1_dB)
semilogx(CEA2_freq, CEA2_B2overB1_dB)
grid()
grid(which='minor', axis='x')
xlim(min(CEA1_freq), max(CEA1_freq))
xlabel('f [Hz]')
ylabel('$b_2/b_1$ [dB]')
title('CEA probe')
a,b = polyfit(CEA1_freq[1000:], CEA1_B2overB1_dB[1000:], deg=1)
print(a,b)
a,b = polyfit(CEA2_freq, CEA2_B2overB1_dB, deg=1)
print(a,b)
plot(CEA1_freq, CEA1_B2overB1_dB)
plot(CEA1_freq[500:], CEA1_B2overB1_dB[500:])
def moving_average(a, n=3) :
    ret = np.cumsum(a, dtype=float)
    ret[n:] = ret[n:] - ret[:-n]
    return ret[n - 1:] / n
import matplotlib.pylab as pylab
pylab.rcParams['figure.figsize'] = 20, 10
f = 62.5e6
G_DUT1 = DUT1_B2overB1_dB + 20*log10(f/DUT1_freq)
G_DUT2 = DUT2_B2overB1_dB + 20*log10(f/DUT2_freq)
G_CEA1 = CEA1_B2overB1_dB + 20*log10(f/CEA1_freq)
G_CEA2 = CEA2_B2overB1_dB + 20*log10(f/CEA2_freq)
subplot(121)
plot(G_CEA1)
plot(G_CEA2)
grid()
ylabel('G [dB]')
xlabel('measurement points')
title('CEA-side probe')
plot(moving_average((G_CEA1+G_CEA2)/2,10), lw=2)
subplot(122)
plot(G_DUT1)
plot(G_DUT2)
grid()
ylabel('G [dB]')
xlabel('measurement points')
title('DUT-side probe')
plot(moving_average((G_DUT1+G_DUT2)/2,10), lw=2)
# average for the last 1000 points
# Standard error calculation to determine the error bar
N=500
CEA_avg = mean(moving_average((G_CEA1[N:]+G_CEA2[N:])/2,10))
DUT_avg = mean(moving_average((G_DUT1[N:]+G_DUT2[N:])/2,10))
DUT_standard_error = std((G_DUT1[N:]+G_DUT2[N:])/2) / sqrt(len(G_DUT1[N:]))
CEA_standard_error = std((G_CEA1[N:]+G_CEA2[N:])/2) / sqrt(len(G_CEA1[N:]))
print('CEA probe average extrapoled gain : G={}+/-{} dB'.format(CEA_avg, CEA_standard_error))
print('DUT probe average extrapoled gain : G={}+/-{} dB'.format(DUT_avg, DUT_standard_error))
Explanation: The gain at 60 MHz extrapolated from low frequency measurement is thus G=-52.99 dB $\pm$ 0.25 dB
Let's look to the welded probe results
CEA-side and DUT-side probes
End of explanation
R = 50
C1 = 1e-15 # unknow
Cs = C1 * (10**(211/20) -1) # fit from the sliding probe linear fit (-211 dB at omega=0)
f = logspace(0, 9, 301)
Z_Cs = 1/(1j*2*pi*f*Cs)
Z_C1 = 1/(1j*2*pi*f*C1)
G = (R + Z_Cs)/(R + Z_Cs + Z_C1)
G_dB = 20*log10(abs(G))
20*log10( 1/Cs / (1/Cs + 1/C1) )
semilogx(f, G_dB, lw=2)
grid()
grid(which='minor', axis='x')
Explanation: A little bit of theory
So, why does the gain evolves that way? Let's look to the probe configuration 
The transfer function is such:
$$
\frac{V_{probe}}{V_{line}}
=
\frac{R + Z_{C_s}}{R + Z_{C_s} + Z_{C_1}}
$$
ie
$$
\frac{V_{probe}}{V_{line}}
=
$$
where $C_s$ is the stray capacitance and $C_1 \ll C_s$.
End of explanation
<END_TASK>