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

% 31.03.2012 - Хомутов С.Ю.
% 24.10.2012 - адаптация для обработки данных MagDas
% 28.10.2012 - с учетом температурной зависимости (в ini_parameters.m)
% 11.03.2014 - учитывается возможность неполных данных
% 13.03.2014 - сделана отбраковка данных о температуре и наклоне
% 25.04.2014 - коррекция времени в данных при неполных файлах
% 13.07.2014 - скачки в данных заполняются NaN
% 29.11.2015 - сделано сравнение с компонентами по FGE
% 03.10.2022 - попытка скорректировать таймер для сек. данных
% 24.03.2024 - данные о коррекции таймера - в файле 
%               ini_parameters/Correct_timer_Magdas.dat
% 07.03.2024 - убираются данные с повторяющимися моментами, возникающие при
%               коррекции времени
% 28.10.2024 - убрана автомат. коррекция сдвига таймера сек. данных (надо
%               через файл Correct_timer_Magdas.dat)
% 07.12.2024 - оценка сдвига таймеров MAGDAS и FGE (строка около 265)
%            - коррекция сдвига таймера
% 08.12.2024 - оценка сдвига таймера MAGDAS делается по FRG
%            - поменял заглавные буквы на прописные в названиях папок 
%               и "\" на "/" - а-ля Linux
% 09.01.2025 - величина сдвига времени Magdas-FRG выводится в файл Delay_Magdas_FRG.dat

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

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

% считываем секундные данные MagDas (в нТл)
file_md=[path_data_magdas date0(7:10) '/S' datestr(T0,'yymmdd') '.ptk'];
disp(['... Считываются данные MagDas (1 sec) ... ' file_md])
[StartPos,EndPos,MagData,Header,NGPS,CORRECT_INF,STATDATA]=readrs1file(file_md,1);
% if StartPos>=0 & EndPos>=0,
disp(CORRECT_INF)

% считываем данные о коррекции времени
dth_corr=0; % инициализируем
file_corr='../ini_parameters/Correct_timer_Magdas.dat';
tmp=load(file_corr);
Tcorr=datenum(tmp(:,1),tmp(:,2),tmp(:,3));
kcorr=find(fix(Tcorr)==T0);
if ~isempty(kcorr)  % нужна коррекция времени
    if length(kcorr)>1
        disp('!!! Больше двух коррекций - пока не умею... Беру по первой...')
        kcorr=kcorr(1);
    end
    thcorr1=tmp(kcorr,4); thcorr2=tmp(kcorr,5);
    dth_corr=dms_deg(tmp(kcorr,6),tmp(kcorr,7),tmp(kcorr,8));
end

% EndPos=86400; % для нештатной обработки
if sum(~isnan(MagData(:,1)))>1,
    thjs=((StartPos:EndPos)'-1)/3600;    % время в часах (для полных суток)
    if dth_corr~=0,
        k=between_my(thjs,thcorr1,thcorr2);
        thjs(k)=thjs(k)+dth_corr;
        fprintf('!!!! Корректирую время на %+7d сек\n',round(dth_corr*3600));
    end
    Tjs=T0+thjs/24;
% EndPos=EndPos-StartPos+1;StartPos=1; % Для обработки 03.02.2017
% tmp=ones(86400,4)*NaN; % для обработки 22.11.2017
    tmp=ones(EndPos-StartPos+1,4)*NaN;
    tmp(StartPos:EndPos,:)=MagData; MagData=tmp; % - " -
    dHjs=MagData(StartPos:EndPos,1);
    dDjs=MagData(StartPos:EndPos,2);
    dZjs=MagData(StartPos:EndPos,3);
    [thjs,m,n]=unique(round(thjs*3600)); thjs=thjs/3600;
    if ~isempty(m), dHjs=dHjs(m); dDjs=dDjs(m); dZjs=dZjs(m); Tjs=T0+thjs/24; end
%     if length(MagData(:,1))<86400,
%         disp('---!!!  Неполный файл секундных данных? Корректирую время !!! ---')
%         Tjs=magdas_timer_correct(Tjs,dHjs); thjs=(Tjs-T0)*24;
%     end
else
    disp('!!!  Проблемы с файлом секундных данных, использую NaN !!!')
    thjs=((1:86400)'-1)/3600;    % время в часах (для полных суток)
    Tjs=T0+thjs/24;  N=length(thjs);
    dHjs=ones(N,1)*NaN; dDjs=dHjs; dZjs=dHjs;
end
clear MagData

% считываем минутные данные MagDas (в нТл)
file_md=[path_data_magdas date0(7:10) '/M' datestr(T0,'yymmdd') '.ptk'];
disp(['... Считываются данные MagDas (1 min) ... ' file_md])
[StartPos,EndPos,MagData,Header_m,NGPS_m,CORRECT_INF_m,STATDATA_m]=readrm1file(file_md,1);
% if StartPos>=0 & EndPos>=0,

% EndPos=1440;
% EndPos=EndPos-StartPos+1;StartPos=1; % Для обраотки 03.02.2017
if sum(~isnan(MagData(:,1)))>1,
    thjm=((StartPos:EndPos)'-0.5)/60;    % время в часах (для полных суток)
    Tjm=T0+thjm/24;
    Tjm=Tjm+dth_corr/24; thjm=thjm+dth_corr;
    tmp=ones(1440,7)*NaN; % для обработки 22.11.2017
    tmp(StartPos:EndPos,:)=MagData; MagData=tmp; % - " -
    dHjm=MagData(StartPos:EndPos,1); 
    dDjm=MagData(StartPos:EndPos,2); dZjm=MagData(StartPos:EndPos,3); % компоненты поля
    EWjm=MagData(StartPos:EndPos,5); NSjm=MagData(StartPos:EndPos,6); 
    TRjm=MagData(StartPos:EndPos,7); % наклоны и температура
    if length(TRjm)<1440,
        disp('---!!!  Неполный файл минутных данных? Корректирую время !!! ---')
        Tjm=magdas_timer_correct(Tjm,TRjm); thjm=rem(Tjm,1)*24;
    end
else
    disp('!!!  Проблемы с файлом минутных данных, использую NaN !!!')
    thjm=((1:1440)'-0.5)/60;
    Tjm=T0+thjm/24; N=length(thjm);
    dHjm=ones(N,1)*NaN; dDjm=dHjm; dZjm=dHjm; % компоненты поля
    EWjm=dHjm; NSjm=dHjm; TRjm=dHjm; % наклоны и температура
end
clear MagData

% принудительная чистка данных при скачках уровня  -  Rev.13.07.2014
dHjs(dHjs<dH_Lim)=NaN; dDjs(dDjs>dD_Lim)=NaN; dZjs(dZjs>dZ_Lim)=NaN;
dHjm(dHjm<dH_Lim)=NaN; dDjm(dDjm>dD_Lim)=NaN; dZjm(dZjm>dZ_Lim)=NaN;
EWjm(EWjm<EW_Lim)=NaN; NSjm(NSjm<EW_Lim)=NaN; TRjm(TRjm<TR_Lim)=NaN;

% Отбраковка задокументированных помех
dHjs=remove_errors('../ini_parameters/Remove_MagDas.dat',T0,thjs,dHjs);
dDjs=remove_errors('../ini_parameters/Remove_MagDas.dat',T0,thjs,dDjs);
dZjs=remove_errors('../ini_parameters/Remove_MagDas.dat',T0,thjs,dZjs);
TRjm=remove_errors('../ini_parameters/Remove_MagDas_TR.dat',T0,thjm,TRjm);
NSjm=remove_errors('../ini_parameters/Remove_MagDas_TR.dat',T0,thjm,NSjm);
EWjm=remove_errors('../ini_parameters/Remove_MagDas_TR.dat',T0,thjm,EWjm);

% Учет температурной зависимости
if sum(~isnan(TRjm))>1,
    if dH_dTR~=0 | dD_dTR~=0 | dZ_dTR~=0,
        disp('... Температурная зависимость учтена, значения температуры < 0.0oC отбракованы ...')
    else
        disp('... Нулевые температурные коэффициенты! ...')
    end
    k=TRjm>0;
    p=0.9999; TRjs=csaps(Tjm(k),TRjm(k),p,Tjs); % аппроксимируем минутную температуру на 1 сек данные
else
    disp('!!! Температурные поправки не учитываются !!!')
    TRjs=ones(length(Tjs),1)*TR0;
end

dHjs=dHjs-dH_dTR*(TRjs-TR0);
dDjs=dDjs-dD_dTR*(TRjs-TR0);  % !!! ВНИМАТЕЛЬНО!!! коэф. в угл.мин или нТл !!!
dZjs=dZjs-dZ_dTR*(TRjs-TR0);

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

% Вычисление полных составляющих поля и модуля
Hjs=sqrt((H0j+dHjs).^2+dDjs.^2);
% H=H0+dHv;
Zjs=Z0j+dZjs;
dDja=r2d(atan(dDjs./Hjs)); Djs=D0j+dDja;
Fjs=total_field(Hjs,Zjs);

HHj=base_to_extract(dHjs); ZZj=base_to_extract(dZjs);

figure(1)
k=isnan(dHjs); dH_NaN=ones(sum(k),1)*nanmean(dHjs);
subplot(311),plot(thjs,dHjs-HHj,thjs(k),dH_NaN-HHj,'r.')
ylabel(sprintf('dH%+d, nT',-HHj))
title([date0 ' ... с учетом температурной зависимости'])
k=isnan(dDjs); dD_NaN=ones(sum(k),1)*nanmean(dDja);
subplot(312),plot(thjs,dDja*60,thjs(k),dD_NaN*60,'r.'), ylabel('dD, arcmin')
k=isnan(dZjs); dZ_NaN=ones(sum(k),1)*nanmean(dZjs);
subplot(313),plot(thjs,dZjs-ZZj,thjs(k),dZ_NaN-ZZj,'r.')
ylabel(sprintf('dZ%+d, nT',-ZZj))
xlabel('Time UT, hours'), %set(gca,'XLim',[0 24],'XTick',0:3:24)
file=['graphs/' datestr(T0,'yyyy') '/' datestr(T0,'yyyymmdd') '_dHdDdZ_TR.png'];
print('-dpng', file)

% считываем данные FGE+GSM (исходные) или FRG
% backup_var='FGE';
backup_var='FRG';
[Tv,H,D,Z,Tf,F]=read_backup_var(date0,'frg');
thv=(Tv-T0)*24; thf=(Tf-T0)*24;

% сохраняем текущие сутки в mat-файл
file_mat=['mat/' num2str(date0(7:10)) '/magdas_' datestr(T0,'yyyymmdd') '_TR'];
save(file_mat,'Tjs','dHjs','dDjs','dZjs','H0j','D0j','Z0j','Tf','F',...
              'Tjm','dHjm','dDjm','dZjm','EWjm','NSjm','TRjm');

figure(2)
HHp=base_to_extract(Hjs); ZZp=base_to_extract(Zjs);
subplot(311),plot(thjs,Hjs-HHp,thv,H-HHp,'r'), axis tight
title([date0 sprintf(';   BLV = %+8.2f, %+9.4f, %+8.2f за ',H0j,D0j,Z0j) datestr(T_bl(kbl))])
ylabel(sprintf('H%+d, nT',-HHp)), legend('Magdas',backup_var,'Location','best')
subplot(312),plot(thjs,Djs*60,thv,D*60,'r'), axis tight
ylabel('D, arcmin')
subplot(313),plot(thjs,Zjs-ZZp,thv,Z-ZZp,'r'), axis tight
ylabel(sprintf('Z%+d, nT',-ZZp)), xlabel('Time UT, hours')
file=['graphs/' num2str(date0(7:10)) '/' datestr(T0,'yyyymmdd') '_HDZ_TR.png'];
print('-dpng', file)

% if max(fv)<50000, return, end   % не рисуем для вариаций

figure(3)
subplot(211),plot(thjs,Fjs-FF,thf,F-FF,'r'), axis tight
title([date0 ...
    sprintf('; dD_dTR=%+6.2f, dD_dTR=%+6.2f, dZ_dTR=%+6.2f',dH_dTR,dD_dTR,dZ_dTR)])
ylabel(sprintf('F%-5d, nT',-FF))
legend('F(MagDas)',['F(' backup_var ')'],'Location','best')
[kj,kf]=ismember_digit(thjs,thf,8);
thdf=thjs(kj); df=Fjs(kj)-F(kf);
subplot(212),plot(thdf,df)
ylabel(['F(MagDas)-F(' backup_var '), nT']), xlabel('UT, hours')
file=['graphs/' num2str(date0(7:10)) '/' datestr(T0,'yyyymmdd') '_F(MagDas)_F(GSM)_TR.png'];
print('-dpng', file)

figure(4)
subplot(311),plot(thjm,TRjm,thjs,TRjs,'r'),title(date0)
ylabel('Temperature, ^oC')
subplot(312),plot(thjm,EWjm*3600), ylabel('Tilt, EW, arcsec')
subplot(313),plot(thjm,NSjm*3600), ylabel('Tilt, NS, arcsec'), xlabel('Time UT, hours')
file=['graphs/' num2str(date0(7:10)) '/' datestr(T0,'yyyymmdd') '_Temperature_Tilt.png'];
print('-dpng', file)

figure(5)
[kj,kv]=ismember_digit(thjs,thv,8);
subplot(311),plot(thjs(kj),Hjs(kj)-H(kv))
title([date0 sprintf(';   BLV = %+8.2f, %+9.4f, %+8.2f за ',H0j,D0j,Z0j) datestr(T_bl(kbl))])
ylabel(['H(Magdas-' backup_var '), nT'])
subplot(312),plot(thjs(kj),(Djs(kj)-D(kv))*60)
ylabel(['D(Magdas-' backup_var '), arcmin'])
subplot(313),plot(thjs(kj),Zjs(kj)-Z(kv))
ylabel(['Z(Magdas-' backup_var '), nT']), xlabel('Time UT, hours')
file=['graphs/' num2str(date0(7:10)) '/' datestr(T0,'yyyymmdd') '_HDZ(MagDas-FGE)_TR.png'];
print('-dpng', file)

% оценка сдвига таймеров MAGDAS и FRG (08.12.2024)
file_mat_frg=['d:/Paratunka/japan_frg/mat/' date0(7:10) ...
          '/FRG_' datestr(T0,'yyyymmdd') '.mat'];
load(file_mat_frg,'Tj','dHj','dDj','dZj')
thj=(Tj-T0)*24;
if sum(~isnan(dHj))>0
    % оценка сдвига по времени
    [C,lag,lag_max_C]=...
        delay_sec_data_fun(thjs,Hjs-nanmean(Hjs),thj,dHj-nanmean(dHj),60,1);
    subplot(211) 
    title(['MAGDAS ' date0 '; Сравнение c FRG'])
    legend('MAGDAS(Japan)','FRG(Japan)'), ylabel('H-scaps(0.1), nT')
    subplot(212)
    title(['MAGDAS ' date0 '; Сдвиг H(magdas) относительно H(frg)), максимум на '...
           sprintf('%+4.1f сек',-lag_max_C)])
end
file_delay='Delay_Magdas_FRG.dat';
disp('Log оценок сдвига - см. файл Delay_Magdas_FRG.dat')
vecT=datevec(T0); tmp=[vecT(1:3) -lag_max_C];
fid=fopen(file_delay,'at');
fprintf(fid, '%4d %02d %02d   %+6.1f\n',tmp');
fclose(fid);
