% Beispiel nichtverschwindendem Drehimpuls

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

% ausgewaehlte Startwerte
tic
y0(1)=0;
y0(3)=0;
y0(6)=0;
y0(8)=0;
y0(2)=rand;
y0(4)=y0(2); 
% 
epsilon=-0.5;   % sinnvolle Werte sind negativ.
% 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);
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)
M = ones(9); M(9,9)=0;
options = odeset('RelTol',1e-11,'AbsTol',1e-11,'MaxStep',0.1);%,'Mass',M);
tmax=200;
[t,y] = ode45(@(t,y) diahphidgl(t,y,epsilon),[0 tmax],y0,options);
toc
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 
