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

% 27.06.2012 - Хомутов С.Ю., Новосибирск
% 24.07.2012 - используется ф-ция dIdD_calc_fun(), в которой вариации dI,dD
%               вычисляются по модулям, а не считываются из файла
% 28.07.2024 - при вычислении dD используется ближайшее базисное значение I0
%              из  Baseline_current_dIdD_KHB.dat
% 30.07.2024 - сделана маркировка удалённых значений Dp,Dm,Ip,Im
% 06.06.2025 - I0 для вычисления dD м.б. взят из заголовка исх. mag-файла
% 20.06.2025 - адаптация под Linux (в путях "\" заменён на "/" и прописные
%               в именах папок)


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

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

disp(['--------------------- Обработка для даты - ' date0 ' -----------'])
disp(['--------------------- dD,dI вычисляются по Dp,Dm,Ip,Im  --------'])

% считываем данные dIdD
file_dIdD=[path_data_dIdD datestr(T0,'yyyy') '/' datestr(T0,'yyyymmdd') '.mag'];
if ~exist(file_dIdD),
    disp('!!! Не могу найти файл с данными dIdD !!! ... Использую NaNs...')
    T=datenum(year0,mon0,day0,0,0,0:dt:86400-1)'; th=(T-T0)*24;
    Fi=ones(length(T),1)*NaN; dIi=Fi; dDi=Fi; SQ_F=Fi;
    F_DI=ones(length(T),4)*NaN; SQ_DI=F_DI;
    head='Нет файла с данными!';
else
    % Считываем базисные на дату (текущая или ближайшая предыдущая)
    tmp=load('../ini_parameters/Baseline_current_dIdD_KHB.dat');
    T_bl=datenum(tmp(:,1),tmp(:,2),tmp(:,3));
    k=max(find(T_bl<=T0));
    D0=tmp(k,4)/60; I0=tmp(k,5)/60;   % D0,I0 - в град.
    fprintf('Текущие базисные ... %+8.2f'', %+8.2f'' (за %s)\n',...
        D0*60,I0*60,datestr(T_bl(k)))
    disp(['... Считываются данные dIdD ... ' file_dIdD])
    [th,F_di,dI_di,dD_di,head_di,SQ_F_di,F_DI_di,SQ_DI_di]=didd_var_read_fun(file_dIdD); % модули до чистки
    k_inc=findstr(head_di{2},'Inclination :')+13; % считываем используемое в файле значение I
    I0=str2num(head_di{2}(k_inc:k_inc+10));
    fprintf('!!! Используется I0=%6.3f из заголовка исходного mag-файла !!!\n',I0)
%     disp('!!! Используется ближайшее I0 из Baseline_current_dIdD_KHB.dat !!!')
    [th,Fi,dIi,dDi,head,SQ_F,F_DI,SQ_DI]=didd_calc_fun(file_dIdD,I0,spikes_lim);
    T=datenum(year0,mon0,day0,th,0,0);
end
% инвертируем данные dD (так получается...)
dDi=-dDi; dD_di=-dD_di;

% Контроль времени в файле данных
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('!!! Нет данных dIdD в начале суточного файла !!!')
end
if max(T)>=T0+1,
    disp('!!! Лишние данные dIdD в конце суточного файла !!!')
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);

% Текущий градиент
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),F_di(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),dD_di(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),dI_di(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)
kdp=isnan(F_DI(:,1)); kdm=isnan(F_DI(:,2)); kip=isnan(F_DI(:,3)); kim=isnan(F_DI(:,4));
plot(th,F_DI(:,1)-FF,'b',th,F_DI(:,2)-FF,'c',...
     th,F_DI(:,3)-FF,'r',th,F_DI(:,4)-FF,'m',...
     th,F_DI_di(:,1)-FF,'b',th(kdp),F_DI_di(kdp,1)-FF,'ko',...
     th,F_DI_di(:,2)-FF,'c',th(kdm),F_DI_di(kdm,2)-FF,'ko',...
     th,F_DI_di(:,3)-FF,'r',th(kip),F_DI_di(kip,3)-FF,'ko',...
     th,F_DI_di(:,4)-FF,'m',th(kim),F_DI_di(kim,4)-FF,'ko',...
     'MarkerSize',3)
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/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'),
    [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
    figure(3)
%     subplot(311), plot(th,F-FF+dFgrad,thp,fp-FF,'r',...
%         th19,F19-FF)
%     legend(sprintf('dIdD%+6.2f',dFgrad),'POS-1','GSM-19W','Location','NorthWest')
    subplot(211), plot(th,F-FF+dFgrad,thp,fp-FF,'r')
    legend(sprintf('dIdD%+6.2f',dFgrad),'POS-1','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(212),plot(th,(F+dFgrad)-fpi,'.-')
    ylabel('dIdD-POS, nT')
%     F19i=interp1(th19,F19,th);    % линейно интерполируем данные GSM-19W
%     subplot(313), plot(th,(F+dFgrad)-F19i,'.-')
%     ylabel('dIdD-GSM19W, nT'), xlabel('Time UT, hours')
    file=['graphs/' datestr(T0,'yyyy') '/' datestr(T0,'yyyymmdd') '_dIdD_POS_GSM19.png'];
    print('-dpng', file)
% end
