function [C,lag,max_delay]=delay_sec_data_fun(tx,x,ty,y,range_t,flag_plot)
% оценка сдвига между двумя рядами секундных данных
%
%   [C,lag]=delay_sec_data_fun(tx,x,ty,y,range_t[,flag_plot]);
%
% tx,x  -   1-й ряд, tx - в часах
% ty,y  -   2-й ряд, ty - в часах
% range_t - диапазон  сдвигов корр. ф-ции (+/-range_t), в сек
% flag_plot - если 1 - рисуем график рядов и корр. ф-ции, если 0 - не рисуем
%
% С     - корр. ф-ция
% lag   - сдвиги (в сек)
% max_delay - сдвиг на максимуме корр. ф-ции

% 20.11.2023 - S. Khomutov, Paratunka, Kamchatka
% 22.11.2023 - сделан вывод графика (при запросе)
% 23.11.2023 - более аккуратно при расчётах с пропусками
%            - вывод сдвига при максимуме корр. ф-ции

% равномерная сетка через 1 сек от мин до макс значений tx,ty
% без NaNs на концах
kx=~isnan(x);ky=~isnan(y);
t=(max([min(tx(kx)) min(ty(ky))])*3600:1:min([max(tx(kx)) max(ty(ky))])*3600)/3600;
t=t(:);
% интерполируем измерен. знач. на равном. сетку
xi=interp1(tx(~isnan(x)),x(~isnan(x)),t);
yi=interp1(ty(~isnan(y)),y(~isnan(y)),t);
% убираем медленные вариации
dx=xi-csaps(t(~isnan(xi)),xi(~isnan(xi)),0.1,t);
dy=yi-csaps(t(~isnan(yi)),yi(~isnan(yi)),0.1,t);
% dx=xi-csaps(t,xi,0.1,t);
% dy=yi-csaps(t,yi,0.1,t);
% смотрим корр. ф-цию
[C,lag]=xcorr(dy,dx,range_t,'coeff');
lag=lag(:);
if flag_plot==1 % рисуем график
    figure(9)
    subplot(211), plot(t,dx,t,dy)
    axis tight, legend('канал 1','канал 2')
    ylabel('dVar, nT'), xlabel('Время UT, час.')
    % ищем положение максимума корр. ф-ции
    pc=0.9999; lag01=min(lag):0.1:max(lag); lag01=lag01(:);
    C01=csaps(lag,C,pc,lag01);
    k01=find(C01==max(C01)); max_delay=lag01(k01);
    subplot(212), plot(lag,C,'o',lag01,C01,'.')
    xlabel('Сдвиги, сек')
    title(sprintf('Сдвиг таймера канала 2 относительно таймера канала 1 = %+4.1f c',lag01(k01)))
end
end

