% Berechnung der Poincareschnitte

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB -
% Beispiele zu Kapitel 1
% Vers. 3.0 2013
%
%% 3. Projekt: MATLAB Funktionen ....
% Eventfunktionen

% ausgewaehlte Startwerte
h1=subplot(2,1,1)  %  Plottet in das zweite Fenster
y0(2)=0;
epsilon=-0.5;
% Zufallswerte fuer die Impulse
for k=1:50
    k
zwi=2*rand;
z(randperm(2))=[zwi;sqrt(4-zwi^2)*rand];
farb=[0,0,0];%rand(1,3);
y0(3)=z(1);
y0(4)=z(2);
% Aus der Energieerhaltung folgt
y0(1)=sqrt((1/2*(y0(3)^2+y0(4)^2)-2)/epsilon);
%options = odeset('RelTol',1e-13,'AbsTol',1e-13,'MaxStep',0.1);
zwei=0.5*(y0(3)*y0(3)+y0(4)*y0(4)) -epsilon*(y0(1)*y0(1)+y0(2)*y0(2))+...
    1/8*(y0(1)*y0(1)*y0(2)*y0(2))*(y0(1)*y0(1)+y0(2)*y0(2))
options = odeset('RelTol',1e-10,'AbsTol',1e-10,'MaxStep',0.1,'Events',@(t,y) poinevent(t,y,epsilon));
tmax=100;
xa=sqrt(-2*epsilon)*y0(1);
plot(xa,y0(3),'o',...
    'MarkerSize',1,'MarkerFaceColor',farb,'MarkerEdgeColor',farb)
hold on
plot(-xa,y0(3),'o',...
    'MarkerSize',1,'MarkerFaceColor',farb,'MarkerEdgeColor',farb)
plot(xa,-y0(3),'o',...
    'MarkerSize',1,'MarkerFaceColor',farb,'MarkerEdgeColor',farb)
plot(-xa,-y0(3),'o',...
    'MarkerSize',1,'MarkerFaceColor',farb,'MarkerEdgeColor',farb)
[t,y,te,ye,ie] = ode45(@(t,y) diahdgl(t,y,epsilon),[0 tmax],y0,options);
% Visualisierung
xa=sqrt(-2*epsilon)*ye(:,1);
ya = ye(:,3);
xa = [xa;-xa;xa;-xa];
ya = [ya;ya;-ya;-ya];
    plot(xa,ya,'.',...
    'MarkerSize',1,'MarkerFaceColor',farb,'MarkerEdgeColor',farb)
%xlim([-2,2]),ylim([-2,2])
%axis equal
end
%hold off
