%
%  tw.m -- ode model
%
%- - - - - - - - - - - - - - - - - - - - - - - - 

% coeffs

global r w d;  d = 0*0.0005;  clf;  figure(1);  clf;  twirl

%  ODE file, IVs, timing, method & accuracy

Xode = 'Xtwirl';  T = 16;  cycles = 8;  method = 'ode15s';  tol = 1e-2;

%  fancy IV input (discretized for repeatability)

disp(' ');  disp(' clik to choose initial point');  disp(' ');

in1 = [];  figure(1);  hold on;  IV = [0 0];
while (isempty(in1))
  in = floor(ginput(1)/(pi/25))*(pi/25); 

  if(isempty(in))
    break
  else
    plot(IV(1),IV(2),'w.')
    disp(['  IV point:  ' num2str(in) '  (rtn to confirm)']); 
    IV = in;  in = [];
    plot(IV(1),IV(2),'k.')
  end
end
plot(IV(1),IV(2),'w.');  plot(IV(1),IV(2),'kx');

in = [];  ID = [0 0];
while (isempty(in))
  in = floor(ginput(1)/(pi/25))*(pi/25); 

  if(isempty(in))
    break
  else
    plot(ID(1),ID(2),'w.')
    disp(['  IV veloc:  ' num2str(in) '  (rtn to confirm)']); 
    ID = in;  in = [];
    plot(ID(1),ID(2),'r.')
  end
end
plot(ID(1),ID(2),'w.');  plot(ID(1),ID(2),'rx');

IVs = [IV(1) ID(1)-IV(1) IV(2) ID(2)-IV(2)];

options = odeset('reltol',tol*1e-3,'abstol',tol*1e-6,'refine',1, ...
		'maxstep',T/20,'stats','off','events','off');

disp(' ');  disp([method ' for ' Xode]);  disp(' ')

%  initial energy

egy = 0.5*(IVs(2)^2 + IVs(4)^2);
for m = 0:1:2
  a = r*(w^m)*exp(i*IVs(1));  za = +1 + a;
  for n = 0:1:2
    b = r*(w^n)*exp(i*IVs(3));  zb = -1 + b;
    egy  = egy + 1./abs(za-zb);
  end
end

disp(['  ' num2str(0) '  ' num2str(egy)]);

%  ODE solve & plot

for C = 1:cycles
  [t y] = feval(method,Xode,[(C-1)*T C*T],IVs,options);
  ang1 = y(:,1);  ang2 = y(:,3);

  IVs = y(end,:);  

  egy = 0.5*(y(end,2)^2 + y(end,4)^2);
  for m = 0:1:2
    a = r*(w^m)*exp(i*y(end,1));  za = +1 + a;

    for n = 0:1:2
      b = r*(w^n)*exp(i*y(end,3));  zb = -1 + b;
      egy  = egy + 1./abs(za-zb);
    end
  end

  disp(['  ' num2str(C) '  ' num2str(egy)])

  figure(1);  clf;  
  contour(th1,th2,vt,10);  colorbar;    hold on
  axis square;  title('\bf billiard trajectory')
  xlabel('\theta_{right}');  ylabel('\theta_{left}')
  plot(mod(ang1(  1,1),2*pi),mod(ang2(  1,1),2*pi),'kx');
  plot(mod(ang1(end,1),2*pi),mod(ang2(end,1),2*pi),'rx');

  for k = 2:size(ang1,1)-1
    p1 = ang1(k-1:k+1)-[1;1;1]*(ang1(k)-mod(ang1(k),2*pi));
    p2 = ang2(k-1:k+1)-[1;1;1]*(ang2(k)-mod(ang2(k),2*pi));
    plot(p1,p2,'k')
  end

  figure(3);
  subplot(2,1,1);
  plot(t,y(:,2),'b');  hold on;  plot(t,y(:,4),'g');
  title('\bf rotation rates (left = blue, right = green)')
  xlabel('time');

  subplot(2,1,2);  hold off
  plot(t,y(:,1)-(y(1,1)-mod(y(1,1),2*pi)),'b');  hold on;
  plot(t,y(:,3)-(y(1,3)-mod(y(1,3),2*pi)),'g');
  title('\bf angle dynamics (left = blue, right = green)')
  xlabel('\bf time');

  pause(3)
end