function [val,W] = TJ(rho,m)
% [val,W] TJ(rho,m) computes the family of coherence measures as defined in Phys. Rev. Lett 121, 010401 (2018)
% rho: dxd density matrix or normalised state vector
% m: non-negative integer coefficient indicating the chosen measure
% the optimal witness W s.t. val=trace(rho*W) is returned only for non-pure rho

dim = size(rho);
d = max(dim);

if m>=d; m = d-1; end

if rank(rho) == 1 % pure state case
    
    if min(dim) ~= 1
        [~,~,v] = svd(rho);
        rho = v(:,1);
    end
    
    m = m+1; % to compute the (m+1)-distillation norm
    psi = sort(abs(rho),'descend');
    
    k = 1;
    minsum = sum(psi(m:d).^2);

    for j=2:m
        jsum = sum(psi(m-j+1:d).^2)/j;
        if jsum < minsum
            minsum = jsum;
            k = j;
        end
    end
    
    val = sum(psi(1:m-k))+k*sqrt(minsum);
    val = val^2-1;
    W = 0;
    
else % general case
    
    rho = (rho+rho')/2; % to avoid numerical problems
    
    cvx_begin quiet sdp
    
        variable W(d,d) hermitian
        -eye(d) <= W <= m*eye(d)
        diag(W) == 0
    
        maximise trace(rho*W);
    
    cvx_end
    
    val = cvx_optval;
    
end

end