"""
    Test script for Q-factor fitting (pure Python version using numpy)
    ==================================================================

      Description
      -----------

         This software is from:
         "Q-factor Measurement by using a Vector Network Analyser",
         A. P. Gregory, National Physical Laboratory Report MAT 58 (2021)

         Fits FL and QL to reflection (S11) data by using the NLQFIT7 algorithm.

         Requires script qfitmod7.py
         Test data is read from file Table6c27.txt (as used in Figure 16 and
         Table 6(c) of MAT 58).

         Tested with Python 2.7 and Python 3.8


      Creative Commons CC0 Public Domain Dedication
      ---------------------------------------------

         This software is released under the Creative Commons CC0
         Public Domain Dedication. You should have received a copy of
         the Legal Code along with this software.
         If not, see <http://creativecommons.org/publicdomain/zero/1.0/>.


     Change log
     ----------

        :30-Oct-2021: Released under CC0.
        :23-Aug-2022: Added code for calculating fitted Trace, and plotting on a
                      polar chart with matplotlib. Uncomment these lines if required.
                      S-parameter data file is now read with loadtxt.
quired.

"""
import numpy as np
from qfitmod7 import *
from math import pi

# Load frequency and S-parameter data from a file.
# File 'Table6c27.txt' contains S21 data for Table 6(c) (27 MHz span) in MAT 58
fn = 'Table6c27.txt';
F, R, I, _Mg, _Ph = np.loadtxt(fn,comments='%',unpack=True)
N = len(F)
D = [complex(R[i],I[i]) for i in range(N)]
#D = R+1j*I # or use numpy array

print('Demonstration of Q-factor fitting (pure Python and numpy implementation)')
print('S11 data is taken from file Table6c27.txt (plotted in Fig. 16 of MAT58)')
print ()

# De-embed cables if they are long (not needed here)

# The convergence algorithm uses a number of steps set by loop_plan, a string of characters as follows:
#   f - fit once without testing for convergence
#   c - repeated fit, iterating until convergence is obtained
#   w - re-calculate weighting factors on basis of previous fit
#       Initially the weighting factors are all unity.
#   r - Exclude 10% of points with largesr residual error from the fit.
loop_plan ='fwfwc'

# Find minimum in |S11| - this is used to give initial value of freq.
Mg = np.absolute(D)
index_min = np.argmin(Mg)
Tol = 1.0E-5
Fseed = F[index_min]

# Set Qseed: An order-of-magnitude estimate for Q-factor
mult = 5.0 # Not critical. A value of around 5 will work well for a wide sweep-range as 
#            shown in Fig. 6(a), and for a narrow (optimised) sweep-range as shown in Fig 6(b).
Qseed = mult*Fseed/(F[-1]-F[0])

quiet = True
if not quiet: print('Initial values for iteration:  Freq=',Fseed,' QL=',Qseed)

# Step 1: Initial unweighted fit --> solution vector
sv = initialfit(F, D, N, Fseed, Qseed)
print('Initial solution (eqn. 17):');
print('  a = %6.4f + j*%6.4f' % (sv[0],sv[1]));
print('  b = %6.4f + j*%6.4f' % (sv[2],sv[3]));
print('  QL = %3.3f' % sv[4]);

print('\nOptimising ....')
# Step 2: Optimised weighted fit --> result vector
mv, weighting_ratio, number_iterations, RMS_Error = optimisefit7(F, D, N, Fseed, sv, loop_plan, Tol, quiet)
m1,m2,m3,m4,QL,FL,m7a = mv

print('Optimised solution (eqn. 26):')
print('  Fr = %10.8f' % FL)
print('  QL = %6.2f  (expect 708)' % QL)
print('  Number of iterations = %i' % number_iterations)
if weighting_ratio is not None: print('  Weighting ratio = %5.3f' % weighting_ratio)
print('  RMS_Error = %10.8f' % RMS_Error)
sqrt_eps = 1.3 # Average for and ptfe sections
print('  Fitted length of uncalibrated line = ',-float(m7a)*2.99792458e2/(4.0*pi*1.3),'mm')
print()

# Now calculate unloaded Q-factor and some other useful quantities.
print('Q-factor of uncoupled one-port resonator (unloaded Q-factor) by Method 1:')
print('Method 1 works when the uncalibrated line is attenuating')
scaling_factor_A = 'Auto'
trmode='reflection_method1'
ifail, p = GetUnloadedData(mv, scaling_factor_A, trmode, quiet)
Qo, cal_diam, cal_gamma_V, cal_gamma_T = p
print('  Qo = %6.2f' % Qo)
print('  Qcircle diameter = %6.4f' % cal_diam)
print('  GammaV = %10.8f +j* %10.8f (S11 detuned)' % (cal_gamma_V.real, cal_gamma_V.imag))
print('  |GammaV| = %10.8f' % abs(cal_gamma_V))
print('  GammaT = %10.8f +j* %10.8f (S11 tuned)' % (cal_gamma_T.real, cal_gamma_T.imag))
print ()

print('Q-factor of uncoupled one-port resonator (unloaded Q-factor) by Method 2:')
print ('Scaling factor A = 1.0 (assume no attenuation in uncalibrated line & the VNA is calibrated)')
scaling_factor_A = 1.0
trmode='reflection_method2'
ifail, p = GetUnloadedData(mv, scaling_factor_A, trmode, quiet)
Qo, cal_diam, cal_gamma_V, cal_gamma_T = p
print('  Qo = %6.2f  (expect 862)' % Qo)
print('  Qcircle diameter = %6.4f' % cal_diam)
print('  GammaV = %10.8f +j* %10.8f (S11 detuned)' % (cal_gamma_V.real, cal_gamma_V.imag))
print('  |GammaV| = %10.8f' % abs(cal_gamma_V))
print('  GammaT = %10.8f +j* %10.8f (S11 tuned)' % (cal_gamma_T.real, cal_gamma_T.imag))


#### Calculate fitted trace and then plot (uncomment if required) ###
"""
import matplotlib.pyplot as plt
from numpy import angle, absolute, empty, vectorize
m6a = QL/FL
vcmplex = vectorize(complex)
y = 1.0/vcmplex(1.0,2*(m6a*F-QL))
expm7 = np.exp(vcmplex(0.0,m7a*(F-FL)))   # m7a = m7/Flwst
Fit = (complex(m1,m2) + complex(m3,m4)*y)*expm7

fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
ax.plot(angle(D), Mg,'-k',label='Measured')
ax.plot(angle(Fit), absolute(Fit),'ro',markevery=5,markersize=2,label='Fit')
ax.set_rmax(1)
ax.set_rticks([0.25, 0.5, 0.75, 1])
ax.set_rlabel_position(-20)  # Move radial labels away from plotted line
ax.grid(True)
ax.legend();

ax.set_title("Measured and fitted Q-circles", va='bottom')
plt.show()

fig2, ax2 = plt.subplots()
ax2.plot(F, Mg,'-k',label='Measured')
ax2.plot(F, absolute(Fit),'ro',markevery=5,markersize=2,label='Fit')
ax2.grid(True)
plt.xlabel('Frequency (GHz)')
plt.ylabel('|S11|')
ax2.legend();

ax2.set_title("Measured and fitted magnitude", va='bottom')
plt.show()
"""

"""
OUTPUT:

Demonstration of Q-factor fitting (pure Python and numpy implementation)
S11 data is taken from file Table6c27.txt (plotted in Fig. 16 of MAT58)

Initial solution (eqn. 17):
  a = 760.9731 + j*67.7804
  b = 0.0609 + j*-0.6432
  QL = 779.068

Optimising ....
Optimised solution (eqn. 26):
  Fr = 3.65293800
  QL = 708.49  (expect 708)
  Number of iterations = 4
  Weighting ratio = 28.317
  RMS_Error = 0.00145957
  Fitted length of uncalibrated line =  57.47056249055141 mm

Q-factor of uncoupled one-port resonator (unloaded Q-factor) by Method 1:
Method 1 works when the uncalibrated line is attenuating
  Qo = 862.58
  Qcircle diameter = 0.3573
  GammaV = 0.09084890 +j* -0.99586469 (S11 detuned)
  |GammaV| = 1.00000000
  GammaT = 0.05878773 +j* -0.64003179 (S11 tuned)

Q-factor of uncoupled one-port resonator (unloaded Q-factor) by Method 2:
Scaling factor A = 1.0 (assume no attenuation in uncalibrated line & the VNA is calibrated)
  Qo = 861.67  (expect 862)
  Qcircle diameter = 0.3538
  GammaV = 0.08996562 +j* -0.98618239 (S11 detuned)
  |GammaV| = 0.99027750
  GammaT = 0.05821616 +j* -0.63380908 (S11 tuned)

"""
