%Check autocorrelation

clear; close all; clc;

format ('compact');
format ('long', 'g');

%--- Include folders with functions ---------------------------------------
addpath include             % The software receiver functions
addpath geoFunctions        % Position calculation related functions


CodePhase =-1:0.001:1; %[chips]
caCode = generateCAcode(1);
%caCode = [caCode caCode caCode];

codeFreq        = 1.023e6;
earlyLateSpc    = 0.5;
Td              = 0.001;

idxSamplingFreq     = 1;
samplingFreqRange   = [4]*1.023e6;
varDLLError         = size(1,length(samplingFreqRange));
maxDLLError         = varDLLError;
minDLLError         =varDLLError;
for samplingFreq    = samplingFreqRange;

    codePhaseStep   = codeFreq / samplingFreq;
    blksize         = Td*samplingFreq;
    tcode           = 0:codePhaseStep : ((blksize-1)*codePhaseStep);
    %tcode2          =ceil(tcode)+1023;
    tcode2          = floor(rem(tcode+1023,1023)) + 1;
    incomingSignal   = caCode(tcode2);

    I_E=size(1,length(CodePhase));
    I_P=I_E;
    I_L=I_E;
    idxCodePhase=1;
    for remCodePhase=CodePhase
        tcode       = (remCodePhase-earlyLateSpc) : ...
                                  codePhaseStep : ...
                                  ((blksize-1)*codePhaseStep+remCodePhase-earlyLateSpc);
        tcode2      = floor(rem(tcode+1023,1023)) + 1;
        %tcode2          =ceil(tcode)+1023;
        earlyCode   = caCode(tcode2);

        % Define index into late code vector
        tcode       = (remCodePhase+earlyLateSpc) : ...
                      codePhaseStep : ...
                      ((blksize-1)*codePhaseStep+remCodePhase+earlyLateSpc);
        tcode2      = floor(rem(tcode+1023,1023)) + 1;
        %tcode2          =ceil(tcode)+1023;
        lateCode    = caCode(tcode2);

        % Define index into prompt code vector
        tcode       = remCodePhase : ...
                      codePhaseStep : ...
                      ((blksize-1)*codePhaseStep+remCodePhase);
        tcode2      = floor(rem(tcode+1023,1023)) + 1;
        %tcode2          =ceil(tcode)+1023;
        promptCode  = caCode(tcode2);

        I_P(idxCodePhase)=sum(promptCode.*incomingSignal)/blksize;
        I_L(idxCodePhase)=sum(lateCode.*incomingSignal)/blksize;
        I_E(idxCodePhase)=sum(earlyCode.*incomingSignal)/blksize;

        idxCodePhase  =idxCodePhase+1;

        %Code Error
    end;
    figure(300);
    plot(CodePhase,I_P,'r','LineWidth',2);grid on;hold on;
    disp(sprintf('Processing Fs =%d',samplingFreq));
    codeError=(I_E-I_L)./(2*(I_E+I_L))- CodePhase;
    codeError2=codeError(find(CodePhase>=-0.5 & CodePhase<=0.5)); 
    minDLLError(idxSamplingFreq)         = min(abs(codeError2));
    maxDLLError(idxSamplingFreq)         = max(codeError2);
    varDLLError(idxSamplingFreq)         = var(abs(codeError2));
    idxSamplingFreq                      = idxSamplingFreq+1;
end;

figure(100);
plot(samplingFreqRange/1.023e6,sqrt(varDLLError)*299792458/codeFreq,'-o');grid on;