% ОБРАБОТКА ДАННЫХ ИЗМЕРЕНИИ С ПОМОЩЬЮ GEM dIdD

% 27.06.2012 - Хомутов С.Ю., Новосибирск
% 17.08.2012 - адаптировано для Паратунки (исходный текст - для Хабаровска)
% 15.02.2013 - зменено расположение папок данных, пути - в dIdD_ini.m
% 15.05.2013 - полная ревизия ПО (С.Ю.Хомутов, Паратунка)
% 31.08.2013 - сделан вывод графиков
% 07.11.2013 - добавлен анализ пропуска данных и вывод сообщения об этом
% 19.06.2014 - обработка данных GSM-19Wc
% 31.07.2015 - сделана возможность коррекции исходных данных до получения
%               минутных
% 16.01.2016 - в качестве независимого измерения F используется POS-1
% 29.06.2016 - ввод даты - через функцию my_date_input()
% 22.08.2017 - добавлено графическое сравнение F(dIdD) с F(GSM-19W)
% 16.07.2019 - контроль лишних файлов в конце файла и недостающих в начале
% 25.05.2022 - обработка идёт при отсутствии данных dIdD - используются
%               NaNs
% 10.04.2024 - на fig3 не рисуется subplot(3) при отсутствии GSM-19
% 20.06.2025 - адаптация под Linux (в путях "\" заменён на "/" и прописные
%               в именах папок)

clear, close all
global date0
run ../ini_parameters/didd_ini_khb  % инициализация постоянных длЯ обработки
run ../ini_parameters/pos_ini_khb.m % настроечные параметры (пути, файлы)
run ../ini_parameters/gsm19w_ini_khb.m 
dt=2.5;         % интервал между измерениями (для создания фильтра Гаусса)

[date0,T0]=my_date_input;
[year0,mon0,day0]=datevec(T0);

disp(['--------------------- Обработка для даты - ' date0 ' -----------'])

% считываем данные dIdD
file_dIdD=[path_data_dIdD datestr(T0,'yyyy') '/' datestr(T0,'yyyymmdd') '.mag'];
if ~exist(file_dIdD),
%     error('!!! Не могу найти файл с данными dIdD !!!')
    th=(0:dt:86400-1)'/3600;
    Fi=ones(length(th),1)*NaN; dIi=Fi; dDi=Fi;
    SQ_F=Fi; F_DI=Fi; SQ_DI=Fi;
    disp(['!!! Не могу найти файл с данными dIdD - использую NaNs ...' file_dIdD])
else
    disp(['!!! Считываю файл ...' file_dIdD])
    [th,Fi,dIi,dDi,head,SQ_F,F_DI,SQ_DI]=didd_var_read_fun(file_dIdD);
end
T=datenum(year0,mon0,day0,th,0,0);
% инвертируем данные dD (так получается...)
dDi=-dDi;

% Контроль времени в файле данных
k=find(diff(th)<=0);
if ~isempty(k), % есть критические сбои таймера
    disp('!!! Есть критические ошибки таймера, см. файл time_error.log !!!')
    [yyy,mmm,ddd,hh,mm,ss]=datevec(T(k-1:k+2));
    fid=fopen('time_error.log','wt');
    fprintf(fid,'%4d %02d %02d  %02d %02d %04.1f\n',[yyy mmm ddd hh mm ss]');
    fclose(fid);
end
dt=round(mean(diff(th))*3600*10)/10; % оценка интервала между измерениями
k=find(diff(th)*3600>dt+0.1);  % пропуски в данных
if length(k)>0,
    disp('!!! Есть пропуски в данных !!! Примерно тут:')
    for kk=1:length(k),
        fprintf('%s ... %s\n',datestr(T(k(kk))),datestr(T(k(kk)+1)))
    end
end

% контроль пропусков в начале файла и лишних данных в конце - 07.01.2015
if min(T)~=T0,
    disp('!!! Нет данных в начале суточного файла !!!')
end
if max(T)>=T0+1,
    disp('!!! Лишние данные в конце суточного файла !!!')
end

% Отбраковка задокументированных помех
F=remove_errors('../ini_parameters/remove_dIdD_KHB.dat',T0,th,Fi,1);
dD=remove_errors('../ini_parameters/remove_dIdD_KHB.dat',T0,th,dDi,2);
dI=remove_errors('../ini_parameters/remove_dIdD_KHB.dat',T0,th,dIi,3); 

% коррекция секундных данных - 31.07.2015
[F,dD,dI]=correction_didd_fun(T,F,dD,dI);

% % Данные GSM-19W
% file_gsm=[path_data_GSM19 date0(7:10) '\gsm' strrep(datestr(T0,26),'/','') '.txt'];
% if exist(file_gsm,'file'),
%     [th19,F19,SQ,head]=gsm19_read_fun_khb(file_gsm);
%     F19=remove_errors('..\ini_parameters\remove_GSM19W_KHB.dat',T0,th19,F19);
% else
% %     disp('!!! Нет данных GSM-19W, использую NaN !!!')
%     th19=(0:23)'; F19=ones(length(th19),1)*NaN;
% end
% 
% % Текущий градиент
% tmp=load('..\ini_parameters\dFgrad_current_dIdD_KHB.dat');
% Tgrad=datenum(tmp(:,1:3));
% k=max(find(Tgrad<=T0));     % ближайшее более ранее значение градиента
% dFgrad=tmp(k,4);            % принятое ближайшее значение градиента

% отбракованные данные
kf=isnan(F); kd=isnan(dD); ki=isnan(dI);

figure(1)
subplot(311), plot(th,F-FF,th(kf),Fi(kf)-FF,'rx','MarkerSize',4)
axis tight, ylabel(sprintf('F%+d, nT',-FF))
title([date0 '  Khabarovsk KHB; GSM-19FD dIdD;'])
subplot(312), plot(th,dD*60,th(kd),dDi(kd)*60,'rx','MarkerSize',4)
axis tight, ylabel('dD, arc.min')
if sum(kd)>0,   % была чистка данных
    legend('Corrected','Removed','Location','best');
end
subplot(313), plot(th,dI*60,th(ki),dIi(ki)*60,'rx','MarkerSize',4)
axis tight, ylabel('dI, arc.min'), xlabel('Time UT, hours')
file=['graphs/' datestr(T0,'yyyy') '/' datestr(T0,'yyyymmdd') '_FdDdI.png'];
print('-dpng', file)

figure(2)
plot(th,F_DI-FF)
title([date0 '  Khabarovsk KHB; GSM-19FD dIdD;'])
axis tight, ylabel(sprintf('F(D,I)%+d, nT',-FF))
legend('Dp','Dm','Ip','Im','Location','northwest')
file=['graphs/' datestr(T0,'yyyy') '/' datestr(T0,'yyyymmdd') '_DpDmIpIm.png'];
print('-dpng', file)

% Текущий градиент
tmp=load('../ini_parameters/dFgrad_current_dIdD_KHB.dat');
Tgrad=datenum(tmp(:,1:3));
k=max(find(Tgrad<=T0));     % ближайшее более ранее значение градиента
dFgrad=tmp(k,4);            % принятое ближайшее значение градиента

% систематич. погрешность скалярного магнитометра
tmp=load('../ini_parameters/dFsys_AbsHut.dat');
TsysAbsHut=datenum(tmp(:,1:3)); 
k=max(find(TsysAbsHut<=T0));
dFsys_Scal=tmp(k,4);
fprintf('Систем. погрешность скалярного магнитометра ... %+6.2f\n',dFsys_Scal)

% считываем данные скалярного магнитометра POS (исходные)
DOY=T0-datenum(str2num(date0(7:10)),1,0);  % день года
file_pos=[path_data_pos sprintf('%4d/pos_%4d.%03d',year0,year0,DOY)];
if exist(file_pos,'file'),
    disp(['... Считываются данные POS-1 ... ' file_pos])
    [fpi,Tp,qmc,err]=pos_read_pm(file_pos,[]);
    thp=rem(Tp,1)*24;
    % чистим брак (недостоверные коды ошибок и выбросы)
    fp=pos_correct_all(fpi,err,thp,dF_spike_lim,accept_code_pos);
    % отбраковка задокументированных помех
    fp=remove_errors('../ini_parameters/remove_pos_khb.dat',T0,thp,fp);
%     % чистка по амплитуде
%     if sum(fp<54500)>0, 
%         disp('!!! чистка ... f<54500 = NaN')
%         fp(fp<54500)=NaN;
%     end
    % контроль рассинхронизации
    k=find(between(rem(Tp*86400,5),1-0.2,1+0.2)); % уход на 1 сек
    if length(k)>0,
        fprintf('Рассинхронизация POS-1 %s --- %s\n', datestr(Tp(min(k))),datestr(Tp(max(k)),13));
    end
    k=find(between(rem(Tp*86400,5),2-0.2,2+0.2)); % уход на 2 сек
    if length(k)>0,
        fprintf('Рассинхронизация POS-1 %s --- %s\n', datestr(Tp(min(k))),datestr(Tp(max(k)),13));
    end
    k=find(between(rem(Tp*86400,5),3-0.2,3+0.2)); % уход на 3 сек
    if length(k)>0,
        fprintf('Рассинхронизация POS-1 %s --- %s\n', datestr(Tp(min(k))),datestr(Tp(max(k)),13));
    end
    k=find(between(rem(Tp*86400,5),4-0.2,4+0.2)); % уход на 4 сек
    if length(k)>0,
        fprintf('Рассинхронизация POS-1 %s --- %s\n', datestr(Tp(min(k))),datestr(Tp(max(k)),13));
    end
    % Текущий градиент POS-1
    tmp=load('../ini_parameters/dFgrad_current_POS_KHB.dat');
    Tgrad_Scal=datenum(tmp(:,1:3));
    k=max(find(Tgrad_Scal<=T0));     % ближайшее более ранее значение градиента
    dFgrad_Scal=tmp(k,4);            % принятое ближайшее значение градиента
    fprintf('Текущий градиент POS-1 ... %+6.2f\n',dFgrad_Scal);
    fp=fp+dFgrad_Scal+dFsys_Scal; fpi=fpi+dFgrad_Scal+dFsys_Scal;
else
    disp('!!! Нет данных POS-1, использую NaN !!!')
    Tp=T0+(0:5:86400-1)'/86400; thp=rem(Tp,1)*24;
    fp=ones(length(Tp),1)*NaN;
end

% Данные GSM-19W
file_gsm=[path_data_GSM19 sprintf('%4d/gsm%s.txt',year0,datestr(T0,'yyyymmdd'))];
if exist(file_gsm,'file'),
    disp(['... Считываются данные GSM-19W ... ' file_gsm])
    [th19,F19,SQ,head]=gsm19_read_fun_khb(file_gsm);
    F19=remove_errors('../ini_parameters/remove_GSM19W_KHB.dat',T0,th19,F19);
    T19=T0+th19/24;
    % Градиент dF в абсол. пав-не - F(DI-pillar)=F(place)+dFgrad_AbsHut
    tmp=load('../ini_parameters/dFgrad_current_GSM19W_KHB.dat');
    Tgrad_Scal=datenum(tmp(:,1:3));
    k=max(find(Tgrad_Scal<=T0));
    dFgrad_Scal=tmp(k,4);
    fprintf('Приведение данных GSM-19W к DI-столбу ... %+6.2f\n',dFgrad_Scal);
    F19=F19+dFgrad_Scal+dFsys_Scal; % учитываем градиент
else
%     disp('!!! Нет данных GSM-19W, использую NaN !!!')
    F19=ones(1440,1)*NaN; T19=T0+(0:1440-1)'/1440; th19=rem(T19,1)*24;
end

% if sum(~isnan(fp))>0
if sum(~isnan(F19))==0 % нет GSM-19
    Nsub=2; % количество subplots
else
    Nsub=3;
end
figure(3)
subplot(Nsub,1,1), plot(th,F-FF+dFgrad,'r',thp,fp-FF,'b',...
    th19,F19-FF,'g')
legend(sprintf('dIdD%+6.2f',dFgrad),'POS-1','GSM-19W','Location','NorthWest')
title([date0 '  Khabarovsk KHB; dIdD vs POS-1 vs GSM-19W;'])
ylabel(sprintf('F%+d, nT',-FF))
fpi=interp1(thp,fp,th);    % линейно интерполируем данные POS-1
subplot(Nsub,1,2),plot(th,(F+dFgrad)-fpi,'.-')
ylabel('dIdD-POS, nT')
if Nsub==3
    F19i=interp1(th19,F19,th);    % линейно интерполируем данные GSM-19W
    subplot(Nsub,1,3), plot(th,(F+dFgrad)-F19i,'.-')
    ylabel('dIdD-GSM19W, nT')
end
xlabel('Time UT, hours')
file=['graphs/' datestr(T0,'yyyy') '/' datestr(T0,'yyyymmdd') '_dIdD_POS_GSM19.png'];
print('-dpng', file)
% end
