function Phi = GaborSteinerETF( M )
% GaborSteinerETF

% Input:
% M, vector of length 1 or 2 of odd integers greater than or equal to 3

% Output:
% Phi, ETF of m^2 vectors in  C^(m(m-1)/2), where m is the product of the
% components of M

% Authors: Emily J. King and Bernhard Bodmann
% Paper: Optimal arrangements of classical and quantum states with limited purity
% https://arxiv.org/abs/1811.11513
% Last Updated: September 13, 2018


% Defining parameters
numFact=length(M);
m=prod(M);
d=m*(m-1)/2;

% Testing input
if ~rem(m,2) || ~prod(M>ones(size(M)))
    error('The input must be a vector of odd integers greater than or equal to 3.')
elseif numFact > 2
    error('The code can currently only handle input of length 2 or less.')
end

% Making the fiducial
E=eye(m);
psi=E(:,1)-E(:,end);
for i=2:(m-1)/2
     psi=vertcat(psi,E(:,i)-E(:,end-i+1));
end

% Building the Gabor-Steiner ETF
Phi=zeros(d,m^2);

if numFact==1
    D=fft(E);
    T=E; T=T([end 1:(end-1)],:);
    
    for j=0:(m-1) % translation
        for k=0:(m-1) % modulation
          Phi(:,m*j+k+1)=kron(eye((m-1)/2),diag(D(k+1,:))*T^j)*psi;
        end
    end
else %numFact==2 by earlier check
    E1=eye(M(1));
    E2=eye(M(2));
    D1=fft(E1);
    D2=fft(E2);
    D=kron(D1,D2);
    T1=E1; T1=T1([end 1:(end-1)],:);
    T2=E2; T2=T2([end 1:(end-1)],:);
    
    for i=0:(M(1)-1) %first component translation
        for j=0:(M(2)-1) %second component translation
            for k=0:(M(1)-1) % first component modululation
                for l=0:(M(2)-1) % second component modululation
                    Phi(:,(M(2)*i+j)*m+(M(2)*k+l+1))=... 
                        kron(eye((m-1)/2), kron(T1^i,T2^j)*diag(D(M(2)*k+l+1,:)))*psi; 
                end
            end
        end
    end
end

end