function l4_phys_mass = READICA_L4_PHYS_MASS(varargin)
%READICA_L4_PHYS_MASS Read RPCICA L4 CORRECTED data for one day.
%   l4_corr = READICA_L4_PHYS_MASS() Reads data for a day. 
%        Skips data with less than 96 e_levels.
%        Asks for rootpath & date
%
%   l4_corr = READICA_L4_PHYS_MASS(stretch) Reads data for a day. 
%        Stretches data with less than 96 e_levels to 96 levels.
%        Asks for rootpath & date
%
%   l4_corr = READICA_L4_PHYS_MASS(stretch, date) Reads data for date given. 
%        If stretch is true, stretches data with less than 96 e_levels to 96 levels.
%        Asks for rootpath
%
%   l4_corr = READICA_L4_PHYS_MASS(stretch, date, rootpath) Reads data for date given.
%        If stretch is true, stretches data with less than 96 e_levels to 96 levels.
%
stretch = false;
data_date = '';
rootpath = '';
if nargin > 2
    if ~ischar(varargin{3})
        help(mfilename)
        return
    end
    rootpath = varargin{3};
end
if nargin > 1
    if ~ischar(varargin{2})
        help(mfilename)
        return
    end
    data_date = varargin{2};
end
if nargin > 0
    if ~islogical(varargin{1})
        help(mfilename)
        return
    end
    stretch = varargin{1};
end

%% Find rootpath of archive where directory DATA is located
if isempty(rootpath) || ~exist(fullfile(rootpath, 'DATA'),'dir')
    for jj=1:5
        if exist(fullfile(rootpath, 'DATA'),'dir')
            break
        else
            rootpath = fullfile('..', rootpath);
        end
    end
end

%% Get path to data
icapath=fullfile(rootpath, 'DATA', 'PHYS_MASS');
while ~exist(icapath, 'dir')
    fprintf(['This script is designed to run inside an archive.\n' ...
        'The DATA/PHYS_MASS directory was not found.\n' ...
        'Please provide the path to the root of the archive i.e.:\n' ...
        '/home/user/psa/RO-C-RPCICA-4-ESC1-PHYS_MASS-V1.0\n' ...
        '\n']);
    rootpath = input('path to the root of the archive:  ','s');
    icapath=fullfile(rootpath, 'DATA', 'PHYS_MASS');
end

fprintf('rootpath: %s\n', rootpath);

if isempty(data_date)
    data_date=input('Date of data to import yyyy-mm-dd:  ','s');
end
dn = datenum(data_date);


cdir=fullfile(icapath, ...
    datestr(dn,'yyyy'), ...
    upper(datestr(dn,'mmm')), ...
    ['D' datestr(dn,'dd')]);
fprintf('cdir: %s\n', cdir);

%% Get files
fcdir=fullfile(cdir, '/RPCICA*.TAB');
filelist=dir(fcdir);

%% initialize matrices
l4_phys_mass.time = [];
l4_phys_mass.delta_t = [];
l4_phys_mass.quality = [];
l4_phys_mass.mode = [];
l4_phys_mass.noise_reduction = [];
l4_phys_mass.pacclevel = [];
l4_phys_mass.azim = [];
l4_phys_mass.elev = [];
l4_phys_mass.swversion = [];
l4_phys_mass.H = [];
l4_phys_mass.He = [];
l4_phys_mass.He2 = [];
l4_phys_mass.Heavy = [];

%% Find function to call
if (exist('IMPORT_L4_CORR') == 2)
    import4_corr = str2func('IMPORT_L4_CORR');
else
    import4_corr = str2func('import_l4_corr');
end

if (exist('IMPORT_L4_PHYS_MASS') == 2)
    import = str2func('IMPORT_L4_PHYS_MASS');
else
    import = str2func('import_l4_phys_mass');
end

%% Read data
e_level_skip = 0;

if ~isempty(filelist)
    for jj=1:length(filelist)
        disp(filelist(jj).name);
        if ~isempty(strfind(filelist(jj).name,'_L4_H.TAB'))
            [table, e_levels] = import(fullfile(cdir, filelist(jj).name));
            lbl = strrep(filelist(jj).name, 'TAB', 'LBL');
            l4_lbl = READ_PDS_LBL(fullfile(cdir, lbl));
            curr_sw = str2double(l4_lbl.ROSETTA_COLON_ICA_SW_VERSION);
            he = strrep(filelist(jj).name, '_L4_H.TAB', '_L4_HE.TAB');
            [hetable, hee_levels] = import(fullfile(cdir, he));
            he2 = strrep(filelist(jj).name, '_L4_H.TAB', '_L4_HE2.TAB');
            [he2table, he2e_levels] = import(fullfile(cdir, he2));
            hvy = strrep(filelist(jj).name, '_L4_H.TAB', '_L4_HVY.TAB');
            [hvytable, hvye_levels] = import(fullfile(cdir, hvy));
        else
            % No support for other files
            continue;
        end
        if e_levels ~= 96 
            if stretch
                % possibly stretch 8 & 32 e_levels to 96 e_levels
                table.flux = [table.flux nan(size(table.time,1),96-e_levels)];
                hetable.flux = [hetable.flux nan(size(table.time,1),96-e_levels)];
                he2table.flux = [he2table.flux nan(size(table.time,1),96-e_levels)];
                hvytable.flux = [hvytable.flux nan(size(table.time,1),96-e_levels)];
            else
                e_level_skip = e_level_skip + 1;
                continue;
            end
        end
        
        if any(table.time ~= hetable.time)
            disp('time is not equal between H & He')
        end
        if any(table.time ~= he2table.time)
            disp('time is not equal between H & He2')
        end
        if any(table.time ~= hvytable.time)
            disp('time is not equal between H & Hvy')
        end
        
        l4_phys_mass.time=[l4_phys_mass.time; table.time];
        l4_phys_mass.delta_t=[l4_phys_mass.delta_t; table.delta_t];
        l4_phys_mass.quality=[l4_phys_mass.quality; table.quality];
        l4_phys_mass.mode= [l4_phys_mass.mode; table.mode];
        l4_phys_mass.noise_reduction=[l4_phys_mass.noise_reduction; table.noise_reduction];
        l4_phys_mass.pacclevel=[l4_phys_mass.pacclevel; table.pacclevel];
        l4_phys_mass.azim=[l4_phys_mass.azim; table.azim];
        l4_phys_mass.elev=[l4_phys_mass.elev; table.elev];
        l4_phys_mass.swversion = [l4_phys_mass.swversion; repmat(curr_sw, length(table.time), 1)];
        l4_phys_mass.H = [l4_phys_mass.H; table.flux]; 
        l4_phys_mass.He = [l4_phys_mass.He; hetable.flux]; 
        l4_phys_mass.He2 = [l4_phys_mass.He2; he2table.flux]; 
        l4_phys_mass.Heavy = [l4_phys_mass.Heavy; hvytable.flux]; 
    end
    save(['ica' data_date '_phys_mass.mat'], 'l4_phys_mass');
    if e_level_skip > 0
        fprintf('Skipped reading %d files due to low number of e_levels\n', e_level_skip);
    end
else
    disp(['Can not find data file: ' fcdir]); 
end
