function acqResults = acquisitionRefineDoppler( fidIn, settings, acqResults )

%==============================================================================%
%                                                                              %
%   acqResults = acquisitionRefineDoppler( fidIn, settings, acqResults )       %
%                                                                              %
%   Refine acquisition Doppler frequency estimate                              %
%   David S. De Lorenzo                                                        %
%   17 July 2008                                                               %
%                                                                              %
%==============================================================================%
%                       GPS Software Receiver in Matlab                        %
%     Copyright (c) 2003-2017 David S. De Lorenzo / 2018-2023 Trimble Inc.     %
%==============================================================================%

%--- Assign temporary variables, making code more readable ---------------------
msToProcessCoherent  = settings.acq.msToProcessCoherent;
coherentBlockStartMs = settings.acq.coherentBlockStartMs;
maxCarrierPhaseError = settings.acq.maxCarrierPhaseError;
dopplerSearchRange   = settings.acq.dopplerSearchRange;
samplesPer1ms        = settings.samplesPer1ms;
chipsPer1ms          = settings.chipsPer1ms;

%--- Samples to process for Doppler frequency refinement -----------------------
samplesToProcess = samplesPer1ms * msToProcessCoherent;

%--- Doppler frequency search spacing in Hz ------------------------------------
dopplerSearchStep = 2*1000*maxCarrierPhaseError / msToProcessCoherent;

%--- Range over which to do frequency refinement -------------------------------
dopplerRefinementRange = min([dopplerSearchRange/2, dopplerSearchStep]);

%--- Number of frequency bins --------------------------------------------------
numberOfDopplerBins = 21;

%--- Nominal Doppler frequency search vector -----------------------------------
dopplerSearchNominal = linspace(-dopplerRefinementRange / 2, ...
                                 dopplerRefinementRange / 2, ...
                                numberOfDopplerBins);

%=== Acquisition - loop over each satellite in the search vector ===============
fprintf('\nRefining Doppler frequency estimates ...\n');
for ii = 1:numel(acqResults)
    
    prn = acqResults(ii).prn;
    
    %--- Start message ---------------------------------------------------------
    fprintf('  --> PRN %d\n', prn);
    
    %--- PRN code chip sequence, values of +1 and -1 ---------------------------
    codeChips = codeGenerateMASTER(settings.signalType, prn);
    
    %--- Repeat the last chip to protect against index roundoff errors ---------
    codeChips = [codeChips codeChips(end)];                          %#ok<AGROW>
    
    %--- Doppler frequency search vector ---------------------------------------
    dopplerSearchVector = acqResults(ii).dopplerFreq + dopplerSearchNominal;
    
    %--- Dimension accumulator output storage matrix ---------------------------
    accumulatorOutput = zeros(length(coherentBlockStartMs), ...
                              numberOfDopplerBins);
    
    %=== Acquisition over (possibly) several time offsets ======================
    for jj = 1:length(coherentBlockStartMs)
        
        %--- Time offset for this iteration ------------------------------------
        msToSkip = settings.msToSkip + coherentBlockStartMs(jj);
        
        %--- Absolute time/phase of each sample point --------------------------
        timePts  = (0:samplesToProcess-1) / settings.samplingFreq;
        phasePts = 2 * pi * timePts;
        phasePts = mod(phasePts, 2*pi);
        
        %--- Seek forward to the specified file position pointer location ------
        placeFilePositionPointer(fidIn, msToSkip, ...
                                 settings.extraSamplesToSkip, settings);
        
        %=== Read data from file ===============================================
        data = readDataFromFile(fidIn, samplesToProcess, settings);
        
        %=== Loop over Doppler frequency =======================================
        for dopplerSearchBin = 1:numberOfDopplerBins
            
            %--- Set signal frequencies ----------------------------------------
            dopplerFreq = dopplerSearchVector(dopplerSearchBin);
            carrierFreq = settings.intermediateFreq + dopplerFreq;
            codeFreq    = settings.codeFreq ...
                          + dopplerFreq * settings.carrierAidingSF;
            
            %--- Remove the carrier from the signal (complex) ------------------
            %    Carrier wipeoff moves the signal to baseband.
            dataBaseband = data .* exp(-1i * carrierFreq * phasePts);
            
            %--- Construct the code replica ------------------------------------
            replicaCodeSequence = ...
                codeChips(floor(mod((codeFreq * timePts ...
                                     - acqResults(ii).codePhase), ...
                                    chipsPer1ms)) + 1);
            
            %--- Code wipeoff and accumulation (complex) -----------------------
            accumulatorOutput(jj, dopplerSearchBin) = ...
                dataBaseband * replicaCodeSequence.';
            
        end % loop over Doppler frequency
        
    end % loop over time offset
    
    %--- Look for the maximum power (complex --> real) -------------------------
    accumulatorOutput    = abs(accumulatorOutput);
    [~,rowNumber]        = max(max(accumulatorOutput, [], 2));
    [~,dopplerResultBin] = max(accumulatorOutput(rowNumber,:));
    
    %--- Store the Doppler frequency estimate ----------------------------------
    acqResults(ii).dopplerFreq = dopplerSearchVector(dopplerResultBin);
    
end % loop over satellite number

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