function [Epoch, Adjust] = NL_fixing_IF(Epoch, Adjust, b_WL, b_NL, elevs, settings)
% This function fixes the NL abiguities for the ionosphere-free linear 
% combination. The WL ambiguities have already been fixed in the function
% WL_fixing.m (using the HMW LC).
% 
% INPUT:
%   Epoch    	epoch-specific data for current epoch [struct]
%   Adjust    	adjustment data and matrices for current epoch [struct]
%   b_WL        WL UPDs for satellites of current epoch, single-differenced
%   b_NL        NL UPDs for satellites of current epoch, single-differenced
%   elevs   	elevations of all satellites of this epoch [°]
%   settings 	struct, settings for processing (from GUI)
% OUTPUT:
%   Epoch     	updated with (integer fixed) NL Ambiguities [struct]
%
% Revision:
%   2023/06/11, MFWG: adding QZSS
%   2024/12/30, MFWG: switching to LAMBDA 4.0
%   2025/03/31, MFWG: convert Q_NN from meters to cycles
%   2025/10/17, MFWG: QZSS fixing
% 
% This function belongs to raPPPid, Copyright (c) 2023, M.F. Glaner
% *************************************************************************



%% Preparations
% get some variables
NO_PARAM = Adjust.NO_PARAM;     % number of estimated parameters
f1 = Epoch.f1;                  % frequency 1st processed frequency
f2 = Epoch.f2;              	% frequency 2nd processed frequency
l1 = Epoch.l1;                  % wavelength 1st processed frequency
l2 = Epoch.l2;                  % wavelength 2nd processed frequency
sats = Epoch.sats;              % satellites of current epoch
no_sats = numel(Epoch.sats);    % number of satellites
% number of GNSS satellites
no_gps  = sum(Epoch.gps);
no_glo  = sum(Epoch.glo);
no_gal  = sum(Epoch.gal);
no_bds  = sum(Epoch.bds);
no_qzs  = sum(Epoch.qzss); 
% reference satellites of GPS, Galileo, BeiDou, QZSS
refSatG_idx = Epoch.refSatGPS_idx;
refSatE_idx = Epoch.refSatGAL_idx;
refSatC_idx = Epoch.refSatBDS_idx;
refSatJ_idx = Epoch.refSatQZS_idx;
% is satellite fixable?
fixable = Epoch.fixable;

% get float ambiguities
idx = (NO_PARAM+1):(NO_PARAM+no_sats);
param_N = Adjust.param(idx);    % float ambiguities
param_N_gps = param_N(Epoch.gps);       % GPS float ambiguities
param_N_gal = param_N(Epoch.gal);       % Galileo float ambiguities
param_N_bds = param_N(Epoch.bds);       % BeiDou float ambiguities
param_N_qzs = param_N(Epoch.qzss);      % QZSS float ambiguities

% get covariance matrix and convert from meters to cycles
Q_NN = Adjust.param_sigma(idx, idx);    % covariance matrix of float ambiguities
lam = (f1.^2 .* l1 - f2.^2 .* l2) ./ (f1.^2 - f2.^2);   % wavelength of IF LC ||| check!
Q_NN = Q_NN ./ lam;      % divide rows by wavelength
Q_NN = Q_NN ./ lam';     % divide columns by wavelength

% remove GLONASS and QZSS (if processed in float solution)
if settings.INPUT.use_GLO
    % GLONASS and QZSS satellites are not needed in the following variables
    is_glo = Epoch.glo;
    sats(is_glo) = [];
    Q_NN(is_glo, :) = [];
    Q_NN(:, is_glo) = [];
    f1(is_glo) = [];
    f2(is_glo) = [];
    l1(is_glo) = [];
    b_WL(is_glo) = [];
    b_NL(is_glo) = [];
    fixable(is_glo) = [];
end


%% Prepare solving the NL ambiguities
% Single-difference covariance matrix
C = eye(no_gps+no_gal+no_bds+no_qzs);
% Create matrix C for covariance propagation of covariance matrix
C = -C;
if settings.INPUT.use_GPS
    C(Epoch.gps,refSatG_idx) = 1;
    C(refSatG_idx,refSatG_idx) = 0;
end
if settings.INPUT.use_GAL
    bool_gal = Epoch.gal;
    bool_gal(Epoch.glo) = [];   % for GLO+GAL processing
    C(bool_gal,refSatE_idx-no_glo) = 1;
    C(refSatE_idx-no_glo,refSatE_idx-no_glo) = 0;
end
if settings.INPUT.use_BDS
    bool_bds = Epoch.bds;
    bool_bds(Epoch.glo) = [];   % for GLO+BDS processing
    C(bool_bds,refSatC_idx-no_glo) = 1;
    C(refSatC_idx-no_glo,refSatC_idx-no_glo) = 0;
end
if settings.INPUT.use_QZSS
    bool_qzs = Epoch.qzss;
    bool_qzs(Epoch.glo) = [];   % for GLO+QZSS processing
    C(bool_qzs,refSatJ_idx-no_glo) = 1;
    C(refSatJ_idx-no_glo,refSatJ_idx-no_glo) = 0;
end
% calculate covariance propagation
Q_NN_SD = C*Q_NN*C';  	% Variance-Covariance-Matrix of SD float ambiguities

% Single difference ambiguities
if ~isempty(refSatG_idx)
    param_N_gps = param_N_gps(refSatG_idx) - param_N_gps;
end
if ~isempty(refSatE_idx)
    param_N_gal = param_N_gal(refSatE_idx-no_gps-no_glo) - param_N_gal;
end
if ~isempty(refSatC_idx)
    param_N_bds = param_N_bds(refSatC_idx-no_gps-no_glo-no_gal) - param_N_bds;
end
if ~isempty(refSatJ_idx)
    param_N_qzs = param_N_qzs(refSatJ_idx-no_gps-no_glo-no_gal-no_bds) - param_N_qzs;
end
param_N_SD = [param_N_gps; param_N_gal; param_N_bds; param_N_qzs];    % SD float ambiguities


% calculate NL ambiguities from the float IF ambiguities
NL_float = param_N_SD./l1 .* (f1+f2)./f1 - Epoch.WL_12(sats) .* (f2 ./ (f1-f2));  % [00]: (4.17)
NL_float = NL_float + b_NL;


% remove reference satellite and satellites which are not fixable
rem = isnan(NL_float);
refSatsIdx = [refSatG_idx, refSatE_idx-no_glo, refSatC_idx-no_glo, refSatJ_idx-no_glo];
rem(refSatsIdx) = true;
rem(~fixable) = true;
% delete float estimation and row&column in covariance matrix
NL_float(rem) = [];
Q_NN_SD(rem, :) = [];
Q_NN_SD(:,rem) = [];


%% Solve all Narrow-Lane ambiguities at once with LAMBDA method
if isempty(NL_float) || numel(NL_float) < 2
    % NL can be fixed because not enough fixable float ambiguities
    Epoch.NL_12(sats) = NaN;
    return
end

% integer fixing with LAMBDA
[afixed, sqnorm] = LAMBDA(NL_float, Q_NN_SD, 5, 3, DEF.AR_THRES_SUCCESS_RATE);


% take best solution
N_best = afixed(:,1);
% take only those ambiguities which are integer
bool_int = (N_best - floor(N_best) == 0);
N_int = NaN(numel(bool_int),1);
N_int(bool_int) = N_best(bool_int);
% save results
NL(~rem) = N_int;
NL(rem) = NaN;
NL(refSatsIdx) = 0;
Epoch.NL_12(sats) = NL;
