function [T,X,Y,Z,D,H,F]=imag_read_fun(code_obs,code_fig);
% Считывание минутных данных из стандартных файлов INTERMAGNET
% длЯ обсерватории с кодом cod_obs
%
%     [T,X,Y,Z,D,H,F]=imag_read_fun(code_obs,code_fig);
%
% code_obs - код обсерватории вида "nvs" или имя файла с путем, если не в
%           текущем подкаталоге)
% code_fig - флаг вывода графики 
%            =1 - рисуем, =0 - нет (по умолчанию =0)
%
% T         - сериальная дата
% X,Y,Z,H,F - в нТл, D - в град.
%
% Минутные файлы должны находитьсЯ в поддиректории tmp\
%

% С.Ю.Хомутов
% 26.05.2010    - адаптировано под сериальную дату и имя файла
% 10.01.2018    - учитывается вариант списка компонент 'XYZS' и 'HDZS'

if nargin<2, code_fig=0; end    % не рисуем по-умолчанию

C1=[];C2=[];X=[];Y=[];Z=[];F=[];t=[];
path_file='tmp\';

if length(code_obs)==3,  % задан код обсерватории, а не имя файла
    D_file=dir([path_file '*.' code_obs]);
    if length(D_file)==0,
        error(['!!! Не могу найти файлы tmp\*.' code_obs ' !!!'])
    end

    file=char(D_file.name); % имена суточных файлов вида "MONDDYY.COD"
    file=imag_file_sort(file);
    % ОпределЯем начальную дату
    tmp=str2num(file(1,6:7));
    if tmp<50
        y_str=num2str(tmp+2000);
    else
        y_str=num2str(tmp+1900);
    end
    tmp_s=[file(1,4:5) '-' file(1,1:3) '-' y_str];
    tmp=datevec(tmp_s);
    T0=cdtmjd(tmp(1),tmp(2),tmp(3));
    disp(['Начальная дата - ' tmp_s])
else
    file=code_obs;
    path_file=[];
end

[ND,tmp]=size(file); N=ND*1440;
C1=zeros(N,1); C2=zeros(N,1); Z=zeros(N,1); F=zeros(N,1); t=zeros(N,1);
C1k=zeros(60,1); C2k=zeros(60,1); Zk=zeros(60,1); Fk=zeros(60,1);
T=zeros(N,1);

k1=1; k2=k1+59;
for kd=1:ND,
    fid=fopen([path_file file(kd,:)],'r');
    disp(file(kd,:))
    for kh=0:23,
       str=fgets(fid);
       if str==-1, break, end
       hh=str2num(str(17:18)); code=str(20:23); tmp=str(5:12);
       for km=1:30,
          str=fgets(fid);
          n1=(km-1)*2+1; n2=n1+1;
          C1k(n1:n2)=[str2num(str(1:7)); str2num(str(33:39))];
          C2k(n1:n2)=[str2num(str(9:15)); str2num(str(41:47))];
          Zk(n1:n2)=[str2num(str(17:23)); str2num(str(49:55))];
          Fk(n1:n2)=[str2num(str(25:30)); str2num(str(57:62))];
       end
       C1(k1:k2)=C1k; C2(k1:k2)=C2k; Z(k1:k2)=Zk; F(k1:k2)=Fk;
       k1=k1+60; k2=k2+60;
    end
    fclose(fid);
    Td=datenum([tmp(4:5) '-' tmp(1:3) '-20' tmp(6:7)]);
    T((kd-1)*1440+1:kd*1440)=Td+(0:1440-1)'/1440;
end

k=C1==999999|C1==9999999; C1(k)=NaN; k=C2==999999|C2==9999999; C2(k)=NaN;
k=Z==999999|Z==9999999; Z(k)=NaN; kf=F==999999; F(kf)=NaN;

switch code
    
    case {'XYZF','XYZS'}
        X=C1/10; Y=C2/10; Z=Z/10; F=F/10; Fv=sqrt(X.^2+Y.^2+Z.^2);
        if mean(Fv(~isnan(Fv)))<10000  % в файле - вариации dX,dY,dZ
            H=zeros(length(Fv),1)*NaN; D=H; Fv=H; % не вычислЯем H, D, Fv
        else
            H=sqrt(X.^2+Y.^2); D=r2d(atan2(Y,X));
        end
        % Базовый уровень вида XX000 нТл, длЯ вычитаниЯ при plot
        XX=base_to_extract(X); YY=base_to_extract(Y); ZZ=base_to_extract(Z); 
        HH=base_to_extract(H); FF=base_to_extract(F);
        DD=base_to_extract(D);
        tit1='Components from data files';
        tit2='Calculated components';

    case {'HDZF','HDZS'}    % D - in arc.min.
        H=C1/10; D=C2/100/60; Z=Z/10; F=F/10; Fv=total_field(H,Z);
        if mean(Fv(~isnan(Fv)))<10000  % в файле - вариации dH,dZ
            X=zeros(length(Fv),1)*NaN; Y=X; Fv=X; % не вычислЯем X, Y, Fv
        else
            X=H.*cos(d2r(D)); Y=H.*sin(d2r(D));
        end
        % Базовый уровень вида XX000 нТл, длЯ вычитаниЯ при plot
        XX=base_to_extract(X); YY=base_to_extract(Y); ZZ=base_to_extract(Z); 
        HH=base_to_extract(H); FF=base_to_extract(F);
        DD=base_to_extract(D);
        tit1='Calculated components';
        tit2='Components from data files';
end

if code_fig==1,
    figure(1)
    subplot(311),plot(T,X-XX), title(tit1), datetick('x')
    ylabel(['X' sprintf('%+d',XX) ', nT']), axis tight
    subplot(312),plot(T,Y-YY), datetick('x')
    ylabel(['Y' sprintf('%+d',YY) ', nT']), axis tight
    subplot(313),plot(T,Z-ZZ), datetick('x'), axis tight
    ylabel(['Z' sprintf('%+d',ZZ) ', nT']), xlabel('Time UT, days')
    set(gcf,'Name','X,Y,Z-components','NumberTitle','off')

    figure(2)
    set(gcf,'Name','H,D,Z-components','NumberTitle','off')
	subplot(311),plot(T,H-HH), title(tit2), datetick('x')
	ylabel(['H' sprintf('%+d',HH) ', nT']), axis tight
	subplot(312),plot(T,(D-DD)*60), datetick('x')
	ylabel(['D' sprintf('%+d',DD) '^o, arcmin']), axis tight
    subplot(313),plot(T,Z-ZZ), datetick('x'), axis tight
    ylabel(['Z' sprintf('%+d',ZZ) ', nT'])
	xlabel('Time UT, days')
	
	figure(3)
	set(gcf,'Name','F(Var-Scal)','NumberTitle','off')
	subplot(211), plot(T,Fv-FF,T,F-FF,'r')
    datetick('x'), axis tight
    ylabel(['F-' sprintf('%d',FF) ', nT'])
    legend('F(Var)','F(Scal)','Location','best')
	subplot(212), plot(T,Fv-F)
    datetick('x'), axis tight
	ylabel('F(Var)-F(Scal), nT'), xlabel('Time UT, days')
end
