% ОБРАБОТКА ДАННЫХ ИЗМЕРЕНИИ С ПОМОЩЬЮ 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 - пути прописаны в инициализа
% 15.05.2013 - полная ревизия ПО (С.Ю.Хомутов, Паратунка)
% 11.12.2013 - сделан вывод минутных данных в формате IAGA2002
% 31.07.2015 - сделана возможность коррекции исходных данных до получения
%               минутных
% 16.01.2016 - в качестве независимого измерения F используется POS-1
% 09.08.2019 - выбор POS-1 или GSM-19W в качестве абсол. магнитометра
% 13.08.2019 - исправлено - не было вывода в mat-файл (после изменений
%               09.08.2019)
% 27.11.2020 - вывод файла iaga2002 сделан в папку года
% 18.05.2022 - исправлена ошибка при генерации F для текущих суток при
%               отсутствии файла (размерность была по предыдущим суткам)
% 20.06.2025 - адаптация под Linux (в путях "\" заменён на "/" и прописные
%               в именах папок)
%            - минутные Fscal вычисляются по Гауссу

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

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

if ~exist('date0','var') | isempty(date0),
    date0=strrep(datestr(now-1,24),'/','.');
end
str=input(['Введите дату DD.MM.YYYY [Enter, если ' date0 '] ........ '],'s');
if ~isempty(str),
    date0=str;
end
T0=my_datenum(date0);   % текущие сутки
[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
% ---------- предыдущие сутки
file_dIdD=[path_data_dIdD datestr(T0_1,'yyyy') '/' datestr(T0_1,'yyyymmdd') '.mag'];
if ~exist(file_dIdD),
    th_1=(0:dt:86400-1)'/3600;
    F_1=ones(length(th_1),1)*NaN; dI_1=F_1; dD_1=F_1;
    disp(['!!! Не могу найти файл с данными dIdD - использую NaNs...' file_dIdD])
else
    disp(['!!! Считываю файл ...' file_dIdD])
    [th_1,F_1,dI_1,dD_1,head]=didd_var_read_fun(file_dIdD);
end
T_1=datenum(year0_1,mon0_1,day0_1,th_1,0,0);
% Отбраковка задокументированных помех
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)
%     error('!!! Не могу найти файл с данными dIdD !!!')
    th=(0:dt:86400-1)'/3600;
    F=ones(length(th),1)*NaN; dI=F; dD=F;
    disp(['!!! Не могу найти файл с данными dIdD - использую NaNs...' file_dIdD])
else
    disp(['!!! Считываю файл ...' file_dIdD])
    [th,F,dI,dD,head]=didd_var_read_fun(file_dIdD);
end
T=datenum(year0,mon0,day0,th,0,0);
% Отбраковка задокументированных помех
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);

% Объединяем сутки
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

% % рассинхронизация - линейно интерполируем на равномерную сетку 2.5 сек
% disp('!!! Рассинхронизация таймера dIdD - интерполирую на равномерную сетку !!!')
% Ti=T; Fi=F; dDi=dD; dIi=dI; thi=th;
% T=datenum(year0,mon0,day0,0,0,0:dt:86400-1); T=T(:);
% th=(T-T0)*24; th=round(th*3600*10)/10/3600;
% F=interp1(Ti,Fi,T,'linear',NaN);
% dD=interp1(Ti,dDi,T,'linear',NaN);
% dI=interp1(Ti,dIi,T,'linear',NaN);

% Сглаживаем данные F
if ~exist(file_dIdD) % нет нормальных измерений, генерируем NaN (1440 знач.)
    Tm=datenum(year0,mon0,day0,0,0:1440-1,0)'; thm=(Tm-T0)*24;
    Fm=ones(1440,1)*NaN;
else
    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);
end

% учитываем текущий градиент 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)); % !!! не совсем корректно, т.к. F - несглажено

% Выбираем данные в 00 сек
[kdi,km]=ismember_digit(T,Tm,9);
Dm=ones(length(Tm),1)*NaN; Im=Dm;
Dm(km)=D(kdi); Im(km)=I(kdi);
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);
% disp('Минутные Fscal - spot в 00 сек')
% интерполируем данные скал. магнитометра на равномерную шкалу
% чтобы убрать рассинхронизацию
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')

% % Вывод по столбцам (включая F)
% [yyy,mmm,ddd,hh,mm,ss]=datevec(T0+tm0/60/24);
% fid=fopen(['Minute\' sprintf('%4d\\%4d%02d%02d',year0,year0,mon0,day0) '_minute.dat'],'wt');
% fprintf(fid,'%% %s\n',datestr(T0));
% fprintf(fid,'%% UT      D('')     H(nT)    Z(nT)\n');
% fprintf(fid,'%02d:%02d %8.3f %9.2f %9.2f\n',[hh round(mm+ss/60) Dm*60 Hm Zm]');
% fclose(fid);

% рисуем графики
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('-------------------------------------------------------------------------------')

