%% Elliptisches Billard  Visualisierung der einzelnen Schritte 
% daher regulaeres Verhalten

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


%% Mittelpunktsgleichung der Ellipse
%clear,clc,close all
a = 2;                 % grosse Halbachse
b = 1;%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
%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);
mstart = mt;
bt = b*0.9*rand;
bstart = bt;
dx = zeros(20,2);
dy = dx;
for k=1:2
    if k==1
        figure%subplot(5,1,1:2)
        plot(x,y,'-k',x,-y,'-k'),axis equal,shg
        hold on
    else
        figure%subplot(5,1,3:4)
        plot(x,y,'-k',x,-y,'-k'),axis equal,shg
        hold on
        bt = bstart*0.999
        vt = vstart;
        mt = mstart;
    end
    [xs,ys] = GerEllSchnitt(a,b,mt,bt,vt,0,0);
    dx(1,k) = xs;
    dy(1,k) = ys;
    xt = linspace(0,xs);
    yt = mt*xt + bt;
    plot(xt,yt,'k')
    
    
    
    % Beispiel
x0 = xs;
%y0 = b*sqrt(1-(x0/a).^2);
y0 = ys;
%xp = linspace(x0-0.75,x0+0.75,10);
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);

 dx(2,k) = xs;
 dy(2,k) = 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)
    
    
    
    n=2;
    while n<20
        n=n+1;
        x0 = xs;
        y0 = ys;
        m  = -(b/a).^2*x0/y0;
%bt  = y0 - m*x0;
        alpha = atan(m);
        beta = atan(mt);
%mt = tan(pi-2*alpha+beta);  % reflektierte Gerade
        mt = tan(-beta+2*alpha);
        bt = ys - mt*xs;
        [xs,ys] = GerEllSchnitt(a,b,mt,bt,vt,x0,y0);
        dx(n,k) = xs;
        dy(n,k) = 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)
    end
end
hold off
%subplot(5,1,1:2), hold off, axis tight
ddx = (dx(:,1)-dx(:,2)).^2;
ddy = (dy(:,1)-dy(:,2)).^2;
da = sqrt(ddx+ddy);
figure%subplot(5,1,5)
plot(da)
