% Demonstration of NLQFIT7 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 and the length of uncalibrated line to S11 data
% by using the NLQFIT7 algorithm.
%
% Load frequency and S-parameter data from a file.
% File 'Table6c27.txt' contains S11 data used in Table 6(c) of MAT 58
% For 27 MHz span with GPC-7 cal.
%
fn = 'Table6c27.txt';
data = load('-ascii', fn);
F = data(:, 1);
DR = data(:, 2);
DI = data(:, 3);
N = length(F);
D = DR + 1j*DI;
DMg = abs(D);

Tol = 1.0E-5;  % for reflection

[minmg, indx_min] = min(DMg);
Fseed = F(indx_min); % Find minimum in |S11|. 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:\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 = OptimiseFit7(F, D, N, Fseed, sv, loop_plan, Tol, quiet);
fprintf('\nOptimised solution:\n');
fprintf('  FL = %8.6f\n',mv(6));
fprintf('  QL = %4.2f\n',mv(5));
L = -2.99724258E2*mv(7)/(4*pi); % for refl.
fprintf('  Len. uncal. line (for ref_index 1.3)= %5.3f (mm)\n\n',L/1.3)

% Calculate GammaV and Qo (reflection case)
%
% Method 1: Coupling has zero series-impedance. A given by eqn. (34) (assumes |A|<1.0
% indicates loss in uncalibrated line).
% D = 2
fprintf('Qo calculations for reflection resonator:\n')
fprintf('Method 1 (coupling has zero series-impedance, attenuating uncalibrated line):\n')
A = 1.0/abs(mv(1)+1j*mv(2));
gv = A*abs(mv(1)+1j*mv(2));
fprintf('  GammaV = %9.6f + j*%9.6f (S11 detuned)\n',mv(1)*A,mv(2)*A);
fprintf('  |GammaV| = %9.6f\n',gv);
d = QCircleDiam(A,mv(1),mv(2),mv(3),mv(4));
d2 = 2.0;
Qo = mv(5)*(1.0+ 1.0/(d2/d - 1.0));
fprintf('  Qo = %4.2f\n',Qo);

% Method 2: Loop coupling. A = 1.0 if uncalibrated line lossless.
%         D (=d2) given by eqn. (35)
fprintf('Method 2 (loop coupling with non-attenuating uncalibrated line):\n')
A = 1.0;
gv = A*abs(mv(1)+1j*mv(2));
fprintf('  GammaV = %9.6f + j*%9.6f (S11 detuned)\n',mv(1)*A,mv(2)*A);
fprintf('  |GammaV| = %9.6f\n',gv);
d = QCircleDiam(A,mv(1),mv(2),mv(3),mv(4));
b = mv(1)+mv(3) + 1j*(mv(2)+mv(4));
gv2 = gv*gv;
mb = abs(b)*A;
cosphi = (gv2 + d*d - mb*mb)/(2.0*gv*d);
d2 = ((1.0-gv2)/(1.0-gv*cosphi)); % d2 is D
fprintf('  D = %6.4f\n',d2);
Qo = mv(5)*(1.0+ 1.0/(d2/d - 1.0));
fprintf('  Qo (loop coupling) = %4.2f\n',Qo);

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

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');
%}
