function result = decc_analyse(ncfile, varargin)

OPT = struct(...
    'radiusFcn', @(ncfile) str2double(regexp(ncfile, '(?<=radius)\d{5}', 'match')),...
    'philimits', [-90 90],...
    'phi_resolution', 1,...
    'tdim', false,...
    't0', false,...
    'zb', false,...
    's', false,...
    'Ero', true,...
    'Sed', true,...
    'Bal', true,...
    'Svsg_cs', true,...
    'Susg_cs', true,...
    'dSvsg_dn', true,...
    'ccg_cs', true,...
    'ceqsg_cs', false,...
    'Fy_mean', false,...
    'zs_mean', true,...
    'u_mean', false,...
    'v_mean', false,...
    'ee_mean', false,...
    'contour', false);

OPT = setproperty(OPT, varargin);

%%
result.ncfile = ncfile;
result.radius = OPT.radiusFcn(ncfile);

%% obtain revision number
revtxt = nc_attget(ncfile, nc_global, 'Build-Revision');
result.revision = str2double(regexp(revtxt, '\d{4}', 'match'));

%%
if OPT.tdim
    tdim = OPT.tdim;
else
    info = nc_info(ncfile);
    tdim = info.Dimension(strcmp({info.Dimension.Name}, 'globaltime')).Length;
end

x = nc_varget(ncfile, 'globalx')';
y = nc_varget(ncfile, 'globaly')';
g = xb_stagger(x, y);

phi = -rad2deg(mean(g.alfaz,1));
yid = phi > min(OPT.philimits) & phi < max(OPT.philimits);
if isnan(OPT.phi_resolution)
    result.phi = phi(yid);
else
    result.phi = min(OPT.philimits):OPT.phi_resolution:max(OPT.philimits);
end
ystart = find(diff(yid) == 1);
ycount = sum(yid);

%% find out when wave energy arrives at waterline near lateral boundary
if OPT.t0
    t = nc_varget(ncfile, 'globaltime');
    H = nc_varget(ncfile, 'H', [0 0 109], [-1 1 1]);
    t0id = find(H/max(H) > .2, 1, 'first');
    result.t0 = roundoff(t(t0id), -2);
end

%% derive cross-shore coordinate
dsz = g.dsz(:,yid);
dsz(1,:) = 0;
s = cumsum(dsz, 1);
if OPT.s
    result.s = s;
end

%%
zb0 = nc_varget(ncfile, 'zb', [0 ystart 0], [1 ycount -1])';
zb1 = nc_varget(ncfile, 'zb', [tdim-1 ystart 0], [1 ycount -1])';
    
if OPT.zb
    result.zb = cat(3, zb0, zb1);
end
    
if OPT.Sed || OPT.Ero || OPT.Bal
    correctionfac = g.dsdnz(:,yid) ./ g.dnz(:,yid);
    Ero = sum(((max(zb1, 5) - max(zb0, 5))) .* correctionfac);
    if OPT.Ero
        result.Ero = -Ero;
    end
    Sed = sum(((min(zb1, 5) - min(zb0, 5))) .* correctionfac);
    if OPT.Sed
        result.Sed = Sed;
    end
    if OPT.Bal
        result.Bal = Sed + Ero;
    end
end

% % % toe of erosion profile
if islogical(OPT.contour) && ~OPT.contour
    % cumvol = mean(cumtrapz(mean(s,2), abs(zb1-zb0), 1), 2);
    cumvol = cumtrapz(mean(s,2), abs(zb1-zb0), 1);
    maxcumvol = repmat(max(cumvol, [], 1), size(cumvol,1), 1);
    relcumvol = cumvol./maxcumvol;
    relcumvol = mean(relcumvol,2);
    imin = find(relcumvol > .025, 1, 'first');
    % [imin imax] = deal(NaN(size(cumvol(1,:))));
    % for i = 1:length(imin)
    %     imin(i) = find(cumvol(:,i)./ maxcumvol(:,i) > .025, 1, 'first');
    %     imax(i) = find(cumvol(:,i)./ maxcumvol(:,i) > .99, 1, 'first');
    % end
else
    id = find(mean(zb0,2) > OPT.contour, 1, 'first');
    id = id-1:id+1;
    imin = round(interp1(mean(zb0(id,:),2), id, OPT.contour));
end
result.contour = mean(zb0(imin,:));
imax = size(x,1);
% %
wmean_dnu = repmat(wmean(g.dnu(imin:imax,yid), g.dsu(imin:imax,yid)), imax-imin+1, 1);
% % Ero = sum(min(zb1(imin:imax,:)-zb0(imin:imax,:),0) .* g.dsu(imin:imax,:) .* g.dnu(imin:imax,:) ./ wmean_dnu);
% % Sed = sum(max(zb1(imin:imax,:)-zb0(imin:imax,:),0) .* g.dsu(imin:imax,:) .* g.dnu(imin:imax,:) ./ wmean_dnu);

% radius = roundoff(diff(minmax(g.yz(110,:)))/2, -2);
% dphi = diff(phi);
% dphi(end+1) = dphi(end);
% dphi_dn = 2*pi*radius/360;

% figure
% subplot(211)
% ph = plot(phi, -Ero, '-r',...
%     phi, Sed, '-b',...
%     phi, Sed+Ero, '-k');
%
% %%
% Svsg = nc_varget(ncfile, 'Svsg', [tdim-1 0 0 0], [1 1 -1 -1])';
if OPT.Svsg_cs
    try
        Svsg = nc_varget(ncfile, 'Svsg_mean');
    catch
        Svsg = NaN(size(x'));
    end
    if length(size(Svsg)) == 3
        Svsg = squeeze(mean(Svsg, 1))';
    else
        Svsg = Svsg';
    end
    result.Svsg_cs = sum(Svsg(imin:imax,yid) .* g.dsu(imin:imax,yid) .* g.dnu(imin:imax,yid) ./ wmean_dnu);
    
    if OPT.dSvsg_dn
        dSvsg = diff(Svsg,1,2);
        dSvsg(:,end+1) = dSvsg(:,end);
        result.dSvsg_dn = sum(dSvsg(imin:imax,yid) ./ g.dnu(imin:imax,yid) .* g.dsu(imin:imax,yid), 1);
    end
end

if OPT.Susg_cs
    try
        Susg = nc_varget(ncfile, 'Susg_mean');%, [0 0 ystart 0], [-1 -1 ycount -1]);
    catch
        Susg = NaN(size(x'));
    end
    idx = repmat({':'}, 1, ndims(Susg));
    idx{size(Susg) == length(yid)} = yid;
    Susg = Susg(idx{:});
    if length(size(Susg)) == 3
        Susg = squeeze(mean(Susg,1))';
    else
        Susg = Susg';
    end
%     result.ccg_cs = sum(ccg(imin:imax,:) .* g.dsu(imin:imax,yid) .* g.dnu(imin:imax,yid) ./ wmean_dnu, 1) ./ sum(g.dsu(imin:imax,yid));
    [tmp ix] = max(Susg,[],1);
    ix = min(ix):max(ix);
%     ix = imin:imax;
    result.Susg_cs = mean(Susg(ix,:)) * 2650;
end

% Svsg_dn = Svsg_cs ./ dphi / dphi_dn;

%
% subplot(212)
% plot(phi, Svsg_cs)
% ylim = get(gca, 'YLim');
%
% ampl = roundoff(min(abs(ylim)) / max(abs(minmax(dSvsg_dn))), -1, @floor);
%
% hold on
% plot(phi, dSvsg_dn*ampl, ':r')
%
% %%
% ve = nc_varget(ncfile, 've', [1 0 0], [1 -1 -1])';
% ue = nc_varget(ncfile, 'ue', [1 0 0], [1 -1 -1])';
% ve_max = max(ve,[],1);
% ue_max = max(ue, [], 1);
% figure; plot(phi, ve_max, phi, ue_max)
%
% %%
% ccg = nc_varget(ncfile, 'ccg', [tdim-1 0 0 0], [1 1 -1 -1])';
if OPT.ccg_cs
    try
        ccg = nc_varget(ncfile, 'ccg_mean');%, [0 0 ystart 0], [-1 -1 ycount -1]);
    catch
        ccg = NaN(size(x'));
    end
    idx = repmat({':'}, 1, ndims(ccg));
    idx{size(ccg) == length(yid)} = yid;
    ccg = ccg(idx{:});
    if length(size(ccg)) == 3
        ccg = squeeze(mean(ccg,1))';
    else
        ccg = ccg';
    end
%     result.ccg_cs = sum(ccg(imin:imax,:) .* g.dsu(imin:imax,yid) .* g.dnu(imin:imax,yid) ./ wmean_dnu, 1) ./ sum(g.dsu(imin:imax,yid));
    [tmp ix] = max(ccg,[],1);
    ix = min(ix):max(ix);
%     ix = imin:imax;
    result.ccg_cs = mean(ccg(ix,:)) * 2650;
end

if OPT.ceqsg_cs
    ceqsg = nc_varget(ncfile, 'ceqsg_mean');%, [0 0 ystart 0], [-1 -1 ycount -1]);
    idx = repmat({':'}, 1, ndims(ceqsg));
    idx{size(ceqsg) == length(yid)} = yid;
    ceqsg = ceqsg(idx{:});
    if length(size(ceqsg)) == 3
        ceqsg = squeeze(mean(ceqsg,1))';
    else
        ceqsg = ceqsg';
    end
%     result.ccg_cs = sum(ccg(imin:imax,:) .* g.dsu(imin:imax,yid) .* g.dnu(imin:imax,yid) ./ wmean_dnu, 1) ./ sum(g.dsu(imin:imax,yid));
    [tmp ix] = max(ccg,[],1);
    ix = min(ix):max(ix);
%     ix = imin:imax;
    result.ceqsg_cs = mean(ceqsg(ix,:)) * 2650;
end

if OPT.zs_mean
    try
        zs_mean = nc_varget(ncfile, 'zs_mean');%, [0 ystart 0], [-1 ycount 1]);
    catch
        zs_mean = NaN(size(x'));
    end
    idx = repmat({':'}, 1, ndims(zs_mean));
    idx{size(zs_mean) == length(yid)} = yid;
    idx{size(zs_mean) == size(x,1)} = 1;
    zs_mean_offshore = zs_mean(idx{:});
    if ~isvector(zs_mean_offshore)
        result.zs_mean_offshore = mean(zs_mean_offshore, 1);
    else
        result.zs_mean_offshore = zs_mean_offshore';
    end
    xstart = 108;%104;
    idx{size(zs_mean) == size(x,1)} = xstart+1;
    zs_mean_nearshore = zs_mean(idx{:});
    if ~isvector(zs_mean_nearshore)
        result.zs_mean_nearshore = mean(zs_mean_nearshore, 1);
    else
        result.zs_mean_nearshore = zs_mean_nearshore';
    end
    result.zs_mean_contour_nearshore = mean(zb0(xstart+1,:));
    dzs_mean = diff(result.zs_mean_nearshore);
    dzs_mean(:,end+1) = dzs_mean(:,end);
    result.dzs_mean_nearshore_dnz = dzs_mean ./ g.dnz(xstart+1,yid);
end

if OPT.Fy_mean
    xstart = 104;
    Fy_mean_nearshore = nc_varget(ncfile, 'Fy_mean', [0 ystart xstart], [-1 ycount 1]);
    if length(size(Fy_mean_nearshore)) == 3
        result.Fy_mean_nearshore = mean(Fy_mean_nearshore, 1);
    else
        result.Fy_mean_nearshore = Fy_mean_nearshore';
    end
    result.Fy_mean_contour_nearshore = mean(zb0(xstart+1,:));
end

if OPT.u_mean
    u_mean = nc_varget(ncfile, 'u_mean');%, [0 ystart 0], [-1 ycount 1]);
    idx = repmat({':'}, 1, ndims(u_mean));
    idx{size(u_mean) == length(yid)} = yid;
    idx{size(u_mean) == size(x,1)} = 1;
    u_mean = u_mean(idx{:});
    if ~isvector(u_mean)
        result.u_mean = mean(u_mean, 1);
    else
        result.u_mean = u_mean';
    end
end

if OPT.v_mean
    v_mean = nc_varget(ncfile, 'v_mean');%, [0 ystart 0], [-1 ycount 1]);
    idx = repmat({':'}, 1, ndims(v_mean));
    idx{size(v_mean) == length(yid)} = yid;
    idx{size(v_mean) == size(x,1)} = 1;
    v_mean = v_mean(idx{:});
    if ~isvector(v_mean)
        result.v_mean = mean(v_mean, 1);
    else
        result.v_mean = v_mean';
    end
end

if OPT.ee_mean
    ee_mean = nc_varget(ncfile, 'ee_mean');%, [0 0 ystart 0], [-1 -1 ycount -1]);
    idx = repmat({':'}, 1, ndims(ee_mean));
    idx{size(ee_mean) == length(yid)} = yid;
    ee_mean = ee_mean(idx{:});
    if length(size(ee_mean)) == 4
        result.ee_mean = squeeze(sum(mean(ee_mean, 1), 4));
    else
        result.ee_mean = squeeze(sum(ee_mean, 3));
    end
    result.ee_mean = mean(result.ee_mean, 2);
    try
        dtheta = nc_attget(ncfile, 'parameter', 'dtheta');
        thetamin = nc_attget(ncfile, 'parameter', 'thetamin');
        thetamax = nc_attget(ncfile, 'parameter', 'thetamax');
        result.theta = ((thetamin+.5*dtheta):dtheta:(thetamax-.5*dtheta))';
    catch
        result.theta = (190:10:350)';
    end
%     result.theta = repmat(result.theta, 1, size(ee_mean,3));
end

result = structfun(@(x) interpphi(phi(yid), result.phi, x), result,...
    'UniformOutput', false);

function new_var =interpphi(phi_orig, phi_new, interpvar)
if ischar(interpvar) || isequal(interpvar, phi_new) || isscalar(interpvar);
    new_var = interpvar;
    return
end
if isvector(interpvar)
    new_var = interp1(phi_orig, interpvar, phi_new);
else
    S = size(interpvar);
    if length(S) == 3
        for i = 1:S(3)
            new_var(:,:,i) = interp2(phi_orig, (1:S(1))', interpvar(:,:,i), phi_new, (1:S(1))');
        end
    else
        new_var = interp2(phi_orig, (1:S(1))', interpvar, phi_new, (1:S(1))');
    end
end