%Es werden Nebenbedingungen f\"ur ein LP-Programm erzeugt, die 
%aus Tangenten an eine Ellipse bestehen. Die Ellipse hat die Form

%      E(t):=(1-rho)z(t)+rho/z(t), rho in [0, 0.5[,  z(t)=exp(it).

%Diese Ellipse hat die Halbachsen a=1 und b=1-2rho. Ist 
%(x_0,y_0), ein Punkt auf der Ellipse, aber kein Punkt auf einer der Achsen,
%so hat die Tangente die Form

%                      x*x_0+y*y_0/(1-2rho)^2=1.

%Diese Tangente durchst\"o\ss t die beiden Achsen an den Stellen

%              (x1,y1) = (1/x_0,0); (x2,y2) = (0,(1-2rho)^2/y_0).

%Zwei aufeinander folgende Nebenbedingungen 
%x*x_0+y*y_0/(1-2rho)^2=1
%x*u_0+y*v_0/(1-2rho)^2=1 schneiden sich im Punkt

%           (xs,ys)= ((y0-v0)/Nenner,(1-2*rho)^2*(u0-x0)/Nenner) 

%                     mit Nenner=u0*y0-x0*v0.
%Das sind Ecken des zul\"assigen Bereichs. F\"ugt man noch die Ecken
%(0, 0), (x1(1), 0), (0, y2(letzte Komponente)) hinzu, erh\"alt man
%alle Ecken.                    Programm: Gerhard Opfer, 22. Januar 1999.
%===========================================================================
%            Vorspann, kann entfallen, dann m\"ussen ell, rho besetzt werden
function []=simplextest(wieviele_Punkte,Form_der_Ellipse);
ell=wieviele_Punkte; if ell<=0, ell=1; end; if ell > 100, ell=100; end;
if nargin==1, rho=0, else
  if Form_der_Ellipse=='Kreis',rho=0;   end;
  if Form_der_Ellipse=='fl_El',rho=0.4; end;
  if Form_der_Ellipse=='rd_El',rho=0.1; end;
end; %if nargin
%===========================================================================
%ell=1;                            %ell ausgew\"ahlte Punkte auf der Ellipse
%rho=0.4;%rho=0;   %rho=0 gibt Einheitskreis, ansonsten Ellipsen bis rho<0.5
T0=linspace(pi/(2*(ell+1)),ell*pi/(2*(ell+1)),ell);  
x0=cos(T0);
y0=sin(T0)*(1-2*rho);  
%===========================================================================
Nenner=ones(1,ell-1);                  %Schnittpunkte benachbarter Tangenten
Nenner=x0(2:ell).*y0(1:ell-1)-x0(1:ell-1).*y0(2:ell);
xs=(y0(1:ell-1)-y0(2:ell))./Nenner;
ys=(1-2*rho)^2*(x0(2:ell)-x0(1:ell-1))./Nenner;
%===========================================================================
                          %Durchsto\ss punkte der Tangenten durch die Achsen
x1=1./x0;y2=(1-2*rho)^2./y0;
%===========================================================================
%Hinzunahme der \"ubrigen Ecken des zul\"assigen Bereichs
xs=[x1(1) xs];ys=[0 ys];xs=[0 xs];ys=[0 ys];xs=[xs 0];ys=[ys y2(ell)];
%xs,ys enthalten jetzt die Koordinaten aller ell+2 Ecken.
%===========================================================================
%       Besetzung des LP-Programms Ax=b, c'x=min, x>=0 und I=Basisindexmenge
%                           A enth\"alt die Koeffizienten der Tangente, s. o.
A=[x0',y0'/(1-2*rho)^2];
A=[A,eye(ell)];I=3:ell+2;
c=zeros(ell+2,1);c(1)=-1;c(2)=-1/(1-2*rho);      %Also max x_1+x_2/(1-2rho)
b=ones(ell,1);
I=3:ell+2;
%===========================================================================
%                                               Ausf\"uhren des LP-Programms
info=0;zeler=0;                      %und Zeichnen der ausgerechneten Punkte
while info==0 & zeler <= ell+3
  zeler=zeler+1;
  if zeler==1
    disp(['Beginn von Simplexiteration Nr. ',int2str(zeler)]); else
    disp(['                            Nr. ',int2str(zeler)]); 
  end; %if
  [wert,x,I,info]=simplex(A,b,c,I);
  plot(x(1),x(2),'+'); hold on
  if zeler==1
    text(x(1),x(2)-y2(ell)/20,'(1)'); 
    else
    text(1.02*x(1),1.02*x(2),['(',int2str(zeler),')']);
  end; %if                              
end; %while
%===========================================================================
t=linspace(0,pi/2);               %Zeichnen der Ellipse im ersten Quadranten
z=exp(i*t);
E=(1-rho)*z+rho./z;
plot(E);                                            %k\"onnte auch entfallen
%===========================================================================
                                         %Zeichnen des zul\"assigen Bereichs
plot(x0,y0,'o',xs,ys,'*',[x1' zeros(ell,1)]',[zeros(ell,1) y2']',...
     [0 max(x1)],[0 0],[0 0],[0 max(y2)]);
pause(5);                              %Bevor der zul. Bereich gefüllt wird, 
                                      %ein kurzer Blick auf die Konstruktion
h=fill(xs,ys,'r');
axis equal      %Ohne diesen Befehl, sind alle Bereiche i. w. Viertelkreise
axis([0  x1(1)  0  y2(ell)]);
legend(h,'Polygon der Nebenbedingungen');
xlabel([int2str(ell),' Tangenten an eine Ellipse mit rho = ',num2str(rho)]);
hold off
%===========================================================================

%Testergebnisse:
%Alle Basisvektoren im Fall ell=2, rho=0.25         
%   not in I        X(I)             x_1, x_2
%-------------------------------------------------------
%    1   2         1         1        (0, 0)      (Nr.1)
%    1   4    0.5774    0.4226        (0, 0.5774) (Nr.4)
%    2   3    1.1547    0.4226   (1.1547, 0)      (Nr 2)
%    3   4    0.7321    0.3660   (0.7321, 0.3660) (Nr 3)









