%
%  chaos.m -- poincare map
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

global par skip

par = [0.005 0.02];  skip = 1;
%par = [0.005 0.0];  skip = 1;

%  plot saved data

figure(3);  hold on;  axis([-pi/2 pi/2 0.2 1.2]);  axis square

if(par(2)==0.02)
  load dataC
elseif(par(2)==0)
  load dataR
end
if(par(2)==0.02|par(2)==0.0)
%  plot(save1(:,1),save1(:,2),'b.')
  plot(save2(:,1),save2(:,2),'r.')
  plot(save3(:,1),save3(:,2),'r.')
  plot(save4(:,1),save4(:,2),'r.')
end

title('\bf poincare section')
xlabel('\bf \theta_- = (\theta_1-\theta_2)');  
ylabel('\bf I_- = (I_1-I_2)/\surd 2')

%  two ODE choices & IVs with 
%  final time, method (alpha string of .m function file) & options

aa = 0.02053/2;
%aa=0.1;
%aa = 0;

Xode = 'Xchaos';  IVs = [0.0 2.0+aa 0.0*pi 1.0-aa];  T = 100;  
method = 'ode113';
options = odeset('reltol',1e-9,'abstol',1e-9, ...
		 'stats','off','events','on');

ham = 2*IVs(2) + IVs(4) + (1/3)*IVs(4)^3            ... 
    + par(1)*4*IVs(2)*IVs(4)*(cos(IVs(1)-IVs(3)))^2 ...
    + par(2)*IVs(2)*sin(IVs(1)) ;
amp = (IVs(2)+IVs(4))/2;

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

save = [];

%  ODE solve

for  j = 1:40;  for jj = 1:20
  [t y] = feval(method,Xode,[0 T],IVs,options);
  IVs = y(end,:);

  tdif = mod((y(end,1)-y(end,3))+pi/2,pi)-pi/2;  
  idif = (y(end,2)-y(end,4))/sqrt(2);
  plot(tdif,idif,'b.');  save = [save ; tdif idif]; 

  con1 = 2*y(end,2) + y(end,4) + (1/3)*y(end,4)^3              ... 
       + par(1)*4*y(end,2)*y(end,4)*(cos(y(end,1)-y(end,3)))^2 ...
       + par(2)*y(end,2)*sin(y(end,1)) - ham;
  con2 = (y(end,2)+y(end,4))/2 - amp;
  disp([num2str((j-1)*20+jj) ' : ' num2str(con1) ' : ' ...
	num2str(con2) ' : ' num2str(cos(y(end,1)/skip))])
end;  pause(0);  end