function tt_ = fetchEurostat_v2(dataset, options)
%FETCHEUROSTAT  Fetch Eurostat JSON-stat data → timetable.
%
%   tt_ = fetchEurostat(dataset, Name, Value)
%
%   Required:
%     dataset   char    - Eurostat dataset ID, e.g. 'namq_10_gdp'
%
%   Name-Value:
%     'params'  cell    - Flat key-value pairs for query string.
%                         Repeated keys allowed for multi-value dims,
%                         e.g. {'nace_r2','TOTAL','nace_r2','A','geo','EA21'}
%     'varName' char    - Column name for single-series output (default: dataset).
%                         For multi-dim output, used as prefix: varName_TOTAL etc.
%     'start'   datetime - Trim observations before this  (default: no trim)
%     'lastN'   double  - lastTimePeriod query value       (default: 1000)

    arguments
        dataset             (1,:) char
        options.params      (1,:) cell     = {}
        options.varName     (1,:) char     = dataset
        options.start       (1,1) datetime = datetime(1,1,1)
        options.lastN       (1,1) double   = 1000
    end

    % --- Build URL -------------------------------------------------------
    baseURL = 'https://ec.europa.eu/eurostat/api/dissemination/statistics/1.0/data/';
    p = options.params;
    assert(mod(numel(p), 2) == 0, 'fetchEurostat: params must be key-value pairs');
    qparts = cell(1, numel(p)/2);
    for k = 1:numel(p)/2
        qparts{k} = [p{2*k-1} '=' p{2*k}];
    end
    url = sprintf('%s%s?lang=en&lastTimePeriod=%d&%s', ...
        baseURL, dataset, options.lastN, strjoin(qparts, '&'));

    % --- Fetch -----------------------------------------------------------
    jsonData = webread(url);

    % --- Dimension metadata ----------------------------------------------
    dimNames = jsonData.id(:)';   % force 1xN
    dimSizes = jsonData.size(:)'; % force 1xN

    timePos  = find(strcmp(dimNames, 'time'));
    nTime    = dimSizes(timePos);
    stride   = prod(dimSizes(timePos+1:end));   % prod([]) = 1

    % --- Time labels → datetime ------------------------------------------
    timeLabels = struct2cell(jsonData.dimension.time.category.label);
    time_      = parseEurostatTime(timeLabels(:));

    % --- Find non-trivial dimensions (size > 1, excluding time) ----------
    multiDims = find(dimSizes > 1 & (1:numel(dimSizes)) ~= timePos);

    % --- Extract values --------------------------------------------------
    vRaw = jsonData.value;

    if numel(multiDims) == 0
        % Single series
        data_ = extractValues(vRaw, nTime, stride, 1, 1);
        colNames = {options.varName};
        dataMatrix = data_;

    elseif isscalar(multiDims)
        % One extra dimension — return one column per slice
        dimPos    = multiDims(1);
        dimName   = dimNames{dimPos};
        dimSize   = dimSizes(dimPos);

        % Get labels for this dimension
        dimLabels = struct2cell( ...
            jsonData.dimension.(dimName).category.label);
        dimLabels = dimLabels(:);

        % Stride for the varying dimension
        % dimStride = prod(dimSizes(dimPos+1:end));

        dataMatrix = nan(nTime, dimSize);
        for s = 1:dimSize
            % Base offset for this slice in the flat index
            % sliceOffset = (s-1) * dimStride * prod(dimSizes(timePos+1:dimPos-1+1));
            % Actually compute per-element using full index decomposition
            dataMatrix(:, s) = extractSlice(vRaw, dimSizes, timePos, dimPos, s-1);
        end
        colNames = strcat(options.varName, '_', dimLabels');

    else
        error('fetchEurostat:multiDim', ...
            'More than one non-time dimension with size>1. Fetch slices separately.');
    end

    % --- Trim & build timetable ------------------------------------------
    mask  = time_ >= options.start;
    tt_   = array2timetable(dataMatrix(mask,:), ...
        'RowTimes', time_(mask), ...
        'VariableNames', colNames);
end


function data_ = extractValues(vRaw, nTime, stride, ~, ~)
%EXTRACTVALUES  Pull a single time series from vRaw.
    data_ = nan(nTime, 1);
    if isnumeric(vRaw)
        v = vRaw(:);
        data_(1:min(numel(v), nTime)) = v(1:min(numel(v), nTime));
    elseif isstruct(vRaw)
        vFields = fieldnames(vRaw);
        for i = 1:numel(vFields)
            flatIdx = str2double(regexprep(vFields{i}, '^x', ''));
            tIdx    = floor(mod(flatIdx, nTime * stride) / stride) + 1;
            if tIdx >= 1 && tIdx <= nTime
                data_(tIdx) = vRaw.(vFields{i});
            end
        end
    end
end


function col = extractSlice(vRaw, dimSizes, timePos, sliceDimPos, sliceIdx)
%EXTRACTSLICE  Extract one slice along sliceDimPos from the full hypercube.
%   sliceIdx is 0-based.
    nTime     = dimSizes(timePos);
    nDims     = numel(dimSizes);
    % totalSize = prod(dimSizes);

    % Precompute strides for each dimension (row-major / C order)
    strides = ones(1, nDims);
    for d = nDims-1:-1:1
        strides(d) = strides(d+1) * dimSizes(d+1);
    end

    col = nan(nTime, 1);

    % Build a flat index for each time step with the slice fixed
    for t = 0:nTime-1
        % Construct multi-index: all singleton dims = 0, time = t, slice = sliceIdx
        multiIdx         = zeros(1, nDims);
        multiIdx(timePos)      = t;
        multiIdx(sliceDimPos)  = sliceIdx;
        flatIdx = sum(multiIdx .* strides);   % 0-based

        val = getValueAt(vRaw, flatIdx);
        if ~isnan(val)
            col(t+1) = val;
        end
    end
end


function val = getValueAt(vRaw, flatIdx)
%GETVALUEAT  Retrieve value at 0-based flatIdx from struct or numeric vRaw.
    if isnumeric(vRaw)
        idx = flatIdx + 1;
        if idx >= 1 && idx <= numel(vRaw)
            val = vRaw(idx);
        else
            val = NaN;
        end
    elseif isstruct(vRaw)
        % Try 'x'-prefixed field name first, then bare number
        fname = sprintf('x%d', flatIdx);
        if isfield(vRaw, fname)
            val = vRaw.(fname);
        else
            fname2 = num2str(flatIdx);
            if isfield(vRaw, fname2)
                val = vRaw.(fname2);
            else
                val = NaN;   % sparse — missing entry means NaN
            end
        end
    else
        val = NaN;
    end
end


function time_ = parseEurostatTime(labels)
%PARSEEUROSTATTIME  Eurostat time strings → datetime.
    ex = labels{1};
    if contains(ex, 'Q')
        time_ = datetime(labels, 'InputFormat', 'yyyy-QQQ');
    elseif numel(ex) == 7 && ex(5) == '-'
        time_ = datetime(labels, 'InputFormat', 'yyyy-MM');
    elseif numel(ex) == 4
        time_ = datetime(labels, 'InputFormat', 'yyyy');
    else
        error('fetchEurostat:parseTime', 'Unrecognised time label format: %s', ex);
    end
    time_ = time_(:);
end