function mv = OptimiseFit7(F, S, N, EstimatedResFreq, sv, ...
                          loop_plan, Tol, quiet)
%
% NLQFIT7 Step (2)  Matlab/GNU Octave.
%   Optimised fit of Q-factor (QL), resonant frequency (FL)
%   and a phase term related to the length of uncalibrated line.
%   Uses 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. 26):
%    mv = [m1, m2, m4, m4, QL, FL, m7/Flwst]
%

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;
m7 = 0.0;

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(7, 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));
            fdn = F(i)/Flwst - m5/m6;
            pj = 1j*m7*fdn;
            expm7 = exp(pj);
            ym = y*expm7;
            X(1,i) = real(expm7);     X(1,i2) = imag(expm7);
            X(2,i) = -imag(expm7);    X(2,i2) = real(expm7);
            X(3,i) = real(ym);        X(3,i2) = imag(ym);
            X(4,i) = -imag(ym);       X(4,i2) = real(ym);
            u = (1j*m3-m4)*y*ym*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);
            fL = Flwst*m5/m6;
            v = ((m1+1j*m2) + (m3+1j*m4)*y)* expm7;
            u = v*fdn;
            X(7,i)=-imag(u);          X(7,i2)=real(u);
            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..m7
        m1 = m1+dm(1);
        m2 = m2+dm(2);
        m3 = m3+dm(3);
        m4 = m4+dm(4);
        m5 = m5+dm(5);
        m6 = m6+dm(6);
        m7 = m7+dm(7);

        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
            fdn = F(i)/Flwst - m5/m6;
            den = 1.0 + 2j*(m6*F(i)/Flwst-m5);
            pj = 1j*m7*fdn;
            E = S(i) - ((m1+1j*m2) + (m3+1j*m4)/den)*exp(pj);
            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, m7/Flwst];
