function mv = OptimiseFit6(F, S, N, EstimatedResFreq, sv, ...
                          loop_plan, Tol, quiet)
%
% NLQFIT6 Step (2) Matlab/GNU Octave.
%   Optimised fit of Q-factor (QL) and resonant frequency (FL)
%   by the gradient-descent method.
%
% 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.
%
% Input data:
%    F                - Frequency array.
%    S                - Complex S-parameter data array.
%    N                - Number of points.
%    EstimatedResFreq - Estimate for resonant frequency.
%    sv               - Solution vector found by initial linear fit.
%    loop_plan        - String of characters which defines order of
%                       steps used by the fitting process.
%                       E.g. 'fwfwc'.
%                       f - fit once without testing for convergence.
%                       c - repeated fit, iterating until convergence is
%                           obtained.
%                       w - re-calculate weighting factors on basis of
%                           previous fit.
%    Tol              - Criterion for the convergence test.
%                       Recommend using 1.0E-5 for reflection
%                       or max(abs(S))*1.0E-5 for transmission.
%    quiet            - Boolean flag controlling output of information to
%                       the console.
%
% Output data (MAT58 eqn. 22):
%    mv = [m1, m2, m4, m4, QL, FL]
%

N2 = N*2;
G = ([1; N2]);
PV = ones(N,1); % column vector
PA = ones(N2,1);
P = eye(N2);
m1 = sv(2)/sv(5);
m2 = -sv(1)/sv(5);
m3 = sv(3)-m1;
m4 = sv(4)-m2;
m5 = sv(5);
Flwst = F(1);   % lowest freq. is a convenient normalisation factor.
m6 = Flwst*m5/EstimatedResFreq;

for k = 1:length(loop_plan)

    op = loop_plan(k);
    if (op == 'w')         % Fr,        QL
        PV = AngularWts(F, Flwst*m5/m6, m5);
        PA = vertcat(PV,PV); % column vector
        P = diag(PA);
        if (~quiet)
                fprintf('Op w,  Calculating weighting factors\n')
        end
        continue
    elseif (op == 'c')
        seek_convergence = 1;
    elseif (op == 'f')
        seek_convergence = 0;
    else
        disp('Unexpected character in loop_plan')
    end

    TerminationConditionMet = 0;
    Valid_RMS_Error = 0;
    Valid_Last_RMS_Error = 0;
    while (~TerminationConditionMet)
        X = zeros(6, N2); % X is transpose of matrix M
        for i = 1:N
            i2 = i+N;
            y = 1.0/(1.0 + 2j*(m6*F(i)/Flwst-m5));
            X(1,i) = 1.0;      X(1,i2) = 0.0;
            X(2,i) = 0.0;      X(2,i2) = 1.0;
            X(3,i) = real(y);  X(3,i2) = imag(y);
            X(4,i) = -imag(y); X(4,i2) = real(y);
            u = (1j*m3-m4)*y*y*2;
            X(5,i) = real(u);  X(5,i2) = imag(u);
            u = -u*F(i)/Flwst;
            X(6,i) = real(u);  X(6,i2) = imag(u);
            v = m1+1j*m2 + (m3+1j*m4)*y;
            r = S(i) - v;              % residual
            G(i) = real(r);    G(i2) = imag(r);
        end
        M = X.';
        T = X*P; % Faster alternative: T = bsxfun(@times,X,PA');
        C = T*M;
        q = T*G;
        dm = C\q;         % solve for changes in m1..m6
        m1 = m1+dm(1);
        m2 = m2+dm(2);
        m3 = m3+dm(3);
        m4 = m4+dm(4);
        m5 = m5+dm(5);
        m6 = m6+dm(6);

        if (Valid_RMS_Error)
            Last_RMS_Error = RMS_Error;
            Valid_Last_RMS_Error = 1;
        end

        SumNum = 0.0;
        SumDen = 0.0;

        for i=1:N
            den = (1.0 + 2j*(m6*F(i)/Flwst-m5));
            E = S(i) - (m1+1j*m2) - (m3+1j*m4)/den;
            ip = PV(i);
            SumNum = SumNum + ip*(real(E)*real(E) + imag(E)*imag(E));
            SumDen = SumDen + ip;
        end
        RMS_Error = sqrt(SumNum/SumDen);
        Valid_RMS_Error = 1;

        if (~quiet)
          fprintf('RMS Error %8.6f, QL %5.2f, FL %9.6f\n', ...
              RMS_Error,m5,Flwst*m5/m6)
        end

        if (seek_convergence)
           if (Valid_Last_RMS_Error)
              delta_S = abs(RMS_Error-Last_RMS_Error);
              TerminationConditionMet = (delta_S < Tol);
           end
        else
           TerminationConditionMet = 1;
        end


     end % while loop
end  % for k loop

% Solution vector.
QL = m5;
FL = m5*Flwst/m6;
mv = [m1, m2, m3, m4, QL, FL];
