function [Th,Hourly,Elements,Td,Daily]=wdc_hourly_read_fun(file,ReqEle)
% Считывание файла часовых значений магнитного поля в формате IAGA (WDC)
%
%  [Th,Hourly[,Elements[,Td,Daily]]]=wdc_hourly_read_fun(file,ReqEle);
%
% file      - файл с данными, с полным путём
% ReqEle    - запрашиваемая компонента ('D','H' и т.п.)
%
% Th,Hourly - сериальная дата и часовые значения (D,I - в град., 
%               X,Y,H,Z,F - в нТл)
% Elements  - список имеющихся в файле составляющих поля
% Td,Daily  - сериальная дата и среднесуточные значения

% 20.05.2016 - S.Khomutov, Paratunka, Kamchatka
% 27.05.2016 - исправлен разбор по часам, вместо sscanf() используется цикл
%               01:24
% 14.06.2018 - заполнение NaN-ами отсутствующей в файле компоненты
% 18.09.2020 - исправлено - не заполнялись NaN-ами суточные средние

% COLUMNS   FORMAT   DESCRIPTION
% 1-3       A3       Observatory 3-letter code
% 4-5       I2       Year. Last 2 digits, 1996 = 96. See also columns 15-16.
% 6-7       I2       Month (01-12)
% 8         Al       Element (D,I,H,X,Y,Z, or F)
% 9-10      I2       Day of month (01-31)
% 11-12     A2       Blanks
% 13        A1       Preliminary or Definitive data.
%                    Preliminary = P , Definitive = D
%                    Intermagnet CD-ROM = I   
% 14        A1       Arbitrary
% 15-16     I2       Century digits
%                    Year = 1889, Century digits = 18
%                    Year = 1996, Century digits = 19
%                    Year = 2014, Century digits = 20
% 17-20     I4       Tabular base.
%                    The values are in degrees for D and I, and in 
%                    hundreds of nanoTeslas for the intensity elements. 
% 21-116    24I4     24 4-digit hourly mean values for the day. 
%                    The values are in tenth-minutes for D and I, and in
%                    nanoTeslas for the intensity elements.
%                    The first hourly mean value represents the mean value 
%                    between 00:00 UT and 01:00 UT, ..., the 24th value 
%                    represents the mean between 23:00 UT and 24:00 UT.
%                    A missing value is identified by 9999.
% 117-120   I4       Daily Mean. 
%                    If any of the hourly mean values for the day are missing 
%                    9999 will appear as the daily mean.
% 121-122            Record end marker.
%                    Two chars 'cr'= 13 and 'nl'= 10.
% The tabular bases, hourly mean values and daily means are right adjusted and 
% signed if negative. 
% Negative values are identified with a minus sign either adjacent to the first 
% significant digit or in the high-order position  of the field. 
% NOTE: A blank digit will not appear between a (-) sign and the first 
% significant digit. For example, a value may appear as -050 or b-50 but not 
% as -b50 (b=blank).
% The 25 values in positions 21-120 will have the range -999 to 9998, with 9999 
% reserved for missing values. 
% The records are sorted according to observatory code, year, month, element,
% day (positions 1-10).
% For the intensity elements we have that
%      hourly value (nanoTeslas) = tab.base*100 + tab.value
% For the angles D and I we have that
%      hourly value (Degrees) = tab.base + tab.value/600

if nargin<2,
    error('!!! Не указана требуемая составляющая поля !!!')
end

if ~exist(file,'file'),
    Th=NaN; Hourly=NaN; Elements=NaN; Td=NaN; Daily=NaN;
    disp(['!!! Не могу найти файл - ' file]);
    return
end

% считываем данные в бинарном виде
fid=fopen(file,'r');
    tmp=fread(fid,Inf,'uint8');
fclose(fid);
Klf=find(tmp==10); % определяем номер символа LF (конец строки)
% превращаем в строковый массив
str=reshape(char(tmp),Klf(1),length(tmp)/Klf(1))';
[Nrec,Nstr]=size(str);

% определяем записанные компоненты
Elements=unique(str(:,8))';
if isempty(strfind(Elements,ReqEle)), % 14.06.2018
    disp(['!!! В файле нет данных для составляющей - ' ReqEle])
    year=str2num(str(1,4:5))+str2num(str(1,15:16))*100;
    T1=datenum(year,1,1,0,30,0); T2=datenum(year,12,31,23,59,59);
    Th=T1:1/24:T2; Th=Th(:);
    Hourly=ones(length(Th),1)*NaN;
    Td=fix(T1):1:T2; Td=Td(:);
    Daily=ones(length(Td),1)*NaN;
    return
end

ObsCod=str(:,1:3);  % код обсерватории
Ele=str(:,8);       % элемент поля
year=str2num(str(:,4:5))+str2num(str(:,15:16))*100; % год
mon=str2num(str(:,6:7));    % месяц
day=str2num(str(:,9:10));   % день
TabBase=str2num(str(:,17:20)); % базовый уровень в таблице
Daily=str2num(str(:,117:120));

aaa=str(:,21:116);
HourlyDay=ones(Nrec,24); tmp=[];
for kr=1:Nrec,
    for kh=1:24,
        k1=20+(kh-1)*4+1; k2=20+kh*4;
        tmp(kh)=str2num(str(kr,k1:k2));
    end
    HourlyDay(kr,:)=tmp;
end

% заполняем пропуск (9999) на NaN
HourlyDay(HourlyDay==9999)=NaN;
Daily(Daily==9999)=NaN;

% Отбор для требуемой компоненты
kc=Ele==ReqEle;
Td=datenum(year(kc),mon(kc),day(kc));
TabBase=TabBase(kc);
Daily=Daily(kc);

if ReqEle=='D' | ReqEle=='I',
    Hourly=repmat(TabBase,1,24)+HourlyDay(kc,:)/600;
    Daily=TabBase+Daily/600;
else
    Hourly=repmat(TabBase*100,1,24)+HourlyDay(kc,:);
    Daily=TabBase*100+Daily;
end
Hourly=Hourly'; Hourly=Hourly(:);

Th=repmat(Td,1,24)+repmat((0.5:23.5)/24,length(Td),1);
Th=Th'; Th=Th(:);
