function [b, p] = field_trend_sen(data, p_level)
% Calculates Sen's slope and Mann-Kendall test for trend analysis on a field
%
% Input:
%   data: Input field with dimensions:
%         - 3D array (lat × lon × time) or
%         - 4D array (lat × lon × month × time)
%   p_level: Significance level for Mann-Kendall test
%
% Output:
%   b: Sen's slope estimator (trend magnitude)
%   p: Mann-Kendall test significance (binary: 1 if significant, 0 if not)
%
% Dependencies:
%   - trend_sen: Function to calculate Sen's slope
%   - mann_kendall: Function to perform Mann-Kendall test

if ndims(data) == 3
    % Process 3D data (lat × lon × time)

    % Create mask for valid (non-NaN) data points
    mask = ~isnan(data(:,:,1));
    % Convert 2D mask to 1D
    mask_1d = reshape(mask, size(mask,1)*size(mask,2), 1);

    % Reshape 3D data to 2D matrix (space × time)
    data_1d = reshape(data, size(data,1)*size(data,2), size(data,3));

    % Initialize output arrays
    b = nan*zeros(size(mask_1d));
    p = nan*zeros(size(mask_1d));

    % Calculate Sen's slope for each valid spatial point
    b(mask_1d) = trend_mk_sen(data_1d(mask_1d,:));
    % Perform Mann-Kendall test along time dimension (dim=2)
    p(mask_1d) = double(mann_kendall(data_1d(mask_1d,:),p_level,'dim',2));

    % Reshape results back to spatial grid
    b = reshape(b, size(data,1), size(data,2));
    p = reshape(p, size(data,1), size(data,2));

elseif ndims(data) == 4
    % Process 4D data (lat × lon × month × time)

    % Create mask for valid (non-NaN) data points
    mask = ~isnan(data(:,:,:,1));
    % Convert 3D mask to 1D
    mask_1d = reshape(mask, size(mask,1)*size(mask,2)*size(mask,3), 1);

    % Reshape 4D data to 2D matrix (space×month × time)
    data_1d = reshape(data, size(data,1)*size(data,2)*size(data,3), size(data,4));

    % Initialize output arrays
    b = nan*zeros(size(mask_1d));
    p = nan*zeros(size(mask_1d));

    % Calculate Sen's slope for each valid spatial/monthly point
    b(mask_1d) = trend_mk_sen(data_1d(mask_1d,:));
    % Perform Mann-Kendall test along time dimension (dim=2)
    p(mask_1d) = double(mann_kendall(data_1d(mask_1d,:),p_level,'dim',2));

    % Reshape results back to spatial/monthly grid
    b = reshape(b, size(data,1), size(data,2), size(data,3));
    p = reshape(p, size(data,1), size(data,2), size(data,3));

end

end