function res_list = ellipseParam(data, k, pcx, pcy, model)
    % model = 1 ---- PCA
    % model = 2 ---- PLSDA
    % model = 3 ---- OPLSDA
    % Check inputs
    if isempty(data)
        warndlg('Seems you forgot to provide data values.','Warning');
        res_list.exit = 1;
        return
    end
    
    if ~istable(data) && ~ismatrix(data)
        warndlg('Data must be a data frame, tbl_df, or a matrix.','Warning');
        res_list.exit = 1;
        return
    end
    
    if pcx == 0 || pcy == 0
        warndlg('No component is provided either in PCx or PCy, or both.','Warning');
        res_list.exit = 1;
        return
    end
    
    if pcx == pcy && model ~= 3
        warndlg('Please provide two different components in PCx and PCy.','Warning');
        res_list.exit = 1;
        return
    end
    
    if k < 2
        warndlg('k must be at least equal to 2.','Warning');
        res_list.exit = 1;
        return
    end
    
    if k > size(data, 2)
        warndlg('k exceeds the number of components in the data.','Warning');
        res_list.exit = 1;
        return
    end
    
    % Convert data to a matrix
    X = data;
    
    % Sample size
    n = size(X, 1);
    
    % Number of principal components
    A = k;
    
    % Squared Mahalanobis distance
%     MDsq = mahalanobis(X, mean(X), cov(X), 'covType', 'full');

    % Assuming X is your data matrix
    % Calculate the covariance matrix
    cov_matrix = cov(X);
    
    % Calculate the mean of each column
    mean_vector = mean(X);
    
    % Calculate the squared Mahalanobis distance
    diff = X - mean_vector;
    inv_cov_matrix = inv(cov_matrix);
    MDsq = sum((diff * inv_cov_matrix) .* diff, 2);
    
    if k == 2
        % 99% and 95% confidence limit
        Tsq_limit1 = (A*(n-1)/(n-A))*finv(0.99, A, n-A);
        Tsq_limit2 = (A*(n-1)/(n-A))*finv(0.95, A, n-A);
        
        % Hotelling ellipse semi-axes
        a_limit1 = sqrt(Tsq_limit1*var(X(:, pcx)));
        a_limit2 = sqrt(Tsq_limit2*var(X(:, pcx)));
        b_limit1 = sqrt(Tsq_limit1*var(X(:, pcy)));
        b_limit2 = sqrt(Tsq_limit2*var(X(:, pcy)));
        
        axis_param = table(...
            a_limit1, ...
            b_limit1, ...
            a_limit2, ...
            b_limit2, ...
            'VariableNames', {'a_99pct', 'b_99pct', 'a_95pct', 'b_95pct'} ...
        );
        
        % Hotelling's T-square
        Tsq = table(((n-2)/(2*(n-1)))*MDsq, 'VariableNames', {'value'});
        
        res_list = struct(...
            'Tsquare', Tsq, ...
            'Ellipse', axis_param, ...
            'cutoff_99pct', Tsq_limit1, ...
            'cutoff_95pct', Tsq_limit2, ...
             exit = 0);
        
        return;
    end

end
