  function T = MakeTransition(theta)
% theta  is a column n-vector that satisfies
%    0<=theta(1)<theta(2)<...<theta(n) < 2*pi.
% T is the nxn 
n = length(theta);
c = cos([theta;theta(1)]);
s = sin([theta;theta(1)]);
% Vector of reciprocal distances...
r = 1./ ((c(1:n)-c(2:n+1)).^2 + (s(1:n)-s(2:n+1)).^2);
% Set up T (scalar-level implementation--no penalty)...
T = zeros(n,n);
for k=1:n
    if k==1
       T(1,1) = 2    /(2+r(1)+r(n));
       T(2,1) = r(1) /(2+r(1)+r(n));
       T(n,1) = r(n) /(2+r(1)+r(n));
    elseif k==n
       T(n,n)   = 2      /(2+r(n-1)+r(n));
       T(1,n)   = r(n)   /(2+r(n-1)+r(n));
       T(n-1,n) = r(n-1) /(2+r(n-1)+r(n));
    else
       T(k,k)   = 2       /(2+r(k-1)+r(k));
       T(k+1,k) = r(k)    /(2+r(k-1)+r(k));
       T(k-1,k) = r(k-1)  /(2+r(k-1)+r(k));
    end
end

% Vectorized set-up....
% The diagonal of D:
d = 2*ones(n,1)+r + [r(n);r(1:n-1)];
% The diagonal of T
diagT = 2./d;
% The subdiagonal of T...
subdiagT = r(1:n-1)./d(1:n-1);
% The superdiagonal of T...
superdiagT = r(1:n-1)./d(2:n);
% Use the diag function to set up the tridiagonal part of T...
T =  diag(subdiagT,-1) + diag(diagT) + diag(superdiagT,1) + diag(r(n),n-1);
% The corner entries..
T(1,n) = r(n)/d(n);
T(n,1) = r(n)/d(1);

% -1 for not exploiting fact that D is diagonal
% -1 for incorrect r

