clc;clear;format longg;

den = [287:300,304:309];

Zpracovani_TBC = [];
Zpracovani_TBC_bez_CTAB_CJIH = [];

%načtení souřadnic
for i = 1:length(den)

Zpracovani_TBC = [Zpracovani_TBC
                  nacti_souradnice(den(i),"Zpracovani_TBC")];
% 
Zpracovani_TBC_bez_CTAB_CJIH = [Zpracovani_TBC_bez_CTAB_CJIH
                  nacti_souradnice(den(i),"Zpracovani_TBC_bez_CTAB_CJIH")];
end

Stanice = ['CJIH';'CTAB';'FRST';'GMUE';'HOR2';'MSTB';'ROHR'];
Stanice2 = ['FRST';'GMUE';'HOR2';'MSTB';'ROHR']; 
%výpočet průměrných hodnot
for i = 1:3
fprintf('%d.iterace\n',i)    
prumery = prumer(Stanice,Zpracovani_TBC);
prumery_bez_CTAB_CJIH = prumer(Stanice,Zpracovani_TBC_bez_CTAB_CJIH);

%výpočet směrodatných odchylek
[smer oprava] = smer_odchylka(Stanice,Zpracovani_TBC,prumery);
[smer_bez_CTAB_CJIH oprava_bez_CTAB_CJIH] = smer_odchylka(Stanice,Zpracovani_TBC_bez_CTAB_CJIH,prumery_bez_CTAB_CJIH);

% testovani odlehlých měření
[Zpracovani_TBC Krit_v]= testovani_odleh_hodnot(smer,oprava,Zpracovani_TBC);
% fprintf('\n')
[Zpracovani_TBC_bez_CTAB_CJIH Krit_v_bez_CTAB_CJIH] = testovani_odleh_hodnot(smer_bez_CTAB_CJIH,oprava_bez_CTAB_CJIH,Zpracovani_TBC_bez_CTAB_CJIH);

end

% rozdíl souřadnic
OFIC_sour = [table(Stanice) array2table([   49.3935539138889        15.5863957277778	 	  576.839  
                                            49.4097967694444        14.6802187194444	      496.233	
                                            48.508733029            14.503356764              621.546
                                            48.771978825            14.986788203              553.17
                                            48.663328819            15.646647458              387.202
                                            48.569907171            16.609795558              248.366
                                            48.5698333              13.995353698              658.065])]
%rozdíl zpracování
Rozdil_a_b = [table(Stanice2) array2table(prumery_bez_CTAB_CJIH{3:end,2:end}-prumery{:,2:end})];
Rozdil_ofic_a = [table(Stanice2) array2table(OFIC_sour{3:end,2:end} - prumery{:,2:end})];
Rozdil_ofic_b = [table(Stanice)  array2table(OFIC_sour{:,2:end} - prumery_bez_CTAB_CJIH{:,2:end})];

% převod fi, lambda na N, E

smer_sec = degrees2sec(smer);
smer_cm = filambda2NE(prumery{:,"PRŮMĚRNÁ ŠÍŘKA"},smer{:,"SMĚR. ODCHYLKA ŠÍŘKY "},smer{:,"SMĚR. ODCHYLKA DÉLKY "});
% 
smer_sec_bez_CTAB_CJIH = degrees2sec(smer_bez_CTAB_CJIH);
smer_cm_bez_CTAB_CJIH = filambda2NE(prumery_bez_CTAB_CJIH{:,"PRŮMĚRNÁ ŠÍŘKA"},smer_bez_CTAB_CJIH{:,"SMĚR. ODCHYLKA ŠÍŘKY "},smer_bez_CTAB_CJIH{:,"SMĚR. ODCHYLKA DÉLKY "});
% % 
oprava_sec = degrees2sec(oprava);

oprava_cm = filambda2NE(rozkopiruj(prumery{:,"PRŮMĚRNÁ ŠÍŘKA"},20),oprava{:,2},oprava{:,3});
% 
oprava_sec_bez_CTAB_CJIH = degrees2sec(oprava_bez_CTAB_CJIH);
oprava_cm_bez_CTAB_CJIH = filambda2NE(rozkopiruj(prumery_bez_CTAB_CJIH{:,"PRŮMĚRNÁ ŠÍŘKA"},20),oprava_bez_CTAB_CJIH{:,2},oprava_bez_CTAB_CJIH{:,3});
% 
Krit_v_sec = degrees2sec(Krit_v);
Krit_v_cm = filambda2NE(prumery{:,"PRŮMĚRNÁ ŠÍŘKA"},Krit_v{:,2},Krit_v{:,3});

Krit_v_sec_bez_CTAB_CJIH = degrees2sec(Krit_v_bez_CTAB_CJIH);
Krit_v_cm_bez_CTAB_CJIH = filambda2NE(prumery_bez_CTAB_CJIH{:,"PRŮMĚRNÁ ŠÍŘKA"},Krit_v_bez_CTAB_CJIH{:,2},Krit_v_bez_CTAB_CJIH{:,3});
%
Rozdil_a_b_sec = degrees2sec(Rozdil_a_b);   
Rozdil_a_b_cm = filambda2NE(prumery{:,"PRŮMĚRNÁ ŠÍŘKA"},Rozdil_a_b{:,2},Rozdil_a_b{:,3});

Rozdil_ofic_a_sec = degrees2sec(Rozdil_ofic_a);
Rozdil_ofic_a_cm = filambda2NE(prumery{:,"PRŮMĚRNÁ ŠÍŘKA"},Rozdil_ofic_a{:,2},Rozdil_ofic_a{:,3})*10;

Rozdil_ofic_b_sec = degrees2sec(Rozdil_ofic_b);
Rozdil_ofic_b_cm = filambda2NE(prumery_bez_CTAB_CJIH{:,"PRŮMĚRNÁ ŠÍŘKA"},Rozdil_ofic_b{:,2},Rozdil_ofic_b{:,3})*10;
function pr = prumer(stanice,tab)
pr = [];

    [s,~] = size(stanice);
    for i = 1:s

    index = all(tab.STANICE == stanice(i,:),2);

        if any(index) == 0 
            continue
        end

    pr = [pr
        table(convertCharsToStrings(stanice(i,:)), ...
        mean(tab.("ŠÍŘKA")(index),'omitnan'), mean(tab.("DÉLKA")(index),'omitnan'),mean(tab.("VÝŠKA")(index),'omitnan'), ...
        'VariableNames',{'STANICE','PRŮMĚRNÁ ŠÍŘKA','PRŮMĚRNÁ DÉLKA','PRŮMĚRNÁ VÝŠKA'});];
    end
end

function [sm opr] = smer_odchylka(stanice,tab,pr)
 [s,~] = size(stanice);

sm  = [];
opr = [];


    for i = 1:s

    index = all(tab.STANICE == stanice(i,:),2);
        if any(index) == 0 
            continue
        end

%     opravy 

    opr_i = [tab.("ŠÍŘKA")(index) - pr.("PRŮMĚRNÁ ŠÍŘKA")(pr.STANICE == stanice(i,:)),...
           tab.("DÉLKA")(index) - pr.("PRŮMĚRNÁ DÉLKA")(pr.STANICE == stanice(i,:)),... 
           tab.("VÝŠKA")(index) - pr.("PRŮMĚRNÁ VÝŠKA")(pr.STANICE == stanice(i,:))
             ];

    sm = [sm 
          table(convertCharsToStrings(stanice(i,:)), sqrt(sum(opr_i(:,1).^2,"omitnan")./(length(opr_i(~isnan(opr_i(:,1))))-1)), sqrt(sum(opr_i(:,2).^2,"omitnan")./(length(opr_i(~isnan(opr_i(:,2))))-1)), sqrt(sum(opr_i(:,3).^2,"omitnan")./(length(opr_i(~isnan(opr_i(:,3))))-1)),...
          'VariableNames',{'STANICE','SMĚR. ODCHYLKA ŠÍŘKY ','SMĚR. ODCHYLKA DÉLKY ','SMĚR. ODCHYLKA VÝŠKY'});];

   opr = [opr
          table(tab.("STANICE")(index,:),'VariableNames',{'STANICE'}),array2table(opr_i,'VariableNames',{'OPRAVA ŠÍŘKY ','OPRAVA DÉLKY ','OPRAVA VÝŠKY'}),table(tab.("DEN")(index),'VariableNames',{'DEN'})];

    end
end

function [Zpracovani v_max_out] = testovani_odleh_hodnot(sm,opr,Zpracovani)

    Krit_hodnoty = [15 2.62
                    16 2.642
                    17 2.664
                    18 2.686
                    19 2.708
                    20 2.73];
    [r,~] = size(sm);
    odlehle_all = [];
    v_max_out = [];
    for i = 1:r
        sm_i = sm(i,:);
        radky = strcmp(opr.STANICE, sm_i.STANICE);
        opr_i = opr{radky,2:end-1};
        r1 = length(opr_i(~isnan(opr_i(:,1))));
        r2 = length(opr_i(~isnan(opr_i(:,2))));
        r3 = length(opr_i(~isnan(opr_i(:,3))));
        r = [r1 r2 r3];
              
           for j = 1:3
               Krit_hodnoty(Krit_hodnoty(:,1) == r(j),2);
               v_max(j) = Krit_hodnoty(Krit_hodnoty(:,1) == r(j),2) * sm_i.(j+1);
               odlehle(:,j) = abs(opr_i(:,j)) > v_max(j);
           end
       odlehle;
       opr_i(odlehle);
       suma_odlehle = sum(odlehle);
       fprintf("%s: %d odlehlých %d odlehlých %d odlehlých \n",sm_i.STANICE,suma_odlehle);
       
       odlehle_all = [odlehle_all
                      odlehle];
v_max_out = [v_max_out
             table(sm_i.STANICE) array2table(v_max)]; 
    end

    Zpracovani = sortrows(Zpracovani,'STANICE');
    Zpracovani{:,2:4}(odlehle_all == 1) = NaN;
    Zpracovani = Zpracovani;
end

function sec = degrees2sec(data)

degrees2dms(data{:,2});
sec(:,1) = ans(:,3);
degrees2dms(data{:,3});
sec(:,2) = ans(:,3);
end

function cm = filambda2NE(fi,delta_fi,delta_lambda)
delta_fi_sec = degrees2dms(delta_fi);
delta_lambda_sec = degrees2dms(delta_lambda);
delta_N = 3094 .* delta_fi_sec(:,3);
delta_E = 3094 .* cosd(fi) .* delta_lambda_sec(:,3);
cm = [delta_N delta_E];
end

function b = rozkopiruj(sm,pocet)
[r,s] = size(sm);
b = [];
for i = 1:r
    b = [b 
         sm(i,:).*ones(pocet,s)];
end
end