% Demonstration of NLQFIT6 Q-factor fitting program for two-port absorption (notch)
% transmission resonator (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.
%            :16-August-2024:  Added experimental correction for pivoting of Q-circle. 
%
% Fits FL and QL to S21 data for notch resonator (MAT 58 Figure 26).
%
% Load frequency and S-parameter data from a file.
% File 'Figure27.txt' contains S21 data for Fig 27 in MAT 58
%
fn = 'Figure27.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 Figure27.txt, which is S21 data shown in Fig 27 of MAT58\n\n')
DMg = abs(D);
[maxmg, indx_max] = max(DMg);
Tol = maxmg*1.0E-5;  % for transmission notch.
[minmg, indx_min] = min(DMg);
Fseed = F(indx_min); % Find minimum in |S21|. Used as initial value of freq.
% 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));
fprintf('Unloaded Q-factor (eqn. 45)\n')
A = 1.0/abs(mv(1)+1j*mv(2)); % Set scale so that Mag(S21)=1.0 at pole.
d = QCircleDiam(A,mv(1),mv(2),mv(3),mv(4));
Qo = mv(5)*(1.0+ 1.0/(1.0/d - 1.0));
fprintf('  Qo = %4.2f\n',Qo);
ST = A*(mv(1)+mv(3) + 1j*(mv(2)+mv(4)));
SV = A*(mv(1) + 1j*mv(2));
rtc = abs(ST+SV)/2+d/2;  % touching circle radius
dc = d/rtc;              % corrected Q-circle diameter
Qo1 = mv(5)*(1.0+ 1.0/(1.0/dc- 1.0));
fprintf('  Qo (with experimental diameter correction) = %4.2f\n\n',Qo1);

%%% 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');
disp('Note: Rotational and scaling factors are applied in MAT 58 Figure 27')
%}
