%% Elliptisches Billard 
% daher regulaeres Verhalten

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB 
% Version 2.3  2016
%% Mittelpunktsgleichung der Ellipse
%clear,clc,close all

% Parameter frei waehlbar
a = 2;                 % grosse Halbachse   
b = 1.5;%2;%2/sqrt(2); % kleine Halbachse
p = b^2/a;             % Halbparameter der Ellipse
%
x = linspace(-a,a);
y = b*sqrt(1-(x/a).^2);
plot(x,y,'k',x,-y,'k'),axis equal,shg    % Darstellung der Ellipse
hold on
%% Teilchenbahn Start
vaus = [-1,1]; % bewegt sich das Teilchen in positive oder negative x-Richtung 
vt = vaus(randperm(2,1));
vstart = vt;
mt = tan(rand*pi);
bt = b*0.9*rand;
bstart = bt;
[xs,ys] = GerEllSchnitt(a,b,mt,bt,vt,0,0);
xt = linspace(0,xs);
yt = mt*xt + bt;
plot(xt,yt,'k')
%%
% Berechnung der Tangente
% Beispiel
x0 = xs;
y0 = ys;
%
m  = -(b/a).^2*x0/y0;
bt  = y0 - m*x0;
%yp = m*xp+bt;
%plot(xp,yp,'k',x0,y0,'mp')
% elastische Streuung
alpha = atan(m);
beta = atan(mt);
%mt = tan(pi-2*alpha+beta);  % reflektierte Gerade
mt = tan(-beta+2*alpha);
bt = ys - mt*xs;
vt = 0;
[xs,ys] = GerEllSchnitt(a,b,mt,bt,vt,xs,ys);
% xp = linspace(x0,x0-0.5*sign(mt))
% yp = mt*xp+bt
% plot(xp,yp,'g'), axis equal
%plot(xs,ys,'kp')
xst = linspace(x0,xs);
yst = mt*xst +bt;
plot(xst,yst,'k')
%%
n=0;
while n<200
    n=n+1;
    x0 = xs;
    y0 = ys;
    m  = -(b/a).^2*x0/y0;

    alpha = atan(m);
    beta = atan(mt);

    mt = tan(-beta+2*alpha);     % reflektierte Gerade
    bt = ys - mt*xs;
    [xs,ys] = GerEllSchnitt(a,b,mt,bt,vt,xs,ys);

    xst = linspace(x0,xs);
    yst = mt*xst +bt;
    plot(xst,yst,'k')
end
hold off
