% ОБРАБОТКА ДАННЫХ ИЗМЕРЕНИИ С ПОМОЩЬЮ GEM dIdD
% Вычислене минутных данных

% 27.06.2012 - Хомутов С.Ю., Новосибирск
% 12.12.2012 - адапитровано из dIdD_Process и CMVS_Process_min
% 13.12.2012 - добавлен анализ сбоя таймера (diff(T)<=0)
% 15.02.2013 - адаптировано для Паратунки (исходный код - для Хабаровска)
% 19.02.2013 - пути прописаны в инициализа
% 23.11.2013 - сделан вывод минутных данных в формате IAGA2002
% 04.12.2014 - паратунская версия адаптирована для Магадана
% 22.05.2016 - автоматическая отбраковка повторяющихся по времени данных
% 14.08.2018 - используется собственное вычисление вариаций dD,dI с помощью
%               didd_calc_fun.m и чисткой импульсных помех
% 04.01.2019 - адаптирован скрипт Магадана
% 19.11.2019 - в качестве скалярного магнитометра используется POS-1 или
%               GSM-19W (на выбор)
% 27.11.2020 - вывод файла iaga2002 сделан в папку года
% 28.07.2024 - при вычислении dD используется ближайшее базисное значение I0
%              из  Baseline_current_dIdD_KHB.dat
% 20.06.2025 - адаптация под Linux (в путях "\" заменён на "/" и прописные
%               в именах папок)

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

dt=2.5;         % интервал между измерениями (для создания фильтра Гаусса)
file_head_XYZ='../ini_parameters/IAGA2002_header_KHB_dIdD_POS.txt';
file_head_HDZ='../ini_parameters/IAGA2002_header_KHB_dIdD_HDZ_POS.txt';
file_head='../ini_parameters/IAGA2002_header_KHB_dIdD_calc_POS_Filter.txt';

% выбираем скалярный магнитометр для абсол. измерений
ScalMag='POS1';
% ScalMag='GSM19W';

[date0,T0]=my_date_input;
[year0,mon0,day0]=datevec(T0);
T0_1=T0-1; date0_1=strrep(datestr(T0_1,24),'/','.');    % Предыдущие сутки
[year0_1,mon0_1,day0_1]=datevec(T0_1);

% считываем данные dIdD
disp(['--------------------- dD,dI вычисляются по Dp,Dm,Ip,Im  --------'])

% ---------- предыдущие сутки
file_dIdD=[path_data_dIdD datestr(T0_1,'yyyy') '/' datestr(T0_1,'yyyymmdd') '.mag'];
if ~exist(file_dIdD),
    disp('!!! Не могу найти файл с данными dIdD - использую NaNs !!!')
    N=86400/dt; % кол-во отсчётов в сутках
    th_1=(0:dt:86400-1)'/3600; F_1=ones(N,1)*NaN; dI_1=F_1; dD_1=F_1;
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_1));
    D0=tmp(k,4)/60; I0=tmp(k,5)/60;   % D0,I0 - в град.
    disp(['... Считываются данные dIdD ... ' file_dIdD])
    [th,F_di,dI_di,dD_di,head_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));
%     disp('!!! Используется ближайшее I0 из Baseline_current_dIdD_KHB.dat !!!')
    fprintf('!!! Используется I0=%6.3f из заголовка исходного mag-файла !!!\n',I0)
%     I0=63; % для тестирования
    [th_1,F_1,dI_1,dD_1,head]=didd_calc_fun(file_dIdD,I0,spikes_lim);
end
T_1=datenum(year0_1,mon0_1,day0_1,th_1,0,0);
% убраем повторяющиеся моменты времени (втрое из двойных)
[T_1,k,l]=unique(T_1); 
if length(k)~=length(l),    % есть повторы времен
    disp('!!! Пытаюсь убрать повторяющиеся значения времени ...')
    th_1=th_1(k); F_1=F_1(k); dI_1=dI_1(k); dD_1=dD_1(k);
end
% Отбраковка задокументированных помех
F_1=remove_errors('../ini_parameters/remove_dIdD_KHB.dat',T0_1,th_1,F_1,1);
dD_1=remove_errors('../ini_parameters/remove_dIdD_KHB.dat',T0_1,th_1,dD_1,2);
dI_1=remove_errors('../ini_parameters/remove_dIdD_KHB.dat',T0_1,th_1,dI_1,3); 
% инвертируем данные dD (так получается...)
dD_1=-dD_1;
% коррекция секундных данных - 31.07.2015
[F_1,dD_1,dI_1]=correction_didd_fun(T_1,F_1,dD_1,dI_1);
% оставляем конец суток (1 мин)
k=th_1>23+59/60;

% текущие сутки
file_dIdD=[path_data_dIdD datestr(T0,'yyyy') '/' datestr(T0,'yyyymmdd') '.mag'];
if ~exist(file_dIdD),
    disp('!!! Не могу найти файл с данными dIdD - использую NaNs !!!')
    N=86400/dt; % кол-во отсчётов в сутках
    th=(0:dt:86400-1)'/3600; F=ones(N,1)*NaN; dI=F; dD=F;
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 - в град.
    disp(['... Считываются данные dIdD ... ' file_dIdD])
    [th,F_di,dI_di,dD_di,head_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));
%     disp('!!! Используется ближайшее I0 из Baseline_current_dIdD_KHB.dat !!!')
    fprintf('!!! Используется I0=%6.3f из заголовка исходного mag-файла !!!\n',I0)
%     I0=63; % для тестирования
    disp('!!! Используется I0 = 63 (из mag-файла)')
    [th,F,dI,dD,head]=didd_calc_fun(file_dIdD,I0,spikes_lim);
end
T=datenum(year0,mon0,day0,th,0,0);
% убраем повторяющиеся моменты времени (втрое из двойных)
[T,k,l]=unique(T); 
if length(k)~=length(l),    % есть повторы времен
    disp('!!! Пытаюсь убрать повторяющиеся значения времени ...')
    th=th(k); F=F(k); dI=dI(k); dD=dD(k);
end
% Отбраковка задокументированных помех
F=remove_errors('../ini_parameters/remove_dIdD_KHB.dat',T0,th,F,1);
dD=remove_errors('../ini_parameters/remove_dIdD_KHB.dat',T0,th,dD,2);
dI=remove_errors('../ini_parameters/remove_dIdD_KHB.dat',T0,th,dI,3); 
% инвертируем данные dD (так получается...)
dD=-dD;
% коррекция секундных данных - 31.07.2015
[F,dD,dI]=correction_didd_fun(T,F,dD,dI);

% Объединяем сутки
k=th_1>23+59/60;    % оставляем конец предыдущих суток (1 мин)
T=[T_1(k);T]; th=[th_1(k)-24;th];
F=[F_1(k);F]; dD=[dD_1(k);dD]; dI=[dI_1(k);dI];

% Контроль времени в файле данных
k=find(diff(th)<=0);
if ~isempty(k), % есть критические сбои таймера
    [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);
    error('!!! Есть критические ошибки таймера, см. файл time_error.log !!!')
end

% Сглаживаем данные F
W=gauss_coeff_fun(-45:dt:45);   % генерация весовых коэффициентов
[Fm,ts,n_good]=one_minute_mean_conv_w(0,86400-dt,th*3600,F,W,dt,0.75);
thm=ts/3600; Tm=datenum(year0,mon0,day0,0,0,ts);

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

% Считываем базисные на дату (текущая или ближайшая предыдущая)
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)))
fprintf('Текущий градиент dFgrad ... %+5.2f нТл\n',dFgrad)
fprintf('Периодичность измерений ... %3.1f сек\n',dt)

% Вычисляем полные компоненты
D=dD+D0; I=dI+I0;
H=F.*cos(d2r(I)); Z=F.*sin(d2r(I));

% Сглаживаем данные D,I
W=gauss_coeff_fun(-45:dt:45);   % генерация весовых коэффициентов
[Dm,ts,n_good]=one_minute_mean_conv_w(0,86400-dt,th*3600,D,W,dt,0.75);
[Im,ts,n_good]=one_minute_mean_conv_w(0,86400-dt,th*3600,I,W,dt,0.75);

Hm=Fm.*cos(d2r(Im)); Zm=Fm.*sin(d2r(Im));
Xm=Hm.*cos(d2r(Dm)); Ym=Hm.*sin(d2r(Dm));

% считываем данные скалярного магнитометра (POS-1 или GSM-19W)
[Tscal,Fscal,dFgrad_Scal]=scalar_mag_select(T0,ScalMag);
% частота измерений конкретного скалярного магнитометра
switch ScalMag
    case 'POS1'
        dt_scal=dt_pos;
    case 'GSM19W'
        dt_scal=dt_gsm19w;
end

% % интерполируем данные скал. магнитометра на начало минуты
% Fscali=interp1(Tscal,Fscal,Tm);
% интерполируем данные скал. магнитометра на равномерную шкалу
% чтобы убрать возможную рассинхронизацию
Tscali=datenum(year0,mon0,day0,0,0,-60:dt_scal:86400-1); Tscali=Tscali(:); % равномерная шкала
Fscali=interp1(Tscal,Fscal,Tscali); thi=(Tscali-T0)*24;
% минутные значения с фильтром Гаусса
Wscal=gauss_coeff_fun(-45:dt_scal:45);
[Fscalm,ts,n_good]=one_minute_mean_conv_w(0,86400-dt_scal,round(thi*3600),Fscali,Wscal,dt_scal,0.75);
disp('Минутные Fscal - по Гауссу...')

% Вывод минутных данных в стандарте Интермагнета (часовыми блоками)
Nm=1440; tm0=(0:Nm-1)';
for kh=1:24,
    hour_str=sprintf('%02d', kh-1);
    k=fix(tm0/60)==kh-1;
    file=hdzf_hour_out_khb_didd(date0,hour_str,Hm(k),Dm(k),Zm(k),Fscalm(k));
end

% Вывод минутных данных в формате IAGA2002
Header=iaga2002_write_fun(Tm,Xm,Ym,Zm,Fscalm,file_head,...
            ['data_IAGA2002/' datestr(T0,'yyyy') '/']);

% вывод mat-файла
k=T>=T0; T=T(k); F=F(k); dD=dD(k); dI=dI(k);
file_didd=[path_all 'dIdD/mat/' datestr(T0,'yyyy') '\dIdD_' datestr(T0,'yyyymmdd') '.mat'];
save(file_didd,'T','F','dD','dI','D0','I0','dFgrad','Fm','Dm','Im','Hm','Zm')

% рисуем графики
figure(1)
set(gcf,'PaperPosition',[1 6 15 10])
subplot(311),plot(tm0/60,Hm-HH,'-','LineWidth',1)
title(sprintf('%s;  Базисные %+8.2f'', %+8.2f''; градиент %+5.2f nT',...
    date0,D0*60,I0*60,dFgrad))
ylabel(sprintf('H%+5d, nT',-HH)), axis tight
set(gca,'XLim',[0 24],'XTick',0:3:24)
subplot(312),plot(tm0/60,Dm*60,'-','LineWidth',1)
axis tight, ylabel('D, arcmin')
set(gca,'XLim',[0 24],'XTick',0:3:24)
subplot(313),plot(tm0/60,Zm-ZZ,'-','LineWidth',1)
axis tight, ylabel(sprintf('Z%+5d, nT',-ZZ))
set(gca,'XLim',[0 24],'XTick',0:3:24)
xlabel('Time UT')
file=['graphs/' datestr(T0,'yyyy') '/' datestr(T0,'yyyymmdd') '_HDZ.png'];
print('-dpng', file)

figure(2)
subplot(411),plot(tm0/60,Xm-XX,'-','LineWidth',1)
title(sprintf('%s;  Базисные %+8.2f'', %+8.2f''; градиент %+5.2f nT',...
    date0,D0*60,I0*60,dFgrad))
ylabel(sprintf('X%+5d, nT',-XX)), axis tight
set(gca,'XLim',[0 24],'XTick',0:3:24)
subplot(412),plot(tm0/60,Ym-YY,'-','LineWidth',1)
ylabel(sprintf('Y%+5d, nT',-YY)), axis tight
set(gca,'XLim',[0 24],'XTick',0:3:24)
subplot(413),plot(tm0/60,Zm-ZZ,'-','LineWidth',1)
axis tight, ylabel(sprintf('Z%+5d, nT',-ZZ))
set(gca,'XLim',[0 24],'XTick',0:3:24)
subplot(414),plot(tm0/60,Fm-FF,'-','LineWidth',1)
axis tight, ylabel(sprintf('F%+5d, nT',-FF))
set(gca,'XLim',[0 24],'XTick',0:3:24)
xlabel('Time UT')

figure(3)
subplot(211),plot(tm0/60,Fm-FF,tm0/60,Fscalm-FF,'r')
title(sprintf('%s;  Базисные %+8.2f'', %+8.2f''; градиенты %+5.2f, %+5.2f nT',...
    date0,D0*60,I0*60,dFgrad,dFgrad_Scal))
axis tight, ylabel(sprintf('F%+d, nT',-FF))
legend('dIdD',ScalMag)
subplot(212),plot(tm0/60,Fm-Fscalm,'.-')
axis tight, ylabel(['F(dIdD)-F(' ScalMag '), nT'])
xlabel('Time UT')

% figure(1)
disp('-------------------------------------------------------------------------------')

