% ñòàò àíàëèç äàííûõ ïî dZ çà 2022

load train_stats_Z_all_days_2022_01_12 thbin SDddZ Tk Kind dth thstat
load schedule_Seyatel_2022.mat LT1 LT2 Naprav train Comments

% ëîêàëüíûå K-èíäåêñû
file_Kindex='d:\Projects\Noise_from_trains\other_data\NVS22k.DKA';
[Tk,Kind]=dka_read_fun(file_Kindex);

k1=Naprav==1; k2=Naprav==2;

figure(1)
plot(thstat,SDddZ,'.',thstat,nanmean(SDddZ),'ko',...
    LT1(k1),zeros(length(LT1(k1)),1),'b^',...
    LT1(k2),zeros(length(LT1(k2)),1),'rv',...
    'MarkerFaceColor','g')
% hold on, 
% plot(LT1(k1),zeros(length(LT1(k1)),1),'b^',LT1(k2),zeros(length(LT1(k2)),1),'rv')
% hold off
title('Novosibirsk   01.01.2022---31.12.2022  SD(diff(Z)) â 15-ìèí èíòåðâàëàõ è îòìåòêè ïîåçäîâ íà Ñåÿòåëå, all days')
ylabel('SD(diff(dZ)), nT'), xlabel('Local time, hours')

% ïðîáóåì ó÷åñòü âîçìóù¸ííîñòü ïîëÿ - äíè ñ áîëüøèì K-èíäåêñàìè èãíîðèðóåì
% ñóììàðíûé çà ñóòêè K-èíäåêñ
tmp=reshape(Kind,8,365); Ksum=sum(tmp); Ksum=Ksum(:);
% óáèðàåì äíè ñ áîëüøèìè K
k10=Ksum>10;
k10=repmat(k10,1,96);

SDddZ_Klt10=SDddZ; SDddZ_Klt10(k10)=NaN;
% stats_dz_mean=nanmean(SDddZ);
% stats_dz_mean_K10=nanmean(SDddZ_Klt10); % ñðåäíèå çà ãîäà äëÿ èíòåðâàëà â ñóòêàõ
stats_dz_mean=nanmedian(SDddZ);
stats_dz_mean_K10=nanmedian(SDddZ_Klt10); % ñðåäíèå çà ãîäà äëÿ èíòåðâàëà â ñóòêàõ

figure(2)
plot(thstat,stats_dz_mean,'ko',thstat,stats_dz_mean_K10,'kp',...
     thstat,SDddZ_Klt10,'.')
legend('All days','Ksum<10')

kSNe=Naprav==1 & train==2; kNSe=Naprav==2 & train==2; % ýëåêòðè÷êè
kSNp=Naprav==1 & train==1; kNSp=Naprav==2 & train==1; % ïàññàæèðñêèå
kpass=train==1;% ïàññàæèðñêèå
figure(3)
plot(thstat,stats_dz_mean,'ko',thstat,stats_dz_mean_K10,'r*','MarkerFaceColor','g')
hold on
plot(LT1(kSNe),zeros(length(LT1(kSNe)),1)+0.051,'b^',...
     LT1(kNSe),zeros(length(LT1(kNSe)),1)+0.049,'rv',...
     LT1(kpass),zeros(length(LT1(kpass)),1)+0.045,'md')
hold off
set (gca,'Xlim',[0 24],'XTick',0:4:24)
legend('Âñå äíè','Ksum<10','ý/ï íà ñåâåð','ý/ï íà þã',...
       'ïàññ.')
title('Íîâîñèáèðñê, 01.01.2022---31.12.2022, SD(diff(Z))')
ylabel('median(SD(diff(dZ))), íÒë'), xlabel('Ìåñòíîå âðåìÿ, ÷àñ')
% ylabel('median(SD(diff(dZ))), nT'), xlabel('Local time, hours')
