%
%  twirl.m -- energy plots
%
%- - - - - - - - - - - - - - - - - - - - - - - - 

%  wheel radius & omega (complex cube-root)

global r w; r = 0.98;  w = exp(i*2*pi/3);

%  angular coordinates & energies

dt = pi/40;  th1 = 0:dt:2*pi;  th2 = 0:dt:2*pi;
[t1,t2] = meshgrid(th1,th2);

vt = zeros(size(t1));  vt1 = zeros(size(t1));  vt2 = zeros(size(t1));

%  indices for rotors (m & n)

for m = 0:1:2
  a = r*(w^m)*exp(i*t1);  za = +1 + a;

  for n = 0:1:2
    b = r*(w^n)*exp(i*t2);  zb = -1 + b;

    v  = 1./abs(za-zb);

    vt  = vt  + v;
    vt1 = vt1 - (v.^3) .* real( i*(za-1).*conj(za-zb));
    vt2 = vt2 - (v.^3) .* real(-i*(zb+1).*conj(za-zb));
  end
end

%  make plots

clf;  figure(1);  colormap(autumn)
contour(th1,th2,vt,10);  colorbar
axis square;  title('\bf billiard trajectory')
xlabel('\bf \theta_{right}');  ylabel('\bf \theta_{left}')

%figure(2);  colormap(autumn)
%mesh(th1,th2,vt);  colorbar
%view(-45,72)
%xlabel('\bf \theta_{right}');  ylabel('\bf \theta_{left}')