% Demonstration of NLQFIT8 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 NLQFIT8 algorithm.
%
% Load frequency and S-parameter data from file Figure23.txt (see Fig. 23 of MAT58).
%
fn = 'Figure23.txt';
data = load('-ascii', fn);
F = data(:, 1);
DR = data(:, 2);
DI = data(:, 3);
N = length(F);
ncable_len = 1.2; % root_eps * cable length in metres
D = (DR + 1j*DI) .* exp(-2j*pi*F*ncable_len*1.0E9/2.99792458E8);
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 = 'fwfwfwc';
quiet = 0;

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

fprintf('\nInitial solution:\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 = OptimiseFit8(F, D, N, Fseed, sv, loop_plan, Tol, quiet);
fprintf('\nOptimised solution:\n');
fprintf('  FL = %8.6f\n',mv(8));
fprintf('  QL = %4.2f\n\n',mv(7));

fprintf('Calculations for transmission resonator (lossless coupling):\n')
% Magnitude |S21| of a 'thru' connection = 0.949
A = 1.0/0.949;
d = QCircleDiam(A,mv(1),mv(2),mv(3),mv(4));
Qo = mv(7)/(1-d);
fprintf('  Qo = %4.2f\n',Qo);
fprintf('  Qcircle diameter = %6.4f\n',d)

%%% Calculate fitted trace and plot (uncomment if required)%%%
%{
QL = mv(7);
FL = mv(8);
m6a = QL/FL;
y = 1.0./(1.0 + 2j.*m6a.*F - 2j.*QL);
t = 2*(F-FL)/FL;
Fit = mv(1)+1j.*mv(2)+ (mv(5)+1j.*mv(6)).*t + (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');
%}
