function result = dewa_analyse(ncfiles, varargin)

% info = nc_info(ncfile);
% 
% tdim = info.Dimension(strcmp({info.Dimension.Name}, 'globaltime')).Length;
% tmdim = 1;
% if any(strcmp({info.Dimension.Name}, 'meantime'))
%     tmdim = info.Dimension(strcmp({info.Dimension.Name}, 'meantime')).Length;
% end
% 
OPT = struct(...
    ...    'tdim', tdim,...
    ...   'tmdim', tmdim,...
    'varnames', {{'zb' 'globalx' 'globaly' 'globaltime' 'Svsg_mean' 'Susg_mean' 'ccg_mean' 've_mean' 'ue_mean' 'zs_mean'}});

OPT = setproperty(OPT, varargin);
% 
% tdim = OPT.tdim;
% tmdim = OPT.tmdim;

%%
% [x y t] = cellfun(@(ncfile) deal(...
%     nc_varget(ncfile, 'globalx'),...
%     nc_varget(ncfile, 'globaly'),...
%     nc_varget(ncfile, 'globaltime')), ncfiles,...
%     'uniformoutput', false); 
% 
% g = cellfun(@xb_stagger, x, y);
% [ny nx] = cellfun(@size, x);
% nt = cellfun(@length, t);
% 


result = cellfun(@(ncfile) readvar(ncfile, OPT.varnames), ncfiles);
% 
% result = mergestructs(result,...
%     struct('ny', num2cell(ny),...
%     'nx', num2cell(nx),...
%     'nt', num2cell(nt)));

% if isfield(result, 'zb')
%     zb0 = cellfun(@(z,nt) repmat(z(1,:,:), [nt 1 1]), {result.zb}, num2cell(nt), 'uniformoutput', false);
%     % A: area per running meter coastline
%     A = cellfun(@(x,y,nt) repmat(shiftdim(x'./y', -1), [nt 1 1]), {g.dsdnz}, {g.dnz}, num2cell(nt), 'uniformoutput', false);
%     Ero = cellfun(@(zb0,zb,A) sum((max(zb, 5) - max(zb0, 5)) .* A,3), zb0, {result.zb}, A, 'uniformoutput', false);
%     Sed = cellfun(@(zb0,zb,A) sum((min(zb, 5) - min(zb0, 5)) .* A,3), zb0, {result.zb}, A, 'uniformoutput', false);
%     result = mergestructs(result,...
%         cell2struct(Ero, 'Ero')',...
%         cell2struct(Sed, 'Sed')');
% end
% 
% dsz = g.dsz;
% dsz(1,:) = 0;
% s = cumsum(dsz, 1);
% 
% zb0 = nc_varget(ncfile, 'zb', [0 0 0], [1 -1 -1])';
% zb1 = nc_varget(ncfile, 'zb', [tdim-1 0 0], [1 -1 -1])';
% 
% Ero = sum((max(zb1, 5) - max(zb0, 5)) .* g.dsdnz ./ g.dnz);
% Sed = sum((min(zb1, 5) - min(zb0, 5)) .* g.dsdnz ./ g.dnz);
% 
% %%
% idymid = floor(size(x,2)/2);
% 
% x = x(:,idymid);
% % y = y(idymid,:);
% 
% zb0 = nc_varget(ncfile, 'zb', [0 0 0], [1 -1 -1]);
% zb0 = select_transect(zb0, idymid, ny);
% zb1 = nc_varget(ncfile, 'zb', [tdim-1 0 0], [1 -1 -1]);
% zb1 = select_transect(zb1, idymid, ny);
% 
% ccg = nc_varget(ncfile, 'ccg_mean');%, [0 0 idymid 0], [tmdim -1 1 -1])';
% ccg = select_transect(ccg, idymid, ny);
% 
% Svsg = nc_varget(ncfile, 'Svsg_mean');%, [0 0 idymid 0], [tmdim -1 1 -1])';
% Svsg = select_transect(Svsg, idymid, ny);
% 
% dSvsg = diff(Svsg([1:end end],:),1,1);
% Svsg_grad = dSvsg ./ g.dnu' .* g.dsu';
% 
% ve = nc_varget(ncfile, 've_mean');%, [0 idymid 0], [tmdim 1 -1])';
% ve = select_transect(ve, idymid, ny);
% 
% zs = nc_varget(ncfile, 'zs_mean');%, [0 idymid 0], [tmdim 1 -1])';
% zs = select_transect(zs, idymid, ny);

function result = readvar(ncfile, varnames)
for varname = varnames
    info = ncinfo(ncfile, varname{1});
    [start count] = deal(zeros(size((info.Dimensions))),...
        ones(size((info.Dimensions)))*-1);
    ydim = strcmp({info.Dimensions.Name}, 'globaly');
    start(ydim) = round(info.Size(ydim) * .3);
    count(ydim) = round(info.Size(ydim) * .4);
    if any(strcmp(varname{1}, {'globalx' 'globaly' 've_mean' 'ue_mean' 'zs_mean'}))
        [start count ydim] = deal(fliplr(start), fliplr(count), fliplr(ydim));
    end
    result.(varname{1}) = nc_varget(ncfile, varname{1},...
        start, count);
    if any(ydim)
        ydim = find(size(result.(varname{1})) == count(ydim));
        result.(varname{1}) = squeeze(mean(result.(varname{1}), ydim));
    end
end
