function sv = InitialFit(F, S, N, EstimatedResFreq, Qseed)
%
% NLQFIT6, NLQFIT7 and NLQFIT8 Step (1). Matlab/GNU Octave.
% Linear least squares Q-factor fit.
% A reasonable estimate for the resonant frequency must be
% supplied. As this is not optimised in this function, the
% solution will only be approximate.
%
% This code listing 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.
%            :12-March-2024:   Tidied equations (now same as latest
%                              Python version. No change to results).
%            :28-June-2024:    Tweaked for Matlab 2024
%
% Input data:
%    F                - Frequency array.
%    S                - Complex S-parameter data array.
%    N                - Number of points.
%    EstimatedResFreq - Estimate for resonant frequency.
%    Qseed            - Estimate for Q-factor QL (it is usually sufficient
%                       to assume Qseed = 1000).
%
% Output data (MAT 58 eqn. 17:
%    sv = [a', a'', b', b'', QL]
%

N2 = 2*N;
X = zeros(5, N2);   % X is the transpose of M
G = zeros(N2, 1);
P = zeros(1, N2);

for i = 1:N
    i2 = i+N;
    t = 2.0*(F(i)/EstimatedResFreq - 1.0);
    y = abs(1.0/(1.0 + 1j*Qseed*t));
    P(i) = y*y;
    P(i2) = P(i);
    X(1,i) = t;  X(1,i2) = 0.0;
    X(2,i) = 0.0; X(2,i2) = t;
    X(3,i) = 1.0;  X(3,i2) = 0.0;
    X(4,i) = 0.0; X(4,i2) = 1.0;
    v1 = S(i);
    G(i) = real(v1);    G(i2) = imag(v1);
    v2 = v1*t;
    X(5,i) = imag(v2);  X(5,i2) = -real(v2);
end

M = X.';
T = X.*P;
C = T*M;
q = T*G;
sv = C\ q;

