function out = aggregateToModelDt(data, spdIn, spdOut, method, area_km2)
% AGGREGATETOMODELDT Da griglia file (es. 10 min) a griglia modello (es. 1 h).
%
% Tabella meteo (method 'meteo' o omesso):
%   PRCP [mm/10 min] sommati -> mm/ora
%   Tmax, Tmin, WSavg, RADavg = media nel blocco
%   RHmin = min nel blocco | RHmax = max nel blocco
%   date = inizio ora di calendario (se colonna date presente)
%
% Vettore:
%   'sum'       — mm nel blocco (pioggia -> mm/ora)
%   'discharge' — media m3/s nel blocco + conversione mm/h (richiede area_km2)

    if nargin < 4 || isempty(method)
        if istable(data)
            method = 'meteo';
        else
            error('aggregateToModelDt: per un vettore specificare method (sum|discharge).');
        end
    end
    method = lower(string(method));

    ratio = round(spdIn / spdOut);

    if istable(data)
        if ratio < 2
            out = data;
            return;
        end
        out = aggregateMeteoTable(data, ratio, spdIn, spdOut);
        return;
    end

    switch method
        case "discharge"
            if nargin < 5 || isempty(area_km2)
                error('aggregateToModelDt: method discharge richiede area_km2.');
            end
            if ratio >= 2
                q_m3s = aggregateVector(data(:), ratio, "mean");
                % fprintf('Portata: 10 min -> orario (media m3/s) -> mm/h\n');
            else
                q_m3s = data(:);
            end
            out = q_m3s * 3600 * 1000 / (area_km2 * 1e6);
        case "sum"
            if ratio < 2
                out = data(:);
            else
                out = aggregateVector(data(:), ratio, "sum");
            end
        otherwise
            error('aggregateToModelDt: method "%s" non supportato.', method);
    end
end

function y_out = aggregateVector(y, ratio, method)
    nUse = floor(numel(y) / ratio) * ratio;
    if nUse < ratio
        error('aggregateToModelDt: serie troppo corta (%d valori).', numel(y));
    end
    y = y(1:nUse);
    nOut = nUse / ratio;
    block = reshape(y, ratio, nOut);

    switch method
        case "sum"
            y_out = sum(block, 1, 'omitnan')';
        case "mean"
            y_out = mean(block, 1, 'omitnan')';
        otherwise
            error('aggregateToModelDt: metodo interno "%s" non valido.', method);
    end
end

function tbl = aggregateMeteoTable(tbl, ratio, spdIn, spdOut)
    cols = tbl.Properties.VariableNames;
    if ismember('date', cols)
        tbl = aggregateMeteoByCalendarHour(tbl, spdIn, spdOut);
        return;
    end

    nRows = height(tbl);
    nUse = floor(nRows / ratio) * ratio;
    if nUse < ratio
        error('aggregateToModelDt: meteo insufficiente (%d righe).', nRows);
    end
    if nUse < nRows
        fprintf('aggregateToModelDt: ignorate %d righe meteo incomplete.\n', nRows - nUse);
    end
    tbl = tbl(1:nUse, :);
    nOut = nUse / ratio;
    cols = tbl.Properties.VariableNames;
    out = table();

    if ismember('PRCP', cols)
        block = reshape(tbl.PRCP, ratio, nOut);
        out.PRCP = sum(block, 1, 'omitnan')';
    end

    meanVars = ["Tmax", "Tmin", "WSavg", "RADavg"];
    for v = meanVars
        if ismember(v, cols)
            block = reshape(tbl.(v), ratio, nOut);
            out.(v) = mean(block, 1, 'omitnan')';
        end
    end

    if ismember('RHmin', cols)
        block = reshape(tbl.RHmin, ratio, nOut);
        out.RHmin = min(block, [], 1, 'omitnan')';
    end
    if ismember('RHmax', cols)
        block = reshape(tbl.RHmax, ratio, nOut);
        out.RHmax = max(block, [], 1, 'omitnan')';
    end
    if ismember('date', cols)
        out.date = tbl.date(1:ratio:nUse);
    end

    tbl = out;
    % fprintf('Meteo %d -> %d passi/gg (%d righe), PRCP in mm/h\n', spdIn, spdOut, nOut);
end

function out = aggregateMeteoByCalendarHour(tbl, spdIn, spdOut)
    dt = datetime(tbl.date, 'InputFormat', 'preserve');
    hourId = dateshift(dt, 'start', 'hour');
    [hours, ~, gid] = unique(hourId, 'stable');

    cols = tbl.Properties.VariableNames;
    out = table();
    out.date = hours;

    if ismember('PRCP', cols)
        out.PRCP = splitapply(@(x) sum(x, 'omitnan'), tbl.PRCP, gid);
    end

    meanVars = ["Tmax", "Tmin", "WSavg", "RADavg"];
    for v = meanVars
        if ismember(v, cols)
            out.(v) = splitapply(@(x) mean(x, 'omitnan'), tbl.(v), gid);
        end
    end

    if ismember('RHmin', cols)
        out.RHmin = splitapply(@(x) min(x, [], 'omitnan'), tbl.RHmin, gid);
    end
    if ismember('RHmax', cols)
        out.RHmax = splitapply(@(x) max(x, [], 'omitnan'), tbl.RHmax, gid);
    end

    % fprintf('Meteo %d -> %d passi/gg (%d righe, ora di calendario), PRCP in mm/h\n', ...
    %     spdIn, spdOut, height(out));
end
