function [T,MinVal,Th,HourVal,Td,DayVal,Tk,Kindex,head]=cdrom_bin_read_fun(year0,mon0,nMon,imo_code,flag_plot)
% считывание бинарных файлов
%
% [T,MinVal,Th,HourVal,Td,DayVal,Tk,Kindex,head]=
%               cdrom_bin_read_fun(year0,mon0,nMon,imo_code,flag_plot);
%
% year0, mon0   - start date (year, month)
% nMon          - number of month
% imo_code      - code of observatory (IMO) или cell({path_of_data
%                   imo_code})
% flag_plot     - if =1 - data are plotted
%
% T,MinVal      - сериальная дата и минутные данные X/H, Y/D, Z, G(dF)
% Th,HourVal    - сериальная дата и часовые данные X/H, Y/D, Z, G(dF)=NaN
% Td,DayVal     - сериальная дата и суточные данные X/H, Y/D, Z, G(dF)=NaN
% Tk,Kindex     - сериальная дата и K-индексы 
% head          - заголовок последних суток

% 23.03.2010 - модифицировано с учетом нового формата IAF v2.00
% 01.05.2011 - учтен формат записи для системы ориентации данных 
%               (для версии IAF = 3) - 'XYZG' или ' XYZ'
% 31.03.2018 - сделана функция из скрипта
% 12.05.2018 - можно задавать только код IMO, а можно - путь к данным и код
%               (с массиве ячеек)
% 10.02.2019 - более корректно учтён знак вычитаемого уровня в подписях на осях ординат
% 10.03.2019 - минимальное F м.б. 20000 нТл
% 23.05.2019 - исправлена логика вывода deltaF на график (~k)

T1=datenum(year0,mon0,1);   % init serial data (1st day of 1st month)
T2=datenum(year0,mon0+nMon,1)-1;  % final serial data (last day of last month)
Nday=T2-T1+1;               % кол-во суток во всем интервале
% mjd1=cdtmjd(year0,mon,1);

if ~iscell(imo_code)  % задан код IMO - путь по умолчанию
    file_path=['CD\' imo_code '\' datestr(T1,'yyyy') '\'];
else    % имеем массив ячеек {dir_of_data imo-code}
    file_path=imo_code{1};
    imo_code=imo_code{2};
end

% инициализация массивов (минутных, часовых, суточных)
E1=zeros(Nday*1440,1)*NaN; E2=E1; Z=E1; dF=E1; T=E1;
E1h=zeros(Nday*24,1)*NaN; E2h=E1h; Zh=E1h; dFh=E1h; Th=E1h;
E1d=zeros(Nday,1)*NaN; E2d=E1d; Zd=E1d; dFd=E1d; Td=E1d;
Kindex=zeros(Nday*8,1)*NaN; Tk=Kindex;

kd_all=1;   % счетчик суток
for kmon=1:nMon,
    Tmon=datenum(year0,mon0+kmon-1,1);   % current serial date
    file=[file_path imo_code datestr(Tmon,'yymmm') '.bin'];
    if ~exist(file),
        error(['!!! No found file - ' file])
    end
    disp(file)
    nd_mon=datenum(year0,mon0+kmon,1)-Tmon;  % кол-во суток в текущем месяце
    fid=fopen(file,'r');
    for kd=1:nd_mon,    % processing for every day
        if feof(fid), break, end
        % Header
        head=cdrom_bin_read_head(fid);
        % One-minute data (1440 minutes x 4 components)
        k1=(kd_all-1)*1440+1; k2=kd_all*1440;
        E1(k1:k2)=fread(fid,1440,'int32');
        E2(k1:k2)=fread(fid,1440,'int32');
        Z(k1:k2)=fread(fid,1440,'int32');
        dF(k1:k2)=fread(fid,1440,'int32');
        T(k1:k2)=T1+kd_all-1+(0:1440-1)'/60/24;
        % Mean hourly data (24 hours x 4 components)
        k1=(kd_all-1)*24+1; k2=kd_all*24;
        E1h(k1:k2)=fread(fid,24,'int32');
        E2h(k1:k2)=fread(fid,24,'int32');
        Zh(k1:k2)=fread(fid,24,'int32');
        dFh(k1:k2)=fread(fid,24,'int32');
        Th(k1:k2)=T1+kd_all-1+(0.5:23.5)'/24;
        % Daily mean data (1 day x 4 components)
        E1d(kd_all)=fread(fid,1,'int32');
        E2d(kd_all)=fread(fid,1,'int32');
        Zd(kd_all)=fread(fid,1,'int32');
        dFd(kd_all)=fread(fid,1,'int32');
        Td(kd_all)=T1+kd_all-1+0.5;
        % 8 digitally derived K-values. Missing value = 999
        k1=(kd_all-1)*8+1; k2=kd_all*8;
        Kindex(k1:k2)=fread(fid,8,'int32');
        Tk(k1:k2)=T1+kd_all-1+(1.5:3:23)'/24;
        tmp=fread(fid,4,'int32');   % reserved

        kd_all=kd_all+1;
    end
    fclose(fid);
end

% T=mjd_ser(MJD); Th=mjd_ser(MJD_h); Td=mjd_ser(MJD_d); Tk=mjd_ser(MJD_i);

E1(E1==999999|E1==888888)=NaN; E2(E2==999999|E2==888888)=NaN;
    Z(Z==999999|Z==888888)=NaN; dF(dF==999999|dF==888888)=NaN;
E1h(E1h==999999)=NaN; E2h(E2h==999999)=NaN;
    Zh(Zh==999999)=NaN; dFh(dFh==999999|dFh==888888)=NaN;
E1d(E1d==999999)=NaN; E2d(E2d==999999)=NaN;
    Zd(Zd==999999)=NaN; dFd(dFd==999999|dFd==888888)=NaN;
Kindex(Kindex==999)=NaN;
    
E1=E1/10;E2=E2/10;Z=Z/10;dF=dF/10;
E1h=E1h/10;E2h=E2h/10;Zh=Zh/10;dFh=dFh/10;
E1d=E1d/10;E2d=E2d/10;Zd=Zd/10;dFd=dFd/10;
Kindex=Kindex/10;

% matrix to out
MinVal=[E1 E2 Z dF]; HourVal=[E1h E2h Zh dFh]; DayVal=[E1d E2d Zd dFd];

fprintf('Измеряемые компоненты - %s\n',head.Orientation)
disp('------ Пропущено минутных данных за этот период:')
disp([sum(isinf(E1)|isnan(E1)) sum(isinf(E2)|isnan(E2))...
        sum(isinf(Z)|isnan(Z)) sum(isinf(dF)|isnan(dF))])

if flag_plot~=1, return, end % выход, если не рисовать
    
% переименовываем в "естественные" компоненты, в зависимости от записи
Orient=head.Orientation;
if Orient(1)==' ',  % без регистрации F(Scal)
    k_shift=1;
else
    k_shift=0;      % имеются данные F(Scal)
end
for k=1:2,
    eval([Orient(k+k_shift) '=E' int2str(k) ';']);
    eval([Orient(k+k_shift) 'h=E' int2str(k) 'h;']);
    eval([Orient(k+k_shift) 'd=E' int2str(k) 'd;']);
end

if exist('D'), 
    D=D/60;
    Fv=total_field(H,Z);
    Fvh=total_field(Hh,Zh);
    Fvd=total_field(Hd,Zd);
else
    Fv=total_field(X,Y,Z);
    Fvh=total_field(Xh,Yh,Zh);
    Fvd=total_field(Xd,Yd,Zd);
end

EE1=base_to_extract(E1);
EE2=base_to_extract(E2);
ZZ=base_to_extract(Z);
if head.IAF_ver<2,
    FF=base_to_extract(dF);
else
    FF=base_to_extract(Fv);
end

figure(1)
subplot(311)
plot(T,E1-EE1,Th,E1h-EE1,'m.',Td,E1d-EE1,'ro','MarkerFaceColor','g')
if nMon==1,
    datetick('x',7)
else
    datetick('x')
end
ylabel([Orient(1+k_shift) sprintf('%+d',-EE1)])
title([upper(imo_code) ' - ' datestr(T1,'yyyy')])

subplot(312)
plot(T,E2-EE2,Th,E2h-EE2,'m.',Td,E2d-EE2,'ro','MarkerFaceColor','g')
if nMon==1,
    datetick('x',7)
else
    datetick('x')
end
ylabel([Orient(2+k_shift) sprintf('%+d',-EE2)])

subplot(313)
plot(T,Z-ZZ,Th,Zh-ZZ,'m.',Td,Zd-ZZ,'ro','MarkerFaceColor','g')
if nMon==1,
    datetick('x',7)
else
    datetick('x')
end
ylabel(['Z' sprintf('%+d',-ZZ)])

figure(2)       % total filed F
if head.IAF_ver<2,
    plot(T,dF-FF,Th,dFh-FF,'m.',Td,dFd-FF,'ro','MarkerFaceColor','g')
    ylabel(['F(scal)' sprintf('%+d',-FF)])
else
    plot(T,Fv-FF,Th,Fvh-FF,'m.',Td,Fvd-FF,'ro','MarkerFaceColor','g')
    ylabel(['F(var)' sprintf('%+d',-FF)])
end
if nMon==1,
    datetick('x',7)
else
    datetick('x')
end
title([upper(imo_code) ' - ' datestr(T1,'yyyy')])

figure(3)   % total field difference dF=F(v)-F(s)
if head.IAF_ver<2,
    plot(T,Fv-dF,Th,Fvh-dFh,'r.',...
        Td,Fvd-dFd,'ko','MarkerFaceColor','y')
else
    k=dF<-20000;        % убираем значения dF, которые = -F(s)
    plot(T(~k),dF(~k),Th,dFh,'r.',Td,dFd,'ko','MarkerFaceColor','y')
end
if nMon==1,
    datetick('x',7)
else
    datetick('x')
end
title([upper(imo_code) ' - ' datestr(T1,'yyyy')])

figure(4)
stem(Tk,Kindex,'.'), ylabel('K-index'), datetick('x')
title(sprintf('K9 = %4d',head.K9))
set(gca,'XLim',[round(min(Tk)) round(max(Tk))])

