% ОБРАБОТКА ДАННЫХ ИЗМЕРЕНИЙ С ПОМОЩЬЮ FRG-601 (КАНАЛЫ H,D,Z)
% Магнитная обсерватория "Паратунка"

% 01.04.2014 - Хомутов С.Ю. (адаптирован текст программы для Magdas)
% 05.05.2014 - убран учет разворота XY
% 23.07.2014 - сделан контроль размера файлов мин. и сек. данных
% 15.04.2016 - вычисление времени в часах - вместо rem(T,1)*24 использую
%               (T-T0)*24 - для получения отрицательных th
% 14.09.2017 - вычисляется корр. ф-ция и делается оценка рассинхронизации
%               таймеров
% 08.10.2018 - в xcorr() используется параметр maxlag=90 - ф-ция считается
%               для сдвигов до +/-90 сек - намного быстрее
% 21.01.2020 - закомментрировано сравнение с FGE (FGE остановлен из-за
%               ФСТЭКа)
%            - сделано сравнение с POS-1
% 07.02.2020 - на Fig1 рисуются пропущенные данные
% 10.05.2022 - обработка при отсутствии данных POS-1
% 10.12.2022 - сделана возможность коррекции сек. данных
% 23.11.2023 - оценка сдвига FGE - ф-цией delay_sec_data_fun
% 06.12.2024 - сделана коррекция таймера (в коде около строки 37)

clear, close all
global date0
% инициализация постоянных длЯ обработки
run ..\ini_parameters\pos1_ini.m
run ..\ini_parameters\frg_ini.m

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

% считываем секундные данные FRG (в нТл)
file_bin=[path_data_frg date0(7:10) '\' date0(7:10) date0(4:5) ...
          '\' strrep(datestr(T0,26),'/','') '.sec'];
disp(['... Считываются данные FRG-601 (1 sec) ... ' file_bin])
[Tj,dHj,dDj,dZj,Hshift,Dshift,Zshift]=read_frg_sec_fun(file_bin);
% ===== коррекция таймера 
[Tj,kc]=correct_time_frg_fun(T0,Tj);
if ~isempty(kc);
    dHj=dHj(kc); dDj=dDj(kc); dZj=dZj(kc);
    Hshift=Hshift(kc); Dshift=Dshift(kc); Zshift=Zshift(kc);
end
% =====
thj=(Tj-T0)*24;    % thj=rem(Tj,1)*24;
dHj=dHj+Hshift; dDj=dDj+Dshift; dZj=dZj+Zshift;

% % заполняем NaN-нами неизмеренные значения (=0 во всех каналах)
% k=(dHj==0 & dDj==0 & dZj==0);
% if sum(k)>0
%     dHj(k)=NaN; dDj(k)=NaN; dZj(k)=NaN;
% end

% проверяем, есть ли пропуски (NaN) в исходных данных
k_NaN=section_limit(Tj,dZj); [m,n]=size(k_NaN);
if ~isempty(k_NaN)
    disp('!!! Есть NaN в исходных данных !!!')
    for ki=1:m 
        fprintf('%s --- %s\n',datestr(Tj(k_NaN(ki,1)),0),datestr(Tj(k_NaN(ki,2)),13));
    end
end

% Отбраковка задокументированных помех
dHj=remove_errors('..\Ini_parameters\Remove_FRG.dat',T0,thj,dHj,1);
dDj=remove_errors('..\Ini_parameters\Remove_FRG.dat',T0,thj,dDj,2);
dZj=remove_errors('..\Ini_parameters\Remove_FRG.dat',T0,thj,dZj,3);

% 10.12.2022 - по образцу файла для 2004 г.
dHji=dHj; dDji=dDj; dZji=dZj;
% коррекция смещения исходных вариаций (dD - в нТл)
[dHj,dDj,dZj]=correction_sec_data_frg_fun(Tj,dHji,dDji,dZji,T0);

% Считываем базисные на дату (текущая или ближайшая предыдущая)
tmp=load('..\Ini_parameters\Baseline_current_FRG.dat');
T_bl=datenum(tmp(:,1),tmp(:,2),tmp(:,3));
k=max(find(abs(T_bl-T0)==min(abs(T_bl-T0)))); % выбираем ближайшие сутки с имеющимися базисными
H0j=tmp(k,4); D0j=tmp(k,5)/60; Z0j=tmp(k,6);   % D0 - в град.
fprintf('Текущие базисные ... %+8.2f nT, %+8.4f'', %+8.2f nT (за %s)\n',...
        H0j,D0j*60,Z0j,datestr(T_bl(k)))

% сохраняем текущие сутки в mat-файл
file_mat=['mat\' num2str(date0(7:10)) '\FRG_' date0(7:10) date0(4:5) date0(1:2)];
save(file_mat,'Tj','dHj','dDj','dZj','H0j','D0j','Z0j');

% Вычисление полных составляющих поля и модуля
Hj=sqrt((H0j+dHj).^2+dDj.^2);
% H=H0+dHv;
Zj=Z0j+dZj;
dDja=r2d(atan(dDj./Hj)); Dj=D0j+dDja;
Fj=total_field(Hj,Zj);

HHj=base_to_extract(dHj); ZZj=base_to_extract(dZj);

figure(1)
knan=isnan(dHj);
subplot(311)
plot(thj,dHji-HHj,'g',thj,dHj-HHj,thj(knan),ones(sum(knan),1)*nanmean(dHj)-HHj,'rx','MarkerSize',4) 
ylabel(sprintf('dH%+d, nT',-HHj)), title(['FRG-601 ' date0])
knan=isnan(dDj);
subplot(312)
plot(thj,dDji,'g',thj,dDj,thj(knan),ones(sum(knan),1)*nanmean(dDj),'rx','MarkerSize',4) 
ylabel('dD, nT')
subplot(313)
knan=isnan(dZj);
subplot(313)
plot(thj,dZji-ZZj,'g',thj,dZj-ZZj,thj(knan),ones(sum(knan),1)*nanmean(dZj)-ZZj,'rx','MarkerSize',4) 
ylabel(sprintf('dZ%+d, nT',-ZZj))
xlabel('Time UT, hours'), %set(gca,'XLim',[0 24],'XTick',0:3:24)
file=['Graphs\' num2str(date0(7:10)) '\FRG_' date0(7:10) date0(4:5) date0(1:2) '_dHdDdZ.png'];
print('-dpng', file)

% Сравниваем с данными FGE
file_fge=['\Paratunka\FGE_GSM\Mat\' date0(7:10)...
          '\fge_gsm_' strrep(datestr(T0,26),'/','') '.mat'];
if exist(file_fge,'file'),
    load(file_fge);
    Tv=Tv-0/86400; % коррекция таймера FGE+GSM
    thv=(Tv-T0)*24;
    % Вычисление полных составляющих поля и модуля
    H=sqrt((H0+dHv).^2+dDv.^2);         % по Линсу
    Z=Z0+dZv;
    dDa=r2d(atan(dDv./H)); D=D0+dDa;    % по Линсу
else
    disp('... Нет данных FGE - генерируются из NaN ...')
    Tv=T0+((1:86400)'-1)/86400; Nv=length(Tv); thv=(Tv-T0)*24;
    H=ones(Nv,1)*NaN; D=H; Z=H;
end

thv=(Tv-T0)*24;    % thv=rem(Tv,1)*24;
[kj,kv]=ismember_digit(thj,thv,8);
figure(2)
subplot(311),plot(thj(kj),Hj(kj)-H(kv))
title(['FRG-601 ' date0 '; Сравнение FRG-601 c FGE'])
ylabel('H(FRG)-H(FGE), nT')
subplot(312),plot(thj(kj),(Dj(kj)-D(kv))*60)
ylabel('D(FRG)-D(FGE), arcmin')
subplot(313),plot(thj(kj),Zj(kj)-Z(kv))
ylabel('Z(FRG)-Z(FGE), nT')
xlabel('Time UT, hours')
file=['Graphs\' num2str(date0(7:10)) '\FRG-FGE_' date0(7:10) date0(4:5) date0(1:2) '.png'];
print('-dpng', file)

figure(3)
subplot(311),plot(thj,Hj-HH,'r',thv,H-HH)
title(['FRG-601 vs FGE ' date0 '; Сравнение FRG-601 c FGE'])
ylabel('H-21000, nT'), legend('FRG-601','FGE')
subplot(312),plot(thj,Dj*60,'r',thv,D*60)
ylabel('D, arcmin')
subplot(313),plot(thj,Zj-ZZ,'r',thv,Z-ZZ)
ylabel('Z-46000, nT')
xlabel('Time UT, hours')
file=['Graphs\' num2str(date0(7:10)) '\FRG_and_FGE_' date0(7:10) date0(4:5) date0(1:2) '.png'];
print('-dpng', file)

if sum(~isnan(H))>0
    % оценка сдвига по времени
    [C,lag,lag_max_C]=delay_sec_data_fun(thj,Hj-nanmean(Hj),thv,H-nanmean(H),60,1);
    subplot(211) 
    title(['FRG-601 ' date0 '; Сравнение c FGE'])
    legend('FRG(Japan)','FGE'), ylabel('H-scaps(0.1), nT')
    subplot(212)
    title(['FRG-601 ' date0 '; Сдвиг по времени с FGE, maximum на '...
           sprintf('%+4.1f сек',lag_max_C)])
end

% данные POS-1 в абсол. пав-не
DOY=T0-datenum(str2num(date0(7:10)),1,0);  % день года
file_pos=[path_data_pos sprintf('%4d\\POS1_%4d_AbsHut.%03d',year0,year0,DOY)];
if exist(file_pos,'file')
    disp(['... Считываются данные POS-1 ... ' file_pos])
    [fpi,Tp,qmc,err]=pos_read_pm(file_pos,[]);
    thp=(Tp-T0)*24;    % thp=rem(Tp,1)*24;
    % % отбраковка задокументированных помех
    fp=remove_errors('..\ini_parameters\remove_pos1.dat',T0,thp,fpi);
    % чистим брак (недостоверные коды ошибок и выбросы)
    fp=pos_correct_all(fp,err,thp,dF_spike_lim_POS,accept_code_pos);
    % принудительная чистка больших выбросов
    k=fp<50500;
    if sum(k)>0
        fp(k)=NaN;
        fprintf('Принудительная чистка с Fp<50500 ... %d',sum(k));
    end
    % Градиент dF в абсол. пав-не - F(DI-pillar)=F(place)+dFgrad_AbsHut
    tmp=load('..\ini_parameters\dFgrad_current_POS1.dat');
    TgradPOS1=datenum(tmp(:,1:3));
    k=max(find(TgradPOS1<=T0));
    dFgrad_POS1=tmp(k,4);
    fprintf('Приведение данных POS-1 к DI-столбу ... %+6.2f\n',dFgrad_POS1);
    % систематич. погрешность абсол. протон. магнитометра
    tmp=load('..\ini_parameters\dFsys_AbsHut.dat');
    TsysAbsHut=datenum(tmp(:,1:3)); 
    k=max(find(TsysAbsHut<=T0));
    dFsys=tmp(k,4);
    fprintf('Система ММП/GSM/POS1 ... %+6.2f\n',dFsys)
    fpi=fpi+dFgrad_POS1+dFsys; fp=fp+dFgrad_POS1+dFsys; % редукция POS-1 
else
    disp('!!! Файла с данными POS-1 нет - использую NaNs...')
    Tp=T0+(0:dt_sec_POS:86400-0.5)/86400'; thp=(Tp-T0)*24;
    fpi=ones(length(Tp),1)*NaN; fp=fpi;
end

% thj=thj+3/3600; Tj=Tj+3/86400; % для теста

figure(5)
subplot(211)
plot(thj,Fj-FF,thp,fp-FF)
title(datestr(T0)), legend('F(FRG)','F(POS)','Location','Best')
ylabel('F-51000, nT')
[kj,kp]=ismember_digit(Tj,Tp,9);
subplot(212), plot(thj(kj),Fj(kj)-fp(kp))
ylabel('F(FRG)-F(POS), nT')

% Контроль размеров суточных файлов
DDD=dir(file_bin);
if DDD.bytes~=Size_sec_data % нестандартный размер
    fprintf('!!! Нестандартный размер файла сек. данных - %d\n',DDD.bytes)
end
file_min=[path_data_frg date0(7:10) '\' date0(7:10) date0(4:5) ...
          '\' strrep(datestr(T0,26),'/','') '.min'];
DDD=dir(file_min);
if DDD.bytes~=Size_min_data % нестандартный размер
    fprintf('!!! Нестандартный размер файла мин. данных - %d\n',DDD.bytes)
end

figure(9)
