% Demonstration of NLQFIT6 Q-factor fitting program (Matlab/GNU Octave).
%
% This software is from:
%   "Q-factor Measurement by using a Vector Network Analyser",
%   A. P. Gregory, National Physical Laboratory Report MAT 58 (2021)
%
% 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-October-2021: Released under CC0.
%            :16-August-2022:  Added code for calculating fitted Trace and plotting on
%                              polar chart. Uncomment if required.
%
% Fits FL and QL to transmission (S21) data for two-port resonator
% by using the NLQFIT6 algorithm.
%
% Load frequency and S-parameter data from a file.
% File 'Figure6b.txt' contains S21 data for Fig 6(b) in MAT 58
%
fn = 'Figure6b.txt';
data = load('-ascii', fn);
F = data(:, 1);
DR = data(:, 2);
DI = data(:, 3);
N = length(F);
D = DR + 1j*DI;

% Note: De-embed cables if they are long (not needed here)

fprintf('Demonstration of Q-factor fitting (Octave/Matlab implementation)\n')
fprintf('Data is taken from file Figure6b.txt, which is S21 data shown in Fig 6(b) of MAT58\n\n')
DMg = abs(D);
[maxmg, indx_max] = max(DMg);
Tol = maxmg*1.0E-5;  % for transmission.
% Find peak in |S21|. Used as initial value of freq.
Fseed = F(indx_max);
% 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(N)-F(1));

loop_plan = 'fwfwc';
quiet = 0;

% Step 1: Initial unweighted fit --> solution vector
sv = InitialFit(F, D, N, Fseed, Qseed);

fprintf('\nInitial solution: (eqn. 17)\n');
fprintf('  a = %6.4f + j*%6.4f\n',sv(1),sv(2));
fprintf('  b = %6.4f + j*%6.4f\n',sv(3),sv(4));
fprintf('  QL = %3.3f\n\n',sv(5));

% Step 2: Optimised weighted fit --> result vector
fprintf('Optimising...\n')
mv = OptimiseFit6(F, D, N, Fseed, sv, loop_plan, Tol, quiet);
fprintf('\nOptimised solution: (eqn. 22)\n');
fprintf('  FL = %8.6f\n',mv(6));
fprintf('  QL = %4.2f\n\n',mv(5));

% Obtain data for uncoupled resonator (scaling factor used for response calibration)
% Scaling factor A is required for calculation of unloaded Q-factor.
% A is the reciprocal of |S21| of a thru in place of resonator
A = 1/0.874;
% A = 1/(0.874*0.92); % Modified version to account for atten. of uncal. lines
d = QCircleDiam(A,mv(1),mv(2),mv(3),mv(4));
Qo = mv(5)/(1-d);
fprintf('Q-factor of uncoupled resonator (scaling factor A applied)')
fprintf('  Qo = %4.2f\n',Qo);
fprintf('Calibrated data (scaling factor A applied)\n')

fprintf('  Qcircle diameter = %6.4f\n',d)
fprintf('  GammaV = %9.6f + j*%9.6f (S21 detuned = leakage vector)\n',mv(1)*A,mv(2)*A);
b = complex(mv(1)+mv(3), mv(2)+mv(4));
fprintf('  GammaT = %9.6f + j*%9.6f (S21 tuned)\n',real(b)*A,imag(b)*A);

%%% Calculate fitted trace and plot (uncomment if required)%%%
%{
QL = mv(5);
FL = mv(6);
m6a = QL/FL;
y = 1.0./(1.0 + 2j.*m6a.*F - 2j.*QL);
Fit = mv(1)+1j.*mv(2)+ (mv(3)+1j.*mv(4)).*y;

figure;
plot(F, abs(D), '-r'); hold on;
plot(F, abs(Fit), '-g');
legend('Meas','Fit');
xlabel('Frequency (GHz)');
ylabel('Scattering Parameter Magnitude');

figure;
Meas_theta = atan2(imag(D),real(D));
Meas_rho = abs(D);
Fit_theta = atan2(imag(Fit),real(Fit));
Fit_rho = abs(Fit);

polar(Meas_theta, Meas_rho); hold on;
polar(Fit_theta, Fit_rho);
legend('Meas','Fit');
%}
