% Poincareschnitt zu nichtverschwindendem Drehimpuls

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB -
% Vers. 1.0 2016


%% ausgewaehlte Startwerte
clear, close all
tic
y0(1)=0;
y0(3)=0;
y0(6)=0;
y0(8)=0;
y0(2)=0.01915;%rand;
y0(4)=y0(2); 

epsilon=-0.5;  % sinnvolle Werte sind negativ.
farb=[0,0,0];%rand(1,3);
neu = 1;
figure, hold on, shg

tphi1 = y0(1)/y0(2);
sphi1 = tphi1/sqrt(1+tphi1^2);
tphi2 = y0(3)/y0(4);
sphi2 = tphi2/sqrt(1+tphi2^2);
mu = y0(1)/sphi1;
nu = y0(3)/sphi2;

xa=sqrt(-2*epsilon)*(mu+nu)/2;
    pmu =  y0(5)*sphi1 + y0(6)*sqrt(1-sphi1^2);
    pnu =  y0(7)*sphi2 + y0(8)*sqrt(1-sphi2^2);
    pxa = (pmu + pnu)/2;
    %zwei=0.5*(y(3)*y(3)+y(4)*y(4)) -epsilon*(y(1)*y(1)+y(2)*y(2))+...
    %1/8*(y(1)*y(1)*y(2)*y(2))*(y(1)*y(1)+y(2)*y(2))
    plot(xa,pxa,'.',...
    'MarkerSize',1,'MarkerFaceColor',farb,'MarkerEdgeColor',farb)
    plot(-xa,pxa,'.',...
    'MarkerSize',1,'MarkerFaceColor',farb,'MarkerEdgeColor',farb)
    plot(xa,-pxa,'.',...
    'MarkerSize',1,'MarkerFaceColor',farb,'MarkerEdgeColor',farb)
    plot(-xa,-pxa,'.',...
    'MarkerSize',1,'MarkerFaceColor',farb,'MarkerEdgeColor',farb)

% Aus der Energieerhaltung folgt
zwi = 4 + 2*epsilon*(y0(2)^2+y0(4)^2)-1/4*y0(2)^2*y0(4)^2*(y0(2)^2 + y0(4)^2);
y0(5) = sqrt(zwi/2);
y0(7) = y0(5);
Lz=y0(2)*y0(5)-y0(1)*y0(6)
%%options = odeset('RelTol',1e-11,'AbsTol',1e-11,'MaxStep',0.1,'events',@(t,y) diaphievent(t,y,epsilon,farb));
zwei = 0.5.*(y0(5).^2 + y0(7).^2) - epsilon*(y0(2).^2 + y0(4).^2) + 1/8*y0(2)^2*y0(4)^2 * (y0(2).^2 + y0(4).^2)
options = odeset('RelTol',1e-11,'AbsTol',1e-11,'MaxStep',0.1,'events',@(t,y) diaphievent(t,y,epsilon));
%options = odeset('events',@(t,y) diaphievent(t,y,epsilon,farb));
tmax=200;
[t,y,te,ye,ie] = ode45(@(t,y) diahphidgl(t,y,epsilon),[0 tmax],y0,options);
toc
ye=y; % per Startwert auf Poincare-Ebene
% Visualisierung
tphi1 = ye(:,1)./ye(:,2);
sphi1 = tphi1./sqrt(1+tphi1.^2);
tphi2 = ye(:,3)./ye(:,4);
sphi2 = tphi2./sqrt(1+tphi2.^2);
mu = ye(:,1)./sphi1;
nu = ye(:,3)./sphi2;
 xa=sqrt(-2*epsilon)*(mu+nu)/2;
 pmu =  ye(:,5).*sphi1 + ye(:,6).*sqrt(1-sphi1.^2);
 pnu =  ye(:,7).*sphi2 + ye(:,8).*sqrt(1-sphi2.^2);
 pxa = (pmu + pnu)/2;
  plot(xa,pxa,'.',...
     'MarkerSize',1,'MarkerFaceColor',farb,'MarkerEdgeColor',farb)
tic
Lz = y(:,2).*y(:,5) - y(:,1).*y(:,6);
Lz2 = y(:,4).*y(:,7) - y(:,3).*y(:,8);
N0 = Lz - Lz2;
pv = y(:,5:8);
pvq = pv.*pv;
Tkin = 0.5*sum(pvq,2);
zzwi = y(:,1:4);
zzwiq = zzwi.^2;
zzwiqsum = sum(zzwiq,2);
zwei = Tkin - epsilon*zzwiqsum + 1/8*(y(:,1).^2 + y(:,2).^2) .* (y(:,3).^2 + y(:,4).^2) .* zzwiqsum;
toc
% Ergebnis testen:
max(N0)
Lz(1)
aLz = 2*(min(Lz)-max(Lz))/(min(Lz)+max(Lz))
azwei = max(abs(zwei-2))/2 
%% ab hier mehrfach ausfuehren -> Schleife oder Run and Advance
farb = [0,0,0];%rand(1,3);
y0 = diaphianfang(y0,epsilon)    % Berechnet neuen Startwert

tphi1 = y0(1)/y0(2);
sphi1 = tphi1/sqrt(1+tphi1^2);
tphi2 = y0(3)/y0(4);
sphi2 = tphi2/sqrt(1+tphi2^2);
mu = y0(1)/sphi1;
nu = y0(3)/sphi2;

    xa=sqrt(-2*epsilon)*(mu+nu)/2
    pmu =  y0(5)*sphi1 + y0(6)*sqrt(1-sphi1^2);
    pnu =  y0(7)*sphi2 + y0(8)*sqrt(1-sphi2^2);
    pxa = (pmu + pnu)/2
%%
%y0 = y(end,:)
clear y
tmax=600;
tphi1 = y0(1)/y0(2);
sphi1 = tphi1/sqrt(1+tphi1^2);
tphi2 = y0(3)/y0(4);
sphi2 = tphi2/sqrt(1+tphi2^2);
mu = y0(1)/sphi1;
nu = y0(3)/sphi2;

xa=sqrt(-2*epsilon)*(mu+nu)/2;
    pmu =  y0(5)*sphi1 + y0(6)*sqrt(1-sphi1^2);
    pnu =  y0(7)*sphi2 + y0(8)*sqrt(1-sphi2^2);
    pxa = (pmu + pnu)/2;
    %zwei=0.5*(y(3)*y(3)+y(4)*y(4)) -epsilon*(y(1)*y(1)+y(2)*y(2))+...
    %1/8*(y(1)*y(1)*y(2)*y(2))*(y(1)*y(1)+y(2)*y(2))
    plot(xa,pxa,'.',...
    'MarkerSize',1,'MarkerFaceColor',farb,'MarkerEdgeColor',farb)
    plot(-xa,pxa,'.',...
    'MarkerSize',1,'MarkerFaceColor',farb,'MarkerEdgeColor',farb)
    plot(xa,-pxa,'.',...
    'MarkerSize',1,'MarkerFaceColor',farb,'MarkerEdgeColor',farb)
    plot(-xa,-pxa,'.',...
    'MarkerSize',1,'MarkerFaceColor',farb,'MarkerEdgeColor',farb)

%


[t,y,te,ye,ie] = ode113(@(t,y) diahphidgl(t,y,epsilon),[0 tmax],y0,options);
toc
% Visualisierung
tphi1 = ye(:,1)./ye(:,2);
sphi1 = tphi1./sqrt(1+tphi1.^2);
tphi2 = ye(:,3)./ye(:,4);
sphi2 = tphi2./sqrt(1+tphi2.^2);
mu = ye(:,1)./sphi1;
nu = ye(:,3)./sphi2;
 xa=sqrt(-2*epsilon)*(mu+nu)/2;
 pmu =  ye(:,5).*sphi1 + ye(:,6).*sqrt(1-sphi1.^2);
 pnu =  ye(:,7).*sphi2 + ye(:,8).*sqrt(1-sphi2.^2);
 pxa = (pmu + pnu)/2;
% test = abs(xa)<0.018;                %Raus oder nicht ??
% xa(test) = [];
% pxa(test) = [];
% test = abs(xa)<0.053 & abs(pxa>1.3175);
% xa(test) = [];
% pxa(test) = [];
 xa = [xa;-xa;xa;-xa];
 pxa = [pxa;pxa;-pxa;-pxa];
  plot(xa,pxa,'.',...
     'MarkerSize',1,'MarkerFaceColor',farb,'MarkerEdgeColor',farb)
tic
Lz = y(:,2).*y(:,5) - y(:,1).*y(:,6);
Lz2 = y(:,4).*y(:,7) - y(:,3).*y(:,8);
N0 = Lz - Lz2;
pv = y(:,5:8);
pvq = pv.*pv;
Tkin = 0.5*sum(pvq,2);
zzwi = y(:,1:4);
zzwiq = zzwi.^2;
zzwiqsum = sum(zzwiq,2);
zwei = Tkin - epsilon*zzwiqsum + 1/8*(y(:,1).^2 + y(:,2).^2) .* (y(:,3).^2 + y(:,4).^2) .* zzwiqsum;
toc
% Ergebnis testen:
max(N0)
aLz = 2*(min(Lz)-max(Lz))/(min(Lz)+max(Lz))
azwei = max(abs(zwei-2))/2 
Lz(1)
shg