% http://geometrie.hs-nb.de/MATLAB-Programme/DAEll.m   Bearbeitungsstand: 12.11.2004
% (C) Martin Nitschke, nitschke@hs-nb.de

% Deklaration der symbolischen Variablen

syms('x0','y0','a','b','alph','t','real');
% Faktor für Umrechnung  ° -> rad
deg=sym(pi)/180;
% Definition A und P
A=sym([5;2]);
P=[x0 + a*cos(alph)*cos(t) - b*sin(alph)*sin(t)
   y0 + a*sin(alph)*cos(t) + b*cos(alph)*sin(t)];
% Einsetzen der Werte a, b usw.
P=subs(P,{a b alph x0 y0},{1.5 1.2 70*deg 2 1});
% Parametrisierung als Inline-Funktionen
Px=inline(char(P(1)),'t'),   Py=inline(char(P(2)),'t')
% zu lösende Gleichung als Inline-Funktion
Gl=inline(char(simple(dot(A-P,diff(P,t)))),'t');
% Aufteilung des Parameterbereichs zur Lösungssuche
tRg=[0:360]*double(deg);
% Initialisierung von AP und d
AP=zeros(2,0);   d=zeros(1,0);
% Schleife über die Teilintervalle; in jedem Intervall
% Vorzeichenwechsel überprüfen und ggf. Lösung berechnen
for k=2:length(tRg)
  if Gl(tRg(k-1))*Gl(tRg(k)) < 0
    tA=fzero(Gl,[tRg(k-1) tRg(k)]);
    AP(:,end+1)=[Px(tA);Py(tA)];
    d(end+1)=norm(AP(:,end)-double(A));
  end
end
[dmin kmin]=min(d);
fprintf('Minimaler Abstand  = %.9f\n',dmin);
fprintf('Zugeh. Ellipsenpkt.=(%.9f,%.9f)\n',AP(:,kmin));
