function srec=multichannel_coils_recon(AMocca,L,N,M,Nc,R)

% input: AMocca Mocca matrix
%        L size of support of coil sensitivities in k-space, L odd
%        (N, M) size of the image m
%        Nc number of coils
%        brain (number of image data taken)
%        R reduction rate
% output: srec of size NxMxNc contains the Nc coils sensitivities

 
 [U,S,V]=svd(AMocca);
 l=floor((L-1)/2);
 anzL=L^2;
 K1=floor(N/2)-l; 
 K2=floor(M/2)-l;

 % illustration of singular values of Mocca matrix 
 sing=diag(S);
 [ns,ns2]=size(sing);
 eigThresh_1=0.01;
 idx = L^2*Nc -max(find(sing >= sing(1)*eigThresh_1)) %number of "small" singular values

 figure('name','Singular values of MOCCA matrix'),
 plot(1:ns,sing,'LineWidth',2)
 set(gca,'FontSize',20)
 colormap(gca,gray(256))
 %print -depsc sing1_L=9.eps;


 num_eig=idx; 
 if num_eig==0 num_eig=1; end
 %num_eig=num_eig
 W=zeros(L,L); W(l+1,l+1)=1;
 w=reshape(W,L^2,1); w=kron(ones(Nc,1),w);
 factors=V(:,end:(-1):end-num_eig+1)'*w;    

   
cvec=V(:,end:(-1):end-num_eig+1)*factors;
crec=zeros(L,L,Nc); srec=zeros(N,M,Nc); crec0=zeros(N,M,Nc); 
 zb=zeros(K1,L); zc=zeros(N+1,K2);
 for j=1:Nc
   cneu=cvec(1+(j-1)*anzL:j*anzL,1); %partion into coil sensitivity vectors
   crec(:,:,j)=reshape(cneu,L,L);
   crec0(K1:K1+L-1,K2:K2+L-1,j)=crec(:,:,j);  % zero padding
   srec(:,:,j)= (N*M)*ifftshift(ifft2(ifftshift(crec0(1:N,1:M,j))));
 end

%Compute d for sos condition
d=zeros(N,M);
for j=1:Nc
    d=d+conj(srec(:,:,j)).*srec(:,:,j);
end


% compute d^+ (componentwise inverse of d if possible)
dinv=ones(N,M)./d;
index=find(abs(dinv)> 10^8|dinv==Inf);
if numel(index)>0 
    dinv= reshape(dinv,N*M,1);
    dinv(index)=dinv(index)./abs(dinv(index));
    maxdinv=max(abs(dinv(:)));
    dinv(index)=dinv(index).*maxdinv;
    dinv=reshape(dinv,N,M);
end
 
% normalize sensitivities
  dinv=sqrt(dinv);
  srec(:,:,:)= dinv.*srec(:,:,:);

%illustrate coils sensitivities

figure('name','coil sensitivities') 
imshow3(abs(srec),[],[2,Nc/2]) 
axis off
set(gca,'FontSize',20)
colormap default
colorbar
set(gca,'FontSize',20);
%saveas(gca,'coilabs2_L=9_Ns76.eps','epsc')
%colormap default

figure('name','coil sensitivities (phase angle)'), 
imshow3(angle(srec),[],[2,Nc/2]) 
axis off
set(gca,'FontSize',20)
colorbar
phasemap
set(gca,'FontSize',20);
%saveas(gca,'coilphase2_L=9_Ns76.eps','epsc')


 

