function [channel, trackResults] = trackingLoops( settings, channel, ...
                                                  trackResults )

%==============================================================================%
%                                                                              %
%   [channel, trackResults] = trackingLoops( settings, channel, trackResults ) %
%                                                                              %
%   Perform code and carrier tracking loops                                    %
%   David S. De Lorenzo                                                        %
%   28 November 2005                                                           %
%                                                                              %
%------------------------------------------------------------------------------%
%                                                                              %
%   The tracking modes are as follows:                                         %
%      carrier-phase tracking:  frequency capture (FLL mode)                   %
%                               phase alignment (pull-in mode)                 %
%                               bit alignnemt (transition mode)                %
%                               steady-state tracking (fine mode)              %
%        code-phase tracking:   initial pull-in (pull-in mode)                 %
%                               coherent refinement (transition mode)          %
%                               steady-state tracking (fine mode)              %
%                                                                              %
%   The discriminators for the various modes are presented by Phil Ward        %
%   in Ch. 5, "Satellite Signal Acquisition, Tracking, and Data                %
%   Demodulation" of Understanding GPS: Principles and Applications,           %
%   2nd ed., 2006.                                                             %
%                                                                              %
%      Frequency capture utilizes a modified arctangent discriminator          %
%      as shown in Table 5.4 on page 171.  Note that the four-quadrant         %
%      arctangent disriminator shown in that table is in error; it             %
%      instead should read "atan2( cross, dot ) / (t2 - t1)".  The FLL         %
%      discriminator used here actually is the two-quadrant arctangent         %
%      (note that sample time (t2-t1) automatically is accounted for in        %
%      computation of the loop filter coefficients):                           %
%          PLL_discrRad = atan( cross / dot )                                  %
%                                                                              %
%      Carrier-phase tracking utilizes a Costas Loop discriminator,            %
%      the two-quadrant arctangent discriminator as shown in Table 5.3         %
%      on page 168:                                                            %
%          PLL_discrRad = atan( Q_prompt / I_prompt )                          %
%                                                                              %
%      Non-coherent code-phase tracking utilizes the normalized                %
%      early-minus-late envelope discriminator as shown in Table 5.5           %
%      on page 174:                                                            %
%          DLL_discr = mean( 1/2 * (E - L) ./ (E + L) )                        %
%                                                                              %
%      Coherent code-phase tracking utilizes the normalized dot product        %
%      discriminator as shown in Table 5.5 on page 174.  Note that the         %
%      discriminator shown in that table is in error; it instead should        %
%      read "1/2 * (I_ES - I_LS) / I_PS":                                      %
%          DLL_discr = 1/2 * (E - L) / P                                       %
%                                                                              %
%==============================================================================%
%                       GPS Software Receiver in Matlab                        %
%     Copyright (c) 2003-2017 David S. De Lorenzo / 2018-2023 Trimble Inc.     %
%==============================================================================%

%--- Assign temporary variables, making code more readable ---------------------
msCount          = channel.msCount;
fllTimeMs        = settings.trk.fllTimeMs;
pullInTimeMs     = settings.trk.pullInTimeMs;
transitionTimeMs = settings.trk.transitionTimeMs;

%--- Store values from the early, prompt, & late code correlators --------------
trackResults.I_early(msCount)  = real(channel.correlators1.Early);
trackResults.I_prompt(msCount) = real(channel.correlators1.Prompt);
trackResults.I_late(msCount)   = real(channel.correlators1.Late);
trackResults.Q_early(msCount)  = imag(channel.correlators1.Early);
trackResults.Q_prompt(msCount) = imag(channel.correlators1.Prompt);
trackResults.Q_late(msCount)   = imag(channel.correlators1.Late);

%--- Store the detailed correlation function -----------------------------------
if (settings.trk.plotCorrelationPeak == true)
    trackResults.correlationFunction(:,msCount) = channel.correlationFunction1;
end

%===============================================================================
%   Set the loop filter coefficients for the current tracking mode.            %
%     The filter formula:                                                      %
%       y(n) = Kp * x(n) + Ki * [x(n) + x(n-1) + ... + x(1) + x(0)]            %
%===============================================================================
if (msCount <= fllTimeMs)
    %--- FLL mode --------------------------------------------------------------
    PLL.Kp              = settings.PLL.F_Kp;
    PLL.Ki              = settings.PLL.F_Ki;
    PLL.integrationTime = settings.PLL.F_integrationTime;
    DLL.Kp              = settings.DLL.P_Kp;
    DLL.Ki              = settings.DLL.P_Ki;
    DLL.integrationTime = settings.DLL.P_integrationTime;
    DLL.discrType       = 'non-coherent';
elseif (msCount <= pullInTimeMs)
    %--- Pull-in mode ----------------------------------------------------------
    PLL.Kp              = settings.PLL.P_Kp;
    PLL.Ki              = settings.PLL.P_Ki;
    PLL.integrationTime = settings.PLL.P_integrationTime;
    DLL.Kp              = settings.DLL.P_Kp;
    DLL.Ki              = settings.DLL.P_Ki;
    DLL.integrationTime = settings.DLL.P_integrationTime;
    DLL.discrType       = 'non-coherent';
elseif (msCount <= transitionTimeMs)
    %--- Transition mode -------------------------------------------------------
    PLL.Kp              = settings.PLL.X_Kp;
    PLL.Ki              = settings.PLL.X_Ki;
    PLL.integrationTime = settings.PLL.X_integrationTime;
    DLL.Kp              = settings.DLL.X_Kp;
    DLL.Ki              = settings.DLL.X_Ki;
    DLL.integrationTime = settings.DLL.X_integrationTime;
    DLL.discrType       = 'coherent';
else
    %--- Fine-tracking mode ----------------------------------------------------
    PLL.Kp              = settings.PLL.T_Kp;
    PLL.Ki              = settings.PLL.T_Ki;
    PLL.integrationTime = settings.PLL.T_integrationTime;
    DLL.Kp              = settings.DLL.T_Kp;
    DLL.Ki              = settings.DLL.T_Ki;
    DLL.integrationTime = settings.DLL.T_integrationTime;
    DLL.discrType       = 'coherent';
end

%--- Find the navigation data bit transition index -----------------------------
%    This happens in the first millisecond of the transition mode.
if (msCount == (pullInTimeMs + 1))
    
    %=== Find the navigation data bit transition index =========================
    %    Pass zeros for blocks during the FLL mode - this is
    %    clearer to read (although trivially slower) than just
    %    passing I_prompt values during pull-in mode and then
    %    adjusting navStartBit to compensate.
    channel.navStartBit = ...
        findNavStartBit(settings.msPerDataBit, settings.secondarySeq, ...
                        [zeros(1, fllTimeMs), ...
                         trackResults.I_prompt(fllTimeMs+1:msCount)]);
    
	%--- Set the data bit counter for the current millisecond ------------------
    channel.databitCounter = ...
        mod(msCount - channel.navStartBit + 1, settings.msPerDataBit) + 1;
	
    %--- Save the index to the tracking results structure ----------------------
    trackResults.navStartBit = channel.navStartBit;
    
end

%--- At the tracking mode transition boundaries, reset the integrated ----------
%    discriminator totals and update the NCO center frequencies.
%    The code generator NCO center frequency is always set so that
%    code doppler is essentially derived from the carrier doppler
%    frequency estimate, and the code tracking loop only needs to
%    track the ionospheric divergence term.
if ((msCount == (fllTimeMs + 1)) ...
 || (msCount == (pullInTimeMs + 1)) ...
 || (msCount == (transitionTimeMs + 1)))
    
    %--- Reset integrated discriminator totals ---------------------------------
    channel.DLL_discrIntegrated = 0;
    channel.PLL_discrIntegrated = 0;
    
    %--- Update the NCO center frequencies -------------------------------------
    channel.carrierFreqNCO = channel.carrierFreq;
    channel.dopplerFreq    = channel.carrierFreq - settings.intermediateFreq;
    channel.codeFreqNCO    = channel.codeFreq ...
                             - (channel.dopplerFreq * settings.carrierAidingSF);
    
    %--- Also set the PLL and DLL counters so that -----------------------------
    %    loop entrance occurs on schedule.
    channel.pllCounter = ...
        mod(channel.databitCounter - 2, PLL.integrationTime) + 1;
    channel.dllCounter = ...
        mod(channel.databitCounter - 2, DLL.integrationTime) + 1;
    
end

%===============================================================================
%                                                                              %
%  Phase-Locked Loop (PLL)                                                     %
%                                                                              %
%===============================================================================

%--- Check if the PLL averaging period is over ---------------------------------
if (channel.pllCounter < PLL.integrationTime)
    
    %--- Increment the PLL counter ---------------------------------------------
    channel.pllCounter = channel.pllCounter + 1;
    
else
    
    %--- Re-set the PLL counter ------------------------------------------------
    channel.pllCounter = 1;
    
    %--- For the first couple of millisecond accumulations ---------------------
    %    don't try frequency/phase correction.
    if (msCount <= 2)
        
        %--- Set the PLL discriminator to zero ---------------------------------
        PLL_discrRad = 0;
        
    %--- This is to compute the FLL discriminator ------------------------------
    elseif (msCount <= fllTimeMs)
        
        %--- Cross term --------------------------------------------------------
        %    c = ( I_1 * Q_2 ) - ( Q_1 * I_2 )
        c = trackResults.I_prompt(msCount-1)*trackResults.Q_prompt(msCount) ...
            - trackResults.Q_prompt(msCount-1)*trackResults.I_prompt(msCount);
        
        %--- Dot term ----------------------------------------------------------
        %    d = ( I_1 * I_2 ) + ( Q_1 * Q_2 )
        d = trackResults.I_prompt(msCount-1)*trackResults.I_prompt(msCount) ...
            + trackResults.Q_prompt(msCount-1)*trackResults.Q_prompt(msCount);
        
        %--- FLL discriminator is the 2-quadrant arctangent --------------------
        PLL_discrRad = atan(c / d);
        
    %--- This is to compute the PLL discriminator ------------------------------
    else
        
        %--- Index of samples to use in the discriminator ----------------------
        idx = msCount + (1-PLL.integrationTime:0);
        
        %--- Accumulate in-phase and quadrature discriminator totals -----------
        I_prompt = sum(trackResults.I_prompt(idx));
        Q_prompt = sum(trackResults.Q_prompt(idx));
        
        %--- Carrier-phase discriminator is the 2-quadrant arctangent ----------
        PLL_discrRad = atan(Q_prompt / I_prompt);
        
    end
    
    %--- Convert from units of radians in the discriminator ouput to cycles ----
    channel.PLL_discr = PLL_discrRad / (2*pi);
    
    %--- Integrate the discriminator output from t=0 to the current epoch ------
    channel.PLL_discrIntegrated = channel.PLL_discrIntegrated ...
                                  + channel.PLL_discr;
    
    %--- PLL filter update -----------------------------------------------------
    %    The filter formula:
    %      y(n) = Kp * x(n) + Ki * [x(n) + x(n-1) + ... + x(1) + x(0)]
    channel.carrierFreqCorrection = PLL.Kp * channel.PLL_discr ...
                                    + PLL.Ki * channel.PLL_discrIntegrated;
    
    %--- Update the carrier frequency ------------------------------------------
    channel.carrierFreq = channel.carrierFreqNCO ...
                          + channel.carrierFreqCorrection;
    
    %--- Update the doppler frequency ------------------------------------------
    channel.dopplerFreq = channel.carrierFreq ...
                          - settings.intermediateFreq;
    
    %--- Calculate the nominal I & Q mixing vector (i.e., zero carrier-phase) --
    %    This mixing vector is valid, once the appropriate phase rotation
    %    is applied in trackingDriver, for the whole PLL integration
    %    period - i.e., the carrier frequency will only be updated
    %    every 20 ms during coherent tracking modes.
    channel.IQ = exp(-1i * (channel.carrierFreq * settings.phasePts));
    
end % PLL averaging period

%===============================================================================
%                                                                              %
%  Delay-Locked Loop (DLL)                                                     %
%                                                                              %
%===============================================================================

%--- Check if the DLL averaging period is over ---------------------------------
if (channel.dllCounter < DLL.integrationTime)
    
    %--- Increment the DLL counter ---------------------------------------------
    channel.dllCounter = channel.dllCounter + 1;
    
else
    
    %--- Re-set the DLL counter ------------------------------------------------
    channel.dllCounter = 1;
    
    %--- Index of samples to use in the discriminator --------------------------
    idx = msCount + (1-DLL.integrationTime:0);
    
    switch DLL.discrType
        
        case 'non-coherent'
            
            %--- Non-coherent code-phase discriminator -------------------------
            E = hypot(trackResults.I_early(idx), ...
                      trackResults.Q_early(idx));
            L = hypot(trackResults.I_late(idx), ...
                      trackResults.Q_late(idx));
            
            %--- Normalized early minus late envelope --------------------------
            channel.DLL_discr = mean(1/2 * (E - L) ./ (E + L));
            
        case 'coherent'
            
            %--- Coherent code-phase discriminator -----------------------------
            E = sum(trackResults.I_early(idx));
            P = sum(trackResults.I_prompt(idx));
            L = sum(trackResults.I_late(idx));
            
            %--- Normalized coherent dot product discriminator -----------------
            channel.DLL_discr = 1/2 * (E - L) / P;
            
        otherwise
            
            error('***  Invalid discriminator type!!  ***');
            
    end
    
    %--- Integrate the discriminator output from t=0 to the current epoch ------
    channel.DLL_discrIntegrated = channel.DLL_discrIntegrated ...
                                  + channel.DLL_discr;
    
    %--- DLL filter update -----------------------------------------------------
    %    The filter formula:
    %      y(n) = Kp * x(n) + Ki * [x(n) + x(n-1) + ... + x(1) + x(0)]
    channel.codeFreqCorrection = DLL.Kp * channel.DLL_discr ...
                                 + DLL.Ki * channel.DLL_discrIntegrated;
    
    %--- Update the code frequency ---------------------------------------------
    %    Include carrier-aiding of the code tracking loop.
    channel.codeFreq = channel.codeFreqNCO ...
                       + channel.codeFreqCorrection ...
                       + (channel.dopplerFreq * settings.carrierAidingSF);
    
    %--- Calculate the nominal code index vector (i.e., zero code-phase) -------
    channel.codePts = channel.codeFreq * settings.timePts;
    
end % DLL averaging period

%===============================================================================
%                                                                              %
%  Save results of the tracking loops                                          %
%                                                                              %
%===============================================================================

%--- Code and carrier parameters -----------------------------------------------
trackResults.msCount                      = msCount;
trackResults.codePhase(msCount)           = channel.codePhase;
trackResults.carrierPhase(msCount)        = channel.carrierPhase;
trackResults.carrierPhaseInt(msCount)     = channel.carrierPhaseInt;
trackResults.codeFreq(msCount)            = channel.codeFreq;
trackResults.carrierFreq(msCount)         = channel.carrierFreq;
trackResults.dopplerFreq(msCount)         = channel.dopplerFreq;

%--- Loop discriminators -------------------------------------------------------
trackResults.DLL_discr(msCount)           = channel.DLL_discr;
trackResults.DLL_discrIntegrated(msCount) = channel.DLL_discrIntegrated;
trackResults.PLL_discr(msCount)           = channel.PLL_discr;
trackResults.PLL_discrIntegrated(msCount) = channel.PLL_discrIntegrated;

return
%===============================================================================