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

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

         Fitting algorithm from:
         "Q-factor Measurement by using a Vector Network Analyser",
         A. P. Gregory, National Physical Laboratory Report MAT 58 (2021).
         Refer to equivalent circuit Figure 8.

         Fits FL and QL to simulated S11 data for one-port resonator
         by using the NLQFIT6 algorithm. This is not the usual case as this
         implies that the uncalibrated line has zero length (for real
         experiments would normally use NLQFIT7).

         Requires script qfitmod6.py

         Tested with 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
     ----------

        :23-Aug-2022: Release under CC0 licence.
"""

import sys
from qfitmod6 import *

Ccav = 3.2E-11
Lcav = 9.0E-11
print ('Lcav =',Lcav, 'Ccav =',Ccav)
FL = 1.0E-9/(2.0*math.pi*math.sqrt(Ccav*Lcav))
print ('FL =',FL, 'GHz')
Ro = 650.0
g = 5.0 # turns ratio (see MAT 58 Figure 8).

# NLQFIT is based on a 3-component equivalent circuit model.
# I.e. assume that series resistance (Rs) and inductance (Ls) of MAT 58 Figure 8 are zero.
# Can experiment with non-zero values, e.g. Ls = 2E-10 and Rs = 1.5.
Ls = 0.0
Rs = 0.0

# Calculate Qo directly.
#   Energy stored = C*V^2  [V is average voltage]
#   Energy dissipated per wave cycle = V^2/(R*w)  [w is angular frequency] 
#   --> Q = w*C*R = R/(w*L)
Qo = Ro/(2.0E9*math.pi*FL*Lcav)
print ('Qo = R/wL = %6.2f\n'%Qo)


# Calculate reflection coefficient (S11) measured by calibrated VNA
def gamma(deltaF, Fr):
    F_GHz = Fr+deltaF
    omega = 2.0E9*math.pi*F_GHz
    Zo = 50.0   # VNA port impedance is 50.0 Ohms.
    Adm = 1.0/Ro + 1.0/complex(0.0, omega*Lcav) + complex(0.0, omega*Ccav)
    Z = 1.0/(Adm*g*g) + complex(0.0,omega*Ls) + Rs
    Gam = (Z-Zo)/(Z+Zo)  # the usual formula for reflection from impedance mismatch.
    return Gam

def GetDataSet():
    Frange = 0.03  # GHz
    Fstep = Frange/200.0
    N = 201
    deltaF = [-Frange/2.0+i*Fstep for i in range(N)]
    F = [FL+deltaF[i] for i in range(N)]    
    Gamma = [gamma(deltaF[i], FL) for i in range(N)]
    return F, Gamma

F,Gamma = GetDataSet()  # Frequency (float) and S11 (complex)

Mg = np.absolute(Gamma)
index_min = np.argmin(Mg) # Use minimum |Gamma| as estimate of resonance frequency

Tol = 1.0E-5
Fseed = F[index_min]

mult = 5.0
Qseed = mult*Fseed/(F[-1]-F[0])
print('Fseed =',Fseed,'Qseed =',Qseed)
N = len(Gamma)
sv = initialfit(F, Gamma, N, Fseed, Qseed)

loop_plan ='fwfwc'
quiet = True

mv, weighting_ratio, number_iterations, RMS_Error = optimisefit6(F, Gamma, N, Fseed, sv, loop_plan, Tol, quiet)
m1,m2,m3,m4,QL,FL2 = mv

print('Optimised solution (eqn. 22):')
print('  FL = %10.8f   (expect %10.8f GHz)'%(FL2,FL))
print('  QL = %6.2f'%QL)
if weighting_ratio is not None: print('  Weighting ratio = %5.3f' % weighting_ratio)
print('  RMS_Error = %10.8f' % RMS_Error)
print()

# Now calculate unloaded Q-factor and some other useful quantities.
# For Method1, use eqn. 32 with D=2. The VNA can be calibrated (preferred) or uncalibrated.
# Would normally supply Scaling factor A 'AUTO', in which case A will be determined automatically
# by using eqn. 34, so attenuation in an uncalibrated line will have no effect on unloaded Q.
scaling_factor_A = 'AUTO' 
trmode='reflection_method1'
ifail, p = GetUnloadedData(mv, scaling_factor_A, trmode, quiet)
if ifail: print(p); sys.exit()
Qo2, cal_diam, cal_gamma_V, cal_gamma_T = p
print('Q-factor of uncoupled one-port resonator (unloaded Q-factor) by METHOD1')
print('  Qo = %6.2f  (expect %0.2f if Ls=Rs=0.0)' % (Qo2,Qo))
print('  Qcircle diameter = %0.4f' % cal_diam)
print('  GammaV = %10.8f +j* %10.8f (S11 detuned)' % (cal_gamma_V.real, cal_gamma_V.imag))
print('  |GammaV| = %10.8f (expect 1.0000000 if AUTO)' % abs(cal_gamma_V))

#### Calculate fitted trace and plot (uncomment if required) ###
"""
import matplotlib.pyplot as plt
from numpy import angle, absolute, empty, array, vectorize
m6a = QL/FL
t = 2*(array(F,dtype=float)-FL)/FL
x = 2*(m6a*F-QL)
vcmplex = vectorize(complex)
y = 1.0/vcmplex(1.0,2*(m6a*F-QL))
Fit = complex(m1,m2) + complex(m3,m4)*y

fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
ax.plot(angle(Gamma), absolute(Gamma),'-k',label='LCR Model')
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("Q-circles", va='bottom')
plt.show()
"""

