function data = IAGA20022Bin (path_in, st_code, path_out)
%    version 30/05/2017

%    data - data array, (days number) x 4

%    path_in  - location of .min files  
%    st_code  - station code (ex., 'aae')
%    path_out - location of output binary files 
   
    namef = [path_in st_code '*.min'];
    files = dir(namef);
    
    if isempty(files)
        return;
    end
    
    data = zeros (1000000, 4);
	miss_val = 99999.00;
    num = 0;
    
    prog_bar = Mbar (1, ['Reading files for station ' st_code]);    
    prog_bar.Bname ('Please wait...', length (files));
    
    tic
    for file_ind = 1: length(files)
        
        f_name = fullfile(path_in, files(file_ind).name);

        fid = fopen (f_name, 'rt');
    
        fseek (fid, 0, 'eof');
        f_size = ftell (fid);

        fseek (fid, 0, 'bof');
    
        str = fgetl(fid);

        while isempty(strfind(str, 'Reported     '))
            str = fgetl(fid);
        end         
        aaa=strtrim(str); aaa=strtrim(aaa(9:length(aaa)));
        orient = aaa(1:4);
        
        while isempty(strfind(str, 'DATE       TIME         DOY'))
            str = fgetl(fid);
        end 

        str = fgetl(fid);
        
        if file_ind == 1
            dt1 = datenum (str(1:16), 'yyyy-mm-dd HH:MM');
        end
        
        num = num + 1;

        data (num, 1) = sscanf (str(28:40), '%f', 1);
        data (num, 2) = sscanf (str(41:50), '%f', 1);
        data (num, 3) = sscanf (str(51:60), '%f', 1);
        data (num, 4) = sscanf (str(61:70), '%f', 1);
        
        while (~feof(fid))
            str = fgetl(fid);
            num = num + 1;
            data (num, 1) = sscanf (str(28:40), '%f', 1);
            data (num, 2) = sscanf (str(41:50), '%f', 1);
            data (num, 3) = sscanf (str(51:60), '%f', 1);
            data (num, 4) = sscanf (str(61:70), '%f', 1);
        end

        if file_ind == length(files)
            dt2 = datenum (str(1:16), 'yyyy-mm-dd HH:MM');
        end

        fclose (fid);        
        
%        data = cat(1, data, dt);

        cur_pos = file_ind;
        prog_bar.Bval (cur_pos);
               
    end
    
    data = data (1:num, 1:4); 
    
    for ch_ind = 1: 4
        channel(ch_ind) = orient (ch_ind:ch_ind);
        
        %--- h->x and d->y conversion ------------
           if (strcmpi(channel(ch_ind), 'H')) % h->x
               data(find(data(:,2)==miss_val),1) = miss_val;
               data(find(data(:,1)~=miss_val),1) = data(find(data(:,1)~=miss_val),1).*cosd(data(find(data(:,1)~=miss_val),2)/60);
               channel(ch_ind) = 'X';
           end
           if (strcmpi(channel(ch_ind), 'D')) % d->y
               data(find(data(:,1)==miss_val),2) = miss_val;
               if (strcmpi(channel(1), 'H'))  % if h->x is not done
                   data(find(data(:,2)~=miss_val),2) = data(find(data(:,2)~=miss_val),1).*sind(data(find(data(:,2)~=miss_val),2)/60);
               else % if h->x is done
                   data(find(data(:,2)~=miss_val),2) = data(find(data(:,2)~=miss_val),1).*tand(data(find(data(:,2)~=miss_val),2)/60);
               end
               channel(ch_ind) = 'Y';
           end
        %-----------------------------------------
        
        %--- x->h conversion ---------------------
%            if (strcmpi(channel(ch_ind), 'X')) % x->h
%                data(find(data(:,2)==miss_val),1) = miss_val;
%                data(find(data(:,1)~=miss_val),1) = sqrt(data(find(data(:,1)~=miss_val),1).^2+data(find(data(:,1)~=miss_val),2).^2);
%                channel(ch_ind) = 'H';
%                disp(st_code);
%            end        
        %-----------------------------------------
        
        %--- linear interpolation of 'miss_val' values --- 
        [keg, beg] = get_potential(find(data (:,ch_ind) == miss_val));
        data(:, ch_ind) = linearInterpol (data(:, ch_ind), keg, beg);
        %-------------------------------------------------         

        %--- conversion of data into binary format --- 
        name_out = [st_code '_' channel(ch_ind) '_' datestr(dt1, 'yymmdd') '_' datestr(dt2, 'yymmdd') '.bin'];
        fid2 = fopen (fullfile (path_out, name_out), 'wb');
        fwrite(fid2, dt1, 'double');
        fwrite(fid2, dt2, 'double');
        fwrite(fid2, data(:, ch_ind), 'float32');
        fclose(fid2);
        %---------------------------------------------
    end
    
	prog_bar.Bclose;
    toc

end

function [keg, beg] = get_potential(indg)

    dl = length(indg);
    if dl<1
        keg=0;
        beg=0;
        return;
    else
        if dl<2
            beg(1,1) = indg;
            beg(1,2) = indg;
        else
            kk = find(indg(2:dl)-indg(1:dl-1)-1);
            yy1 = cat(1, indg(1), indg(kk+1));
            yy2 = cat(1, indg(kk), indg(dl));
            beg = cat(2, yy1, yy2);
        end
        [keg rw] = size(beg);
    end

end

function data_out = linearInterpol (data, keg, beg)
    data_out = data;
    for i = 1: keg
        i1 = beg(i, 1);
        i2 = beg(i, 2);
        if i1 == 1
            if (i2 ~= length(data_out))
                data_out(i1:i2) = data_out(i2+1);
            end
        elseif i2 == length(data_out)
            data_out(i1:i2) = data_out(i1-1);
        else
            i1 = i1-1;
            i2 = i2+1;
            data_out(i1:i2) = linspace(data_out(i1), data_out(i2), i2-i1+1);
        end
    end
end