function [ C ] = Noise_Covariance_Matrix( I,J,K,mode_1,mode_2,sigma_2 )

if mode_1==mode_2
    C=eye(I*J*K);
else
    C=zeros(I*J*K);
    for i=1:I
        for j=1:J
            for k=1:K
                
                if mode_1==1
                    index_1=j+(k-1)*J+(i-1)*J*K;
                elseif mode_1==2
                    index_1=j+(k-1)*J+(i-1)*J*K;
                elseif mode_1==3
                    index_1=i+(k-1)*I+(j-1)*I*K;
                elseif mode_1==4
                    index_1=i+(j-1)*I+(k-1)*I*J;
                elseif mode_1==5
                    index_1=i+(j-1)*I+(k-1)*I*J;
                end
                
                if mode_2==1
                    index_2=j+(k-1)*J+(i-1)*J*K;
                elseif mode_2==2
                    index_2=j+(k-1)*J+(i-1)*J*K;
                elseif mode_2==3
                    index_2=i+(k-1)*I+(j-1)*I*K;
                elseif mode_2==4
                    index_2=i+(j-1)*I+(k-1)*I*J;
                elseif mode_2==5
                    index_2=i+(j-1)*I+(k-1)*I*J;                    
                end
                    
                C(index_1,index_2)=sigma_2;
                
            end
        end
    end
end
end

