function [U_p, V_p, pc_U, pc_V, scf, s2_sum] = mca(data1, data2, n_mca, method, lat1, lon1, lat2, lon2)

% Maximum Covariance Analysis (MCA) between two spatiotemporal fields
%
% Input:
%   data1(lat,lon,time): First climate field (left field for SVD), detrended and demeaned
%   data2(lat,lon,time): Second climate field (right field for SVD), detrended and demeaned
%   n_mca: Number of MCA modes to retain
%   method: 'aw' for area-weighted analysis
%   lat1, lon1: Latitude and longitude arrays for data1
%   lat2, lon2: Latitude and longitude arrays for data2
%
% Output:
%   U_p: Spatial patterns for left field (lat × lon × n_mca)
%   V_p: Spatial patterns for right field (lat × lon × n_mca)
%   pc_U: Time series (expansion coefficients) for left field (time × n_mca)
%   pc_V: Time series (expansion coefficients) for right field (time × n_mca)
%   scf: Squared covariance fraction for each mode
%   s2_sum: Total squared covariance
%
% Example usage:
%   mca(data1, data2);
%   mca(data1, data2, 'aw', lat1, lon1, lat2, lon2);

% Apply area weighting if method is 'aw'
if nargin > 3
    if strcmp(method,'aw')         % area weighting
        % Weight data by square root of cosine of latitude
        [~,LAT1] = meshgrid(lon1, lat1);
        data1 = data1.*(cosd(LAT1)).^0.5;

        [~,LAT2] = meshgrid(lon2, lat2);
        data2 = data2.*(cosd(LAT2)).^0.5;
    end
end

% Reshape 3D data arrays to 2D matrices (space × time)
% Remove NaN points (typically land points in ocean data)
data1_1d = reshape(data1, size(data1,1)*size(data1,2), size(data1,3));
mask1 = ~isnan(data1(:,:,1));
mask1_1d = mask1(:);
data1_1d = data1_1d(mask1_1d,:);

data2_1d = reshape(data2, size(data2,1)*size(data2,2), size(data2,3));
mask2 = ~isnan(data2(:,:,1));
mask2_1d = mask2(:);
data2_1d = data2_1d(mask2_1d,:);

% Compute cross-covariance matrix and perform SVD
Cxy = data1_1d*data2_1d'/(size(data1,3)-1);  % Cross-covariance matrix
[U,Sigma,V] = svd(Cxy,0);  % Singular Value Decomposition
s = diag(Sigma);  % Extract singular values
s2 = s.^2;  % Square of singular values
s2_sum = sum(s2);  % Total squared covariance
scf = s2/sum(s2);  % Squared covariance fraction

% Initialize output arrays for patterns and time series
U_p = nan(size(data1,1),size(data1,2),n_mca);
pc_U = nan(size(data1,3),n_mca);

V_p = nan(size(data2,1),size(data2,2),n_mca);
pc_V = nan(size(data2,3),n_mca);

% Calculate spatial patterns and time series for each mode
for mode = 1:n_mca
    % Left field (data1)
    Uall = nan(size(data1,1)*size(data1,2),1);
    Uall(mask1) = U(:,mode); 
    U_p(:,:,mode) = reshape(Uall,size(data1,1),size(data1,2));  % Spatial pattern
    pc_U(:,mode) = (U(:,mode))'*data1_1d;  % Time series
    
    % Right field (data2)
    Vall = nan(size(data2,1)*size(data2,2),1);
    Vall(mask2) = V(:,mode)'; 
    V_p(:,:,mode) = reshape(Vall,size(data2,1),size(data2,2));  % Spatial pattern
    pc_V(:,mode) = (V(:,mode))'*data2_1d;  % Time series
end

% Keep only the requested number of modes for SCF
scf = scf(1:n_mca);

end