function [b, p, r] = field_regress(data, index)
% Performs linear regression of a spatial field against a time series
%
% Input:
%   data: Input field with dimensions:
%         - 3D array (lat × lon × time) or
%         - 4D array (lat × lon × month × time)
%   index: Time series to regress against (predictor variable)
%
% Output:
%   b: Regression coefficients (slope)
%   p: P-values from correlation analysis with autocorrelation consideration
%   r: Residuals from regression
%
% Dependencies:
%   - corr_auto: Function to compute correlation with autocorrelation adjustment

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

    % Reshape 3D data to 2D matrix (space × time) and remove NaN points
    data_1d = reshape(data, size(data,1)*size(data,2), size(data,3));
    mask = ~isnan(data(:,:,1));
    mask_1d = mask(:);
    data_1d = data_1d(mask_1d,:);

    if nargout > 1
        % If p-values are requested, compute full statistics
        b_1d = nan*zeros(length(data_1d),1);
        p_1d = nan*zeros(length(data_1d),1);
        r_1d = nan*zeros(size(data_1d));
        
        % Loop through each spatial point
        for i=1:size(data_1d,1)
            % Perform regression with constant term
            [bx,~,rx] = regress(data_1d(i,:)',[index; ones(1,length(index))]');
            b_1d(i) = bx(1);  % Store regression coefficient
            r_1d(i,:) = rx;   % Store residuals
            % Compute correlation p-value with autocorrelation consideration
            [~, p_1d(i)] = corr_auto(data_1d(i,:)',index);
        end

    else
        % If only regression coefficients are needed
        b_1d = nan*zeros(length(data_1d),1);
        r_1d = nan*zeros(size(data_1d));
        for i=1:size(data_1d,1)
            [bx,~,rx] = regress(data_1d(i,:)',[index; ones(1,length(index))]');
            b_1d(i) = bx(1);
            r_1d(i,:) = rx;
        end    
    end    

    % Reshape results back to spatial grid
    ball = nan(size(data,1)*size(data,2),1);
    ball(mask) = b_1d; 
    b = reshape(ball,size(data,1),size(data,2));
    
    % Reshape residuals
    rall = nan(size(data,1)*size(data,2),size(data,3));
    mask_r = repmat(mask,1,1,size(data,3));
    rall(mask_r) = r_1d; 
    r = reshape(rall,size(data,1),size(data,2),size(data,3));

    if nargout > 1
        % Reshape p-values if computed
        pall = nan(size(data,1)*size(data,2),1);
        pall(mask) = p_1d; 
        p = reshape(pall,size(data,1),size(data,2));
    end

elseif ndims(data) == 4
    % Process 4D data (lat × lon × month × time)
    % Similar process as above but with additional dimension

    % 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));
    mask = ~isnan(data(:,:,:,1));
    mask_1d = mask(:);
    data_1d = data_1d(mask_1d,:);

    if nargout > 1
        % Full statistics computation
        b_1d = nan*zeros(length(data_1d),1);
        p_1d = nan*zeros(length(data_1d),1);
        r_1d = nan*zeros(size(data_1d));
        for i=1:size(data_1d,1)
            [bx,~,rx] = regress(data_1d(i,:)',[index; ones(1,length(index))]');
            b_1d(i) = bx(1);
            r_1d(i,:) = rx;
            [~, p_1d(i)] = corr_auto(data_1d(i,:)',index);
        end

    else
        % Regression coefficients only
        b_1d = nan*zeros(length(data_1d),1);
        r_1d = nan*zeros(size(data_1d));
        for i=1:size(data_1d,1)
            [bx,~,rx] = regress(data_1d(i,:)',[index; ones(1,length(index))]');
            b_1d(i) = bx(1);
            r_1d(i,:) = rx;
        end    
    end    

    % Reshape results back to spatial/monthly grid
    ball = nan(size(data,1)*size(data,2)*size(data,3),1);
    ball(mask) = b_1d; 
    b = reshape(ball,size(data,1),size(data,2),size(data,3));

    rall = nan(size(data,1)*size(data,2)*size(data,3),size(data,4));
    mask_r = repmat(mask,1,1,1,size(data,4));
    rall(mask_r) = r_1d; 
    r = reshape(rall,size(data,1),size(data,2),size(data,3),size(data,4));
    
    if nargout > 1
        pall = nan(size(data,1)*size(data,2)*size(data,3),1);
        pall(mask) = p_1d; 
        p = reshape(pall,size(data,1),size(data,2),size(data,3));
    end
end

end