function [area_sum_mat] = ...
    map_binary_wet(folder_path, polygon_250m, part_num)
% Function takes a folder path, loads all map files, determines how many
% time steps of map files there are, analyzes each of the *map.nc files for
% inundation at each element. Returns a matrix with inundated area totals
% for each 250 meter segment

fnames = struct2cell(dir(strcat(folder_path,'*map.nc')))'; 
[fnum,~] = size(fnames); % number of files

if (fnum / part_num) > 1
    
    file_datetime = cell(fnum,1);
    
    for j = 1:fnum
        [~,name,~] = fileparts(fnames{j,1});
        file_datetime(j,1) = {name(end-18:end-4)};
    end
    
    time_stamps = unique(file_datetime(:,1));
    [run_num,~] = size(time_stamps);
    idxywA_cell = cell(run_num,1);
    
    for l = 1:run_num
        % Retrieve files only having correct time stamp
        fnames_time = struct2cell(dir(strcat(folder_path,'*',...
            time_stamps{l,1},'_map.nc')))'; 
        [fnum_time,~] = size(fnames_time); % number of files
        
        for m = 1:fnum_time
            
            fname_m = strcat(folder_path,fnames_time{m,1});
        
            if m == 1
                id = double(ncread(fname_m,'FlowElemGlobalNr'));
                x = ncread(fname_m,'FlowElem_xcc');
                y = ncread(fname_m,'FlowElem_ycc');
                d = ncread(fname_m,'waterdepth');
                wet_yn = any(d > 0,2);
                A = ncread(fname_m,'FlowElem_bac');
            else
                id = ([id; double(ncread(fname_m,'FlowElemGlobalNr'))]);
                x = ([x; ncread(fname_m,'FlowElem_xcc')]);
                y = ([y; ncread(fname_m,'FlowElem_ycc')]);
                d = ncread(fname_m,'waterdepth');
                wet_yn = ([wet_yn; any(d > 0, 2)]);
                A = ([A; ncread(fname_m,'FlowElem_bac')]);
            end

        end
        
        % Combine all extracted data
        idxywA_raw = ([id x y double(wet_yn) A]);
        [~,unique_r,~] = unique(idxywA_raw(:,2:3),'rows','stable');
        idxywA_cell(l,1) = {idxywA_raw(unique_r,:)};
        
        clear id x y wet_yn A idxywA_raw idxywA_unique
        
    end
    
    idxywA = cell2mat(idxywA_cell(1,1));
    idxywA_one = idxywA;
    for n = 2:run_num
        idxywA_hold = cell2mat(idxywA_cell(n,1));
        idxywA(:,4) = max(idxywA(:,4),idxywA_hold(:,4));
        change_wet = sum(idxywA(:,4)~=idxywA_hold(:,4));
    end
    
else
    
    for k = 1:fnum
        
        fname_k = strcat(folder_path,fnames{k,1});
        
        if k == 1
            id = double(ncread(fname_k,'FlowElemGlobalNr'));
            x = ncread(fname_k,'FlowElem_xcc');
            y = ncread(fname_k,'FlowElem_ycc');
            d = ncread(fname_k,'waterdepth');
            wet_yn = any(d > 0,2);
            A = ncread(fname_k,'FlowElem_bac');
        else
            id = ([id; double(ncread(fname_k,'FlowElemGlobalNr'))]);
            x = ([x; ncread(fname_k,'FlowElem_xcc')]);
            y = ([y; ncread(fname_k,'FlowElem_ycc')]);
            d = ncread(fname_k,'waterdepth');
            wet_yn = ([wet_yn; any(d > 0, 2)]);
            A = ([A; ncread(fname_k,'FlowElem_bac')]);
        end
        
    end
    
    % Combine all extracted data
    idxywA_raw = ([id x y double(wet_yn) A]);
    [~,unique_r,~] = unique(idxywA_raw(:,2:3),'rows','stable');
    idxywA = idxywA_raw(unique_r,:);
    
end

% Determine maximum cross-sectional area
xs_poly_num = max(polygon_250m(:,3));

% Create empty matrix for area information
area_sum_mat = zeros(xs_poly_num,3);

% For each cross-sectional area check what mesh points are within the
% bounds and calculate areas
for xs_num = 1:xs_poly_num

    in_poly_yn = inpolygon(idxywA(:,2),idxywA(:,3),...
        polygon_250m(polygon_250m(:,3)==xs_num,1),...
        polygon_250m(polygon_250m(:,3)==xs_num,2));

    area_sum_mat(xs_num,1) = sum(idxywA(in_poly_yn==1,5));
    area_sum_mat(xs_num,2) = ...
        sum(idxywA(in_poly_yn==1&idxywA(:,4)==1,5));
    area_sum_mat(xs_num,3) = ...
        sum(idxywA(in_poly_yn==1&idxywA(:,4)==0,5));

end

end