function y=diaphianfang(y0,epsilon)

% Berechnung der Anfangswerte der einzelnen Bahnen zu nicht-verschwindendem Drehimpuls

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

Lz = y0(2).*y0(5) - y0(1).*y0(6);
Lz2 = y0(4).*y0(7) - y0(3).*y0(8);
% Poincare-Ebene
popmu = [1,0,-8*epsilon,-8,4*Lz^2];
popmu=roots(popmu);
n = 0;
for k=1:4
    if isreal(popmu(k))
        n=n+1;
        muqi(n) = popmu(k);
    end
end
muqi=sort(muqi);
muq = muqi(1)+(muqi(end)-muqi(1))*rand; % mu quadrat
nuq = muq;
y(1) = sqrt(muq)*rand;                  % hier Vorzeichen variieren
y(2) = sqrt(muq - y(1)^2);              % y(1) ... y(4)
y(3) = sqrt(nuq)*rand;
y(4) = sqrt(nuq - y(3)^2);
% Kinetische Energie
zzwi = y(1:4);
zzwiq = zzwi.^2;
zzwiqsum = sum(zzwiq);
Tkina = 2 + epsilon*zzwiqsum - 1/8*(y(1).^2 + y(2).^2) .* (y(3).^2 + y(4).^2) .* zzwiqsum;
%
Tkinmunu = Tkina - Lz^2/muq;
pmu = sqrt(2*Tkinmunu)*rand;
pnu = sqrt(2*Tkinmunu-pmu^2);
Lztest1 = @(phi) y(2).*(pmu*sin(phi) + Lz/sqrt(muq)*cos(phi)) - y(1).*(pmu*cos(phi) - Lz/sqrt(muq)*sin(phi))-Lz;
phi1 = fzero(Lztest1,0);
Lztest2 = @(phi) y(4).*(pnu*sin(phi) + Lz/sqrt(nuq)*cos(phi)) - y(3).*(pnu*cos(phi) - Lz/sqrt(nuq)*sin(phi))-Lz;
phi2 = fzero(Lztest2,0);
y(5) = pmu*sin(phi1)+Lz/sqrt(muq)*cos(phi1);
y(6) = pmu*cos(phi1)-Lz/sqrt(muq)*sin(phi1);
y(7) = pnu*sin(phi2)+Lz/sqrt(nuq)*cos(phi2);
y(8) = pnu*cos(phi2)-Lz/sqrt(nuq)*sin(phi2);
%
%Testen
% Lz,  y(2).*y(5) - y(1).*y(6)
% Lz2, y(4).*y(7) - y(3).*y(8)
% pv = y(5:8);
% pvq = pv.*pv;
% Tkina,Tkin = 0.5*sum(pvq)
% zwei = Tkin - epsilon*zzwiqsum + 1/8*(y(:,1).^2 + y(:,2).^2) .* (y(:,3).^2 + y(:,4).^2) .* zzwiqsum