function wdc_write_fun(Th,Hh,Dh,Zh,Obs_Code)
% запись часовых данных посуточно в формате WDC
%
%       wdc_write_fun(Th,Hh,Dh,Zh,Obs_Code)
%
% Th       - сериальная дата
% Hh,Dh,Zh - часовые данные (H,Z - в нТл, D - в угл.мин.)
%           все вектора д.б. одинаковой длины. Если данные отсутствуют,
%           то д.б. заполнены NaN-ами
% Obs-Code - IAGA-код обсерватории
%
% результаты выводятся в файл wdc_temporary.dat

% 25.08.2018 - S. Khomutov, Paratunka, Kamchatka, Russia
% 18.09.2020 - сделан вывод по месяцам (D,H,Z)

% формат строки с данными
format_string='%s%02d%02d%c%02d    %02d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d%4d\n';

Th=unique(fix(Th)); % сериальная дата суточного блока (строки)
Nd=length(Th);

Obs=repmat(upper(Obs_Code),Nd,1); % вектор с кодами обсерватории
[year,mon,day]=datevec(Th); century=fix(year/100);
Ele_H=repmat('H',Nd,1); Ele_D=repmat('D',Nd,1); Ele_Z=repmat('Z',Nd,1);

Hd=reshape(Hh,24,Nd)'; % разбираем вектор по суткам/часам [Nd,24]
BaseH=repmat(fix(nanmean(Hh)/100),Nd,1); % базовый уровень H=BaseH*100+dH;
dH=round(Hd-repmat(BaseH*100,1,24)); % остаток для файла
Hmean=round(mean(dH')');       % среднесуточное значение
dH(isnan(dH))=9999; Hmean(isnan(Hmean))=9999;

Dd=reshape(Dh,24,Nd)'; % разбираем вектор по суткам/часам [Nd,24]
BaseD=repmat(fix(nanmean(Dh)/60),Nd,1); % базовый уровень D(deg) = BaseD + dD/10/60
dD=round((Dd-repmat(BaseD*60,1,24))*10); % остаток для файла
Dmean=round(mean(dD')');       % среднесуточное значение
dD(isnan(dD))=9999; Dmean(isnan(Dmean))=9999;

Zd=reshape(Zh,24,Nd)'; % разбираем вектор по суткам/часам [Nd,24]
BaseZ=repmat(fix(nanmean(Zh)/100),Nd,1); % базовый уровень H=BaseH*100+dH;
dZ=round(Zd-repmat(BaseZ*100,1,24)); % остаток для файла
Zmean=round(mean(dZ')');       % среднесуточное значение
dZ(isnan(dZ))=9999; Zmean(isnan(Zmean))=9999;

[yyy,mmm,ddd]=datevec(Th);
fid=fopen('wdc_temporary.dat','wt');
for km=min(mmm):max(mmm)    % перебор по месяцам
    days_in_month=find(km==mmm);    % номера дней для месяца km
    for k=1:length(days_in_month)   % перебор по суткам
        kd=days_in_month(k);
        fprintf(fid,format_string,...
            Obs(kd,:),year(kd)-century(kd)*100,mon(kd),Ele_D(kd),day(kd),...
            century(kd),BaseD(kd),dD(kd,:),Dmean(kd));
    end
    for k=1:length(days_in_month)   % перебор по суткам
        kd=days_in_month(k);
        fprintf(fid,format_string,...
            Obs(kd,:),year(kd)-century(kd)*100,mon(kd),Ele_H(kd),day(kd),...
            century(kd),BaseH(kd),dH(kd,:),Hmean(kd));
    end
    for k=1:length(days_in_month)   % перебор по суткам
        kd=days_in_month(k);
        fprintf(fid,format_string,...
            Obs(kd,:),year(kd)-century(kd)*100,mon(kd),Ele_Z(kd),day(kd),...
            century(kd),BaseZ(kd),dZ(kd,:),Zmean(kd));
    end
end
fclose(fid);
