%% Chaotisches Billard 

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% Version 2.7  2016


%% Festlegung der Berandung und Startwerte
R0 = 1;
epsilon = .25;
phi = linspace(0,2*pi,200);
x = R0*(1+epsilon*cos(phi)).*cos(phi);
y = R0*(1+epsilon*cos(phi)).*sin(phi);
plot(x,y,'k'),shg
hold on
% Teilchenbahn Start
vaus = [-1,1]; % bewegt sich das Teilchen in positive oder negative x-Richtung 
vt = vaus(randperm(2,1));
vstart = vt;
phi0 = rand*pi;
mt = tan(phi0);
bt = R0*0.9*rand;
bstart = bt;
%% Berechnung des Einfallwinkels
hphi = @(phi)  R0*(1+epsilon*cos(phi)).*sin(phi) - mt * R0 * (1+epsilon*cos(phi)).*cos(phi)-bt;
phi0 = fzero(hphi,phi0);
%% Plotten der Starttrajektorie
xs = R0*(1+epsilon*cos(phi0)).*cos(phi0);
x0 = xs;
xt = linspace(0,xs);
yt = mt*xt + bt;
ys = yt(end);
%plot(xt(1),yt(1),'mp')
plot(xt,yt,'k')
%plot(xt(end),yt(end),'gp')
%%  Berechnen der Tangente
dx = -epsilon*sin(phi0)*cos(phi0) - (1+epsilon*cos(phi0))*sin(phi0);
dy = -epsilon*sin(phi0)*sin(phi0) + (1+epsilon*cos(phi0))*cos(phi0);
mtan = dy/dx;
phit = atan(mtan);
%%  Visualisieren der Tangente
mtan = dy/dx;
btn = R0*(1+epsilon*cos(phi0)).*sin(phi0) - mtan*xs;
xt = linspace(xs-0.2,xs+.2);
yt = mtan*xt + btn;
%plot(xt,yt,'r'),shg
axis equal
%%  Berechnen der Normale fuer Visualisierungen
h1 = @(x) (1 + x.*mtan)./sqrt((1+x.^2).*(1+mtan.^2));
msen = fzero(h1,0);
%%  dient der Visualisierung
x = linspace(xs-0.05,xs+0.05);
y = msen*x + ys - msen*xs; 
hold on
%plot(x,y,'r'),shg
%% Berechnung der reflektierten Bahn
tanp = (mtan - mt)/(1+mtan*mt);
mt2 = (mtan+tanp)/(1-tanp*mtan);
%mt
%% dient der Visualisierung
x = linspace(xs,xs+0.5);
y = mt2*x + ys - mt2*xs;
%plot(x,y,'k'),shg
%%  Berechnung des naechsten Schnittpunkts
bt2 = ys - mt2*xs;
hphi2 = @(phi)  R0*(1+epsilon*cos(phi)).*sin(phi) - mt2 * R0 * (1+epsilon*cos(phi)).*cos(phi)-bt2;
phi0 = fzero(hphi2,-pi+1);
xs2 = R0*(1+epsilon*cos(phi0)).*cos(phi0);
% xs,xs2
%% Plotten der naechsten Bahn der 
xt = linspace(xs,xs2);
yt = mt2*xt + bt2;
ys2 = yt(end);
plot(xt,yt,'k')
%% und los gehts
for k=1:100
    xs = xs2;
    [xs2,mt2,bt2,phi0] = GerGESchnitt(epsilon,R0,phi0,mt2,xs2,ys2);
    %
    xt = linspace(xs,xs2);
    yt = mt2*xt + bt2;
    ys2 = yt(end);
    plot(xt,yt,'k')
end

