function AMocca = Moccamatrix(YP,M1,L)
% computes the Moccamatrix of size (Nc*M1^2) x (Nc*L^2) from the (M1+L-1) x (M1+L-1) centered 
% (completely acquired k-space data
%
% G. Plonka 2024

[N,M,Nc]=size(YP);

%extract the centered calibration data of size M1+L-1 x M1+L-1 for each coil
K1=floor(N/2)-floor((M1+L-1)/2); K2=floor(M/2)-floor((M1+L-1)/2);
Y1=YP(K1+1:K1+M1+L-1,K2+1:K2+M1+L-1,:);

%build partial MOCCA matrix from shifts of calibration data

l=floor((L-1)/2);
anzL=L^2;
A1=[];
for j=1:Nc
Yb=Y1(:,:,j);
    for l1=-l:l
        for l2=-l:l
            Ybc=circshift(Yb,[l2,l1]);
            Ybcpart=Ybc(l+1:M1+l,l+1:M1+l); %partial calibration matrix blocks
            Ybv=reshape(Ybcpart,M1^2,1);   %calibration matrix blocks of size M1^2 x L^2
            A1=[A1 Ybv];
        end
    end
end

%Build the Mocca matrix of size (M1^2 * Nc) x (L^2 * Nc)
AMocca=[];
ysum=zeros(M1^2,anzL);
for j=1:Nc
    ysum=ysum+A1(:,1+(j-1)*anzL:j*anzL);
end
for k=1:Nc
    Z=[];
    for j=1:Nc
        if j==k
            Z=[Z A1(:,1+(k-1)*anzL:k*anzL)-ysum];
        else
        Z=[Z A1(:,1+(k-1)*anzL:k*anzL)];
        end
    end
    AMocca = [AMocca;Z];
end
end