 close all 
clear all
clc

%definování zdrojového adresáøe pro èerpání dat
s = dir ('./Archiv/*.dat');
pocet_files=size(s,1);

%naètení vstupních dat v požadovaném formátu
i=1:pocet_files;
nazev = s(i).name;
nazev=nazev(1,1:(end-4));
filename = ['./Archiv/',nazev,'.dat'];

delimiter = ';'; %definice znaku rozdìlujícího vstupní data do jednotlivých èástí (segmentù)

formatSpec = '%f%f%s%s%[^\n\r]';

fileID = fopen(filename,'r');
dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'EmptyValue' ,NaN, 'ReturnOnError', false);
fclose(fileID);

%rozdìlení zprávy ve vstupním formátu na 4 segmenty
Unix_time = dataArray{:, 1};
ReceiverID = dataArray{:, 2};
Time = dataArray{:, 3};
Message = dataArray{:, 4};

%definování poètu uvažovaných zpráv ze zdrojového souboru pro další analýzy 
n=320062; %111 000 zpráv cca 1 minuta záznamu ve 3:00 pro 4 pøijímaèe ve špièce 300000 více zpráv 1 minuta, zpìtnì lze definovat podle promìnné T
%definice souøadnic jednotlivých pøijímaèù v systému
ReceiverStrahov=[50.080512805,14.395710655,338.33];
ReceiverHorska=[50.067873582,14.423431681,208.62];
ReceiverTower=[50.050384622,14.436212947,377.21];
ReceiverLKLT=[50.129189,14.525771,285.0];

DF=NaN(n,1);
c=299792458; %rychlost elektromagnetického záøení ve vakuu
r=6371000; %polomìr Zemì
Adress=cell(n,1);
DF17_error=zeros(n,1); 

%nadefinování sledované ICAO adresy odpovídaèe
ID='' % ID vybrané pro analýzu záznamu 1706190300.dat
% ID='3C6488'; % ID vybrané pro analýzu záznamu 1709170640.dat

i=1;

%cyklus pro zjištìní DF, MODE_S_ID a vytvoøení matice zpráv pøijatých od vybraného ID


for j=1:n
 
    m=length(Message{j,1});
    
    if m==14||m==28
    
    MessageBin{j,1}=hex2bin(Message{j,1});
        
    DF(j,1)=polyval(MessageBin{j,1}(1,1:5),2);
    
    if DF(j,1)==17    
       MODE_S_ID(j,1)=string(bin2hex(MessageBin{j,1}(1,9:32)));
       [Adress_DF17HEX]=crc_division4(MessageBin{j,1}); %funkce pro dekódování adresnì paritního pole pro kontrolu pøijatých zpráv typu DF17
       DF17_paritycheck(j,1)=string(Adress_DF17HEX);
        
       if DF17_paritycheck(j,1)~='000000' 
          DF17_error(j,1)=1;
       end
       
    end         
     
    if DF(j,1)==11
       MODE_S_ID(j,1)=string(bin2hex(MessageBin{j,1}(1,9:32)));
    end
    
    if DF(j,1)==4||DF(j,1)==5||DF(j,1)==20||DF(j,1)==21||DF(j,1)==0||DF(j,1)==16
       [AdressHEX]=crc_division4(MessageBin{j,1}); %funkce pro dekódování MODE_S_ID z adresnì paritního pole
       MODE_S_ID(j,1)=string(AdressHEX);
    end
    
    if DF(j,1)==1||DF(j,1)==2||DF(j,1)==3||DF(j,1)==6||DF(j,1)==7||DF(j,1)==8||DF(j,1)==9||DF(j,1)==10||DF(j,1)==12||DF(j,1)==13||DF(j,1)==14||DF(j,1)==15||DF(j,1)==18||DF(j,1)==22||DF(j,1)==23||DF(j,1)==24
       MODE_S_ID(j,1)='AAAAAA'; %jiné DF kódy
    end
    
    if DF(j,1)==19
       MODE_S_ID(j,1)='MMMMMM'; %Military ES
    end
    
    end
    
    E(j,1)=double(isnan(DF(j,1)));
         
    if E(j,1)
       MODE_S_ID(j,1)='EEEEEE'; %chybné zprávy
    end

    if string(MODE_S_ID(j,1))==ID&&DF17_error(j,1)==0 %vytváøení matice Message_ID obsahující zprávy pouze od vybraného letadla a dekódování èasu pøíjmu zprávy pøijímaèem
      %if ReceiverID(j)==1||ReceiverID(j)==2||ReceiverID(j)==3 %podmínka pro vyøazení zpráv od jednotlivých pøijímaèù, vzorové nastavení vyøazuje pøijímaè è.4
          Message_ID{i,2}=Message{j};
          Message_ID{i,3}=DF(j);
          Message_ID{i,4}=ReceiverID(j);
          MessageBin_ID{i,1}=hex2bin(Message_ID{i,2});
       
          TimeBin_ID=hex2bin(Time{j}); %pøevedení promìnné Time do binárního tvaru
    
          TimeBin_sec=(TimeBin_ID(1:18)); %separace èásti promìnné TimeBin, která nese informaci o sekundách
          TimeBin_nano=(TimeBin_ID(19:48));%separace èásti promìnné TimeBin, která nese informaci o nanosekundách
   
          TimeDec_sec=polyval(TimeBin_sec,2);   %pøevod sekund do desítkové soustav
          TimeDec_nano=polyval(TimeBin_nano,2); %pøevod nano do desítkové soustavy
    
          Message_ID{i,5}=TimeDec_sec+(TimeDec_nano.*(10^-9)); %výsledný èas pøíjmu zprávy na konkrétním pøijímaèi
       
          i=i+1;
       %end
           
     end
      
    end
  

n1=length(Message_ID);
j=1;
Poloha=zeros(n1,2);
TimeReceiver=NaN(n1,1);
TimeReceiverAircraft=NaN(n1,1);
TimeAircraft=NaN(n1,1);
Receiver1Time=zeros(n1,4);
Receiver2Time=zeros(n1,4);
Receiver3Time=zeros(n1,4);
Receiver4Time=zeros(n1,4);
DistanceReceiverAircraft=NaN(n1,1);
Poloha_L_ECEF=zeros(n1,3);
r1=1;
r2=1;
r3=1;
r4=1;

%transformace souøadnic pøijímaèù do Kartézského ECEF souøadného systému pomocí funkce WGS2ECEF
[ReceiverStrahov(1,1),ReceiverStrahov(1,2),ReceiverStrahov(1,3)]=WGS2ECEF(50.080512805,14.395710655,338.33);
[ReceiverHorska(1,1),ReceiverHorska(1,2),ReceiverHorska(1,3)]= WGS2ECEF(50.067873582,14.423431681,208.62);  
[ReceiverTower(1,1),ReceiverTower(1,2),ReceiverTower(1,3)]=WGS2ECEF(50.050384622,14.436212947,377.21);
[ReceiverLKLT(1,1),ReceiverLKLT(1,2),ReceiverLKLT(1,3)]=WGS2ECEF(50.129189,14.525771,285.0);

%cyklus pro dekódování polohy, výšky a dopoèítání èasu vyslání zprávy odpovídaèem u zpráv se známou polohou

for i=1:n1
    
 if Message_ID{i,3}==17
    Typecode(i,1)=polyval(MessageBin_ID{i,1}(1,33:37),2);
           
       if Typecode(i,1)>8&&Typecode(i,1)<19 % pokud je to polohový squitter (airborn possition message)
           [Lat,Lon]=DF17_poloha(MessageBin_ID{i,1},50.1152,14.1993); % funkce pro lokální dekódování polohy
           Poloha(i,1)=Lat;
           Poloha(i,2)=Lon;
           
           [Altitude]=DF17_TC9_18_vyska(MessageBin_ID{i,1}); % funkce pro dekódování výšky
           Poloha(i,3)=Altitude*0.3048; %pøepoèet Altitude na metry do promìnné Poloha
            
           [Poloha_L_ECEF(i,1),Poloha_L_ECEF(i,2),Poloha_L_ECEF(i,3)]=WGS2ECEF(Poloha(i,1),Poloha(i,2),Poloha(i,3)); %transformace souøadnic letadla do kartézského souøadného systému ECEF pomocí funkce WGS2ECEF
       
       end
 end
           % rozdìlení zpráv dle pøijímaèe pøíjmu zprávy a výpoèet èasu vyslání zprávy u zpráv se známou polohou
           if Message_ID{i,4}==1
              Receiver1Time(r1,1)=i;
              Receiver1Time(r1,2)=Message_ID{i,5}; % èas pøijetí zprávy
               
               if Poloha_L_ECEF(i,1)~=0
                  DistanceReceiverAircraft(i,1)=sqrt((ReceiverStrahov(1,1)-Poloha_L_ECEF(i,1)).^2+(ReceiverStrahov(1,2)-Poloha_L_ECEF(i,2)).^2+(ReceiverStrahov(1,3)-Poloha_L_ECEF(i,3)).^2); % výpoèet vzdálenosti dvou bodù v prostoru
                  TimeReceiverAircraft(i,1)=DistanceReceiverAircraft(i,1)/c; % èas šíøení zprávy od letadla k jednotlivému pøijímaèi
                  TimeAircraft(i,1)=(Message_ID{i,5})-(TimeReceiverAircraft(i,1)); %výsledný èas vyslání zprávy letadlem
                  Receiver1Time(r1,3)=TimeAircraft(i);
               else
                  Receiver1Time(r1,3)=NaN;
               end
               
               r1=r1+1;

           elseif Message_ID{i,4}==2
                  Receiver2Time(r2,1)=i;
                  Receiver2Time(r2,2)=Message_ID{i,5};
               
               if Poloha_L_ECEF(i,1)~=0
                  DistanceReceiverAircraft(i,1)=sqrt((ReceiverHorska(1,1)-Poloha_L_ECEF(i,1)).^2+(ReceiverHorska(1,2)-Poloha_L_ECEF(i,2)).^2+(ReceiverHorska(1,3)-Poloha_L_ECEF(i,3)).^2); % výpoèet vzdálenosti dvou bodù v prostoru
                  TimeReceiverAircraft(i,1)=DistanceReceiverAircraft(i,1)/c; %èas šíøení zprávy od letadla k jednotlivému pøijímaèi
                  TimeAircraft(i,1)=(Message_ID{i,5})-(TimeReceiverAircraft(i,1)); %výsledný èas vyslání zprávy letadlem
                  Receiver2Time(r2,3)=TimeAircraft(i);
               else
                  Receiver2Time(r2,3)=NaN;
               end
               
               r2=r2+1;

           elseif Message_ID{i,4}==3
                  Receiver3Time(r3,1)=i;
                  Receiver3Time(r3,2)=Message_ID{i,5};
               
               if Poloha_L_ECEF(i,1)~=0
                  DistanceReceiverAircraft(i,1)=sqrt((ReceiverTower(1,1)-Poloha_L_ECEF(i,1)).^2+(ReceiverTower(1,2)-Poloha_L_ECEF(i,2)).^2+(ReceiverTower(1,3)-Poloha_L_ECEF(i,3)).^2); % výpoèet vzdálenosti dvou bodù v prostoru
                  TimeReceiverAircraft(i,1)=DistanceReceiverAircraft(i,1)/c; %èas šíøení zprávy od letadla k jednotlivému pøijímaèi
                  TimeAircraft(i,1)=(Message_ID{i,5})-(TimeReceiverAircraft(i,1)); %výsledný èas vyslání zprávy letadlem
                  Receiver3Time(r3,3)=TimeAircraft(i);
               else
                  Receiver3Time(r3,3)=NaN;
               end
               
               r3=r3+1;
       
           elseif Message_ID{i,4}==4
                  Receiver4Time(r4,1)=i;
                  Receiver4Time(r4,2)=Message_ID{i,5};
               
               if Poloha_L_ECEF(i,1)~=0
                  DistanceReceiverAircraft(i,1)=sqrt((ReceiverLKLT(1,1)-Poloha_L_ECEF(i,1)).^2+(ReceiverLKLT(1,2)-Poloha_L_ECEF(i,2)).^2+(ReceiverLKLT(1,3)-Poloha_L_ECEF(i,3)).^2); % výpoèet vzdálenosti dvou bodù v prostoru
                  TimeReceiverAircraft(i,1)=DistanceReceiverAircraft(i,1)/c; %èas šíøení zprávy od letadla k jednotlivému pøijímaèi
                  TimeAircraft(i,1)=(Message_ID{i,5})-(TimeReceiverAircraft(i,1)); %výsledný èas vyslání zprávy letadlem
                  Receiver4Time(r4,3)=TimeAircraft(i);
               else
                  Receiver4Time(r4,3)=NaN;
               end
               
               r4=r4+1;
           end 
       end 
          
%pøiøazení èasu vyslání zprávám bez znalosti polohy na pøijímaèi 1

j=1;
    while isnan(Receiver1Time(j,3)) %zjištìní indexu zprávy se známým èasem vyslání
    j=j+1;
    end
    
T=Receiver1Time(j,3);
I=j;
    for i=1:j-1
    Receiver1Time(i,3)=T-Receiver1Time(I,2)+Receiver1Time(i,2); %doplnìní èasù vyslání zprávám pøed zprávou se známou polohou
    end

    for i=j:r1-1 %%protože r1 je +1 oproti matici Receiver1Time z tohoto dùvodu -1
       if ~isnan(Receiver1Time(i,3))
           T=Receiver1Time(i,3);
           I=i; %"zafixování" I (index zprávy se známou polohou - známým èasem vyslání zprávy)
       else
       Receiver1Time(i,3)=T+Receiver1Time(i,2)-Receiver1Time(I,2); %pøièítání èasu rozdílu pøijmu dvou zpráv ke zprávì u které znám èas vyslání
       end
    end

%pøiøazení èasu vyslání zprávám bez znalosti polohy na pøijímaèi 2

j=1;
while isnan(Receiver2Time(j,3)) %zjištìní indexu zprávy se známým èasem vyslání
    j=j+1;
end
    
T=Receiver2Time(j,3);
I=j;
for i=1:j-1
    Receiver2Time(i,3)=T-Receiver2Time(I,2)+Receiver2Time(i,2); %doplnìní èasù vyslání zprávám pøed zprávou se známou polohou
end

for i=j:r2-1 %protože r2 je +1 oproti matici Receiver2Time z tohoto dùvodu -1
       if ~isnan(Receiver2Time(i,3))
           T=Receiver2Time(i,3);
           I=i; %"zafixování" I (index zprávy se známou polohou - známým èasem vyslání zprávy)
       else
       Receiver2Time(i,3)=T+Receiver2Time(i,2)-Receiver2Time(I,2); %pøièítání èasu rozdílu pøijmu dvou zpráv ke zprávì u které znám èas vyslání
       end
end

%pøiøazení èasu vyslání zprávám bez znalosti polohy na pøijímaèi 3

j=1;
while isnan(Receiver3Time(j,3)) %zjištìní indexu zprávy se známým èasem vyslání
    j=j+1;
end
    
T=Receiver3Time(j,3);
I=j;
for i=1:j-1
    Receiver3Time(i,3)=T-Receiver3Time(I,2)+Receiver3Time(i,2); %doplnìní èasù vyslání zprávám pøed zprávou se známou polohou
end

for i=j:r3-1 %protože r3 je +1 oproti matici Receiver3Time z tohoto dùvodu -1
       if ~isnan(Receiver3Time(i,3))
           T=Receiver3Time(i,3);
           I=i; %"zafixování" I (index zprávy se známou polohou - známým èasem vyslání zprávy)
       else
       Receiver3Time(i,3)=T+Receiver3Time(i,2)-Receiver3Time(I,2); %pøièítání èasu rozdílu pøijmu dvou zpráv ke zprávì u které znám èas vyslání
       end
end

%pøiøazení èasu vyslání zprávám bez znalosti polohy na pøijímaèi 4
   
j=1;
while isnan(Receiver4Time(j,3)) %zjištìní indexu zprávy se známým èasem vyslání
    j=j+1;
end
    
T=Receiver4Time(j,3);
I=j;
for i=1:j-1
    Receiver4Time(i,3)=T-Receiver4Time(I,2)+Receiver4Time(i,2); %doplnìní èasù vyslání zprávám pøed zprávou se známou polohou
end

for i=j:r4-1 %protože r4 je +1 oproti matici Receiver4Time z tohoto dùvodu -1
       if ~isnan(Receiver4Time(i,3))
           T=Receiver4Time(i,3);
           I=i; %"zafixování" I (index zprávy se známou polohou - známým èasem vyslání zprávy)
       else
       Receiver4Time(i,3)=T+Receiver4Time(i,2)-Receiver4Time(I,2); %pøièítání èasu rozdílu pøijmu dvou zpráv ke zprávì u které znám èas vyslání
       end
end

%cyklus pro vyhledávání odražených zpráv na pøijímaèi 1
x=10; %nadefinování parametru urèujícího poèet porovnávaných zpráv
y=r1-1;
Receiver1Time(:,4)=0;

for i=1:y-1
     m=length(Message_ID{Receiver1Time(i,1),2});
            for j=i+1:(i+x)*(i+x<y)+y*(i+x>=y) %logická podmínka
                if Receiver1Time(j,4)~=1
                    if m==length(Message_ID{Receiver1Time(j,1),2})&&abs(Receiver1Time(i,3)-Receiver1Time(j,3))<(m==14)*56*(10^-6)+(m==28)*112*(10^-6) %logická podmínka
                       Receiver1Time(j,4)=1;
                    end
                end
            end
end

Odraz1=sum(Receiver1Time(:,4));
Receiver1=r1-1-Odraz1;

%cyklus pro vyhledávání odražených zpráv na pøijímaèi 2
x=10; %nadefinování parametru urèujícího poèet porovnávaných zpráv
y=r2-1;
Receiver2Time(:,4)=0;

for i=1:y-1
     m=length(Message_ID{Receiver2Time(i,1),2});
            for j=i+1:(i+x)*(i+x<y)+y*(i+x>=y) %logická podmínka
                if Receiver2Time(j,4)~=1
                    if m==length(Message_ID{Receiver2Time(j,1),2})&&abs(Receiver2Time(i,3)-Receiver2Time(j,3))<(m==14)*56*(10^-6)+(m==28)*112*(10^-6) %logická podmínka
                       Receiver2Time(j,4)=1;
                    end
                end
            end
end

Odraz2=sum(Receiver2Time(:,4));
Receiver2=r2-1-Odraz2;

%cyklus pro vyhledávání odražených zpráv na pøijímaèi 3         
x=10; %nadefinování parametru urèujícího poèet porovnávaných zpráv
y=r3-1;
Receiver3Time(:,4)=0;

for i=1:y-1
     m=length(Message_ID{Receiver3Time(i,1),2});
            for j=i+1:(i+x)*(i+x<y)+y*(i+x>=y) %logická podmínka
                if Receiver3Time(j,4)~=1
                    if m==length(Message_ID{Receiver3Time(j,1),2})&&abs(Receiver3Time(i,3)-Receiver3Time(j,3))<(m==14)*56*(10^-6)+(m==28)*112*(10^-6) %logická podmínka
                       Receiver3Time(j,4)=1;
                    end
                end
            end
end

Odraz3=sum(Receiver3Time(:,4));
Receiver3=r3-1-Odraz3;

%cyklus pro vyhledávání odražených zpráv na pøijímaèi 4         
x=10; %nadefinování parametru urèujícího poèet porovnávaných zpráv
y=r4-1;
Receiver4Time(:,4)=0;

for i=1:y-1
     m=length(Message_ID{Receiver4Time(i,1),2});
            for j=i+1:(i+x)*(i+x<y)+y*(i+x>=y) %logická podmínka
                if Receiver4Time(j,4)~=1
                    if m==length(Message_ID{Receiver4Time(j,1),2})&&abs(Receiver4Time(i,3)-Receiver4Time(j,3))<(m==14)*56*(10^-6)+(m==28)*112*(10^-6) %logická podmínka
                       Receiver4Time(j,4)=1;
                    end
                end
            end
end
 
Odraz4=sum(Receiver4Time(:,4));
Receiver4=r4-1-Odraz4;
         
% slouèení matic jednotlivých pøijímaèù do jedné spoleèné matice

ReceiversTime=zeros(n1,4);

for i=1:n1
       d1=find(Receiver1Time(:,1)==i); %hledá index i v recieveru pokud ano do d1 zapíše na jakém øádku pro i=1 najde 1 nakopíruje zprávu do promìnné ReceiversTime na pozici 1
       d2=find(Receiver2Time(:,1)==i);
       d3=find(Receiver3Time(:,1)==i);
       d4=find(Receiver4Time(:,1)==i);
       if ~isempty(d1)
          ReceiversTime(i,:)=Receiver1Time(d1,:);
          elseif ~isempty(d2)
          ReceiversTime(i,:)=Receiver2Time(d2,:);
          elseif ~isempty(d3)
          ReceiversTime(i,:)=Receiver3Time(d3,:);
          elseif ~isempty(d4)
          ReceiversTime(i,:)=Receiver4Time(d4,:);
       end
end

Odrazy=sum(ReceiversTime(:,4));

%vyhledávání duplicitních zpráv ve sjednocené matici ReceiversTime (napøíè všemi pøijímaèi)

x=500; %nadefinování parametru urèujícího poèet porovnávaných zpráv 
y=n1;

for i=1:y-1
     m=length(Message_ID{i,2}); %hledám délku zprávy,která je v Message_ID na pozici i
            for j=i+1:(i+x)*(i+x<y)+y*(i+x>=y) %logická podmínka
                if ReceiversTime(j,4)~=1
                    if m==length(Message_ID{j,2})&&abs(ReceiversTime(i,3)-ReceiversTime(j,3))<(m==14)*56*(10^-6)+(m==28)*112*(10^-6) %logická podmínka
                       ReceiversTime(j,4)=1;
                    end
                end
            end
end

DuplicityTotal=sum(ReceiversTime(:,4));

MessagesID=length(Message_ID)-sum(ReceiversTime(:,4))

DF0=0;
DF4=0;
DF5=0;
DF11=0;
DF16=0;
DF17=0;
TC4=0;
TC18=0;
TC19=0;
DF20=0;
DF21=0;

for k=1:n1
    
    if Message_ID{k,3}==0&&ReceiversTime(k,4)~=1 %celkový poèet ID pøijatých zpráv DF 0 od daného letadla
    DF0=DF0+1;
    end
    
    if Message_ID{k,3}==4&&ReceiversTime(k,4)~=1 %celkový poèet ID pøijatých zpráv DF 4 od daného letadla
    DF4=DF4+1;
    end
    
    if Message_ID{k,3}==5&&ReceiversTime(k,4)~=1 %celkový poèet ID pøijatých zpráv DF 5 od daného letadla
    DF5=DF5+1;
    end
    
    if Message_ID{k,3}==11&&ReceiversTime(k,4)~=1 %celkový poèet ID pøijatých zpráv DF 11 od daného letadla
    DF11=DF11+1;
    end
    
    if Message_ID{k,3}==16&&ReceiversTime(k,4)~=1 %celkový poèet ID pøijatých zpráv DF 16 od daného letadla
    DF16=DF16+1;
    end
    
    if Message_ID{k,3}==17&&ReceiversTime(k,4)~=1 %celkový poèet ID pøijatých zpráv DF 17 od daného letadla bez duplicitních zpráv
       DF17=DF17+1;
    
      if Typecode(k,1)>0&&Typecode(k,1)<5
          TC4=TC4+1;
      end
        
      if Typecode(k,1)>8&&Typecode(k,1)<19
           TC18=TC18+1;
      end
        
      if Typecode(k,1)==19
           TC19=TC19+1;
      end
      
    end
    
    if Message_ID{k,3}==20&&ReceiversTime(k,4)~=1 %celkový poèet ID pøijatých zpráv DF 20 od daného letadla
    DF20=DF20+1;
    end
    
    if Message_ID{k,3}==21&&ReceiversTime(k,4)~=1 %celkový poèet ID pøijatých zpráv DF 21 od daného letadla
    DF21=DF21+1;
    end

end

%finální matice zpráv vyslaných odpovídaèem vybrané ICAO adresy (ID)
j=1;
for i=1:n1
    if ReceiversTime(i,4)==0
        DFinTime(j,1)=ReceiversTime(i,3);
        DFinTime(j,2)=Message_ID{i,4};
        DFinTime(j,3)=Message_ID{i,3};
        j=j+1;
    end
end

T=DFinTime(j-1,1)-DFinTime(1,1);


