%
%   Beispiel (9.7) - Anpassung einer Ellipse an Messdaten
%
%   Messdaten :
%
datname='bahn.dat'; 	% Name der Datei mit den Daten 
load(datname);		% Laden der Daten
%
%	Da der Name der Datei variabel ist, werden die Daten in einer Datei mit
%	veraenderlichem Namen gespeichrt. Die folgenden zwei Zeilen speichern sie
%	in der Matrix werte
% 
matrixname= datname(1:findstr(datname,'.')-1);
werte=eval(matrixname);
%
xi=werte(:,1);
eta=werte(:,2);
%
%   Loesung des linearen Ausgleichsproblems :
%
m = length(xi);
format long e
A=[eta.*eta xi.*eta xi eta ones(m,1)];
b=xi.*xi;
'Die Parameter sind'
x=A\b
%
%   Normalform der Quadrik (Eigenwerte u. Eigenvektoren) :
%
D  = [1 -x(2)/2;-x(2)/2 -x(1)];
[W,E] = eig(D);
la = diag(E);
if min(la)<=0, error('keine Ellipse'),end;
d  = [x(3); x(4)];
d  = W'*d;
d  = d./(2*la);
d1 = sqrt(la).*d;
e  = x(5) + d1'*d1;
%   Halbachsen
al = sqrt(e/la(1));
be = sqrt(e/la(2));
%   Brennpunkte
brenn1 = [sqrt(al*al-be*be);0]+d;
brenn1 = W*brenn1;
brenn2 = [-sqrt(al*al-be*be);0]+d;
brenn2 = W*brenn2;
%   Mittelpunkt (Ursprung)
mitte = W*d;
%
%   Plot der Ellipse :
%
t  = [0:0.01:6.3]';
y  = [al*cos(t)+d(1) be*sin(t)+d(2)];
y  = (W*y')';
y1 = y(:,1);
y2 = y(:,2);
plot(y1,y2)
%
u1 = min(y1);
o1 = max(y1);
u2 = min(y2);
o2 = max(y2);
l1 = u1 - 0.1*(o1-u1);
r1 = o1 + 0.1*(o1-u1);
l2 = u2 - 0.1*(o2-u2);
r2 = o2 + 0.1*(o2-u2);
 axis([l1 r1 l2 r2])
hold on
plot(xi,eta,'o')
plot(brenn1(1),brenn1(2),'o')
plot(brenn2(1),brenn2(2),'o')
plot(mitte(1),mitte(2),'+')
hold off




