unit d2basic;
{Copyright 2008-20011 Manfred Rosendahl
 Private Nutzung erlaubt
 Bei kommerzieller Nutzung bitte Mail an: ros at uni-koblenz.de
}
{
Datenrepräsentation
Geraden,Linien  Hesse Form a,b,c  a*x+b*y=c
a cos der Normalen b sin der Normalen  c Abstand vom Nullpunkt(gerichtet)
geraden mit Normalen Richtung alpha und Abstand d in Richtung der Normalen;
Richtung der Geraden(Strahl) ist alpha-90;0<=alpha<360; d=c

Segmente x1,y1,x2,y2
Kreis   x,y,r
Winkel in Grad
}
interface
type real=double;
const errorvalue:double=0.0;
var delta:double=0.000001; // Geraden Parallel
var tol:double=0.000001;
var eps:double=1e-6;
var tolaaff:double=0.0000001; // für xxcirclecc0 wenn eine Lösung eine gerade ergibt
const maxpunkte=50;
const radiant=pi/180;
type tpunkte=array[1..maxpunkte]of record
 x,y:double
end;
type tmatrix=array[1..2,1..2]of double;
 tvector=array[1..2]of double;

function sign(x:double):integer;

function eq(a,b:double):boolean;

function eq3(a1,b1,c1,a2,b2,c2:double):boolean;

var xloesung,yloesung,rloesung:array[0..15]of double;
 loesung:integer;

function newloesung(x,y,r:double):integer;
{newloesung>0 neue loesung
           =-1 alte loesung
}

function atan2(y,x:double):double;

procedure rotate90p(xp,yp,xc,yc:double;var x,y:double);
{xp,yp wird um xc,yc um 90 Grad gedreht}

procedure rotateminus90p(xp,yp,xc,yc:double;var x,y:double);
{xp,yp wird um xc,yc um -90 Grad gedreht}

procedure rotatep(xp,yp,xc,yc,alpha:double;var x,y:double);
{xp,yp wird um xc,yc um alpha Grad gedreht}

procedure line2p(x1,y1,x2,y2:double;var a,b,c:double;var fail:boolean);
{Umwandlung Linie mit 2 Punkten in Hesse Normalform}
procedure gerade(d,alpha:double;var a,b,c:double);
// gerade in a,b,c
procedure hesse2gerade(a,b,c:double;var d,alpha:double);
// a,b,c nach gerade
function winkell(a,b:double):double;
{Winkel aus sin,cos}

function parallel2l(a1,b1,c1,a2,b2,c2:double;var gleich:boolean):boolean;
// 2 Linien sind parallel
// gleich: die beiden Linien sind gleich

procedure pointonline(x0,y0,x1,y1,d:double;var x,y:double);
// x,y liegt auf de Linie x0,y0-x1,y1 im Abstand d von x0,y0

procedure line2gerade(x1,y1,x2,y2:double;var c,alpha:double;var fail:boolean);
{Umwandlung Linie mit 2 Punkten in Hesse Normalform}

procedure pointcl(x2,y2,r2,a1,b1,c1:double;mode:integer;var x,y:double;
                  var fail:boolean);
{Schnittpunkt linie a1,b1,c1 und Kreis x2,y2,r2}

procedure pointcl2p(x2,y2,r2,x1,y1,xx1,yy1:double;mode:integer;var x,y:double;
                  var fail:boolean);
{Schnittpunkt linie x1,y1,xx2,yy2 und Kreis x2,y2,r2}

procedure pointcg(x2,y2,r2,c,alpha:double;mode:integer;var x,y:double;
                  var fail:boolean);
{Schnittpunkt gerade c,alpha und Kreis x2,y2,r2}

procedure point2l(a1,b1,c1,a2,b2,c2:double;var x,y:double;var fail:boolean);
{Schnittpunkt Linien a1,b1,c1 und a2,b2,c2}

procedure pointl2pl2p(x1,y1,x2,y2,x3,y3,x4,y4:double;var x,y:double;var fail:boolean);
{Schnittpunkt Linien x1,y1,x2,y2 und x3,y3,x4,y4}

procedure point2g(c1,alpha1,c2,alpha2:double;var x,y:double;var fail:boolean);
{Schnittpunkt gerade c1,alpha1 und c2,alpha2}

procedure point2c(x1,y1,r1,x2,y2,r2:double;mode:integer;var x,y:double;
                   var fail:boolean);
{Schnittpunkt Kreis x1,y1,r1 und Kreis x2,y2,r2}

procedure tancc(x1,y1,r1,x2,y2,r2:double;ialter:integer;var x,y:double;var fail:boolean);
{get one tangent point of two circles
 ialter=0,1,2,3 --> point on the second circle;
 ialter=4,5,6,7 --> point on the first circle
}

procedure tangentecc(xx1,yy1,rr1,xx2,yy2,rr2:double;si:integer;
                  var x1,y1,x2,y2:double;var fail:boolean);
{Tangentiale Linie zwischen 2 Kreisen}

procedure tangentepc(x1,y1,x2,y2,r2:double;ialter:integer;var x,y:double;
                  var fail:boolean) ;
{tangent point x,y from a given point x1,y1 onto a circle x2,y2,r2}
procedure tangentlinepc(x1,y1,x2,y2,r2:double;ialter:integer;var a,b,c:double;
                  var fail:boolean) ;
{tangent linie a,b,c from a given point x1,y1 onto a circle x2,y2,r2}

procedure llotlc(a,b,c,x,y,r:double;mode:integer;var aa,bb,cc:double);
// Linie als Lot von Kreis x,y,r an Linie a,b,c

procedure lotlc(xx2,yy2,xx3,yy3,xx1,yy1,r:double;alter:integer;var x,y:double);
{ get a tangent point of circle generated by a vertical line of the given line
  alter=0,1 --> point is on the circle
  alter=2,3 --> point is on the line
}

procedure lotpc(x1,y1,x2,y2,r2:double;mode:integer;var x,y:double);
{lot line from a point onto a circle }

procedure parallelepl(xx0,yy0,xx1,yy1,xx2,yy2:double;var x1,y1,x2,y2:double);
{ get parallel line of the given line.this parallel line is located by the
  given point }

procedure parallelelc(lx1,ly1,lx2,ly2,cx,cy,r:double;alter:integer;var x1,y1,x2,y2:double);
{ a tangent line is generated which is also parallel to the given line }

procedure Lotpl2p(const Px,Py,x1,y1,x2,y2:double;var x,y:double;var fail:boolean);
//procedure lotpl2p(x1,y1,x2,y2,x3,y3:double;var x,y:double;var fail:boolean);
{Lot vom punkt x1,y1 auf die Linie x2,y2,x3,y3}
//var a,b,c:double;

procedure linelc(x1,y1,x2,y2,x3,y3,r3:double;lalter:integer;
                  var xx1,yy1,xx2,yy2:double);
{ get a tangent line of the circle. this line is vertical or parallel
  to the given line .
  lalter=0,1 --> vertical to the given line.
  lalter=2,3 --> parallel to the given line.
  lalter=4   --> vertical line between line and circle.  }

procedure linepc(x1,y1,x4,y4,r:double;lalter:integer;
                  var xx1,yy1,xx2,yy2:double;var fail:boolean);
{tangent or lot line from a point onto a circle }

procedure linell(l1x1,l1y1,l1x2,l1y2,l2x1,l2y1,l2x2,l2y2:double;lalter:integer;
                 var x1,y1,x2,y2:double);
{ get a line between two lines }

procedure winkelhalbll(l1x1,l1y1,l1x2,l1y2,l2x1,l2y1,l2x2,l2y2:double;lalter:integer;
                 var x1,y1,x2,y2:double;var fail:boolean);
{ winkelhalbierende }

procedure mittelsenkrechte(x1,y1,x2,y2:double;var a,b,c:double);

function getphi(xa,ya,xi,yi:double):double;
{winkel von xa,ya nach xi,yi in grad}

function diffphi(phi2,phi1:real):real;
// diffphi >0 Drehung um diffphi nach links führt von phi1 nach phi2
// diffphi <0 analog nach rechts
function drehung(links:boolean;phi1,phi2:real):real;
//drehwinkel um links oder rechtsrum von phi1 nach phi2 zu kommen
function normphi(phi:real):real;
// phi in denBereich 0 .. 360 bringen
function gegenwinkel(phi:real):real;
// gegenwinkel 90-> -270 oder -180 -> 180

procedure parallel2p(xx1,yy1,xx2,yy2,d:double;lalter:integer;var x1,y1,x2,y2:double);
{get a parallel line of the given line with a given distance d}

function distance2p(x1,y1,x2,y2:double):double;

function distancepl(x1,y1,a2,b2,c2:double):double;
{Abstand zwischen Punkt x1,y1 und der linie a2,b2,c2 }

function distancepl2p(x1,y1,x2,y2,x3,y3:double):double;
{Abstand zwischen Punkt x1,y1 und der linie x2,y2,x3,y3 }

procedure lotpl(x1,y1,a2,b2,c2:double;var x,y:double);
{Lot vom punkt x1,y1 auf die Linie a2,b2,c2}

procedure rotatel(al,bl,cl,xc,yc,alpha:double;var a,b,c:double);
{al,bl,cl wird um xc,yc mit alpha gedreht}

procedure winkelhalb(a1,b1,c1,a2,b2,c2:double;mode:integer;
                     var a,b,c:double;var fail:boolean);
{Gerade a,b,c die Winkelhalbierende von Geraden a1,b1,c1 und a2,b2,c2}

function distancepc(x1,y1,x2,y2,r2:double;s:integer):double;
{Abstand zwischen Punkt x1,y1 und dem Kreis x2,y2,r2 }

procedure circle2lr(a1,b1,c1,a2,b2,c2,r:double;mode:integer;
                     var x,y:double;var fail:boolean);
{Mittelpunkt eines Kreises tangential zu 2 linien und mit radius r}

procedure circle2cr(x1,y1,r1,x2,y2,r2,r:double;mode:integer;
                     var x,y:double;var fail:boolean);
{Mittelpunkt eines Kreises tangential zu 2 Kreisen und mit radius r}

procedure circlelcr(a1,b1,c1,x2,y2,r2,r:double;mode:integer;
                     var x,y:double;var fail:boolean);
{Mittelpunkt eines Kreises tangential zu Linie und Kreis und mit radius r}

procedure circle3p(x1,y1,x2,y2,x3,y3:double;var x,y,r:double;var fail:boolean);
{Kreis aus 3 Punkten}

procedure circle3s(xa1,ya1,x1,y1,xa2,ya2,x2,y2,xa3,ya3,x3,y3:double;mode:integer;
                    var x,y,r:double;var fail:boolean);
{Kreis Tangential an 3 Segmente}

procedure circle3l(a1,b1,c1,a2,b2,c2,a3,b3,c3:double;mode:integer;
                    var x,y,r:double;var fail:boolean);
{Kreis Tangential an 3 Geraden}

procedure change(var x,y:double);

function flaeche(var punkte:tpunkte;n:integer):double;

procedure schwerpunkt(var punkte:tpunkte;n:integer;
                      var flaeche:double;var xschwer,yschwer:double);
{xschwer,yschwer schwerpunkt}

procedure getkoord(x0,y0,xu,yu,xv,yv,xp,yp:double;var u,v,unorm,vnorm:double);
{von den vektoren 0-u und 0-v wird ein Koordinatensystem aufgespannt;
 xp,yp ist ein punkt im kartesischen koordinaten
 u,v sind die koordinaten des Punktes im schiefen Koordinaten
 es ist die l”sung des Gleichungssystems
 u*vekxu+v*vekxv=distx
 u*vekyu+v*vekyv=disty
 unorm,vnorm sind dieselben u,v Kordinaten mit der l„nge der 0-u,0-v
 Vektoren
}

procedure orthokoord(x0,y0,xu,yu,xp,yp:double;var u,v:double);
{Es wird ein Koordinatensystem mit dem Ursprung x0,y0 und der x-achse xu,yu
 aufgespannt;
 xp,yp ist ein punkt im kartesischen koordinaten
 u,v sind die koordinaten des Punktes im neuen Koordinaten
 es ist die l”sung des Gleichungssystems
 u*vekxu+v*vekxv=distx
 u*vekyu+v*vekyv=disty
 unorm,vnorm sind dieselben u,v Kordinaten mit der l„nge der 0-u,0-v
 Vektoren
}

procedure getpolar(xp,yp,xc,yc:double;var alpha,r:double);
{xp,yp ist Polarpunkt von xc,yc mit abstand r and Winkel alpha}

procedure polarpoint(xc,yc,alpha,d:double;var x,y:double);
{xp,yp ist Polarpunkt von xc,yc mit abstand r and Winkel alpha}

procedure linepw(x,y,alpha:double;var a,b,c:double);
{Linie aus Punkt und Richtungswinkel}

function alpha3p(x0,y0,x1,y1,x2,y2:double):double;
{Kleinster Winkel an x0,y0 nach x1,y1 und x2,y2}

function alpha3pdir(x0,y0,x1,y1,x2,y2:double):double;
{Gerichteter Winkel an x0,y0 von x1,y1 nach x2,y2}

procedure ctransform(isline:boolean;x,y,r,xc,yc,rc:double;
                     var istline:boolean;var xt,yt,rt:double);
{Kreistransformation Linie/Kreis x,y,r wird am Kreis xc,yc,rc transformiert
 und ergibt xt,yt,rt
}

procedure xcirclecpar(x1,y1,r1,a2,b2,c2,c3:double;mode:integer;var x,y,r:double;
                     var fail:boolean);
{Hilfroutine:
 x,y,r Kreis tangential an die beiden parallelen Geraden a2,b2,c2 und
 a2,b2,c3 und den Kreis x1,y1,r1
}

procedure xcirclec2l(alpha,x1,y1,r1:double;mode:integer;var x,y,r:double;
                     var fail:boolean);
{ Bestimmt die Kreise an 2 symmetrische Linien durch den Nullpunkt,
 die einen Winkel von +/- alpha zur X-Achse haben
Die Lösungen 0 bis 3 werden erzielt durch folgendes Gleichungssystem:
(x-x1)^2+y1^2=(r+/-r1)^2
r=x*sin(alpha)
Sie liegen auf der X Achse
Die Lösungen 4 bis 7 durch das Gleichungssystem:
(y-y1)^2+x1^2=(r+/-r1)^2
r=y*cos(alpha)
Sie liegen auf der y Achse
}
procedure xcirclel2p(a2,b2,c2,x1,y1,x2,y2:double;mode:integer;
          var x,y,r:double;var fail:boolean);
// Kreis Linie und 2 Punkte

procedure circlec2l(x1,y1,r1,a2,b2,c2,a3,b3,c3:double;mode:integer;
                    var x,y,r:double;var fail:boolean);
{Kreis Tangential an Kreis und 2 Linien}

procedure circle2lp(a1,b1,c1,a2,b2,c2,x3,y3:double;mode:integer;
                    var x,y,r:double;var fail:boolean);
{Kreis Tangential an Punkt und 2 Linien}

procedure circlel2p(a1,b1,c1,x2,y2,x3,y3:double;var x,y,r:double;var fail:boolean);
{Kreis Tangential an 2 Punkte und 1 Linie}

procedure xcircleclp(x1,y1,r1,a2,b2,c2,x3,y3:double;mode:integer;
          var x,y,r:double;var fail:boolean);
{Hilfroutine f�r Kreis Tangential an Kreis,Linie und Punkt}

procedure circleclp(x1,y1,r1,a2,b2,c2,x3,y3:double;mode:integer;
          var x,y,r:double;var fail:boolean);
{Kreis Tangential an Kreis,Linie und Punkt}

procedure xcirclecc0(px1,py1,pr1,px2,py2,pr2:double;mode:integer;
                   var x,y,r:double;var fail:boolean);
{Hilfsroutine
 Emittelt L”sungen tangential zu den Kreisen px1,py1,pr1 und px2,py2,pr2
 und zu dem Punkt 0,0.
 Entwickelt nach einem Hinweis von Ron Capelli(IBM)
}

function pointinc(x1,y1,x2,y2,r2:double):boolean;
{Punkt innerhalb des Kreises}

procedure senkrechtedurchPunkt(a,b,c,x,y:double;var aa,bb,cc:double);

procedure circleptx(x1,y1,a2,b2,c2,x3:double;mode:integer;var x,y,r:double;var fail:boolean);
// Kreis durch Punkt x1,y1 und Tangential an a2,b2,c2 und Mittelpunkt X-Koordinate x3
procedure circlepty(x1,y1,a2,b2,c2,y3:double;mode:integer;var x,y,r:double;var fail:boolean);
// Kreis durch Punkt x1,y1 und Tangential an a2,b2,c2 und Mittelpunkt Y-Koordinate y3
procedure circlectx(x1,y1,r1,a2,b2,c2,x3:double;mode:integer;var x,y,r:double;var fail:boolean);
// Kreis Tangential zu Kreis x1,y1,r1 und Tangential an a2,b2,c2 und Mittelpunkt x-Koordinate x3
procedure circlecty(x1,y1,r1,a2,b2,c2,y3:double;mode:integer;var x,y,r:double;var fail:boolean);
// Kreis Tangential zu Kreis x1,y1,r1 und Tangential an a2,b2,c2 und Mittelpunkt y-Koordinate y3
procedure circleptl(x1,y1,a1,b1,c1,a2,b2,c2:double;mode:integer;var x,y,r:double;var fail:boolean);
// Kreis durch Punkt x1,y1 und Tangential an a1,b1,c1 und Mittelpunkt durch Linie a2,y2,c2
procedure circle0yc(y1,x2,y2,r2:double;mode:integer;var x,y,r:double;var fail:boolean);
// Kreis durch Nullpunkt, Tangential an y=y1 und Mittelpunkt auf Kreis x2,y2,r2
procedure circle0tc(a1,b1,c1,x2,y2,r2:double;mode:integer;var x,y,r:double;var fail:boolean);
// Kreis durch Nullpunkt tangential an a1,b1,c1 und Mittelpunkt durch Kresi x2,y2,r2
procedure circleptc(x1,y1,a1,b1,c1,x2,y2,r2:double;mode:integer;var x,y,r:double;var fail:boolean);
// Kreis durch x1,y1 tangential an a1,b1,c1 und Mittelpunkt auf Kreis x2,y2,r2
procedure circle3c(px1,py1,pr1,px2,py2,pr2,px3,py3,pr3:double;pmode:integer;
                  var x,y,r:double;var fail:boolean);
{Kreis tangential an 3 Kreise}

procedure circlec2p(x1,y1,r1,x2,y2,x3,y3:double;pmode:integer;
                  var x,y,r:double;var fail:boolean);
{Kreis tangential an Kreise und 2 Punkten}

procedure circle2cp(x1,y1,r1,x2,y2,r2,x3,y3:double;pmode:integer;
                    var x,y,r:double;var fail:boolean);
{Kreis tangential an 2 Kreise und 1 Punkt}

procedure circle2cl(x1,y1,r1,x2,y2,r2,a3,b3,c3:double;mode:integer;
                    var x,y,r:double;var fail:boolean);
{Kreis tangential an 2 Kreise und 1 Linie}

procedure circleclmc(x1,y1,r1,a2,b2,c2,x3,y3,r3:double;mode:integer;
 var x,y,r:double;var fail:boolean);
procedure circleclml(x1,y1,r1,a2,b2,c2,a3,b3,c3:double;mode:integer;
 var x,y,r:double;var fail:boolean);
{Kreis tangential an Kreis,Linie und Mittelpunkt auf Kreis bzw. Linie}

procedure circle2cmc(x1,y1,r1,x2,y2,r2,x3,y3,r3,rold:double;mode:integer;
 var x,y,r:double;var fail:boolean);
procedure circle2cml(x1,y1,r1,x2,y2,r2,a3,b3,c3:double;mode:integer;
 var x,y,r:double;var fail:boolean);
{Kreis tangential an 2 Kreise und Mittelpunkt auf Kreis bzw. Linie}

procedure circlecpmx(x1,y1,r1,x2,y2,x3:double;mode:integer;var x,y,r:double;var fail:boolean);
{Kreis tangential an Kreis, Punkt und x Mittelpunkt}
procedure circlecpmy(x1,y1,r1,x2,y2,y3:double;mode:integer;var x,y,r:double;var fail:boolean);
{Kreis tangential an Kreis, Punkt und x Mittelpunkt}
procedure circle2cmx(x1,y1,r1,x2,y2,r2,x:double;mode:integer;
 var y,r:double;var fail:boolean);
{Kreis tangential an 2 Kreise und Mittelpunkt auf Y Achse}
procedure circle2cmy(x1,y1,r1,x2,y2,r2,y:double;mode:integer;
 var x,r:double;var fail:boolean);
{Kreis tangential an 2 Kreise und Mittelpunkt auf X Achse}

function min(a,b:double):double;
{ liefert minimum zweier Zahlen}
function max(x,y:double):double;
{ liefert maximum zweier Zahlen}
procedure extenarc(xm,ym,r,phimin,phimax:double;
          var xmin,ymin,xmax,ymax:double);

procedure fillet(x1,y1,x2,y2,x3,y3,x4,y4,r:double;
                 var xm,ym,xp1,yp1,xp2,yp2:double;var fail,vertauscht:boolean);
// Verundung von Linie x1,y1-x2,y2 und Linie x3,y3,x4,y4 mit Radius r;
// Ergebnis Bogen mit mittelpunkt xm,ym und von xp1,yp1 nach xp2,yp2

function posauflinie(x1,y1,x2,y2,x,y:double;var posi:double):boolean;
// Wenn x,y auf der linie von x1,y1 nach x2,y2 liegt wird mit posi die Position
// angegeben (posi=0 x,y=x1,y1  posi=1 x,y=x2,y2)
// Wenn x,y nicht auf der geraden der Linie liegt ist result=false;
function test1(a,b,c:double):boolean;

procedure tangente(x,y,xm,ym:double;var a,b,c:double);
// Tangente an Kreis mit Mittelpunkt xm,ym im Punkt x,y
procedure parallelep(a,b,c,x,y:double;var aa,bb,cc:double);
// Parallele zu a,b,c durch den Punkt x,y
procedure punktauflinie(x1,y1,x2,y2,d:double;var xx,yy:double);
// Punkt auf x1,y1-x2,y2 mit abstand d von x1,y2
procedure pointmittecl(x1,y1,r1,a2,b2,c2:double;var x,y:double);
// Mittenpunkt zwischen Kreis und Linie
procedure matrixmul(a,b:tmatrix;var c:tmatrix);
procedure matvecmul(b:tmatrix;a:tvector;var c:tvector);
procedure setmat(var mat:tmatrix;x11,x12,x21,x22:double);

implementation
uses math,sysutils;

procedure drawkreis(x,y,r:double;c:integer);
begin
end;
procedure drawgerade(a,b,c:double;col:integer);
begin
end;
procedure tracea(a:array of const);
begin
end;
procedure trace(a:string);
begin
end;

function sign(x:double):integer;
begin
if x<0 then sign:=-1
else if x>0 then sign:=1
else sign:=0;
end;

var epseq:double=0.00001;
function eq(a,b:double):boolean;
begin
eq:=abs(a-b)<epseq;
end;

function eq3(a1,b1,c1,a2,b2,c2:double):boolean;
begin
epseq:=(abs(a1)+abs(b1)+abs(c1)+abs(a2)+abs(b2)+abs(c2))/1e8;
eq3:=eq(a1,a2)and eq(b1,b2)and eq(c1,c2);
end;

function newloesung(x,y,r:double):integer;
{newloesung>0 neue loesung
           =-1 alte loesung
}
var i:integer;
 sx,sy,sr,d,xl,yl,rl:double;
begin
for i:=0 to loesung-1 do begin
 xl:=xloesung[i];yl:=yloesung[i];rl:=rloesung[i];
 sx:=abs(xl)+abs(x);
 sy:=abs(yl)+abs(y);
 sr:=abs(rl)+abs(r);
 if sx<tol then sx:=1;
 if sy<tol then sy:=1;
 if sr<tol then sr:=1;
 d:=abs(xl-x)/sx+abs(yl-y)/sy+abs(rl-r)/sr;
 if d<0.001 then begin
  newloesung:=-1;exit;
 end;
end;
xloesung[loesung]:=x;yloesung[loesung]:=y;rloesung[loesung]:=r;newloesung:=loesung;
inc(loesung);
end;

function atan2(y,x:double):double;
begin
if x=0 then begin
 if y>=0 then result:=pi/2 else result:=-pi/2
end else result:=arctan2(y,x);
end;

function arccos(x:double):double;
{arccos von x ergebnis in bogenwinkel
-1<=x<=1
}
begin
arccos:=atan2(sqrt(1-x*x),x);
end;

function amax(x,y:double):double;
begin
if x>=y then amax:=x else amax:=y;
end;

function amin(x,y:double):double;
begin
if x<=y then amin:=x else amin:=y;
end;

procedure rotate90p(xp,yp,xc,yc:double;var x,y:double);
{xp,yp wird um xc,yc um 90 Grad gedreht}
var dx,dy:double;
begin
dx:=xp-xc;dy:=yp-yc;
x:=xc-dy;y:=yc+dx;
end;

procedure rotateminus90p(xp,yp,xc,yc:double;var x,y:double);
{xp,yp wird um xc,yc um -90 Grad gedreht}
var dx,dy:double;
begin
dx:=xp-xc;dy:=yp-yc;
x:=xc+dy;y:=yc-dx;
end;

procedure rotatep(xp,yp,xc,yc,alpha:double;var x,y:double);
{xp,yp wird um xc,yc um alpha Grad gedreht}
var sina,cosa,dx,dy:double;
begin
alpha:=alpha*radiant;
sina:=sin(alpha);
cosa:=cos(alpha);
dx:=xp-xc;
dy:=yp-yc;
x:=xc+dx*cosa-dy*sina;
y:=yc+dx*sina+dy*cosa;
end;
(*
procedure line2p(x1,y1,x2,y2:double;var a,b,c:double;var fail:boolean);
{Umwandlung Linie mit 2 Punkten in Hesse Normalform
 Es wird die Richtung der geraden berücksichtigt
 a cos der Normalen bzw. Richtung+90
 b sin der Normalen bzw. Richtung+90
}
var dx,dy,d,a1,b1:double;
begin
  fail:=false;dx:=x2-x1;dy:=y2-y1;d:=dx*dx+dy*dy;
  if d<tol then begin
   fail:=true; a:=0; b:=0; c:=0; exit;
  end;
  d:=sqrt(d);a1:=dy/d;b1:=dx/d;
  a:=-a1;b:=b1;c:=-(x1*a1)+(y1*b1);
end;
*)

procedure line2p(x1,y1,x2,y2:double;var a,b,c:double;var fail:boolean);
{Umwandlung Linie mit 2 Punkten in Hesse Normalform
 Es wird die Richtung der geraden berücksichtigt
 a cos der Normalen bzw. Richtung+90
 b sin der Normalen bzw. Richtung+90
}
var dx,dy,d,a1,b1:double;
begin
fail:=false;dx:=x2-x1;dy:=y2-y1;d:=dx*dx+dy*dy;
if d<tol then begin
 fail:=true; a:=0; b:=0; c:=0; exit;
end;
d:=sqrt(d);a1:=dy/d;b1:=dx/d;
a:=-a1;b:=b1;c:=-(x1*a1)+(y1*b1);
if (c<0)or((c=0)and (a<0))or((c=0)and(a=0)and(b<0))then begin
 a:=-a;b:=-b;c:=-c;
end;
end;

function winkell(a,b:double):double;
begin
winkell:=atan2(-a,b)/radiant;
end;

procedure pointcl(x2,y2,r2,a1,b1,c1:double;mode:integer;var x,y:double;
                  var fail:boolean);
{Schnittpunkt linie a1,b1,c1 und Kreis x2,y2,r2}
var d,t,u:double;
 s3:double;
label 99;
begin
  if mode=0 then s3:=1 else s3:=-1;
  x:=errorvalue; y:=errorvalue;
  fail:=false;t:=c1-x2*a1-y2*b1;u:=r2;
  d:=(u-t)*(u+t);
  (*
  if abs(d)>4*(abs(u)+abs(t)){-tol}then begin
   if d<0 then begin
    fail:=true;exit
   end else d:=0;
  end else;
  *)
  if (abs(d)<tol)and(mode>0)then goto 99;
  if d<tol then begin
   if d>-tol then d:=0
   else begin
    99:fail:=true;exit;
   end;
  end;
  d:=s3*sqrt(d);
  x:=x2+a1*t+b1*d;
  y:=y2+b1*t-a1*d;
end;

procedure pointcl2p(x2,y2,r2,x1,y1,xx1,yy1:double;mode:integer;var x,y:double;
                  var fail:boolean);
{Schnittpunkt linie x1,y1,xx2,yy2 und Kreis x2,y2,r2}
var a1,b1,c1:double;
begin
x:=errorvalue; y:=errorvalue;
line2p(x1,y1,xx1,yy1,a1,b1,c1,fail);
if not fail then pointcl(x2,y2,r2,a1,b1,c1,mode,x,y,fail);
end;

procedure pointcg(x2,y2,r2,c,alpha:double;mode:integer;var x,y:double;
                  var fail:boolean);
{Schnittpunkt gerade c,alpha und Kreis x2,y2,r2}
var a1,b1,c1:double;
begin
gerade(c,alpha,a1,b1,c1);
pointcl(x2,y2,r2,a1,b1,c1,mode,x,y,fail);
end;

procedure point2l(a1,b1,c1,a2,b2,c2:double;var x,y:double;var fail:boolean);
{Schnittpunkt Linien a1,b1,c1 und a2,b2,c2}
var d:double;
begin
fail:=false;d:=a1*b2-a2*b1;
if abs(d)<0.0001 then begin fail:=true; exit;end;
x:=(b2*c1-b1*c2)/d;y:=(c2*a1-c1*a2)/d;
end;

procedure pointl2pl2p(x1,y1,x2,y2,x3,y3,x4,y4:double;var x,y:double;var fail:boolean);
{Schnittpunkt Linien x1,y1,x2,y2 und x3,y3,x4,y4}
var a1,b1,c1,a2,b2,c2:double;
begin
line2p(x1,y1,x2,y2,a1,b1,c1,fail);
if fail then exit;
line2p(x3,y3,x4,y4,a2,b2,c2,fail);
if fail then exit;
point2l(a1,b1,c1,a2,b2,c2,x,y,fail);
end;

procedure point2g(c1,alpha1,c2,alpha2:double;var x,y:double;var fail:boolean);
{Schnittpunkt gerade c1,alpha1 und c2,alpha2}
var a1,b1,cc1,a2,b2,cc2:double;
begin
gerade(c1,alpha1,a1,b1,cc1);
gerade(c2,alpha2,a2,b2,cc2);
point2l(a1,b1,cc1,a2,b2,cc2,x,y,fail);
end;

procedure point2c(x1,y1,r1,x2,y2,r2:double;mode:integer;var x,y:double;
                   var fail:boolean);
{Schnittpunkt Kreis x1,y1,r1 und Kreis x2,y2,r2}
label 99;
var a1,b1,d,d1,t,u:double;
 s3:integer;
 tol:double;
begin
  if mode=0 then s3:=1 else s3:=-1;
  fail:=false;a1:=x2-x1;b1:=y2-y1;d:=sqrt(a1*a1+b1*b1);
  tol:=(d+r1+r2)/10000;
  if d<=tol then goto 99;
  if (abs(d-(r1+r2))<tol)and(mode>0)then goto 99;       // berühren
  if (abs(d-abs(r1-r2))<tol)and(mode>0)then goto 99;
  d1:=((x2+x1)*a1+(y2+y1)*b1-(r2+r1)*(r2-r1))/(2*d);
  a1:=a1/d;b1:=b1/d;t:=d1-x2*a1-y2*b1;u:=r2;
  d:=(u-t)*(u+t);
  if d<-tol*tol then begin
   99: fail:=true; exit;
  end;
  if d<0 then d:=0;
  d:=s3*sqrt(d);x:=x2+a1*t+b1*d;y:=y2+b1*t-a1*d;
end;

procedure angle2l(a1,b1,c1,a2,b2,c2:double;var w1,w2:double);
{Winkel zwischen 2 Geraden in Grad}
var x:double;
begin
x:=arccos(a1*a2+b1*b2);
w1:=arccos(amin(1,amax(0,x)))*180/pi;
w2:=180-w1;
end;

procedure tancc(x1,y1,r1,x2,y2,r2:double;ialter:integer;var x,y:double;var fail:boolean);
{get one tangent point of two circles
 ialter=0,1,2,3 --> point on the second circle;
 ialter=4,5,6,7 --> point on the first circle
}
//const tol=0.001;
var
si,s1,s2:integer;
xx1,xx2,yy1,yy2,xx3,yy3,rr1,rr2,dx,dy,d,dd,dr,t,tt1,a1,b1,a,b,c:double;
change:boolean;
fact,my,d1,d2:double;
begin
fail:=false; si:=ialter mod 8; x:=errorvalue; y:=errorvalue;
if si<4 then begin
 xx1:=x1;yy1:=y1;rr1:=r1;
 xx2:=x2;yy2:=y2;rr2:=r2;
end else begin
 xx1:=x2;yy1:=y2;rr1:=r2;
 xx2:=x1;yy2:=y1;rr2:=r1;
end;
dx:=xx2-xx1;dy:=yy2-yy1; dd:=dx*dx+dy*dy;
if dd<tol then begin fail:=true;{keinschnitt:=true;}exit; end;
dd:=sqrt(dd);a1:=dy/dd;b1:=-dx/dd;
s1:=((si and 2)shr 1)*2-1;s2:=(si and 1)*2-1;
dr:=s2*rr2-s1*rr1;t:=dr/dd;d:=(1-t)*(1+t);
if d<-tol then begin fail:=true; exit;end;
if d<0 then d:=0; d:=sqrt(d);
a:=a1*d+b1*t;b:=b1*d-a1*t;c:=xx2*a+yy2*b+r2*s2;
tt1:=c-xx1*a-yy1*b;
x:=xx1+a*tt1;y:=yy1+b*tt1;
end;

procedure tangentecc(xx1,yy1,rr1,xx2,yy2,rr2:double;si:integer;
                  var x1,y1,x2,y2:double;var fail:boolean);
{Tangentiale Linie zwischen 2 Kreisen}
var fail1,fail2:boolean;
begin
fail:=false;fail1:=false;fail2:=false;
if si = 0 then begin
 tancc(xx1,yy1,rr1,xx2,yy2,rr2,0,x1,y1,fail1);
 tancc(xx2,yy2,rr2,xx1,yy1,rr1,3,x2,y2,fail2);
end else if si = 1 then begin
 tancc(xx1,yy1,rr1,xx2,yy2,rr2,3,x1,y1,fail1);
 tancc(xx2,yy2,rr2,xx1,yy1,rr1,0,x2,y2,fail2);
end else if si = 2 then begin
 tancc(xx1,yy1,rr1,xx2,yy2,rr2,1,x1,y1,fail1);
 tancc(xx2,yy2,rr2,xx1,yy1,rr1,1,x2,y2,fail2);
end else if si = 3 then begin
 tancc(xx1,yy1,rr1,xx2,yy2,rr2,2,x1,y1,fail1);
 tancc(xx2,yy2,rr2,xx1,yy1,rr1,2,x2,y2,fail2);
end;
fail:=fail1 or fail2;
end;

procedure tangentepc(x1,y1,x2,y2,r2:double;ialter:integer;var x,y:double;
                  var fail:boolean) ;
{tangent point x,y from a given point x1,y1 onto a circle x2,y2,r2}
begin
 fail:=false;
 x:=errorvalue; y:=errorvalue;
 tancc(x2,y2,r2,x1,y1,0,ialter*2,x,y,fail);
end;
procedure tangentlinepc(x1,y1,x2,y2,r2:double;ialter:integer;var a,b,c:double;
                  var fail:boolean) ;
var x,y:double;
begin
tangentepc(x1,y1,x2,y2,r2,ialter,x,y,fail);
if not fail then line2p(x1,y1,x,y,a,b,c,fail);
end;

function getphi(xa,ya,xi,yi:double):double;
{winkel von xa,ya nach xi,yi in grad}
var kx,ky:double;
    phi:double;
begin
kx:=xi-xa;ky:=yi-ya;
if ky=0.0 then  phi:=0.0
else if kx=0.0 then begin
 if ky>0.0 then phi:=90.0 else phi:=270.0
end else PHI:=ARCTAN(ky/kx)/RADIANT;
if kx<0.0 then phi:=180.0+phi;
if phi<0.0 then phi:=360.0+phi;
getphi:=phi;
end;

procedure parallel2p(xx1,yy1,xx2,yy2,d:double;lalter:integer;var x1,y1,x2,y2:double);
{get a parallel line of the given line with a given distance d}
var
phi,a:double;
si:integer;
begin
phi:=getphi(xx1,yy1,xx2,yy2);
if phi >180 then phi:=phi-180;
a:=pi*phi/180;
if lalter=0 then si:=-1 else si:=1;
x1:=xx1-si*d*sin(a);y1:=yy1+si*d*cos(a);
x2:=xx2-si*d*sin(a);y2:=yy2+si*d*cos(a);
end;

function distance2p(x1,y1,x2,y2:double):double;
var d,dx,dy:double;
begin
dx:=x1-x2;dy:=y2-y1;
d:=sqrt(dx*dx+dy*dy);
distance2p:=d;
end;

function distancepl(x1,y1,a2,b2,c2:double):double;
{Abstand zwischen Punkt x1,y1 und der linie a2,b2,c2
 mit Vorzeichen
 f�r normierte a2,b2 ist dies einfach a2*x1+b2*y1+c2
}
var q,r:double;
begin
q:=a2*a2+b2*b2;
if q<eps then distancepl:=0
else begin
 r:=a2*x1+b2*y1-c2;
 distancepl:=r/sqrt(q);
end;
end;
function distancepl2p(x1,y1,x2,y2,x3,y3:double):double;
{Abstand zwischen Punkt x1,y1 und der linie x2,y2,x3,y3 }
var a,b,c:double;fail:boolean;
begin
line2p(x2,y2,x3,y3,a,b,c,fail);
result:=distancepl(x1,y1,a,b,c);
end;
procedure lotpl(x1,y1,a2,b2,c2:double;var x,y:double);
{Lot vom punkt x1,y1 auf die Linie a2,b2,c2
gel”st wird a2*x+b2*y+c2=0 und b2*x-a2*y+(a2*y1-b2*x1)=0
}
var x2,y2,d:double;
begin
{
d:=a2*a2+b2*b2;
x:=-(a2*c2 - b2*b2*x1 + a2*b2*y1)/d;
y:=-(b2*c2 + a2*b2*x1 - a2*a2*y1)/d;
}
{alt wohl nur richtig}
x2:=a2*c2;y2:=b2*c2;
d:=(x1-x2)*a2+(y1-y2)*b2;
x:=x1-d*a2;y:=y1-d*b2;

end;

procedure rotatel(al,bl,cl,xc,yc,alpha:double;var a,b,c:double);
{al,bl,cl wird um xc,yc mit alpha gedreht}
var w,xx,yy:double;
begin
w:=winkell(al,bl);
lotpl(xc,yc,al,bl,cl,xx,yy);
w:=w+alpha;
rotatep(xx,yy,xc,yc,alpha,xx,yy);
a:=sin(w*radiant);b:=cos(w*radiant);
c:=xx*a+yy*b;
end;

procedure winkelhalb(a1,b1,c1,a2,b2,c2:double;mode:integer;
                     var a,b,c:double;var fail:boolean);
var alpha,x,y:double;
begin
point2l(a1,b1,c1,a2,b2,c2,x,y,fail);
if fail then exit;
alpha:=(winkell(a1,b1)+winkell(a2,b2))/2;
if mode=1 then alpha:=alpha+90;
if alpha>=360 then alpha:=alpha-360;
a:=-sin(alpha*radiant);
b:=cos(alpha*radiant);
c:=a*x+b*y;
end;

function distancepc(x1,y1,x2,y2,r2:double;s:integer):double;
{Abstand zwischen Punkt x1,y1 und dem Kreis x2,y2,r2 }
var d:double;
var dx,dy:double;
begin
dx:=x1-x2;dy:=y2-y1;
d:=sqrt(dx*dx+dy*dy);
if s<>0 then distancepc:=d+r2
else distancepc:=abs(d-r2);
end;

procedure lotpc(x1,y1,x2,y2,r2:double;mode:integer;var x,y:double);
{Lot vom punkt x1,y1 auf den Kreis x2,y2,r2}
var fail:boolean;
begin
pointcl2p(x2,y2,r2,x1,y1,x2,y2,mode,x,y,fail);
end;

procedure circle2lr(a1,b1,c1,a2,b2,c2,r:double;mode:integer;
                     var x,y:double;var fail:boolean);
{Mittelpunkt eines Kreises tangential zu 2 Linien und mit Radius r}
var d,d1,d2:double;
 si,s1,s2:integer;
begin
loesung:=0;
d:=a1*b2-a2*b1;
if abs(d)<delta*delta then fail:=true
else begin
 fail:=false;
 for si:=0 to 3 do begin
  s1:=(si shr 1)mod 2*2-1;s2:=si mod 2*2-1;
  d1:=c1+s1*r;d2:=c2+s2*r;
  x:=(b2*d1-b1*d2)/d;y:=(d2*a1-d1*a2)/d;
  if newloesung(x,y,r)=mode then exit;
 end;
end;
end;

procedure circle2cr(x1,y1,r1,x2,y2,r2,r:double;mode:integer;
                     var x,y:double;var fail:boolean);
{Mittelpunkt eines Kreises tangential zu 2 Kreisen und mit radius r}
var ddd,aa1,bb1,d,a1,b1,d1,u,dd,t:double;
 si,s1,s2,s3:integer;
begin
fail:=false;
aa1:=x2-x1;bb1:=y2-y1;d:=sqrt(aa1*aa1+bb1*bb1);
if d<tol then fail:=true else begin
 a1:=aa1/d;b1:=bb1/d;
 loesung:=0;
 for si:=0 to 7 do begin
  s1:=(si shr 2)mod 2*2-1;s2:=(si shr 1)mod 2*2-1;s3:=si mod 2*2-1;
  d1:=((x2+x1)*aa1+(y2+y1)*bb1-(2*r+s1*r1+s2*r2)*(s2*r2-s1*r1))/(2*d);
  t:=d1-x2*a1-y2*b1;u:=r2+s2*r;dd:=(u-t)*(u+t);
  if abs(dd)>0.001*(abs(u)+abs(t))then begin
   if dd<0 then continue;
  end else dd:=0;
  ddd:=s3*sqrt(dd);
  x:=x2+a1*t+b1*ddd;y:=y2+b1*t-a1*ddd;
  if newloesung(x,y,0)=mode then exit;
 end;
end;
fail:=true;
end;

procedure xcircle2cr(x1,y1,r1,x2,y2,r2,r:double;mode:integer;
                     var x,y:double;var fail:boolean);
{Mittelpunkt eines Kreises tangential zu 2 Kreisen und mit radius r
 als Hilfsroutinen beim Nachiterieren von xcirclexx0 wenn ein Term durch den
 dividiert wird Null wird.
 Es wird dann tol eingesetzt und das Ergebnis hat einen kleinen Fehler
 Der Radius wird dann solange iteriert bis der Fehler verschwindet
}
var ddd,aa1,bb1,d,a1,b1,d1,u,dd,t:double;
 si,s1,s2,s3:integer;
begin
fail:=true;
aa1:=x2-x1;bb1:=y2-y1;d:=sqrt(aa1*aa1+bb1*bb1);
if d<tol then exit else begin
 a1:=aa1/d;b1:=bb1/d;
 si:=mode;
 s1:=(si shr 2)mod 2*2-1;s2:=(si shr 1)mod 2*2-1;s3:=si mod 2*2-1;
 d1:=((x2+x1)*aa1+(y2+y1)*bb1-(2*r+s1*r1+s2*r2)*(s2*r2-s1*r1))/(2*d);
 t:=d1-x2*a1-y2*b1;u:=r2+s2*r;dd:=(u-t)*(u+t);
 if abs(dd)>0.001*(abs(u)+abs(t))then begin
  if dd<0 then exit;
 end else dd:=0;
 ddd:=s3*sqrt(dd);
 x:=x2+a1*t+b1*ddd;y:=y2+b1*t-a1*ddd;
 fail:=false;
end;
end;

procedure circleclr(x1,y1,r1,a2,b2,c2,r:double;mode:integer;
                     var x,y:double;var fail:boolean);
{Mittelpunkt eines Kreises tangential zu Linie und Kreis und mit radius r}
var si,s1,s2,s3:integer;
 d1,t,u,d:double;
begin
loesung:=0;fail:=false;
for si:=0 to 7 do begin
 s1:=(si shr 2)mod 2*2-1;s2:=(si shr 1)mod 2*2-1;s3:=si mod 2*2-1;
 d1:=c2+s1*r;t:=d1-x1*a2-y1*b2; u:=r1+s2*r;d:=(u-t)*(u+t);
 if d<-abs(u+t) then continue;
 if d<0 then d:=0;d:=s3*sqrt(d);
 x:=x1+a2*t+b2*d;y:=y1+b2*t-a2*d;
 if newloesung(x,y,r)=mode then exit;
end;
fail:=true;
end;

procedure circlelcr(a1,b1,c1,x2,y2,r2,r:double;mode:integer;
                     var x,y:double;var fail:boolean);
{Mittelpunkt eines Kreises tangential zu Linie und Kreis und mit radius r}
var si,s1,s2,s3:integer;
 d1,t,u,d:double;
begin
loesung:=0;fail:=false;
for si:=0 to 7 do begin
 s1:=(si shr 2)mod 2*2-1;s2:=(si shr 1)mod 2*2-1;s3:=si mod 2*2-1;
 d1:=c1+s1*r;t:=d1-x2*a1-y2*b1; u:=r2+s2*r;d:=(u-t)*(u+t);
 if d<-abs(u+t) then continue;
 if d<0 then d:=0;d:=s3*sqrt(d);
 x:=x2+a1*t+b1*d;y:=y2+b1*t-a1*d;
 if newloesung(x,y,r)=mode then exit;
end;
fail:=true;
end;
procedure circle3p(x1,y1,x2,y2,x3,y3:double;var x,y,r:double;var fail:boolean);
var mx1,my1,mx2,my2,dx1,dx2,dy1,dy2:double;
 dx,dy:double;
begin
mx1:=(x1+x2)/2;my1:=(y1+y2)/2;mx2:=(x2+x3)/2;my2:=(y2+y3)/2;
dx1:=x2-x1;dy1:=y2-y1;dx2:=x3-x2;dy2:=y3-y2;
pointl2pl2p(mx1,my1,mx1-dy1,my1+dx1,mx2,my2,mx2-dy2,my2+dx2,x,y,fail);
if not fail then begin
 dx:=x-x1;dy:=y-y1;
 r:=sqrt(dx*dx+dy*dy);
 loesung:=1;
 xloesung[0]:=x;yloesung[0]:=y;rloesung[0]:=r;
end else begin
 r:=errorvalue;
 loesung:=0;
end;
end;
procedure xcircle3p(x1,y1,x2,y2,x3,y3:double;var x,y,r:double;var fail:boolean);
// wie oben aber ohne newloesung
var mx1,my1,mx2,my2,dx1,dx2,dy1,dy2:double;
dx,dy:double;
begin
mx1:=(x1+x2)/2;my1:=(y1+y2)/2;mx2:=(x2+x3)/2;my2:=(y2+y3)/2;
dx1:=x2-x1;dy1:=y2-y1;dx2:=x3-x2;dy2:=y3-y2;
pointl2pl2p(mx1,my1,mx1-dy1,my1+dx1,mx2,my2,mx2-dy2,my2+dx2,x,y,fail);
if not fail then begin
 dx:=x-x1;dy:=y-y1;
 r:=sqrt(dx*dx+dy*dy);
end else r:=errorvalue;
end;


procedure circle3s(xa1,ya1,x1,y1,xa2,ya2,x2,y2,xa3,ya3,x3,y3:double;mode:integer;
                    var x,y,r:double;var fail:boolean);
var f:double;
    a,b,c:array[0..2] of double;
    w:array[0..1,0..1,1..3] of double;
    si,i,j,k,s:integer;
    fail0,fail1,fail2:boolean;
begin
loesung:=0;
fail:=true; fail0:=false; fail1:=false; fail2:=false;
line2p(xa1,ya1,x1,y1,a[0],b[0],c[0],fail0);
line2p(xa2,ya2,x2,y2,a[1],b[1],c[1],fail1);
line2p(xa3,ya3,x3,y3,a[2],b[2],c[2],fail2);
if fail0 or fail1 or fail2 then exit;
for i:=0 to 1 do begin
 j:=(i+2) mod 3;
 for k:=0 to 1 do begin
  s:=2*k-1;
  f:=sqrt(sqr(a[i]+s*a[j])+sqr(b[i]+s*b[j]));
  w[i,k,1]:=(a[i]+s*a[j])/f;
  w[i,k,2]:=(b[i]+s*b[j])/f;
  w[i,k,3]:=(c[i]+s*c[j])/f;
 end;
end;
for si:=0 to 3 do begin
 i:=si and 1; j:=si shr 1;
 point2l(w[1,i,1],w[1,i,2],w[1,i,3],w[0,j,1],w[0,j,2],w[0,j,3],x,y,fail);
 if not fail then begin
  r:=abs(a[0]*x+b[0]*y-c[0]);
  if newloesung(x,y,r)=mode then exit;
 end;
end;
fail:=true;
end;

procedure circle3l(a1,b1,c1,a2,b2,c2,a3,b3,c3:double;mode:integer;
                    var x,y,r:double;var fail:boolean);
var f:double;
    a,b,c:array[0..2] of double;
    w:array[0..1,0..1,1..3] of double;
    si,i,j,k,s:integer;
begin
loesung:=0;
fail:=true;
a[0]:=a1;b[0]:=b1;c[0]:=c1;
a[1]:=a2;b[1]:=b2;c[1]:=c2;
a[2]:=a3;b[2]:=b3;c[2]:=c3;
for i:=0 to 1 do begin
 j:=(i+2) mod 3;
 for k:=0 to 1 do begin
  s:=2*k-1;
  f:=sqrt(sqr(a[i]+s*a[j])+sqr(b[i]+s*b[j]));
  w[i,k,1]:=(a[i]+s*a[j])/f;
  w[i,k,2]:=(b[i]+s*b[j])/f;
  w[i,k,3]:=(c[i]+s*c[j])/f;
 end;
end;
for si:=0 to 3 do begin
 i:=si and 1; j:=si shr 1;
 point2l(w[1,i,1],w[1,i,2],w[1,i,3],w[0,j,1],w[0,j,2],w[0,j,3],x,y,fail);
 if not fail then begin
  r:=abs(a[0]*x+b[0]*y-c[0]);
  if newloesung(x,y,r)=mode then exit;
 end;
end;
fail:=true;
end;

procedure change(var x,y:double);
var z:double;
begin
z:=x;x:=y;y:=z;
end;

function flaeche(var punkte:tpunkte;n:integer):double;
var i,j:integer;
 f:double;
begin
f:=0;
for i:=1 to n do begin
 j:=i+1;if j>n then j:=1;
 f:=f+(punkte[j].x-punkte[i].x)*(punkte[j].y+punkte[i].y)/2;
end;
flaeche:=f;
end;

procedure schwerpunkt(var punkte:tpunkte;n:integer;
                      var flaeche:double;var xschwer,yschwer:double);
{xschwer,yschwer schwerpunkt}
var i,j:integer;
 f,fe,mx,my:double;
begin
f:=0;mx:=0;
for i:=1 to n do begin
 j:=i+1;if j>n then j:=1;
 fe:=(punkte[j].x-punkte[i].x)*(punkte[j].y+punkte[i].y)/2;
 f:=f+fe;
 mx:=mx+fe*((punkte[j].x+punkte[i].x)/2);
end;
flaeche:=abs(f);
xschwer:=mx/f;
f:=0;my:=0;
for i:=1 to n do begin
 j:=i+1;if j>n then j:=1;
 fe:=(punkte[j].y-punkte[i].y)*(punkte[j].x+punkte[i].x)/2;
 f:=f+fe;
 my:=my+fe*((punkte[j].y+punkte[i].y)/2);
end;
yschwer:=my/f;
end;

procedure getkoord(x0,y0,xu,yu,xv,yv,xp,yp:double;var u,v,unorm,vnorm:double);
{von den vektoren 0-u und 0-v wird ein Koordinatensystem aufgespannt;
 xp,yp ist ein punkt im kartesischen koordinaten
 u,v sind die koordinaten des Punktes im schiefen Koordinaten
 es ist die l”sung des Gleichungssystems
 u*vekxu+v*vekxv=distx
 u*vekyu+v*vekyv=disty
 unorm,vnorm sind dieselben u,v Kordinaten mit der l„nge der 0-u,0-v
 Vektoren
}
var vekxu,vekyu,lenu,vekxv,vekyv,lenv:double;
 det,distx,disty:double;
begin
 vekxu:=xu-x0;vekxv:=xv-x0;
 vekyu:=yu-y0;vekyv:=yv-y0;
 distx:=xp-x0;
 disty:=yp-y0;
 det:=vekxu*vekyv-vekxv*vekyu;
 u:=(distx*vekyv-vekxv*disty)/det;
 v:=(vekxu*disty-distx*vekyu)/det;
 lenu:=sqrt(sqr(vekxu)+sqr(vekyu));
 lenv:=sqrt(sqr(vekxv)+sqr(vekyv));
 unorm:=u*lenu;
 vnorm:=v*lenv;
end;

procedure orthokoord(x0,y0,xu,yu,xp,yp:double;var u,v:double);
{Es wird ein Koordinatensystem mit dem Ursprung x0,y0 und der x-achse xu,yu
 augespannt;
 xp,yp ist ein punkt im kartesischen koordinaten
 u,v sind die koordinaten des Punktes im neuen Koordinaten
 es ist die l”sung des Gleichungssystems
 u*vekxu+v*vekxv=distx
 u*vekyu+v*vekyv=disty
 unorm,vnorm sind dieselben u,v Kordinaten mit der l„nge der 0-u,0-v
 Vektoren
}
var vekxu,vekyu,l,vekxv,vekyv:double;
 det,distx,disty:double;
begin
 vekxu:=xu-x0;
 vekyu:=yu-y0;
 vekxv:=-vekyu;
 vekyv:=vekxu;
 distx:=xp-x0;
 disty:=yp-y0;
 det:=vekxu*vekyv-vekxv*vekyu;
 if det<eps then begin
  det:=det;
 end;
 u:=(distx*vekyv-vekxv*disty)/det;
 v:=(vekxu*disty-distx*vekyu)/det;
 l:=sqrt(sqr(vekxu)+sqr(vekyu));
 u:=u*l;
 v:=v*l;
end;

procedure getpolar(xp,yp,xc,yc:double;var alpha,r:double);
begin
r:=distance2p(xc,yc,xp,yp);
alpha:=getphi(xc,yc,xp,yp);
end;

procedure polarpoint(xc,yc,alpha,d:double;var x,y:double);
begin
x:=xc+d*cos(alpha*radiant);
y:=yc+d*sin(alpha*radiant);
end;

function alpha3p(x0,y0,x1,y1,x2,y2:double):double;
var xx,yy,a:double;
begin
orthokoord(x0,y0,x1,y1,x2,y2,xx,yy);
a:=abs(getphi(0,0,xx,yy));
if a>180 then a:=360-a;
if a=180 then
a:=a;
alpha3p:=a;
end;

function alpha3pdir(x0,y0,x1,y1,x2,y2:double):double;
var xx,yy,a:double;
begin
orthokoord(x0,y0,x1,y1,x2,y2,xx,yy);
a:=getphi(0,0,xx,yy);
alpha3pdir:=a;
end;

function test1(a,b,c:double):boolean;
{Test der Lösung a-b oder a-c müssen 0 sein}
begin
result:=(abs(a-b)<10*tol)or(abs(a-c)<10*tol);
if not result then begin
 a:=a;
end;
end;
function test2(a,b,c:double):boolean;
{Test der Lösung a-b oder a-c müssen 0 sein}
begin
result:=(abs(a-b)<sqr(10*tol))or(abs(a-c)<sqr(10*tol));
if not result then begin
 a:=a;
end;
end;

procedure ctransform(isline:boolean;x,y,r,xc,yc,rc:double;
                     var istline:boolean;var xt,yt,rt:double);
var a:double absolute x;
 b:double absolute y;
 c:double absolute r;
 at:double absolute xt;
 bt:double absolute yt;
 ct:double absolute rt;
 x1,y1,x2,y2,x3,y3,alpha,d:double;
 fail:boolean;

 procedure transp(x,y:double;var xt,yt:double);
 {Circle point transfomation}
 var d,alpha:double;
 begin
 getpolar(x,y,xc,yc,alpha,d);
 polarpoint(xc,yc,alpha,rc*rc/d,xt,yt);
 end;

begin
if isline then begin
 if abs(distancepl(xc,yc,a,b,c))<tol then begin
  {Line stays as it is}
  at:=a;bt:=b;ct:=c;istline:=true;
 end else begin
  {Line becomes circle}
  lotpl(xc,yc,a,b,c,x1,y1);
  getpolar(x1,y1,xc,yc,alpha,d);
  rt:=sqr(rc)/d/2;
  polarpoint(xc,yc,alpha,rt,xt,yt);istline:=false;
 end;
end else begin
 if abs(distance2p(x,y,xc,yc)-r)<tol then begin
  {circle becomes line}
  alpha:=getphi(xc,yc,x,y);
  polarpoint(xc,yc,alpha,sqr(rc)/(2*r),x1,y1);
  polarpoint(x1,y1,alpha+90,r,x2,y2);
  line2p(x1,y1,x2,y2,at,bt,ct,fail);
  istline:=true;
 end else begin
  {circle to circle}
  polarpoint(x,y,0,r,x1,y1);
  polarpoint(x,y,90,r,x2,y2);
  polarpoint(x,y,180,r,x3,y3);
  transp(x1,y1,x1,y1);
  transp(x2,y2,x2,y2);
  transp(x3,y3,x3,y3);
  xcircle3p(x1,y1,x2,y2,x3,y3,xt,yt,rt,fail);
  istline:=false;
 end;
end;
end;

procedure xcirclecpar(x1,y1,r1,a2,b2,c2,c3:double;mode:integer;var x,y,r:double;
                     var fail:boolean);
var rr,cmitte:double;
 mode2:integer;
begin
fail:=true;
cmitte:=(c2+c3)/2;
mode2:=mode div 2;
r:=abs(c3-c2)/2;

if not odd(mode) then rr:=r+r1 else rr:=abs(r-r1);
if rr<0 then exit;
if cmitte<0 then begin
 a2:=-a2;b2:=-b2;cmitte:=-cmitte;
end;
pointcl(x1,y1,rr,a2,b2,cmitte,mode2,x,y,fail);
end;

procedure xcirclec2l(alpha,x1,y1,r1:double;mode:integer;var x,y,r:double;
                     var fail:boolean);
{ Bestimmt die Kreise an 2 symmetrische Linien durch den Nullpunkt,
 die einen Winkel von +/- alpha zur X-Achse haben
Die Lösungen 0 bis 3 werden erzielt durch folgendes Gleichungssystem:
(x-x1)^2+y1^2=(r+/-r1)^2
r=x*sin(alpha)
Sie liegen auf der X Achse
Die Lösungen 4 bis 7 durch das Gleichungssystem:
(y-y1)^2+x1^2=(r+/-r1)^2
r=y*cos(alpha)
Sie liegen auf der y Achse
}
var ff,a,b,c,sina,sina2,cosa,cosa2,x12,y12,r12:double;
 s1,s2:integer;
begin
alpha:=alpha*radiant;
sina:=sin(alpha);
sina2:=sqr(sina);
cosa:=cos(alpha);
cosa2:=sqr(cosa);
r12:=sqr(r1);
x12:=sqr(x1);
y12:=sqr(y1);
{
drawgerade(sina,cosa,0,7);
drawgerade(-sina,cosa,0,7);
drawkreis(x1,y1,r1,7);
}
if mode<4 then begin
 a:=(x12+y12-r12)/(sina2-1);
 if mode<=1 then b:=(sina*r1+x1)/(sina2-1)
 else b:=(-sina*r1+x1)/(sina2-1);
 c:=a+sqr(b);
 if c<0 then exit;
 case mode of
  0,2:x:=sqrt(c)-b;
  1,3:x:=-sqrt(c)-b;
  {
  2:x:=sqrt(c)+b;
  3:x:=-sqrt(c)+b
  }
  else exit;
 end;
 y:=0;
 r:=x*sina;
 {
 s1:=(mode mod 2)*2-1;
 s2:=(mode div 2)*2-1;
 fail:=true;
 ff:=2*x1*r1*s1*sina-y12+r12+sina2*x12+sina2*y12;
 if ff<0 then exit;
 x:=0.5*(2*x1+2*r1*s1*sina+s2*2*sqrt(ff))/(1-sina2);
 r:=abs(x)*sina;
 }
 // testen
 fail:=not test2(sqr(x-x1)+sqr(y-y1),sqr(r+r1),sqr(r-r1));
 //drawkreis(x,y,r,4+ord(fail));
end else begin
 mode:=mode-4;
 a:=(x12+y12-r12)/(cosa2-1);
 if mode<=1 then b:=(cosa*r1+y1)/(cosa2-1)
 else b:=(-cosa*r1+y1)/(cosa2-1);
 c:=a+sqr(b);
 if c<0 then exit;
 case mode of
  0,2:y:=sqrt(c)-b;
  1,3:y:=-sqrt(c)-b;
  {
  2:y:=sqrt(c)+b;
  3:y:=-sqrt(c)+b
  }
  else exit;
 end;
 r:=y*cosa;
 x:=0;
 {
 s1:=(mode mod 2)*2-1;
 s2:=(mode div 2)*2-1;
 fail:=true;
 ff:=2*x1*r1*s1*sina-y12+r12+sina2*x12+sina2*y12;
 if ff<0 then exit;
 x:=0.5*(2*x1+2*r1*s1*sina+s2*2*sqrt(ff))/(1-sina2);
 r:=abs(x)*sina;
 }
 // testen
 fail:=not test2(sqr(x-x1)+sqr(y-y1),sqr(r+r1),sqr(r-r1));
 //drawkreis(x,y,r,4+ord(fail));
end;
end;

procedure circlec2l(x1,y1,r1,a2,b2,c2,a3,b3,c3:double;mode:integer;
                    var x,y,r:double;var fail:boolean);
var mode1,mode2:integer;
 d,xx,yy,xx1,yy1,xx2,yy2,xm,ym,xs,ys,alpha,al,xxx,yyy:double;
 x11,y11:double;
 i,k,m:integer;
 phi1,phi2:double;
label 1,2;
begin
loesung:=0;
if abs(a2-a3)+abs(b2-b3)<delta then begin //Parallel
 for i:=0 to 3 do begin
  xcirclecpar(x1,y1,r1,a2,b2,c2,c3,i,x,y,r,fail);
  if fail then continue;
  if newloesung(x,y,r)=mode then exit;
 end;
 fail:=true;
end else begin
 point2l(a2,b2,c2,a3,b3,c3,xx,yy,fail);
 if fail then exit;
 phi1:=winkell(a2,b2);
 phi2:=winkell(a3,b3);
 alpha:=abs(phi2-phi1)/2;
 al:=(phi2+phi1)/2;
 x11:=x1-xx;y11:=y1-yy; // Schnittpunkt in den Nullpunkt
 rotatep(x11,y11,0,0,-al,x11,y11);
 //drawkreis(x11,y11,r1,7);
 for k:=0 to 7 do begin
  xcirclec2l(alpha,x11,y11,r1,k,xxx,yyy,r,fail);
  if fail then continue;
  //drawkreis(xxx,yyy,r,5);
  rotatep(xxx,yyy,0,0,al,x,y);    //Ergebnis drehen
  x:=x+xx;y:=y+yy; // und verschieben
  //drawkreis(x,y,r,6);
  if newloesung(x,y,r)= mode then exit
 end;
 fail:=true;
end;
end;

procedure circle2lp(a1,b1,c1,a2,b2,c2,x3,y3:double;mode:integer;
                    var x,y,r:double;var fail:boolean);
{Kreis Tangential an Punkt und 2 Linien}
var aa1,bb1,cc1,aa2,bb2,cc2,xx,yy,xx1,yy1,xx2,yy2,rr,dd,d,f1,f2,f3:double;
begin
//Schnittpunkt
point2l(a1,b1,c1,a2,b2,c2,xx,yy,fail);
//Winkelhalbierende
winkelhalb(a1,b1,c1,a2,b2,c2,0,aa1,bb1,cc1,fail);
winkelhalb(a1,b1,c1,a2,b2,c2,1,aa2,bb2,cc2,fail);
if abs(distancepl(x3,y3,aa2,bb2,cc2))<abs(distancepl(x3,y3,aa1,bb1,cc1))then begin
 aa1:=aa2;bb1:=bb2;cc1:=cc2;
end;
// Testen ob Punkt auf einer der Linien
if (abs(distancepl(x3,y3,a1,b1,c1))<0.1)or(abs(distancepl(x3,y3,a2,b2,c2))<0.1) then begin
 rotate90p(xx,yy,x3,y3,xx1,yy1);
 line2p(x3,y3,xx1,yy1,aa2,bb2,cc2,fail);
 point2l(aa1,bb1,cc1,aa2,bb2,cc2,x,y,fail);
end else begin
 // xx1,yy1 Punkt auf Winkelhalbierender mit abstand von x3,y3
 d:=distance2p(xx,yy,x3,y3);
 pointcl(xx,yy,d,aa1,bb1,cc1,0,xx1,yy1,fail);
 pointcl(xx,yy,d,aa1,bb1,cc1,1,xx2,yy2,fail);
 if distance2p(xx,yy,xx2,yy2)<distance2p(xx,yy,xx1,yy1)then begin
  xx1:=xx2;yy1:=yy2;
 end;
 // Fehler: Radius-Abstand
 if mode=0 then begin // Mittelpunkt zwischen xx1,yy1 und schnittpunkt
  // Schnittpunkt
  f1:=-d;
 end else begin
  rr:=abs(distancepl(xx1,yy1,a1,b1,c1));dd:=distance2p(x3,y3,xx1,yy1);f2:=rr-dd;
  xx:=xx1+xx1-xx;yy:=yy1+yy1-yy;
  rr:=abs(distancepl(xx,yy,a1,b1,c1));dd:=distance2p(x3,y3,xx,yy);f1:=rr-dd;
 end;
 while true do begin
  xx2:=(xx+xx1)/2;yy2:=(yy+xx1)/2;
  rr:=abs(distancepl(xx2,yy2,a1,b1,c1));dd:=distance2p(x3,y3,xx2,yy2);f3:=rr-dd;
  if f3<0.1 then break;
  if f1*f3>0 then begin
   xx:=xx2;yy:=yy2;f1:=f3;
  end else begin
   xx1:=xx2;yy1:=yy2;f2:=f3;
  end;
 end;
 x:=xx2;y:=yy2;
end;
r:=distance2p(x,y,x3,y3);
end;
{
procedure circle2lp(a1,b1,c1,a2,b2,c2,x3,y3:double;mode:integer;
                    var x,y,r:double;var fail:boolean);
begin
circlec2l(x3,y3,0,a1,b1,c1,a2,b2,c2,mode,x,y,r,fail);
end;
}
procedure circlel2p(a1,b1,c1,x2,y2,x3,y3:double;var x,y,r:double;var fail:boolean);
var x0,a,b,c,xx,phi,u,v:double;
 waschange:boolean;
 mode,i:integer;
begin
fail:=true;
loesung:=0;
for i:=0 to 1 do begin
 xcirclel2p(a1,b1,c1,x2,y2,x3,y3,i,x,y,r,fail);
 if not fail then newloesung(x,y,r);
end;
(*
exit;
mode:=pmode;
waschange:=false;
//if mode >1 then exit;
if abs(a1)>=abs(b1)then begin
 x0:=c1/a1;phi:=arctan(b1/-a1)/radiant;
end else begin
 waschange:=true;
 change(a1,b1);change(x2,y2);change(x3,y3);
 if a1<0 then begin
  a1:=-a1;b1:=-b1;c1:=-c1;
 end;
 x0:=c1/a1;phi:=arctan(b1/-a1)/radiant+90;
end;
rotatep(x2,y2,x0,0,phi,x2,y2);
rotatep(x3,y3,x0,0,phi,x3,y3);
if abs(y2)<eps then y2:=0;
if abs(y3)<eps then y3:=0;
if (y3*y2<-eps)or(y3+y2=0)then begin // Verschiedene Seiten, beide auf Linie
 fail:=true;exit;
end;
xx:=x2;
x2:=0;x3:=x3-xx;
{Linie= X Achse Punkt 2 auf y Achse}
if abs(y2-y3)>tol then begin
 a:=-2*x3;b:=-2*y3+2*y2;c:=sqr(x3)+sqr(y3)-sqr(y2);
 v:=-y2*(-y2*sqr(a)+2*c*b+y2*sqr(b));
 if v<-eps then exit
 else if v<0 then v:=0;
 u:=sqrt(v)/b;
 for mode:=0 to 1 do begin
  if mode=0 then x:=-y2*a/b+u else x:=-y2*a/b-u;
  y:=(-c-a*x)/b;r:=y;
  x:=x+xx;
  rotatep(x,y,x0,0,-phi,x,y);
  if waschange then change(x,y);
  if newloesung(x,y,r)=pmode then begin
   fail:=false;exit;
  end;
 end;
end else begin
 r:=(sqr(x2-x3)/4+sqr(y2))/(2*y2);
 x:=(x2+x3)/2;y:=r;
 x:=x+xx;
 rotatep(x,y,x0,0,-phi,x,y);
 if waschange then change(x,y);
 newloesung(x,y,r);
 fail:=pmode<>0;exit;
end;
fail:=true;
*)
end;

procedure xcircleclp(x1,y1,r1,a2,b2,c2,x3,y3:double;mode:integer;
          var x,y,r:double;var fail:boolean);
var d,d1,d2,ff,f1,f2,x0,a,xx,phi:double;
 waschange,wasminus:boolean;
 x32,y32,r12,y12:double;
 z1,z2,z3:double;
begin
fail:=true;waschange:=false;wasminus:=false;
if abs(a2)<abs(b2) then begin
 waschange:=true;
 change(a2,b2);change(x1,y1);change(x3,y3);
end;
x0:=c2/a2;
phi:=atan2(a2,b2)/radiant;
rotatep(x1,y1,x0,0,phi,x1,y1);
rotatep(x3,y3,x0,0,phi,x3,y3);
xx:=x1;
x1:=0;x3:=x3-xx;
{Linie= X Achse Mittelpunkt Kreis auf y Achse}
x32:=sqr(x3);y32:=sqr(y3);r12:=sqr(r1);y12:=sqr(y1);
if abs(y3)<10*tol then begin
 y3:=0;
 x:=x3;
 y:=(r12-y12-x32)/-(2*(y1-r1));
 r:=abs(y);
end else begin
 if y3<0 then begin
  wasminus:=true;
  y1:=-y1;y3:=-y3;
 end;
 case mode of
  0:begin
   if abs(r1 + y1 - y3)<tol then exit;
   if r1 + y1<0 then exit;
   if -Power(r1,2) + Power(x3,2) + Power(y1,2) - 2*y1*y3 + Power(y3,2)<0 then exit;
   x:=(r1*x3 + x3*y1 + Sqrt(r1 + y1)*Sqrt(y3)*Sqrt(-Power(r1,2) + Power(x3,2) + Power(y1,2) - 2*y1*y3 + Power(y3,2)))/(r1 + y1 - y3);
   y:=(Power(x,2) - 2*x*x3 + Power(x3,2) + Power(y3,2))/(2.*y3);
  end;
  1:begin
   if abs(r1 + y1 - y3)<tol then exit;
   if r1 + y1<0 then exit;
   if -Power(r1,2) + Power(x3,2) + Power(y1,2) - 2*y1*y3 + Power(y3,2)<0 then exit;
   if abs(r1 + y1 - y3)<tol then exit;
   x:=(r1*x3 + x3*y1 - Sqrt(r1 + y1)*Sqrt(y3)*Sqrt(-Power(r1,2) + Power(x3,2) + Power(y1,2) - 2*y1*y3 + Power(y3,2)))/(r1 + y1 - y3);
   y:=(Power(x,2) - 2*x*x3 + Power(x3,2) + Power(y3,2))/(2.*y3);
  end;
  2:begin
   if abs(r1 - y1 + y3)<tol then exit;
   if Power(r1,3)*y3 - r1*Power(x3,2)*y3 - Power(r1,2)*y1*y3 + Power(x3,2)*y1*y3 - r1*Power(y1,2)*y3 + Power(y1,3)*y3 + 2*r1*y1*Power(y3,2) -
          2*Power(y1,2)*Power(y3,2) - r1*Power(y3,3) + y1*Power(y3,3)<0 then exit;
   x:=(r1*x3 - x3*y1 + Sqrt(Power(r1,3)*y3 - r1*Power(x3,2)*y3 - Power(r1,2)*y1*y3 + Power(x3,2)*y1*y3 - r1*Power(y1,2)*y3 + Power(y1,3)*y3 + 2*r1*y1*Power(y3,2) -
          2*Power(y1,2)*Power(y3,2) - r1*Power(y3,3) + y1*Power(y3,3)))/(r1 - y1 + y3);
   y:=(Power(x,2) - 2*x*x3 + Power(x3,2) + Power(y3,2))/(2.*y3);
  end;
  3:begin
   if abs(r1 - y1 + y3)<tol then exit;
   if Power(r1,3)*y3 - r1*Power(x3,2)*y3 - Power(r1,2)*y1*y3 + Power(x3,2)*y1*y3 - r1*Power(y1,2)*y3 + Power(y1,3)*y3 + 2*r1*y1*Power(y3,2) -
          2*Power(y1,2)*Power(y3,2) - r1*Power(y3,3) + y1*Power(y3,3)<0 then exit;
   x:=(r1*x3 - x3*y1 - Sqrt(Power(r1,3)*y3 - r1*Power(x3,2)*y3 - Power(r1,2)*y1*y3 + Power(x3,2)*y1*y3 - r1*Power(y1,2)*y3 + Power(y1,3)*y3 + 2*r1*y1*Power(y3,2) -
          2*Power(y1,2)*Power(y3,2) - r1*Power(y3,3) + y1*Power(y3,3)))/(r1 - y1 + y3);
   y:=(Power(x,2) - 2*x*x3 + Power(x3,2) + Power(y3,2))/(2.*y3);
  end;
  else exit;
 end;
 (*
 if mode<=1 then
 ff:=x32*sqr(-2*r1 + 2*y1) -
      4*(r1 - y1 + y3)*(r1*x32 - x32*y1 - r12*y3 + y12*y3 +
              r1*y32 - y1*y32)
 else
 ff:=x32*sqr(-2*r1 - 2*y1) -
      4*(r1 + y1 - y3)*(r1*x32 + x32*y1 + r12*y3 - y12*y3 +
              r1*y32 + y1*y32);
 if ff<0 then ff:=0;//exit;
 f1:=r1 - y1 + y3;
 f2:=r1 + y1 - y3;
 case mode of
 0:
 if abs(f1)<tol then exit else
 x:=(-(x3*(-2*r1 + 2*y1)) +
     sqrt(ff)
    )/(2*(r1 - y1 + y3));
 1:
 if abs(f1)<tol then exit else
 x:=(-(x3*(-2*r1 + 2*y1)) -
     sqrt(ff)
    )/(2*(r1 - y1 + y3));
 2:
 if abs(f2)<tol then exit else
 x:=(-(x3*(-2*r1 - 2*y1)) +
     sqrt(ff)
    )/(2*(r1 + y1 - y3));
 3:
 if abs(f2)<tol then exit else
 x:=(-(x3*(-2*r1 - 2*y1)) -
     sqrt(ff)
    )/(2*(r1 + y1 - y3));
 end;
 a:=sqr(x-x3);
 y:=(a+y32)/(2*y3);
 r:=abs(y);
 *)
 if y<0 then exit;
 r:=y;
 // Testen
 if r>1e9 then exit;
 if mode<=1 then begin
  d1:=sqr(x)+sqr(y1-y)-sqr(r1+r);
  if (abs(d1)>500*tol) then exit;
 end else begin
  d2:=sqr(x)+sqr(y1-y)-sqr(r1-r);
  if (abs(d2)>500*tol) then exit;
 end;
 //if (abs(d1)>500*tol)and(abs(d2)>500*tol) then exit;
end;
if wasminus then y:=-y;
x:=x+xx;
rotatep(x,y,x0,0,-phi,x,y);
if waschange then change(x,y);
fail:=false;
end;
procedure xcirclel2p(a2,b2,c2,x1,y1,x2,y2:double;mode:integer;
          var x,y,r:double;var fail:boolean);
// Kreis Linie und 2 Punkte
var d,d1,d2,ff,f1,f2,x0,a,xx,phi:double;
 waschange,wasminus:boolean;
 x22,y22,y12:double;
 z1,z2,z3:double;
begin
fail:=true;waschange:=false;wasminus:=false;
if abs(a2)<abs(b2) then begin
 waschange:=true;
 change(a2,b2);change(x1,y1);change(x2,y2);
end;
x0:=c2/a2;
phi:=atan2(a2,b2)/radiant;
rotatep(x1,y1,x0,0,phi,x1,y1);
rotatep(x2,y2,x0,0,phi,x2,y2);
xx:=x1;
x1:=0;x2:=x2-xx;
{Linie= X Achse Mittelpunkt Kreis auf y Achse}
x22:=sqr(x2);y22:=sqr(y2);y12:=sqr(y1);
if abs(y2)<10*tol then begin
 y2:=0;
 x:=x2;
 y:=(-y12-x22)/-(2*(y1));
 r:=abs(y);
end else begin
 if y2<0 then begin
  wasminus:=true;
  y1:=-y1;y2:=-y2;
 end;
 case mode mod 2 of
 0:begin
  y:=(Power(x2,2) + (Power(x2,2)*Power(y1,2))/Power(y1 - y2,2) - (2*Power(x2,2)*y1)/(y1 - y2) + Power(y2,2) +
        (2*x2*Power(y1,1.5)*Sqrt(y2)*Sqrt(Power(x2,2) + Power(y1,2) - 2*y1*y2 + Power(y2,2)))/Power(y1 - y2,2) -
        (2*x2*Sqrt(y1)*Sqrt(y2)*Sqrt(Power(x2,2) + Power(y1,2) - 2*y1*y2 + Power(y2,2)))/(y1 - y2) +
        (y1*y2*(Power(x2,2) + Power(y1,2) - 2*y1*y2 + Power(y2,2)))/Power(y1 - y2,2))/(2.*y2);
 x:=(x2*y1 + Sqrt(y1)*Sqrt(y2)*Sqrt(Power(x2,2) + Power(y1,2) - 2*y1*y2 + Power(y2,2)))/(y1 - y2);
 end;
 1:begin
  y:=(Power(x2,2) + (Power(x2,2)*Power(y1,2))/Power(y1 - y2,2) - (2*Power(x2,2)*y1)/(y1 - y2) + Power(y2,2) -
        (2*x2*Power(y1,1.5)*Sqrt(y2)*Sqrt(Power(x2,2) + Power(y1,2) - 2*y1*y2 + Power(y2,2)))/Power(y1 - y2,2) +
        (2*x2*Sqrt(y1)*Sqrt(y2)*Sqrt(Power(x2,2) + Power(y1,2) - 2*y1*y2 + Power(y2,2)))/(y1 - y2) +
        (y1*y2*(Power(x2,2) + Power(y1,2) - 2*y1*y2 + Power(y2,2)))/Power(y1 - y2,2))/(2.*y2);
  x:=(x2*y1 - Sqrt(y1)*Sqrt(y2)*Sqrt(Power(x2,2) + Power(y1,2) - 2*y1*y2 + Power(y2,2)))/(y1 - y2);
 end;
 end;
 if y<0 then exit;
 r:=y;
 // Testen
 if r>1e9 then exit;
 d1:=sqr(x)+sqr(y1-y)-sqr(r);
 if (abs(d1)>500*tol) then exit;
end;
if wasminus then y:=-y;
x:=x+xx;
rotatep(x,y,x0,0,-phi,x,y);
if waschange then change(x,y);
fail:=false;
end;

procedure xcircleclp2(x1,y1,r1,a2,b2,c2,x3,y3:double;mode:integer;
          var x,y,r:double;var fail:boolean);
// Spezielversion, wenn x1,y1 und x2,y3 den gleichen abstand von a2,b2,c2 haben
// und die Radien gleich sin
var d,d1,d2,ff,f1,f2,x0,a,xx,phi:double;
 waschange,wasminus:boolean;
 x32,y32,r12,y12:double;
 z1,z2,z3:double;
begin
fail:=true;waschange:=false;wasminus:=false;
if abs(a2)<abs(b2) then begin
 waschange:=true;
 change(a2,b2);change(x1,y1);change(x3,y3);
end;
x0:=c2/a2;
phi:=atan2(a2,b2)/radiant;
rotatep(x1,y1,x0,0,phi,x1,y1);
rotatep(x3,y3,x0,0,phi,x3,y3);
if abs(y1-y3)>tol then exit;
xx:=x1;
x1:=0;x3:=x3-xx;
{Linie= X Achse Mittelpunkt Kreis auf y Achse}
x32:=sqr(x3);y32:=sqr(y3);r12:=sqr(r1);y12:=sqr(y1);
if y3<0 then begin
 wasminus:=true;
 y1:=-y1;y3:=-y3;
end;
case mode of
 0:begin
  x:=x3/2;
  y:=(sqr(r1)-sqr(y1)-sqr(x))/-(2*(y1+r1));
 end;
 1:begin
  x:=x3/2;
  y:=(sqr(r1)-sqr(y1)-sqr(x))/-(2*(y1-r1));
 end;
 else exit;
end;
r:=y;
if wasminus then y:=-y;
x:=x+xx;
rotatep(x,y,x0,0,-phi,x,y);
if waschange then change(x,y);
fail:=false;
end;

procedure circleclp(x1,y1,r1,a2,b2,c2,x3,y3:double;mode:integer;
          var x,y,r:double;var fail:boolean);
var i:integer;
begin
loesung:=0;
for i:=0 to 3 do begin
 xcircleclp(x1,y1,r1,a2,b2,c2,x3,y3,i,x,y,r,fail);
 if fail then continue;
 if newloesung(x,y,r)=mode then exit;
end;
fail:=true;
end;
procedure xmcircle(px1,py1,pr1,px2,py2,pr2,px3,py3,pr3:double;
                  var x,y,r:double;var fail:boolean);forward;
{spezial Prozedur um einen Kreis der innerhalb von 3 Kreisen ist zu finden}

procedure xcirclecc0(px1,py1,pr1,px2,py2,pr2:double;mode:integer;
                   var x,y,r:double;var fail:boolean);
{Hilfsroutine
 Emittelt L”sungen tangential zu den Kreisen px1,py1,pr1 und px2,py2,pr2
 und zu dem Punkt 0,0.
 Entwickelt nach einem Hinweis von
}
var a,b,c,u,w,x1,y1,r1,x2,y2,r2,k1,k2,k3,f,g,h1,h2,hh1,hh2,gg:double;
 xx,yy,rr,ff:double;
 mode2,i,k:integer;
 a1,a2,b1,b2,c1,c2,d1,d2,r12,r22,x12,x22,y12,y22,aaff:double;
 iterieren:boolean;
//var tol,tolaff:double;
label 1;
begin
{
tolaaff:=(abs(px1)+abs(py1)+abs(pr1)+abs(px2)+abs(py2)+abs(pr2))/1e9;
tol:=tolaaff*10;
}
fail:=true;iterieren:=false;
mode2:=mode div 2;
x1:=px1;y1:=py1;r1:=pr1;x2:=px2;y2:=py2;r2:=pr2;
if (r1<0)or(r2<0) then exit;
if abs(r1)+abs(r2)<tol then begin
 xcircle3p(x1,y1,x2,y2,0,0,x,y,r,fail);exit;
end;
if eq3(x1,y1,0,0,0,0)then begin
 {x1,y1 Nullpunkt}
 r:=r1/2;
 point2c(0,0,r,x2,y2,r+r2,mode,x,y,fail);
 exit;
end else if eq3(x2,y2,0,0,0,0) then begin
 r:=r2/2;
 point2c(0,0,r,x1,y1,r+r1,mode,x,y,fail);
 exit;
end;
r12:=r1*r1;r22:=r2*r2;
x12:=x1*x1;x22:=x2*x2;
y12:=y1*y1;y22:=y2*y2;
a1:=-2*x1*r2;b1:=-2*y1*r2;
a2:=-2*x2*r1;b2:=-2*y2*r1;
c1:=r2*x12+r2*y12;
c2:=r1*x22+r1*y22;
d1:=r12*r2;
d2:=r22*r1;
if mode2=0 then begin
 k1:=a1-a2;
 k2:=b1-b2;
 k3:=c1-c2-d1+d2;
end else begin
 k1:=a1+a2;
 k2:=b1+b2;
 k3:=c1+c2-d1-d2;
end;
if abs(k2)>abs(k1) then begin
 f:=-k1/k2;g:=-k3/k2;
 {
 dline(-4,-4*f+g,4,4*f+g);
 }
 h1:=g-y1;
 h2:=g-y2;
 hh1:=h1*h1;hh2:=h2*h2;gg:=g*g;
 if r1<>0 then begin
  a:=(-2*x1-2*f*y1)/(2*r1);
  b:=(x1*x1-2*g*y1+y1*y1-r1*r1)/(2*r1);
  aaff:=a*a-1-f*f;
  if abs(aaff)<tolaaff then begin
   aaff:=tolaaff;iterieren:=true;
  end;
  u:=(2*a*b-2*f*g)/(aaff);
  w:=(g*g-b*b)/(aaff);
  c:=w+(u/2)*(u/2);
  if c<-eps then exit
  else if c<0 then c:=0;
  if mode mod 2 =0 then x:=sqrt(c)-(u/2)
  else x:=-sqrt(c)-(u/2);
 end else if r2<>0 then begin
  a:=(-2*x2-2*f*y2)/(2*r2);
  b:=(x2*x2-2*g*y2+y2*y2-r2*r2)/(2*r2);
  aaff:=a*a-1-f*f;
  if abs(aaff)<tolaaff then begin
   aaff:=tolaaff;;iterieren:=true;
  end;
  u:=(2*a*b-2*f*g)/(aaff);
  w:=(g*g-b*b)/(aaff);
  c:=w+(u/2)*(u/2);
  if c<-eps then exit
  else if c<0 then c:=0;
  if mode mod 2 =0 then x:=sqrt(c)-(u/2)
  else x:=-sqrt(c)-(u/2);
 end else begin
  x:=(2*g*y1-x1*x1-y1*y1)/(-2*x1-2*f*y1);
 end;
 y:=f*x+g;
end else begin
 if k1=0 then begin
  k1:=k1;
 end;
 if abs(k1)<eps then exit;
 f:=-k2/k1;g:=-k3/k1;
 {
 dline(-4*f+g,-4,4*f+g,4);
 }
 h1:=g-y1;
 h2:=g-y2;
 hh1:=h1*h1;hh2:=h2*h2;gg:=g*g;

 if r1<>0 then begin
  a:=(-2*y1-2*f*x1)/(2*r1);
  b:=(y1*y1-2*g*x1+x1*x1-r1*r1)/(2*r1);
  aaff:=a*a-1-f*f;
  if abs(aaff)<tolaaff then begin
   aaff:=tolaaff;iterieren:=true;
  end;
  u:=(2*a*b-2*f*g)/(aaff);
  w:=(g*g-b*b)/(aaff);
  c:=w+(u/2)*(u/2);
  if c<-tol then exit
  else if c<0 then c:=0;
  if mode mod 2 =0 then y:=sqrt(c)-(u/2)
  else y:=-sqrt(c)-(u/2);
 end else if r2<>0 then begin
  a:=(-2*y2-2*f*x2)/(2*r2);
  b:=(y2*y2-2*g*x2+x2*x2-r2*r2)/(2*r2);
  aaff:=a*a-1-f*f;
  if abs(aaff)<tolaaff then begin
   aaff:=tolaaff;iterieren:=true;
  end;
  u:=(2*a*b-2*f*g)/(aaff);
  w:=(g*g-b*b)/(aaff);

  c:=w+(u/2)*(u/2);
  if c<-tol then exit
  else if c<0 then c:=0;
  if mode mod 2 =0 then y:=sqrt(c)-(u/2)
  else y:=-sqrt(c)-(u/2);
 end else begin
  y:=(2*g*x1-y1*y1-x1*x1)/(-2*y1-2*f*x1);
 end;
 x:=f*y+g;
end;
r:=sqrt(x*x+y*y);
r1:=r-tol*1000;r2:=r+tol*1000;
ff:=sqr(r)-sqr(x)-sqr(y);
{
fail:=false;
exit; // kein iterieren
}
if iterieren then begin
 if round(r)=56 then begin
  r:=r;
 end;
 ff:=sqr(r)-sqr(x)-sqr(y);
 if abs(ff)>100*tol then begin
  fail:=true;exit;
 end;
 {
 if abs(ff)<tol*0.001 then begin
  fail:=false;exit;
 end;
 } 
 for k:=1 to 20 do begin
  for i:=0 to 7 do begin
   xcircle2cr(px1,py1,pr1,px2,py2,pr2,r,i,x,y,fail);
   rr:=r;
   if not fail then begin
    ff:=sqr(r)-sqr(x)-sqr(y);
    if abs(ff)>tol*10000 then continue;
    if abs(ff)<tol*0.001 then begin
     r:=r;fail:=false;
     exit;
    end;
    if ff>0 then r2:=r //r zu groß
    else r1:=r;
    r:=(r1+r2)/2;
    if abs(r1-r2)<tol*0.001 then begin
     fail:=false;exit;
    end; 
    goto 1;
   end;
  end;
  1:
 end;
 if abs(ff)>tol then begin
  fail:=true;exit;
 end;
end;

{
setcolor(12);
dcircle(x1,y1,r1);
dcircle(x2,y2,r2);
dpoint(0,0);
setcolor(11);
dcircle(x,y,r);

readln;
}
fail:=false;
end;

function pointinc(x1,y1,x2,y2,r2:double):boolean;
begin
pointinc:=distance2p(x1,y1,x2,y2)<r2;
end;

procedure xmcircle(px1,py1,pr1,px2,py2,pr2,px3,py3,pr3:double;
                  var x,y,r:double;var fail:boolean);
{spezial Prozedur um einen Kreis der innerhalb von 3 Kreisen ist zu finden}
var mode,si,s,i,j,k,rmin,mode0,mode1,mode2:integer;
 f,x1,y1,r1,x2,y2,r2,x3,y3,r3:double;
 ok,fail1,fail2,fail3:boolean;
 a1,a2,a3,b1,b2,b3,c1,c2,c3,d1,d2,d3,xx,yy,rr:double;
 xl1,yl1,xl2,yl2,rr1,xx1,yy1,xx2,yy2,xxx,yyy,rrr:double;
 x4,y4,x5,y5,x6,y6:double;
 istline:boolean;
 alpha,xs,ys,x11,y11,al,d,cmitte,xm,ym:double;
 a,b,c:array[0..2]of double;
 w:array[0..1,0..1,1..3]of double;
 rrold,xxxx,yyyy,rrrr:double;
begin
x1:=px1;y1:=py1;r1:=pr1;x2:=px2;y2:=py2;r2:=pr2;x3:=px3;y3:=py3;r3:=pr3;
point2c(x1,y1,r1,x2,y2,r2,0,xx1,yy1,fail);
if fail then exit;
point2c(x1,y1,r1,x2,y2,r2,1,xx2,yy2,fail);
if fail then exit;
if pointinc(xx2,yy2,x3,y3,r3) then begin
 change(xx1,xx2);change(yy1,yy2);
end;
rr1:=distance2p(xx1,yy1,xx2,yy2);

{Kreis an dem transformiert wird}
{
setcolor(13);
dcircle(xx1,yy1,rr1);
}
point2c(x1,y1,r1,xx1,yy1,rr1,0,xl1,yl1,fail);
if fail then exit;
if distance2p(xl1,yl1,xx2,yy2)<tol then
point2c(x1,y1,r1,xx1,yy1,rr1,1,xl1,yl1,fail);
if fail then exit;
point2c(x2,y2,r2,xx1,yy1,rr1,0,xl2,yl2,fail);
if fail then exit;
if distance2p(xl2,yl2,xx2,yy2)<tol then
point2c(x2,y2,r2,xx1,yy1,rr1,1,xl2,yl2,fail);
if fail then exit;

line2p(xx2,yy2,xl1,yl1,a1,b1,c1,fail);{trans. c1}
if fail then exit;
line2p(xx2,yy2,xl2,yl2,a2,b2,c2,fail);{trans. c2}
if fail then exit;
{in Linien Tranformierte Kreise}
{
setcolor(12);
dhline(a1,b1,c1);
dhline(a2,b2,c2);
}
(*
polarpoint(x3,y3,0,r3,x4,y4);
polarpoint(x3,y3,90,r3,x5,y5);
polarpoint(x3,y3,180,r3,x6,y6);
ctransform(x4,y4,xx1,yy1,rr1,x4,y4,fail);
ctransform(x5,y5,xx1,yy1,rr1,x5,y5,fail);
ctransform(x6,y6,xx1,yy1,rr1,x6,y6,fail);

xcircle3p(x4,y4,x5,y5,x6,y6,xx,yy,rr,fail);{trans. c3}
*)
ctransform(false,x3,y3,r3,xx1,yy1,rr1,istline,xx,yy,rr);
{In kreis Tranformierter 3.Kreis}
{
dcircle(xx,yy,rr);
}
if istline then begin
 {circle3l(xx,yy,rr,a1,b1,c1,a2,b2,c2,i,xxx,yyy,rrr,fail);}
 a[0]:=xx;b[0]:=yy;c[0]:=rr;
 a[1]:=a1;b[1]:=b1;c[1]:=c1;
 a[2]:=a2;b[2]:=b2;c[2]:=c2;
 for i:=0 to 1 do begin
  j:=(i+2) mod 3;
  for k:=0 to 1 do begin
   s:=2*k-1;
   f:=sqrt(sqr(a[i]+s*a[j])+sqr(b[i]+s*b[j]));
   w[i,k,1]:=(a[i]+s*a[j])/f;
   w[i,k,2]:=(b[i]+s*b[j])/f;
   w[i,k,3]:=(c[i]+s*c[j])/f;
  end;
 end;
 for si:=0 to 3 do begin
  i:=si and 1; j:=si shr 1;
  point2l(w[1,i,1],w[1,i,2],w[1,i,3],w[0,j,1],w[0,j,2],w[0,j,3],xxx,yyy,fail);
  if not fail then begin
   rrr:=abs(a[0]*x+b[0]*y-c[0]);
   ctransform(false,xxx,yyy,rrr,xx1,yy1,rr1,istline,x,y,r);
   fail:=istline;
   if (not fail) and(r3>=r)and pointinc(x,y,x3,y3,r3) and pointinc(x,y,x2,y2,r2)
   and pointinc(x,y,x1,y1,r1) then exit;
  end;
 end;
 fail:=true;exit;
end else begin
 if abs(a1-a2)+abs(b1-b2)<delta then begin
  rrold:=rr;
  for mode:=0 to 3 do begin
   cmitte:=(c1+c2)/2;
   mode2:=mode div 2;
   rrr:=abs(c2-c1)/2;
   if not odd(mode) then rr:=rrr+rrold else rr:=abs(rrr-rrold);
   if rr<0 then continue;
   if cmitte<0 then begin
    a1:=-a1;b1:=-b1;cmitte:=-cmitte;
   end;
   pointcl(xx,yy,rr,a1,b1,cmitte,mode2,xxx,yyy,fail);
   if fail then continue;
   ctransform(false,xxx,yyy,rrr,xx1,yy1,rr1,istline,x,y,r);
   fail:=istline;
   if (not fail) and(r3>=r)and pointinc(x,y,x3,y3,r3) and pointinc(x,y,x2,y2,r2)
   and pointinc(x,y,x1,y1,r1) then exit;
  end;
  fail:=true;
 end else begin
  point2l(a1,b1,c1,a2,b2,c2,xxx,yyy,fail);
  d:=1{distance2p(xx,yy,x1,y1)};
  for i:=0 to 3 do begin
   mode1:=i mod 2;mode2:=i div 2;
   pointcl(xxx,yyy,d,a1,b1,c1,mode1,xx1,yy1,fail);
   if fail then continue;
   pointcl(xxx,yyy,d,a2,b2,c2,1-mode2,xx2,yy2,fail);
   {
   setcolor(13);
   dhline(a2,b2,c2);dcircle(xxx,yyy,d);
   }
   if fail then continue;
   xm:=(xx1+xx2)/2;ym:=(yy1+yy2)/2;
   al:=getphi(xxx,yyy,xm,ym);
   orthokoord(xxx,yyy,xm,ym,x1,y1,x11,y11);
   orthokoord(xxx,yyy,xm,ym,xx2,yy2,xs,ys);
   alpha:=getphi(0,0,xs,ys);
   if alpha>180 then alpha:=360-alpha;
   for k:=0 to 3 do begin
    xcirclec2l(alpha,x11,y11,r1,k,xxxx,yyyy,rrrr,fail);
    if fail then continue;
    polarpoint(xxx,yyy,al,xxxx,xxxx,yyyy);
    ctransform(false,xxxx,yyy,rrrr,xx1,yy1,rr1,istline,x,y,r);
    fail:=istline;
    if (not fail) and(r3>=r)and pointinc(x,y,x3,y3,r3) and pointinc(x,y,x2,y2,r2)
    and pointinc(x,y,x1,y1,r1) then exit;
   end;
  end;
  fail:=true;
 end;
end;
end;


procedure circle3c(px1,py1,pr1,px2,py2,pr2,px3,py3,pr3:double;pmode:integer;
                  var x,y,r:double;var fail:boolean);
var mode0,mode1,mode2:integer;
 x1,y1,r1,x2,y2,r2,x3,y3,r3,rr:double;
label 1;
begin
fail:=true;
loesung:=0;
if eq3(px1,py1,pr1,px2,py2,pr2)then exit;
if eq3(px1,py1,pr1,px3,py3,pr3)then exit;
if eq3(px2,py2,pr2,px3,py3,pr3)then exit;
x1:=px1;y1:=py1;r1:=pr1;x2:=px2;y2:=py2;r2:=pr2;x3:=px3;y3:=py3;r3:=pr3;
if r2<r3 then begin
 change(x2,x3);change(y2,y3);change(r2,r3);
end;
if r1<r3 then begin
 change(x1,x3);change(y1,y3);change(r1,r3);
end;
x1:=x1-x3;y1:=y1-y3;
x2:=x2-x3;y2:=y2-y3;
for mode0:=0 to 3 do begin
 for mode1:=0 to 3 do begin
  case mode1 of
   0:xcirclecc0(x1,y1,r1+r3,x2,y2,r2+r3,mode0,x,y,r,fail);
   1:xcirclecc0(x1,y1,r1-r3,x2,y2,r2+r3,mode0,x,y,r,fail);
   2:xcirclecc0(x1,y1,r1-r3,x2,y2,r2-r3,mode0,x,y,r,fail);
   3:xcirclecc0(x1,y1,r1+r3,x2,y2,r2-r3,mode0,x,y,r,fail);
  end;
  if not fail then begin
   x:=x+x3;y:=y+y3;rr:=r;
   for mode2:=0 to 1 do begin
    if not odd(mode2) then r:=rr+r3 else begin
     r:=rr-r3;
    end;
    r:=abs(r);
    {
    if not fail then fail:=not test1(distance2p(px1,py1,x,y),pr1+r,abs(r-pr1));
    if not fail then fail:=not test1(distance2p(px2,py2,x,y),pr2+r,abs(r-pr2));
    if not fail then fail:=not test1(distance2p(px3,py3,x,y),pr3+r,abs(r-pr3));
    }
    if not test1(distance2p(px1,py1,x,y),pr1+r,abs(r-pr1)) then continue;
    if not test1(distance2p(px2,py2,x,y),pr2+r,abs(r-pr2)) then continue;
    if not test1(distance2p(px3,py3,x,y),pr3+r,abs(r-pr3)) then continue;
    if newloesung(x,y,r)=pmode then exit;
   end;
   fail:=true;
  end;
 end;
end;
{
xmcircle(px1,py1,pr1,px2,py2,pr2,px3,py3,pr3,x,y,r,fail);
if fail or (not newloesung(y,y,r)=pmode) then exit;
}
end;

(*
procedure circle3c(px1,py1,pr1,px2,py2,pr2,px3,py3,pr3:double;pmode:integer;
                  var x,y,r:double;var fail:boolean);
var i,j,mode,rmin,mode0,mode1,mode2:integer;
 x1,y1,r1,x2,y2,r2,x3,y3,r3:double;
 ok,fail1,fail2,fail3:boolean;
 xl,yl,rl:array[0..7]of double;
label 1;
begin
fail:=true;
if eq3(px1,py1,pr1,px2,py2,pr2)then exit;
if eq3(px1,py1,pr1,px3,py3,pr3)then exit;
if eq3(px2,py2,pr2,px3,py3,pr3)then exit;
ok:=false;
{test ob alle 3 �berlappen}
{
if (distance2p(px1,py1,px2,py2)<pr1+pr2)
and (distance2p(px1,py1,px3,py3)<pr1+pr3)
and (distance2p(px2,py2,px3,py3)<pr2+pr3)then begin
if fail
 if not fail then exit;
 end else dec(pmode);
end;
}
loesung:=0;
for i:=-1 to 95 do begin
 if i<0 then begin
  xmcircle(px1,py1,pr1,px2,py2,pr2,px3,py3,pr3,x,y,r,fail);
  if fail or (not newloesung(y,y,r)=pmode) then continue
  else exit;
 end;
 x1:=px1;y1:=py1;r1:=pr1;x2:=px2;y2:=py2;r2:=pr2;x3:=px3;y3:=py3;r3:=pr3;
 mode0:=i mod 4;mode:=i div 4;
 mode1:=mode mod 4;mode:=mode div 4;
 mode2:=mode mod 2;mode:=mode div 2;

 if r2<r1 then begin
  if r3<r2 then rmin:=3 else rmin:=2;
 end else if r3<r1 then rmin:=3 else rmin:=1;
 if mode=1 then begin
  change(x1,x3);change(y1,y3);change(r1,r3);
 end else if mode=2 then begin
  change(x2,x3);change(y2,y3);change(r2,r3);
 end;
 {
 if (mode1=2)and(mode=0)and(mode2=0) then begin
  setcolor(12);
  dcircle(x1,y1,r1-r3);
  dcircle(x2,y2,r2-r3);
  dpoint(x3,y3);
 end;
 }
 x1:=x1-x3;y1:=y1-y3;
 x2:=x2-x3;y2:=y2-y3;
 case mode1 of
  0:xcirclecc0(x1,y1,r1+r3,x2,y2,r2+r3,mode0,x,y,r,fail);
  1:xcirclecc0(x1,y1,r1-r3,x2,y2,r2+r3,mode0,x,y,r,fail);
  2:xcirclecc0(x1,y1,r1-r3,x2,y2,r2-r3,mode0,x,y,r,fail);
  3:xcirclecc0(x1,y1,r1+r3,x2,y2,r2-r3,mode0,x,y,r,fail);
 end;
 if not fail then begin
  x:=x+x3;y:=y+y3;
  {
  if (mode1=2)and(mode=0)and(mode2=0) then begin
   setcolor(11);
   dcircle(x,y,r);
  end;
  }
  if odd(mode2) then r:=r+r3 else r:=r-r3;
  {fail:=(r<eps);}
  r:=abs(r);
  {
  if (mode0=1)and(mode1=2)and(mode=0)and(mode2=0) then begin
   setcolor(13);
   dcircle(x,y,r);
   x:=x;
  end;
  }
  if not fail then fail:=not test1(distance2p(px1,py1,x,y),pr1+r,abs(r-pr1));
  if not fail then fail:=not test1(distance2p(px2,py2,x,y),pr2+r,abs(r-pr2));
  if not fail then fail:=not test1(distance2p(px3,py3,x,y),pr3+r,abs(r-pr3));
 end;
 if not fail then begin
  if newloesung(x,y,r)=pmode then exit;
  1:
 end;
end;
fail:=true;
end;
*)
procedure circlec2p(x1,y1,r1,x2,y2,x3,y3:double;pmode:integer;
                  var x,y,r:double;var fail:boolean);
var mode:integer;
begin
fail:=true;loesung:=0;
if eq3(x2,y2,0,x3,y3,0)then exit;
x1:=x1-x3;y1:=y1-y3;
x2:=x2-x3;y2:=y2-y3;
for mode:=0 to 0 do begin
 xcirclecc0(x1,y1,r1,x2,y2,0,mode,x,y,r,fail);
 if not fail then begin
  x:=x+x3;y:=y+y3;
  if newloesung(x,y,r)=pmode then begin
   exit;
  end;
 end;
end;
fail:=true;
end;

procedure circle2cp(x1,y1,r1,x2,y2,r2,x3,y3:double;pmode:integer;
                    var x,y,r:double;var fail:boolean);
var mode:integer;
begin
fail:=true;loesung:=0;
for mode:=0 to 3 do begin
 xcirclecc0(x1-x3,y1-y3,r1,x2-x3,y2-y3,r2,mode,x,y,r,fail);
 //if r>(abs(x1-x2)+abs(y2-y1))*1000 then fail:=true;
 if not fail then begin
  x:=x+x3;y:=y+y3;
  if newloesung(x,y,r)=pmode then begin
   exit;
  end;
 end;
end;
fail:=true;
end;

procedure circle2cl(x1,y1,r1,x2,y2,r2,a3,b3,c3:double;mode:integer;
                    var x,y,r:double;var fail:boolean);
var xl,yl,xc,yc,rc:double;
 a1,b1,c1,a2,b2,c2,x3,y3,r3:double;
 xxx,yyy,alpha,xs,ys,x11,y11,al,xm,ym,xx1,yy1,xx2,yy2,d,xx,yy:double;
 phi1,phi2:double;
 isline,fail1,fail2:boolean;
 i,mode1,mode2,k,j,kk:integer;
 rr,drr,dr,dc,px1,py1,pr1,px2,py2,pr2,pa3,pb3,pc3,d1,d2,dd,dd2:double;
 test,second:boolean;
label 1,2,3,10;
begin
loesung:=0;
{second:=false;}
test:=false;
point2c(x1,y1,r1,x2,y2,r2,1,xc,yc,fail);
if not fail then begin //Kreise scheiden sich echt
 lotpl(xc,yc,a3,b3,c3,xl,yl);
 rc:=distance2p(xc,yc,xl,yl);
 if rc<r1 then rc:=r1;
 if rc<r2 then rc:=r2;
 ctransform(false,x1,y1,r1,xc,yc,rc,isline,a1,b1,c1);
 ctransform(false,x2,y2,r2,xc,yc,rc,isline,a2,b2,c2);
 ctransform(true,a3,b3,c3,xc,yc,rc,isline,x3,y3,r3);
 10:
 drawgerade(a1,b1,c1,3);
 drawgerade(a2,b2,c2,3);
 drawkreis(x3,y3,r3,3);
 drawkreis(xc,yc,rc,5);
 if abs(a1-a2)+abs(b1-b2)<delta then begin // Transformierte Kreise geben parallel Geraden
  for i:=0 to 3 do begin
   xcirclecpar(x3,y3,r3,a1,b1,c1,c2,i,x,y,r,fail);
   if fail then continue;
   ctransform(false,x,y,r,xc,yc,rc,isline,x,y,r);
   fail:=isline;
   if fail then continue;
   if test and (abs(distancepl(x,y,a3,b3,c3)-r)>tol) then continue;
   if newloesung(x,y,r)= mode then exit;
  end;
 end else begin
  point2l(a1,b1,c1,a2,b2,c2,xx,yy,fail);
  if fail then exit;
  dd:=distance2p(xx,yy,x3,y3);
  phi1:=winkell(a1,b1);
  phi2:=winkell(a2,b2);
  al:=(phi2+phi1)/2;
  alpha:=abs(phi2-phi1)/2;
  x11:=x3-xx;y11:=y3-yy; // Schnittpunkt in den Nullpunkt
  rotatep(x11,y11,0,0,-al,x11,y11);
  drawkreis(x11,y11,r3,7);
  for k:=0 to 7 do begin
   {dd2:=distance2p(0,0,x11,y11);}
   xcirclec2l(alpha,x11,y11,r3,k,xxx,yyy,r,fail);
   if fail then goto 2;
   rotatep(xxx,yyy,0,0,al,x,y);    //Ergebnis drehen
   x:=x+xx;y:=y+yy; // und verschieben
   drawkreis(x,y,r,3);
   ctransform(false,x,y,r,xc,yc,rc,isline,x,y,r);
   if isline then goto 2;
   drawkreis(x,y,r,1);
   if (abs(abs(distancepl(x,y,a3,b3,c3))-r)>tol) then goto 2;
   if newloesung(x,y,r)= mode then begin
    exit;
   end;
   2:
  end;
 end;
 fail:=true;
end else begin
 pointcl(x1,y1,r1,a3,b3,c3,1,xc,yc,fail);
 if not fail then begin // Kreis 1 und Linie schneiden sich
  1:
  rc:=2*r1;
  ctransform(false,x1,y1,r1,xc,yc,rc,isline,a1,b1,c1);
  ctransform(false,x2,y2,r2,xc,yc,rc,isline,x3,y3,r3);
  a2:=a3;b2:=b3;c2:=c3;
  goto 10;
 end else begin
  pointcl(x2,y2,r2,a3,b3,c3,1,xc,yc,fail); // Kreis 2 und Lnie schneiden sich
  if not fail then begin
   change(x1,x2);change(y1,y2);change(r1,r2);goto 1;
  end else begin
   {alle 3 Elemente sind disjunkt}
   if r1<r2 then begin
    change(x1,x2);change(y1,y2);change(r1,r2);
   end;
   px1:=x1;py1:=y1;pr1:=r1;
   px2:=x2;py2:=y2;pr2:=r2;
   pa3:=a3;pb3:=b3;pc3:=c3;
   dr:=r2;dc:=r2;drr:=r2;
   for i:=0 to 1 do begin
    for j:=0 to 1 do begin
     r1:=pr1+dr;c3:=pc3+dc;
     for k:=0 to 3 do begin
      xcircleclp(x1,y1,r1,a3,b3,c3,x2,y2,k,x,y,rr,fail);
      if fail then continue;
      for kk:=0 to 1 do begin
       r:=rr-drr;
       drr:=-drr;
       if r<0 then continue;
       if not test1(abs(distance2p(x1,y1,x,y)),pr1+r,r-pr1) then continue;
       if abs(abs(distancepl(x,y,pa3,pb3,pc3))-r)>tol then continue;
       drawkreis(x,y,r,4);
       if newloesung(x,y,r)=mode then exit;
      end;
     end;
     dr:=-dr;
    end;
    dc:=-dc;
   end;
   if abs(pr1-pr2)<tol then
   for i:=0 to 1 do begin
    xcircleclp2(px1,py1,pr1,pa3,pb3,pc3,px2,py2,i,x,y,r,fail);
    if fail then continue;
    drawkreis(x,y,r,4);
    if newloesung(x,y,r)=mode then exit;
   end;
   fail:=true;
  end;
 end;
end;
end;

procedure circle2cx(px1,py1,pr1,px2,py2,pr2,px:double;mode:integer;var y,r:double);
// Kreis Tangential an 2 Kreise mit gegebener x Koordinate des Mittelpunktes
var dd,s,s2:double;
 a,b,c,d,e:double;
 x,x1,y1,x2,y2,r1,r2,rr:double;
i:integer;
const epsfaktor=0.00001;
var eps:double;
begin
eps:=(pr1+pr2)*epsfaktor;
// Keineren Kreis auf Radius null setzen und größeren Kreis in den Nullpunkt
if pr2>pr1 then begin
 d2basic.change(px1,px2);
 d2basic.change(py1,py2);
 d2basic.change(pr1,pr2);
end;
x1:=px1;
x2:=px2-px1;y2:=py2-py1;r1:=pr1-pr2;x:=px-px1;r2:=pr2;
if r1>eps then begin         // Kreise unterschiedlich groß
 // x²+y²=(r1+r)²
 a:=(-2*y2)/(-2*r1);
 b:=(-2*x*x2+sqr(x2)+sqr(y2)+sqr(r1))/(-2*r1);
 if abs(1-sqr(a))>eps then begin
  c:=(a*r1+a*b)/(1-sqr(a));
  d:=(-sqr(x)+sqr(r1)+2*B*r1+sqr(b))/(1-sqr(a));
  dd:=d+sqr(c);
  if dd>-eps then begin
   if dd<0 then dd:=0;
   e:=sqrt(dd);
   y:=e+c;
   r:=a*y+b;
   if r-r2>eps then begin
    y:=y+py1;r:=r+r2;
    if newloesung(px,y,r)=mode then exit;
   end;
   if e<>0 then begin
    y:=-e+c;
    r:=a*y+b;
    if r-r2>eps then begin
     y:=y+py1;r:=r-r2;
     if newloesung(px,y,r)=mode then exit;
    end;
   end;
  end;
 end;
 // x²+y²=(r1-r)²
 a:=(-2*y2)/(+2*r1);
 b:=(-2*x*x2+sqr(x2)+sqr(y2)+sqr(r1))/(+2*r1);
 if abs(1-sqr(a))>eps then begin
  c:=(-a*r1+a*b)/(1-sqr(a));
  d:=(-sqr(x)+sqr(r1)-2*B*r1+sqr(b))/(1-sqr(a));
  dd:=d+sqr(c);
  if dd>-eps then begin
   if dd<0 then dd:=0;
   e:=sqrt(dd);
   y:=e+c;
   r:=a*y+b;
   if r+r2>eps then begin
    y:=y+py1;r:=r+r2;
    newloesung(px,y,r);
   end;
   if e<>0 then begin
    y:=-e+c;
    r:=a*y+b;
    if r+r2>eps then begin
     y:=y+py1;r:=r+r2;
     if newloesung(px,y,r)=mode then exit;
    end;
   end;
  end;
 end;
end else begin // Kreise gleich groß
 if abs(y2)>eps then begin
  y:=(-2*x*x2+sqr(x2)+sqr(y2))/(2*y2);
  r:=sqrt(sqr(x)+sqr(y));
  y:=y+py1;rr:=r;
  r:=rr+r2;
  if newloesung(px,y,r)=mode then exit;
  r:=rr-r2;
  if newloesung(px,y,r)=mode then exit;
 end;
end;
end;
procedure llotlc(a,b,c,x,y,r:double;mode:integer;var aa,bb,cc:double);
// Linie als Lot von Kreis x,y,r an Linie a,b,c
var d:double;
begin
aa:=b;bb:=-a;
d:=distancepl(x,y,aa,bb,0);
if mode=0 then cc:=d+r else cc:=d-r;
end;

procedure lotlc(xx2,yy2,xx3,yy3,xx1,yy1,r:double;alter:integer;var x,y:double);
{ get a tangent point of circle generated by a vertical line of the given line
  alter=0,1 --> point is on the circle
  alter=2,3 --> point is on the line
}
const tol=0.001;
var
si,s1:integer;
dx,dy,fact,my,d1,d2:real;
begin
alter:=alter mod 4;
if alter>1 then si:=alter-2 else si:=alter;
s1:=2*(si and 1)-1;
dx:=xx3-xx2;dy:=yy3-yy2;
if dx=0 then yy2:=yy1 else begin
  fact:=dy/dx;d1:=dy;d2:=-dx;
  my:=(yy2-xx2*fact-yy1+xx1*fact)/(d2-d1*fact);
  xx2:=xx1+(my*d1);yy2:=yy1+(my*d2);
end;
fact:=r/sqrt(sqr(dx)+sqr(dy)); dx:=s1*dx*fact;dy:=s1*dy*fact;
if alter>1 then begin x:=xx2+dx; y:=yy2+dy;end
           else begin x:=xx1+dx; y:=yy1+dy;end;
end;
procedure parallelepl(xx0,yy0,xx1,yy1,xx2,yy2:double;var x1,y1,x2,y2:double);
{ get parallel line of the given line( xx1,yy1,xx2,yy2).this parallel line is located by the
  given point(xx0,yy0) }
var
 phi,k,dx,dy,a0,b0,c0,a1,b1,c1,a2,b2,c2:double;
 fail:boolean;
begin
dx:=xx2-xx1;dy:=yy2-yy1;
if dy = 0 then begin x1:=xx1;y1:=yy0;x2:=xx2;y2:=yy0;end
else if dx = 0 then begin x1:=xx0;y1:=yy1;x2:=xx0;y2:=yy2;end
else begin
  k:= dy/dx;
  a0:= k; b0:= -1; c0:= yy0-k*xx0;
  a1:= 1; b1:= k;  c1:= -xx1-k*yy1;
  a2:= 1; b2:= k;  c2:= -xx2-k*yy2;
  point2l(a0,b0,c0,a1,b1,c1,x1,y1,fail);
  point2l(a0,b0,c0,a2,b2,c2,x2,y2,fail);
end;
end;

procedure parallelelc(lx1,ly1,lx2,ly2,cx,cy,r:double;alter:integer;var x1,y1,x2,y2:double);
{ a tangent line is generated which is also parallel to the given line }
var
k,dx,dy,phi,x,y,f,d:real;
s1,s2:integer;
begin
dx:=lx2-lx1; dy:=ly2-ly1;
if dx=0 then k:=0 else k:=dy/dx;
phi:=getphi(lx1,ly1,lx2,ly2);
f:=cy-k*cx-ly1+k*lx1;
s1:=0;s2:=0;
if k>=0  then begin
  if f>0 then begin s1:=(2*alter-1)*(-1); s2:=-s1;end
         else begin s1:=2*alter-1; s2:=-s1;end;
end else if k<0 then begin
  if f>0 then begin s1:=2*alter-1; s2:=s1;end
         else begin s1:=(2*alter-1)*(-1); s2:=s1;end;
end;
x:=cx+s1*r*abs(sin(phi*radiant)); y:=cy+s2*r*abs(cos(phi*radiant));
parallelepl(x,y,lx1,ly1,lx2,ly2,x1,y1,x2,y2);
end;
(*
procedure lotpl2p(x1,y1,x2,y2,x3,y3:double;var x,y:double;var fail:boolean);
{Lot vom punkt x1,y1 auf die Linie x2,y2,x3,y3}
var a,b,c:double;
begin
line2p(x2,y2,x3,y3,a,b,c,fail);
if not fail then lotpl(x1,y1,a,b,c,x,y);
end;
*)
procedure Lotpl2p(const Px,Py,x1,y1,x2,y2:double;var x,y:double;var fail:boolean);
var
  Vx    : double;
  Vy    : double;
  Wx    : double;
  Wy    : double;
  c1    : double;
  c2    : double;
  Ratio : double;
begin
fail:=true;
Vx := x2 - x1;
Vy := y2 - y1;
Wx := Px - x1;
Wy := Py - y1;
c1 := Vx * Wx + Vy * Wy;
c2 := Vx * Vx + Vy * Vy;
if c2=0 then exit;
Ratio := c1 / c2;
x := x1 + Ratio * Vx;
y := y1 + Ratio * Vy;
fail:=false;
end;


procedure linelc(x1,y1,x2,y2,x3,y3,r3:double;lalter:integer;
                  var xx1,yy1,xx2,yy2:double);
{ get a tangent line of the circle. this line is vertical or parallel
  to the given line .
  lalter=0,1 --> vertical to the given line.
  lalter=2,3 --> parallel to the given line.
  lalter=4   --> vertical line between line and circle.  }
var f:real;
 s1:integer;
 fail,change:boolean;
begin
s1:= lalter mod 5;
if (s1=0) or (s1=1) then begin
 lotlc(x1,y1,x2,y2,x3,y3,r3,s1,xx2,yy2);
 lotlc(x1,y1,x2,y2,x3,y3,r3,s1+2,xx1,yy1);
end else if (s1=2) or (s1=3) then begin
 s1:=s1-2; parallelelc(x1,y1,x2,y2,x3,y3,r3,s1,xx1,yy1,xx2,yy2);
end else begin
 lotpl2p(x3,y3,x1,y1,x2,y2,xx1,yy1,fail);
 f:=(distance2p(xx1,yy1,x2,y3)-r3)/r3;
 xx2:=(xx1+f*x3)/(1+f); yy2:=(yy1+f*y3)/(1+f);
end;
end;

procedure linepc(x1,y1,x4,y4,r:double;lalter:integer;
                  var xx1,yy1,xx2,yy2:double;var fail:boolean);
{tangent or lot line from a point onto a circle }
var f:real;
 si:integer;
begin
xx1:=x1;yy1:=y1;
si:=lalter mod 3;
if si<>2 then tangentepc(x1,y1,x4,y4,r,si,xx2,yy2,fail) else begin
 f:=(distance2p(x4,y4,x1,y1)-r)/r;
 xx2:=(x1+f*x4)/(1+f); yy2:=(y1+f*y4)/(1+f);
end;
end;
(*
procedure lotpc(x1,y1,x4,y4,r:double;var xx1,yy1,xx2,yy2:double;var fail:boolean);
{lot line from a point onto a circle }
var f:real;
 si:integer;
begin
xx1:=x1;yy1:=y1;
f:=(distance2p(x4,y4,x1,y1)-r)/r;
xx2:=(x1+f*x4)/(1+f); yy2:=(y1+f*y4)/(1+f);
end;
*)
procedure winkelhalbll(l1x1,l1y1,l1x2,l1y2,l2x1,l2y1,l2x2,l2y2:double;lalter:integer;
                 var x1,y1,x2,y2:double;var fail:boolean);
{ winkelhalbierende }
var a1,b1,c1,a2,b2,c2,d,dmax,a,b,c:double;
begin
line2p(l1x1,l1y1,l1x2,l1y2,a1,b1,c1,fail);
if fail then exit;
line2p(l2x1,l2y1,l2x2,l2y2,a2,b2,c2,fail);
if fail then exit;
point2l(a1,b1,c1,a2,b2,c2,x1,y1,fail);
if fail then exit;
dmax:=distance2p(x1,y1,l1x1,l1y1);
d:=distance2p(x1,y1,l1x2,l1y2);if d>dmax then dmax:=d;
d:=distance2p(x1,y1,l2x1,l2y1);if d>dmax then dmax:=d;
d:=distance2p(x1,y1,l2x2,l2y2);if d>dmax then dmax:=d;
winkelhalb(a1,b1,c1,a2,b2,c2,0,a,b,c,fail);
case lalter mod 4 of
 0:begin x2:=x1+dmax*a;y2:=y1+dmax*b;end;
 1:begin x2:=x1-dmax*b;y2:=y1+dmax*a;end;
 2:begin x2:=x1-dmax*a;y2:=y1-dmax*b;end;
 3:begin x2:=x1+dmax*b;y2:=y1-dmax*a;end;
end;
end;
procedure linell(l1x1,l1y1,l1x2,l1y2,l2x1,l2y1,l2x2,l2y2:double;lalter:integer;
                 var x1,y1,x2,y2:double);
{ get a line between two lines }
var
si:integer;
begin
si:=lalter mod 9;
if si = 0 then begin
    x1:=l1x1; y1:=l1y1; x2:=l2x1; y2:=l2y1;
end else if si = 1 then begin
    x1:=l1x1; y1:=l1y1; x2:=l2x2; y2:=l2y2;
end else if si = 2 then begin
    x1:=l1x2; y1:=l1y2; x2:=l2x1; y2:=l2y1;
end else if si = 3 then begin
    x1:=l1x2; y1:=l1y2; x2:=l2x2; y2:=l2y2; //0..3:Endpunkte verbunden
end else if si=4 then begin
    x1:=(l1x1+l1x2)/2; y1:=(l1y1+l1y2)/2;
    x2:=(l2x1+l2x2)/2; y2:=(l2y1+l2y2)/2; // Mittelpunkte
end else if si=5 then begin
    x1:=l1x1; y1:=l1y1;
    x2:=(l2x1+l2x2)/2; y2:=(l2y1+l2y2)/2;
end else if si=6 then begin
    x1:=l1x2; y1:=l1y2;
    x2:=(l2x1+l2x2)/2; y2:=(l2y1+l2y2)/2;
end else if si=7 then begin
    x1:=l2x1; y1:=l2y1;
    x2:=(l1x1+l1x2)/2; y2:=(l1y1+l1y2)/2;
end else if si=8 then begin
    x1:=l2x2; y1:=l2y2;
    x2:=(l1x1+l1x2)/2; y2:=(l1y1+l1y2)/2; // Mittel mit Endpunkt
end;
end;

function min(a,b:double):double;
{ liefert minimum zweier Zahlen}
begin
 if a<b then
 begin
  min:=a;
  exit;
 end else
  min:=b
end;

function max(x,y:double):double;
{ liefert maximum zweier Zahlen}
begin
 if x > y then max := x else max := y;
end;

function between(w,wmin,wmax:double):boolean;
begin
between:=(wmax<=wmin)xor((w-wmin)*(wmax-w)>=0.0);
end;

procedure extenarc(xm,ym,r,phimin,phimax:double;
          var xmin,ymin,xmax,ymax:double);
var sinphimin,sinphimax,cosphimin,cosphimax,xx:double;
begin
sinphimin:=(r*sin(phimin*radiant));
sinphimax:=(r*sin(phimax*radiant));
if sinphimin>sinphimax then begin
 xx:=sinphimin;sinphimin:=sinphimax;sinphimax:=xx;
end;
cosphimin:=(r*cos(phimin*radiant));
cosphimax:=(r*cos(phimax*radiant));
if cosphimin>cosphimax then begin
 xx:=cosphimin;cosphimin:=cosphimax;cosphimax:=xx;
end;
if between(90,phimin,phimax) then ymax:=ym+r
else ymax:=ym+sinphimax;
if between(270,phimin,phimax) then ymin:=ym-r
else ymin:=ym+sinphimin;
if between(0,phimin,phimax) then xmax:=xm+r
else xmax:=xm+cosphimax;
if between(180,phimin,phimax) then xmin:=xm-r
else xmin:=xm+cosphimin;
end;

function posauflinie(x1,y1,x2,y2,x,y:double;var posi:double):boolean;
// Wenn x,y auf der linie von x1,y1 nach x2,y2 liegt wird mit posi die Position
// angegeben (posi=0 x,y=x1,y1  posi=1 x,y=x2,y2)
// Wenn x,y nicht auf der geraden der Linie liegt ist result=false;
var a,b,c,d,dx,dy:double;
 fail:boolean;
begin
result:=false;
line2p(x1,y1,x2,y2,a,b,c,fail);
tracea(['posauflinie x1,y1,x2,y2=',x1,y1,x2,y2,' a,b,c=',a,b,c]);
if fail then exit;
d:=distancepl(x,y,a,b,c);
tracea(['posauflinie d=',d]);
if abs(d)<=tol then begin
 dx:=x2-x1;dy:=y2-y1;
 if abs(dx)>=abs(dy) then posi:=(x-x1)/dx else posi:=(y-y1)/dy;
 result:=true;
end else posi:=99;
end;

procedure fillet(x1,y1,x2,y2,x3,y3,x4,y4,r:double;
                 var xm,ym,xp1,yp1,xp2,yp2:double;var fail,vertauscht:boolean);
// Verundung von Linie x1,y1-x2,y2 und Linie x3,y3,x4,y4 mit Radius r;
// Ergebnis Bogen mit mittelpunkt xm,ym und von xp1,yp1 nach xp2,yp2
var a1,b1,c1,a2,b2,c2,dmin,d1,d2,z,posi1,posi2:double;
 xx,yy,xx1,yy1,xx2,yy2:array[0..3]of double;
 afail:array[0..3]of boolean;
 imin,i:integer;
begin
trace('fillet');
tracea(['line2p x1,y1,x2,y2=',x1,y1,x2,y2]);
line2p(x1,y1,x2,y2,a1,b1,c1,fail);
if fail then exit;
tracea(['line2p x3,y3,x4,y4=',x3,y3,x4,y4]);
line2p(x3,y3,x4,y4,a2,b2,c2,fail);
if fail then exit;
for i:=0 to 3 do begin
 circle2lr(a1,b1,c1,a2,b2,c2,r,i,xx[i],yy[i],afail[i]);
 trace('circle2lr i='+inttostr(i)+'fail='+booltostr(afail[i]));
 if not afail[i] then begin
  pointcl(xx[i],yy[i],r,a1,b1,c1,0,xx1[i],yy1[i],afail[i]);
  tracea(['point 1.Linie i=',i,' fail=',afail[i],'xx1=',xx1[i],' yy1=',yy1[i]]);
  pointcl(xx[i],yy[i],r,a2,b2,c2,0,xx2[i],yy2[i],afail[i]);
  tracea(['point 2.Linie i=',i,' fail=',afail[i],'xx2=',xx2[i],' yy2=',yy2[i]]);
 end;
end;
//Prüfen welche lösung auf den Linien liegt
dmin:=maxlongint;imin:=0;
for i:=0 to 3 do begin
 posauflinie(x1,y1,x2,y2,xx1[i],yy1[i],posi1);
 posauflinie(x3,y3,x4,y4,xx2[i],yy2[i],posi2);
 tracea(['i=',i,' posi1=',posi1,' posi2=',posi2]);
 if between(posi1,0,1)then d1:=0 else if posi1>1 then d1:=posi1-1
 else d1:=-posi1;
 if between(posi2,0,1)then d2:=0 else if posi2>1 then d2:=posi2-1
 else d2:=-posi2;
 tracea(['d1=',d1,' d2=',d2]);
 if dmin>d1+d2 then begin
  dmin:=d1+d2;imin:=i;
  tracea(['new dmin=',dmin]);
 end;
end;
xm:=xx[imin];ym:=yy[imin];
xp1:=xx1[imin];yp1:=yy1[imin];xp2:=xx2[imin];yp2:=yy2[imin];
vertauscht:=false;
if alpha3pdir(xm,ym,xp1,yp1,xp2,yp2)>180 then begin
 z:=xp1;xp1:=xp2;xp2:=z;z:=yp1;yp1:=yp2;yp2:=z;vertauscht:=true;
end;
end;

function diffphi(phi2,phi1:real):real;
// diffphi >0 Drehung um diffphi nach links führt von phi1 nach phi2
// diffphi <0 analog nach rechts
begin
result:=phi2-phi1;
if result>180 then result:=result-360
else if result<-180 then result:=result+360;
end;

function drehung(links:boolean;phi1,phi2:real):real;
//drehwinkel um links oder rechtsrum von phi1 nach phi2 zu kommen
var dphi:real;
begin
dphi:=phi2-phi1;
if links then begin
 if dphi<0 then dphi:=360+dphi;
end else begin
 if dphi>0 then dphi:=-360+dphi;
 dphi:=-dphi;
end;
result:=dphi;
end;

function normphi(phi:real):real;
// phi in denBereich 0 .. 360 bringen
begin
while phi>=360 do phi:=phi-360;
while phi<0 do phi:=phi+360;
result:=phi;
end;

function gegenwinkel(phi:real):real;
// gegenwinkel 90-> -270 oder -180 -> 180
begin
result:=-(360-abs(phi))*sign(round(phi))
end;

procedure gerade(d,alpha:double;var a,b,c:double);
// gerade in a,b,c
begin
c:=d;
alpha:=alpha*radiant;
a:=cos(alpha);b:=sin(alpha);
{
if (a<0)or((a=0)and(b<0)) then begin
 a:=-a;b:=-b;c:=-c;
end;
}
end;

procedure hesse2gerade(a,b,c:double;var d,alpha:double);
// a,b,c nach gerade
begin
d:=c;
alpha:=atan2(b,a)/radiant;
if alpha<0 then alpha:=360+alpha;
end;
procedure line2gerade(x1,y1,x2,y2:double;var c,alpha:double;var fail:boolean);
{Umwandlung Linie mit 2 Punkten in Punkt Gerade}
var a,b:double;
begin
line2p(x1,y1,x2,y2,a,b,c,fail);
if not fail then begin
end;
 alpha:=atan2(b,a)/radiant;
 if alpha<0 then alpha:=360+alpha;
end;
procedure mittelsenkrechte(x1,y1,x2,y2:double;var a,b,c:double);
var xm,ym,x3,y3:double;
 fail:boolean;
begin
xm:=(x1+x2)/2;ym:=(y1+y2)/2;
rotate90p(x2,y2,xm,ym,x3,y3);
line2p(xm,ym,x3,y3,a,b,c,fail);
end;

procedure tangente(x,y,xm,ym:double;var a,b,c:double);
// Tangente an Kreis mit Mittelpunkt xm,ym im Punkt x,y
var xx,yy:double;
 fail:boolean;
begin
rotate90p(xm,ym,x,y,xx,yy);
line2p(x,y,xx,yy,a,b,c,fail);
end;

procedure parallelep(a,b,c,x,y:double;var aa,bb,cc:double);
// Parallele zu a,b,c durch den Punkt x,y
var d:double;
begin
d:=distancepl(x,y,a,b,c);
aa:=a;bb:=b;cc:=c+d;
end;

procedure pointmittecl(x1,y1,r1,a2,b2,c2:double;var x,y:double);
var xx,yy,d,d2,f,dx,dy:double;
begin
lotpl(x1,y1,a2,b2,c2,xx,yy);
d:=distance2p(x1,y1,xx,yy);
dx:=xx-x1;dy:=yy-y1;
d2:=(d+d-r1)/2; //Abstand des mittenpunktes
f:=d2/d;
x:=x1+dx*f;y:=y1+dy*f;
end;

procedure punktauflinie(x1,y1,x2,y2,d:double;var xx,yy:double);
// Punkt auf x1,y1-x2,y2 mit abstand d von x1,y2
var d2,dx,dy,f:double;
begin
dx:=x2-x1;dy:=y2-y1;
d2:=sqrt(sqr(x1-x2)+sqr(y1-y2));
f:=d/d2;
xx:=x1+dx*f;yy:=y1+dy*f;
end;

procedure matrixmul(a,b:tmatrix;var c:tmatrix);
var i,j,k:integer;
begin
for i:=1 to 2 do for j:=1 to 2 do begin
 c[i,j]:=0;
 for k:=1 to 2 do c[i,j]:=c[i,j]+a[i,k]*b[k,j];
end;
end;
procedure matvecmul(b:tmatrix;a:tvector;var c:tvector);
var i,k:integer;
begin
for i:=1 to 2 do begin
 c[i]:=0;
 for k:=1 to 2 do c[i]:=c[i]+a[k]*b[i,k];
end;
end;
procedure setmat(var mat:tmatrix;x11,x12,x21,x22:double);
begin
mat[1,1]:=x11;mat[1,2]:=x12;mat[2,1]:=x21;mat[2,2]:=x22;
end;

procedure linepw(x,y,alpha:double;var a,b,c:double);
{Linie aus Punkt und Richtungswinkel}
var sina,cosa:double;
begin
sina:=sin(alpha*radiant);cosa:=cos(alpha*radiant);
a:=sina;b:=-cosa;c:=x*sina-y*cosa;
end;

function parallel2l(a1,b1,c1,a2,b2,c2:double;var gleich:boolean):boolean;
// 2 Linien sind parallel
// gleich: die beiden Linien sind gleich
begin
gleich:=false;
result:=abs(a1*b2-a2*b1)<0.01;
if result then begin
 if abs(a1)>abs(b1) then begin
  if a1<0 then c1:=-c1;
  if a2<0 then c2:=-c2;
 end else begin
  if b1<0 then c1:=-c1;
  if b2<0 then c2:=-c2;
 end;
 gleich:=abs(c1-c2)<0.5;
end;
end;

procedure pointonline(x0,y0,x1,y1,d:double;var x,y:double);
// x,y liegt auf de Linie x0,y0-x1,y1 im Abstand d von x0,y0
var f,dx,dy:double;
begin
dx:=x1-x0;dy:=y1-y0;
f:=d/sqrt(dx*dx+dy*dy);
dx:=dx*f;dy:=dy*f;
x:=x0+dx;y:=y0+dy;
end;

procedure senkrechtedurchPunkt(a,b,c,x,y:double;var aa,bb,cc:double);
var xx,yy:double;
begin
lotpl(x,y,a,b,c,xx,yy);
aa:=b;bb:=a;cc:=xx*b+yy*a;
end;

procedure circleptx(x1,y1,a2,b2,c2,x3:double;mode:integer;var x,y,r:double;var fail:boolean);
var a,b,c,aa,bb,cc,xx1,yy,yy1,yy2,xx2:double;
 d,t1,t2,t3:double;
begin
fail:=false;
if a2<0 then begin
 a2:=-a2;b2:=-b2;c2:=-c2;
end;
x:=x3;
d:=-distancepl(x1,y1,a2,b2,c2);
if abs(d)<tol then begin
 // Punkt liegt auf Linie
 senkrechtedurchPunkt(a2,b2,c2,x1,y1,aa,bb,cc);
 point2l(1,0,x3,aa,bb,cc,x,y,fail);
 r:=distance2p(x,y,x1,y1);
end else begin
 // p in Nullpunkt
 x3:=x3-x1;
 a:=a2;b:=b2;c:=d;
 if abs(b)<>1 then begin
  t1:=c - (Power(b,2)*c)/(-1 + Power(b,2)) - a*x3 + (a*Power(b,2)*x3)/(-1 + Power(b,2));
  t2:=Power(c,2) - 2*a*c*x3 - Power(x3,2) + Power(a,2)*Power(x3,2) + Power(b,2)*Power(x3,2);
  t3:=b*sqrt(t2);
  if t2<-0.01 then begin
   fail:=true;exit;
  end else if t2<0 then t2:=0;
  if mode=0 then begin
   r:=abs(t1-t3/(b*b-1));
   yy:=(b*c - a*b*x3 + Sqrt(t2))/(-1 + Power(b,2));
  end else begin
   r:=abs(t1+t3/(b*b-1));
   yy:=(b*c - a*b*x3 - Sqrt(t2))/(-1 + Power(b,2));
  end;
  y:=yy+y1;
 end else begin
  r:=abs((Power(c,2) + Power(x3,2))/(2.*c));
  if b=1 then yy:=c-r else yy:=c+r;
  y:=yy+y1;
 end;
end;
end;
procedure circlepty(x1,y1,a2,b2,c2,y3:double;mode:integer;var x,y,r:double;var fail:boolean);
var a,b,c,aa,bb,cc,xx:double;
 d,t1,t2,t3:double;
begin
fail:=false;
if a2<0 then begin
 a2:=-a2;b2:=-b2;c2:=-c2;
end;
y:=y3;
d:=-distancepl(x1,y1,a2,b2,c2);
if abs(d)<tol then begin
 // Punkt liegt auf Linie
 senkrechtedurchPunkt(a2,b2,c2,x1,y1,aa,bb,cc);
 point2l(0,1,y3,aa,bb,cc,x,y,fail);
 r:=distance2p(x,y,x1,y1);
end else begin
 // p in Nullpunkt
 y3:=y3-y1;
 a:=a2;b:=b2;c:=d;
 if abs(a)<>1 then begin
  t1:=c - (a*a*c)/(-1 + a*a) - b*y3 + (a*a*b*y3)/(-1 + a*a);
  t2:=c*c - 2*b*c*y3 - y3*y3 + a*a*y3*y3 + b*b*y3*y3;
  if t2<-0.01 then begin
   fail:=true;exit;
  end else if t2<0 then t2:=0;
  t3:=a*sqrt(t2);
  if mode=0 then begin
   r:=abs(t1-t3/(a*a-1));
   xx:=(a*c - a*b*y3 + Sqrt(t2))/(-1 + a*a);
  end else begin
   r:=abs(t1+t3/(a*a-1));
   xx:=(a*c - a*b*y3 - Sqrt(t2))/(-1 + a*a);
  end;
 end else begin
  r:=abs((Power(c,2) + Power(y3,2))/(2.*c));
  if a=1 then xx:=c-r else xx:=c+r;
 end;
 x:=xx+x1;
end;
end;

procedure circlectx(x1,y1,r1,a2,b2,c2,x3:double;mode:integer;var x,y,r:double;var fail:boolean);
var a,b,c,aa,bb,cc:double;
 d,rr:double;
begin
fail:=true;
if a2<0 then begin
 a2:=-a2;b2:=-b2;c2:=-c2;
end;
x:=x3;
d:=-distancepl(x1,y1,a2,b2,c2);
if abs(d)<tol then begin
 // mittelPunkt liegt auf Linie
 senkrechtedurchPunkt(a2,b2,c2,x1,y1,aa,bb,cc);
 point2l(1,0,x3,aa,bb,cc,x,y,fail);
 r:=distance2p(x,y,x1,y1)-r1;
 if r<0 then exit;
 fail:=false;exit;
end else begin
 if mode<2 then begin
  if d>0 then c2:=c2+r1 else c2:=c2-r1;
  circleptx(x1,y1,a2,b2,c2,x3,mode mod 2,x,y,rr,fail);
  if fail then exit;
  r:=abs(rr)-r1;
 end else begin
  if d<0 then c2:=c2+r1 else c2:=c2-r1;
  circleptx(x1,y1,a2,b2,c2,x3,mode mod 2,x,y,rr,fail);
  if fail then exit;
  r:=abs(rr)+r1;
 end;
 {
 if c2>=0 then c2:=c2+r1 else c2:=c2-r1;
 circleptx(x1,y1,a2,b2,c2,x3,mode,x,y,rr,fail);
 if fail then exit;
 r:=abs(rr)-r1;
 fail:=false;
 }
end;
end;

procedure circlecty(x1,y1,r1,a2,b2,c2,y3:double;mode:integer;var x,y,r:double;var fail:boolean);
var a,b,c,d,aa,bb,cc,xx1,yy1,yy2,xx2:double;
 rr:double;
begin
fail:=true;r:=0;
if a2<0 then begin
 a2:=-a2;b2:=-b2;c2:=-c2;
end;
y:=y3;
d:=distancepl(x1,y1,a2,b2,c2);
if abs(d)<tol then begin
 // KreismittelPunkt liegt auf Linie
 senkrechtedurchPunkt(a2,b2,c2,x1,y1,aa,bb,cc);
 point2l(0,1,y3,aa,bb,cc,x,y,fail);
 r:=distance2p(x,y,x1,y1)-r1;
 if r<0 then exit;
 fail:=false;exit;
end else begin
 if mode<2 then begin
  if d<0 then c2:=c2+r1 else c2:=c2-r1;
  circlepty(x1,y1,a2,b2,c2,y3,mode mod 2,x,y,rr,fail);
  if fail then exit;
  r:=abs(rr)-r1;
 end else begin
  if d>0 then c2:=c2+r1 else c2:=c2-r1;
  circlepty(x1,y1,a2,b2,c2,y3,mode mod 2,x,y,rr,fail);
  if fail then exit;
  r:=abs(rr)+r1;
 end;
end;
end;
procedure circleptl(x1,y1,a1,b1,c1,a2,b2,c2:double;mode:integer;var x,y,r:double;var fail:boolean);
var d1,d2,aa,bb,cc:double;
begin
fail:=false;
if a1<0 then begin
 a1:=-a1;b1:=-b1;c1:=-c1;
end;
if a2<0 then begin
 a2:=-a2;b2:=-b2;c2:=-c2;
end;
d1:=-distancepl(x1,y1,a1,b1,c1);
d2:=-distancepl(x1,y1,a2,b2,c2);
if abs(d1)<tol then begin
 // Punkt liegt auf Linie
 senkrechtedurchPunkt(a1,b1,c1,x1,y1,aa,bb,cc);
 point2l(a2,b2,c2,aa,bb,cc,x,y,fail);
 r:=distance2p(x,y,x1,y1);
end else begin
 // p in Nullpunkt
 c1:=d1;
 c2:=d2;
 if abs(b1)<>1 then begin
  if mode=0 then begin
   x:=(-((Power(a2,2)*b1*b2*c1)/(-Power(a2,2) + Power(a2,2)*Power(b1,2) - 2*a1*a2*b1*b2 - Power(b2,2) + Power(a1,2)*Power(b2,2))) +
     (a1*a2*Power(b2,2)*c1)/(-Power(a2,2) + Power(a2,2)*Power(b1,2) - 2*a1*a2*b1*b2 - Power(b2,2) + Power(a1,2)*Power(b2,2)) + c2 +
     (a1*a2*b1*b2*c2)/(-Power(a2,2) + Power(a2,2)*Power(b1,2) - 2*a1*a2*b1*b2 - Power(b2,2) + Power(a1,2)*Power(b2,2)) +
     (Power(b2,2)*c2)/(-Power(a2,2) + Power(a2,2)*Power(b1,2) - 2*a1*a2*b1*b2 - Power(b2,2) + Power(a1,2)*Power(b2,2)) -
     (Power(a1,2)*Power(b2,2)*c2)/(-Power(a2,2) + Power(a2,2)*Power(b1,2) - 2*a1*a2*b1*b2 - Power(b2,2) + Power(a1,2)*Power(b2,2)) -
     (a2*b2*Sqrt(Power(a2,2)*Power(c1,2) + Power(b2,2)*Power(c1,2) - 2*a1*a2*c1*c2 - 2*b1*b2*c1*c2 - Power(c2,2) + Power(a1,2)*Power(c2,2) +
      Power(b1,2)*Power(c2,2)))/(-Power(a2,2) + Power(a2,2)*Power(b1,2) - 2*a1*a2*b1*b2 - Power(b2,2) + Power(a1,2)*Power(b2,2)))/a2;
   y:=(Power(a2,2)*b1*c1 - a1*a2*b2*c1 - a1*a2*b1*c2 - b2*c2 + Power(a1,2)*b2*c2 +
        a2*Sqrt(Power(a2,2)*Power(c1,2) + Power(b2,2)*Power(c1,2) - 2*a1*a2*c1*c2 - 2*b1*b2*c1*c2 - Power(c2,2) + Power(a1,2)*Power(c2,2) + Power(b1,2)*Power(c2,2)))/
      (-Power(a2,2) + Power(a2,2)*Power(b1,2) - 2*a1*a2*b1*b2 - Power(b2,2) + Power(a1,2)*Power(b2,2));
  end else begin
   x:=(-((Power(a2,2)*b1*b2*c1)/(-Power(a2,2) + Power(a2,2)*Power(b1,2) - 2*a1*a2*b1*b2 - Power(b2,2) + Power(a1,2)*Power(b2,2))) +
      (a1*a2*Power(b2,2)*c1)/(-Power(a2,2) + Power(a2,2)*Power(b1,2) - 2*a1*a2*b1*b2 - Power(b2,2) + Power(a1,2)*Power(b2,2)) + c2 +
      (a1*a2*b1*b2*c2)/(-Power(a2,2) + Power(a2,2)*Power(b1,2) - 2*a1*a2*b1*b2 - Power(b2,2) + Power(a1,2)*Power(b2,2)) +
      (Power(b2,2)*c2)/(-Power(a2,2) + Power(a2,2)*Power(b1,2) - 2*a1*a2*b1*b2 - Power(b2,2) + Power(a1,2)*Power(b2,2)) -
      (Power(a1,2)*Power(b2,2)*c2)/(-Power(a2,2) + Power(a2,2)*Power(b1,2) - 2*a1*a2*b1*b2 - Power(b2,2) + Power(a1,2)*Power(b2,2)) +
      (a2*b2*Sqrt(Power(a2,2)*Power(c1,2) + Power(b2,2)*Power(c1,2) - 2*a1*a2*c1*c2 - 2*b1*b2*c1*c2 - Power(c2,2) + Power(a1,2)*Power(c2,2) +
      Power(b1,2)*Power(c2,2)))/(-Power(a2,2) + Power(a2,2)*Power(b1,2) - 2*a1*a2*b1*b2 - Power(b2,2) + Power(a1,2)*Power(b2,2)))/a2;
   y:=(Power(a2,2)*b1*c1 - a1*a2*b2*c1 - a1*a2*b1*c2 - b2*c2 + Power(a1,2)*b2*c2 -
      a2*Sqrt(Power(a2,2)*Power(c1,2) + Power(b2,2)*Power(c1,2) - 2*a1*a2*c1*c2 - 2*b1*b2*c1*c2 - Power(c2,2) + Power(a1,2)*Power(c2,2) + Power(b1,2)*Power(c2,2)))/
      (-Power(a2,2) + Power(a2,2)*Power(b1,2) - 2*a1*a2*b1*b2 - Power(b2,2) + Power(a1,2)*Power(b2,2));
  end;
  x:=x+x1;y:=y+y1;
  r:=distance2p(x1,y1,x,y);
 end else begin
  // b1=1 d.h. Waggrechte mit y=b2*c2;
  circlepty(x1,y1,a1,b1,c1,b2*c2,mode,x,y,r,fail);
 end;
end;
end;

procedure circle0yc(y1,x2,y2,r2:double;mode:integer;var x,y,r:double;var fail:boolean);
begin
try
case mode of
 0:y:=y1 + y2 - Sqrt(4*Power(y1 + y2,2) + (4*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)))/3. +
         (16*Power(2,0.3333333333333333)*(Power(r2,4) - 2*Power(r2,2)*Power(x2,2) + Power(x2,4) + 3*Power(x2,2)*Power(y1,2) + 4*Power(r2,2)*y1*y2 +
              2*Power(x2,2)*y1*y2 + 4*Power(y1,2)*Power(y2,2)))/
          (3.*Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) +
              1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 +
              1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) +
              Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) -
                442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) +
                884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 -
                442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 +
                1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) -
                442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) -
                442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) +
                1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)) +
         Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) +
            1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 +
            1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) +
            Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) -
              442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) +
              884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 -
              442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 +
              1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) -
              442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) -
              442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) +
              1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)/(3.*Power(2,0.3333333333333333)))/2. +
      Sqrt(8*Power(y1 + y2,2) + (8*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)))/3. -
         (16*Power(2,0.3333333333333333)*(Power(r2,4) - 2*Power(r2,2)*Power(x2,2) + Power(x2,4) + 3*Power(x2,2)*Power(y1,2) + 4*Power(r2,2)*y1*y2 +
              2*Power(x2,2)*y1*y2 + 4*Power(y1,2)*Power(y2,2)))/
          (3.*Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) +
              1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 +
              1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) +
              Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) -
                442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) +
                884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 -
                442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 +
                1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) -
                442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) -
                442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) +
                1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)) -
         Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) +
            1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 +
            1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) +
            Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) -
              442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) +
              884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 -
              442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 +
              1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) -
              442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) -
              442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) +
              1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)/(3.*Power(2,0.3333333333333333)) -
         (64*Power(y1 + y2,3) + 32*(y1 + y2)*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)) -
            32*(Power(r2,2)*y1 + Power(x2,2)*y1 - Power(y1,3) + Power(r2,2)*y2 - Power(x2,2)*y2 - Power(y1,2)*y2 - y1*Power(y2,2) - Power(y2,3)))/
          (4.*Sqrt(4*Power(y1 + y2,2) + (4*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)))/3. +
              (16*Power(2,0.3333333333333333)*(Power(r2,4) - 2*Power(r2,2)*Power(x2,2) + Power(x2,4) + 3*Power(x2,2)*Power(y1,2) + 4*Power(r2,2)*y1*y2 +
                   2*Power(x2,2)*y1*y2 + 4*Power(y1,2)*Power(y2,2)))/
               (3.*Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) +
                   1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 +
                   1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) +
                   Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) -
                     442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) +
                     884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 -
                     442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 +
                     1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) +
                     884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 442368*Power(x2,8)*Power(y1,2)*Power(y2,2) +
                     3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) -
                     442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) -
                     1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) + 1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),
                  0.3333333333333333)) + Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) +
                 576*Power(r2,2)*Power(x2,2)*Power(y1,2) + 1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 -
                 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 + 1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) +
                 1024*Power(y1,3)*Power(y2,3) + Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) +
                   1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) - 442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) +
                   2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) + 884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) +
                   2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 -
                   442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 + 1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) +
                   884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 442368*Power(x2,8)*Power(y1,2)*Power(y2,2) +
                   3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) - 442368*Power(x2,4)*Power(y1,6)*Power(y2,2) -
                   1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) +
                   1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)/(3.*Power(2,0.3333333333333333)))))/2.
      ;
1:y:=y1 + y2 - Sqrt(4*Power(y1 + y2,2) + (4*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)))/3. +
         (16*Power(2,0.3333333333333333)*(Power(r2,4) - 2*Power(r2,2)*Power(x2,2) + Power(x2,4) + 3*Power(x2,2)*Power(y1,2) + 4*Power(r2,2)*y1*y2 +
              2*Power(x2,2)*y1*y2 + 4*Power(y1,2)*Power(y2,2)))/
          (3.*Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) +
              1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 +
              1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) +
              Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) -
                442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) +
                884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 -
                442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 +
                1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) -
                442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) -
                442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) +
                1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)) +
         Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) +
            1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 +
            1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) +
            Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) -
              442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) +
              884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 -
              442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 +
              1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) -
              442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) -
              442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) +
              1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)/(3.*Power(2,0.3333333333333333)))/2. -
      Sqrt(8*Power(y1 + y2,2) + (8*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)))/3. -
         (16*Power(2,0.3333333333333333)*(Power(r2,4) - 2*Power(r2,2)*Power(x2,2) + Power(x2,4) + 3*Power(x2,2)*Power(y1,2) + 4*Power(r2,2)*y1*y2 +
              2*Power(x2,2)*y1*y2 + 4*Power(y1,2)*Power(y2,2)))/
          (3.*Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) +
              1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 +
              1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) +
              Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) -
                442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) +
                884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 -
                442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 +
                1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) -
                442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) -
                442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) +
                1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)) -
         Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) +
            1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 +
            1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) +
            Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) -
              442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) +
              884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 -
              442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 +
              1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) -
              442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) -
              442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) +
              1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)/(3.*Power(2,0.3333333333333333)) -
         (64*Power(y1 + y2,3) + 32*(y1 + y2)*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)) -
            32*(Power(r2,2)*y1 + Power(x2,2)*y1 - Power(y1,3) + Power(r2,2)*y2 - Power(x2,2)*y2 - Power(y1,2)*y2 - y1*Power(y2,2) - Power(y2,3)))/
          (4.*Sqrt(4*Power(y1 + y2,2) + (4*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)))/3. +
              (16*Power(2,0.3333333333333333)*(Power(r2,4) - 2*Power(r2,2)*Power(x2,2) + Power(x2,4) + 3*Power(x2,2)*Power(y1,2) + 4*Power(r2,2)*y1*y2 +
                   2*Power(x2,2)*y1*y2 + 4*Power(y1,2)*Power(y2,2)))/
               (3.*Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) +
                   1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 +
                   1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) +
                   Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) -
                     442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) +
                     884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 -
                     442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 +
                     1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) +
                     884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 442368*Power(x2,8)*Power(y1,2)*Power(y2,2) +
                     3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) -
                     442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) -
                     1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) + 1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),
                  0.3333333333333333)) + Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) +
                 576*Power(r2,2)*Power(x2,2)*Power(y1,2) + 1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 -
                 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 + 1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) +
                 1024*Power(y1,3)*Power(y2,3) + Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) +
                   1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) - 442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) +
                   2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) + 884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) +
                   2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 -
                   442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 + 1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) +
                   884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 442368*Power(x2,8)*Power(y1,2)*Power(y2,2) +
                   3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) - 442368*Power(x2,4)*Power(y1,6)*Power(y2,2) -
                   1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) +
                   1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)/(3.*Power(2,0.3333333333333333)))))/2.;
2:y:=y1 + y2 + Sqrt(4*Power(y1 + y2,2) + (4*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)))/3. +
         (16*Power(2,0.3333333333333333)*(Power(r2,4) - 2*Power(r2,2)*Power(x2,2) + Power(x2,4) + 3*Power(x2,2)*Power(y1,2) + 4*Power(r2,2)*y1*y2 +
              2*Power(x2,2)*y1*y2 + 4*Power(y1,2)*Power(y2,2)))/
          (3.*Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) + 
              1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 + 
              1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) + 
              Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) - 
                442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) + 
                884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 - 
                442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 + 
                1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 
                442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) - 
                442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) + 
                1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)) + 
         Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) + 
            1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 + 
            1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) + 
            Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) - 
              442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) + 
              884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 - 
              442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 + 
              1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 
              442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) - 
              442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) + 
              1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)/(3.*Power(2,0.3333333333333333)))/2. + 
      Sqrt(8*Power(y1 + y2,2) + (8*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)))/3. - 
         (16*Power(2,0.3333333333333333)*(Power(r2,4) - 2*Power(r2,2)*Power(x2,2) + Power(x2,4) + 3*Power(x2,2)*Power(y1,2) + 4*Power(r2,2)*y1*y2 + 
              2*Power(x2,2)*y1*y2 + 4*Power(y1,2)*Power(y2,2)))/
          (3.*Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) + 
              1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 + 
              1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) + 
              Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) - 
                442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) + 
                884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 -
                442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 +
                1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 
                442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) - 
                442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) +
                1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)) - 
         Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) + 
            1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 + 
            1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) + 
            Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) -
              442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) + 
              884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 - 
              442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 + 
              1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 
              442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) - 
              442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) + 
              1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)/(3.*Power(2,0.3333333333333333)) + 
         (64*Power(y1 + y2,3) + 32*(y1 + y2)*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)) - 
            32*(Power(r2,2)*y1 + Power(x2,2)*y1 - Power(y1,3) + Power(r2,2)*y2 - Power(x2,2)*y2 - Power(y1,2)*y2 - y1*Power(y2,2) - Power(y2,3)))/
          (4.*Sqrt(4*Power(y1 + y2,2) + (4*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)))/3. + 
              (16*Power(2,0.3333333333333333)*(Power(r2,4) - 2*Power(r2,2)*Power(x2,2) + Power(x2,4) + 3*Power(x2,2)*Power(y1,2) + 4*Power(r2,2)*y1*y2 + 
                   2*Power(x2,2)*y1*y2 + 4*Power(y1,2)*Power(y2,2)))/
               (3.*Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) + 
                   1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 + 
                   1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) + 
                   Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) - 
                     442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) + 
                     884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 - 
                     442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 + 
                     1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 
                     884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 
                     3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) - 
                     442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 
                     1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) + 1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),
                  0.3333333333333333)) + Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 
                 576*Power(r2,2)*Power(x2,2)*Power(y1,2) + 1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 
                 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 + 1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 
                 1024*Power(y1,3)*Power(y2,3) + Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 
                   1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) - 442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 
                   2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) + 884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 
                   2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 
                   442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 + 1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 
                   884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 442368*Power(x2,8)*Power(y1,2)*Power(y2,2) +
                   3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) - 442368*Power(x2,4)*Power(y1,6)*Power(y2,2) -
                   1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) + 
                   1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)/(3.*Power(2,0.3333333333333333)))))/2.
      ;
3:  y:=y1 + y2 + Sqrt(4*Power(y1 + y2,2) + (4*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)))/3. +
         (16*Power(2,0.3333333333333333)*(Power(r2,4) - 2*Power(r2,2)*Power(x2,2) + Power(x2,4) + 3*Power(x2,2)*Power(y1,2) + 4*Power(r2,2)*y1*y2 + 
              2*Power(x2,2)*y1*y2 + 4*Power(y1,2)*Power(y2,2)))/
          (3.*Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) + 
              1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 + 
              1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) + 
              Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) - 
                442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) + 
                884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 - 
                442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 + 
                1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 
                442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) - 
                442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) + 
                1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)) + 
         Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) + 
            1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 + 
            1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) + 
            Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) - 
              442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) + 
              884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 - 
              442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 + 
              1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 
              442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) - 
              442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) + 
              1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)/(3.*Power(2,0.3333333333333333)))/2. - 
      Sqrt(8*Power(y1 + y2,2) + (8*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)))/3. - 
         (16*Power(2,0.3333333333333333)*(Power(r2,4) - 2*Power(r2,2)*Power(x2,2) + Power(x2,4) + 3*Power(x2,2)*Power(y1,2) + 4*Power(r2,2)*y1*y2 + 
              2*Power(x2,2)*y1*y2 + 4*Power(y1,2)*Power(y2,2)))/
          (3.*Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) + 
              1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 + 
              1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) + 
              Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) - 
                442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) + 
                884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 - 
                442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 + 
                1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 
                442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) - 
                442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) +
                1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)) -
         Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) + 
            1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 + 
            1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) + 
            Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) -
              442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) + 
              884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 - 
              442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 + 
              1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 
              442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) - 
              442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) + 
              1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)/(3.*Power(2,0.3333333333333333)) + 
         (64*Power(y1 + y2,3) + 32*(y1 + y2)*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)) - 
            32*(Power(r2,2)*y1 + Power(x2,2)*y1 - Power(y1,3) + Power(r2,2)*y2 - Power(x2,2)*y2 - Power(y1,2)*y2 - y1*Power(y2,2) - Power(y2,3)))/
          (4.*Sqrt(4*Power(y1 + y2,2) + (4*(Power(r2,2) - Power(x2,2) - 3*Power(y1,2) - 4*y1*y2 - 3*Power(y2,2)))/3. + 
              (16*Power(2,0.3333333333333333)*(Power(r2,4) - 2*Power(r2,2)*Power(x2,2) + Power(x2,4) + 3*Power(x2,2)*Power(y1,2) + 4*Power(r2,2)*y1*y2 + 
                   2*Power(x2,2)*y1*y2 + 4*Power(y1,2)*Power(y2,2)))/
               (3.*Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 576*Power(r2,2)*Power(x2,2)*Power(y1,2) + 
                   1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 + 
                   1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 1024*Power(y1,3)*Power(y2,3) + 
                   Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) - 
                     442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) + 
                     884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 - 
                     442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 + 
                     1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 
                     884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 
                     3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) - 
                     442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 
                     1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) + 1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),
                  0.3333333333333333)) + Power(128*Power(r2,6) - 384*Power(r2,4)*Power(x2,2) + 384*Power(r2,2)*Power(x2,4) - 128*Power(x2,6) + 
                 576*Power(r2,2)*Power(x2,2)*Power(y1,2) + 1152*Power(x2,4)*Power(y1,2) + 768*Power(r2,4)*y1*y2 - 384*Power(r2,2)*Power(x2,2)*y1*y2 - 
                 384*Power(x2,4)*y1*y2 + 1152*Power(x2,2)*Power(y1,3)*y2 + 1536*Power(r2,2)*Power(y1,2)*Power(y2,2) + 768*Power(x2,2)*Power(y1,2)*Power(y2,2) + 
                 1024*Power(y1,3)*Power(y2,3) + Sqrt(442368*Power(r2,6)*Power(x2,4)*Power(y1,2) - 1327104*Power(r2,4)*Power(x2,6)*Power(y1,2) + 
                   1327104*Power(r2,2)*Power(x2,8)*Power(y1,2) - 442368*Power(x2,10)*Power(y1,2) - 110592*Power(r2,4)*Power(x2,4)*Power(y1,4) + 
                   2211840*Power(r2,2)*Power(x2,6)*Power(y1,4) + 884736*Power(x2,8)*Power(y1,4) - 442368*Power(x2,6)*Power(y1,6) + 
                   2211840*Power(r2,4)*Power(x2,4)*Power(y1,3)*y2 - 442368*Power(r2,2)*Power(x2,6)*Power(y1,3)*y2 - 1769472*Power(x2,8)*Power(y1,3)*y2 - 
                   442368*Power(r2,2)*Power(x2,4)*Power(y1,5)*y2 + 1769472*Power(x2,6)*Power(y1,5)*y2 - 442368*Power(r2,4)*Power(x2,4)*Power(y1,2)*Power(y2,2) + 
                   884736*Power(r2,2)*Power(x2,6)*Power(y1,2)*Power(y2,2) - 442368*Power(x2,8)*Power(y1,2)*Power(y2,2) + 
                   3538944*Power(r2,2)*Power(x2,4)*Power(y1,4)*Power(y2,2) - 884736*Power(x2,6)*Power(y1,4)*Power(y2,2) - 442368*Power(x2,4)*Power(y1,6)*Power(y2,2) - 
                   1769472*Power(r2,2)*Power(x2,4)*Power(y1,3)*Power(y2,3) - 1769472*Power(x2,6)*Power(y1,3)*Power(y2,3) + 
                   1769472*Power(x2,4)*Power(y1,5)*Power(y2,3) - 1769472*Power(x2,4)*Power(y1,4)*Power(y2,4)),0.3333333333333333)/(3.*Power(2,0.3333333333333333)))))/2.
      ;
 end;
except
fail:=true;exit;
end;
r:=y1-y;
x:=sqrt(r*r-y*y);
fail:=abs(sqr(r2)-sqr(x-x2)-sqr(y-y2))>0.1;
end;
procedure circle0tc(a1,b1,c1,x2,y2,r2:double;mode:integer;var x,y,r:double;var fail:boolean);
var d,xx1,yy1,xx2,yy2,xx,yy,phi:double;
begin
d:=distancepl(0,0,a1,b1,c1);
lotpl(0,0,a1,b1,c1,xx1,yy1);
phi:=getphi(0,0,xx1,yy1);
rotatep(x2,y2,0,0,90-phi,xx2,yy2);
circle0yc(-d,xx2,yy2,r2,mode,xx,yy,r,fail);
rotatep(xx,yy,0,0,phi-90,x,y);
end;

procedure circleptc(x1,y1,a1,b1,c1,x2,y2,r2:double;mode:integer;var x,y,r:double;var fail:boolean);
begin
// x1,y1 in Nullpunkt;
x2:=x2-x1;y2:=y2-y1;
c1:=-distancepl(x1,y1,a1,b1,c1);
circle0tc(a1,b1,c1,x2,y2,r2,mode,x,y,r,fail);
x:=x+x1;y:=y+y1;
end;
procedure circlecpmx(x1,y1,r1,x2,y2,x3:double;mode:integer;var x,y,r:double;var fail:boolean);
{Kreis tangential an Kreis, Punkt und x Mittelpunkt}
begin
if mode=0 then begin
y:=(-4*Power(r1,2)*y1 - 8*x*x1*y1 + 4*Power(x1,2)*y1 + 8*x*x2*y1 - 4*Power(x2,2)*y1 + 4*Power(y1,3) - 4*Power(r1,2)*y2 + 8*x*x1*y2 -
 4*Power(x1,2)*y2 - 8*x*x2*y2 + 4*Power(x2,2)*y2 - 4*Power(y1,2)*y2 - 4*y1*Power(y2,2) + 4*Power(y2,3) -
 4*r1*Sqrt(Power(r1,2) - 4*Power(x,2) + 4*x*x1 - Power(x1,2) + 4*x*x2 - 2*x1*x2 - Power(x2,2) - Power(y1,2) + 2*y1*y2 - Power(y2,2))*
 Sqrt(Power(r1,2) - Power(x1,2) + 2*x1*x2 - Power(x2,2) - Power(y1,2) + 2*y1*y2 - Power(y2,2)))
 /(2.*(-4*Power(r1,2) + 4*Power(y1,2) - 8*y1*y2 + 4*Power(y2,2)));
end else begin
y:=(-4*Power(r1,2)*y1 - 8*x*x1*y1 + 4*Power(x1,2)*y1 + 8*x*x2*y1 - 4*Power(x2,2)*y1 + 4*Power(y1,3) - 4*Power(r1,2)*y2 + 8*x*x1*y2 -
 4*Power(x1,2)*y2 - 8*x*x2*y2 + 4*Power(x2,2)*y2 - 4*Power(y1,2)*y2 - 4*y1*Power(y2,2) + 4*Power(y2,3) +
 4*r1*Sqrt(Power(r1,2) - 4*Power(x,2) + 4*x*x1 - Power(x1,2) + 4*x*x2 - 2*x1*x2 - Power(x2,2) - Power(y1,2) + 2*y1*y2 - Power(y2,2))*
 Sqrt(Power(r1,2) - Power(x1,2) + 2*x1*x2 - Power(x2,2) - Power(y1,2) + 2*y1*y2 - Power(y2,2)))
 /(2.*(-4*Power(r1,2) + 4*Power(y1,2) - 8*y1*y2 + 4*Power(y2,2)));
end;
x:=x3;
r:=distance2p(x,y,x2,y2);
end;
procedure circlecpmy(x1,y1,r1,x2,y2,y3:double;mode:integer;var x,y,r:double;var fail:boolean);
{Kreis tangential an Kreis, Punkt und x Mittelpunkt}
begin
circlecpmx(y1,x1,r1,y2,x2,y3,mode,y,x,r,fail);
end;
procedure circleclmc(x1,y1,r1,a2,b2,c2,x3,y3,r3:double;mode:integer;
 var x,y,r:double;var fail:boolean);
begin
end;
procedure circleclml(x1,y1,r1,a2,b2,c2,a3,b3,c3:double;mode:integer;
 var x,y,r:double;var fail:boolean);
begin
end;

procedure circle2cmc(x1,y1,r1,x2,y2,r2,x3,y3,r3,rold:double;mode:integer;
 var x,y,r:double;var fail:boolean);
{Kreis tangential an 2 Kreise und Mittelpunkt auf Kreis}
var rmin,rmax,rmitte,xx1,yy1,xx2,yy2,xx3,yy3,d1,d2,d3:double;
begin
rmin:=distance2p(x1,y1,x2,y2)-r1-r2;
if rmin<0 then rmin:=0;
rmax:=rold*2;
repeat
 xcircle2cr(x1,y1,r1,x2,y2,r2,rmin,mode,xx1,yy1,fail);
 //dcircle(xx1,yy1,rmin);
 d1:=distance2p(xx1,yy1,x3,y3)-r3;
 xcircle2cr(x1,y1,r1,x2,y2,r2,rmax,mode,xx2,yy2,fail);
 //dcircle(xx2,yy2,rmax);
 d2:=distance2p(xx2,yy2,x3,y3)-r3;
 rmitte:=(rmin+rmax)/2;
 xcircle2cr(x1,y1,r1,x2,y2,r2,rmitte,mode,xx3,yy3,fail);
 d3:=distance2p(xx3,yy3,x3,y3)-r3;
 if abs(d3)<0.1 then break;
 if d1*d3<0 then rmax:=rmitte // zwischen rmin und rmitte
 else if d2*d3<0 then rmin:=rmitte
 else begin
  fail:=true;exit;
 end;
until false;
r:=rmitte;x:=xx3;y:=y3;
end;
procedure circle2cml(x1,y1,r1,x2,y2,r2,a3,b3,c3:double;mode:integer;
 var x,y,r:double;var fail:boolean);
{Kreis tangential an 2 Kreise und Mittelpunkt auf Linie}
var xx1,yy1,xx2,yy2,xx,phi:double;
begin
if a3=0 then begin // waagrechte linie
 circle2cmy(x1,y1,r1,x2,y2,r2,c3,mode,x,r,fail);
 y:=c3;
end else begin
 xx:=c3/a3;
 phi:=atan2(-a3,b3)/radiant;
 rotatep(x1,y1,xx,0,-phi,xx1,yy1);rotatep(x2,y2,xx,0,-phi,xx2,yy2);
 circle2cmy(xx1,yy1,r1,xx2,yy2,r2,0,mode,x,r,fail);
 rotatep(x,0,xx,0,phi,x,y);
end;
end;
procedure circle2cm0(x1,y1,r1,x2,y2,r2:double;mode:integer;
 var y,r:double;var fail:boolean);
{Kreis tangential an 2 Kreise und Mittelpunkt auf Y Achse}
begin
if mode=0 then
y:=(-4*Power(r1,2)*y1 + 8*r1*r2*y1 - 4*Power(r2,2)*y1 + 4*Power(x1,2)*y1 - 4*Power(x2,2)*y1 + 4*Power(y1,3) - 4*Power(r1,2)*y2 + 8*r1*r2*y2 -
        4*Power(r2,2)*y2 - 4*Power(x1,2)*y2 + 4*Power(x2,2)*y2 - 4*Power(y1,2)*y2 - 4*y1*Power(y2,2) + 4*Power(y2,3) +
        4*(r1 - r2)*Sqrt(Power(r1,4) - 4*Power(r1,3)*r2 + 6*Power(r1,2)*Power(r2,2) - 4*r1*Power(r2,3) + Power(r2,4) - 2*Power(r1,2)*Power(x1,2) +
           4*r1*r2*Power(x1,2) - 2*Power(r2,2)*Power(x1,2) + Power(x1,4) - 2*Power(r1,2)*Power(x2,2) + 4*r1*r2*Power(x2,2) - 2*Power(r2,2)*Power(x2,2) -
           2*Power(x1,2)*Power(x2,2) + Power(x2,4) - 2*Power(r1,2)*Power(y1,2) + 4*r1*r2*Power(y1,2) - 2*Power(r2,2)*Power(y1,2) + 2*Power(x1,2)*Power(y1,2) +
           2*Power(x2,2)*Power(y1,2) + Power(y1,4) + 4*Power(r1,2)*y1*y2 - 8*r1*r2*y1*y2 + 4*Power(r2,2)*y1*y2 - 4*Power(x1,2)*y1*y2 - 4*Power(x2,2)*y1*y2 -
           4*Power(y1,3)*y2 - 2*Power(r1,2)*Power(y2,2) + 4*r1*r2*Power(y2,2) - 2*Power(r2,2)*Power(y2,2) + 2*Power(x1,2)*Power(y2,2) + 2*Power(x2,2)*Power(y2,2) +
           6*Power(y1,2)*Power(y2,2) - 4*y1*Power(y2,3) + Power(y2,4)))/(2.*(-4*Power(r1,2) + 8*r1*r2 - 4*Power(r2,2) + 4*Power(y1,2) - 8*y1*y2 + 4*Power(y2,2)))
else
y:=(-4*Power(r1,2)*y1 + 8*r1*r2*y1 - 4*Power(r2,2)*y1 + 4*Power(x1,2)*y1 - 4*Power(x2,2)*y1 + 4*Power(y1,3) - 4*Power(r1,2)*y2 + 8*r1*r2*y2 -
        4*Power(r2,2)*y2 - 4*Power(x1,2)*y2 + 4*Power(x2,2)*y2 - 4*Power(y1,2)*y2 - 4*y1*Power(y2,2) + 4*Power(y2,3) -
        4*(r1 - r2)*Sqrt(Power(r1,4) - 4*Power(r1,3)*r2 + 6*Power(r1,2)*Power(r2,2) - 4*r1*Power(r2,3) + Power(r2,4) - 2*Power(r1,2)*Power(x1,2) +
           4*r1*r2*Power(x1,2) - 2*Power(r2,2)*Power(x1,2) + Power(x1,4) - 2*Power(r1,2)*Power(x2,2) + 4*r1*r2*Power(x2,2) - 2*Power(r2,2)*Power(x2,2) -
           2*Power(x1,2)*Power(x2,2) + Power(x2,4) - 2*Power(r1,2)*Power(y1,2) + 4*r1*r2*Power(y1,2) - 2*Power(r2,2)*Power(y1,2) + 2*Power(x1,2)*Power(y1,2) +
           2*Power(x2,2)*Power(y1,2) + Power(y1,4) + 4*Power(r1,2)*y1*y2 - 8*r1*r2*y1*y2 + 4*Power(r2,2)*y1*y2 - 4*Power(x1,2)*y1*y2 - 4*Power(x2,2)*y1*y2 -
           4*Power(y1,3)*y2 - 2*Power(r1,2)*Power(y2,2) + 4*r1*r2*Power(y2,2) - 2*Power(r2,2)*Power(y2,2) + 2*Power(x1,2)*Power(y2,2) + 2*Power(x2,2)*Power(y2,2) +
           6*Power(y1,2)*Power(y2,2) - 4*y1*Power(y2,3) + Power(y2,4)))/(2.*(-4*Power(r1,2) + 8*r1*r2 - 4*Power(r2,2) + 4*Power(y1,2) - 8*y1*y2 + 4*Power(y2,2)));
r:=distance2p(0,y,x1,y1)-r1;
end;
procedure circle2cmx(x1,y1,r1,x2,y2,r2,x:double;mode:integer;
 var y,r:double;var fail:boolean);
{Kreis tangential an 2 Kreise und Mittelpunkt auf Y Achse}
begin
x1:=x1-x;x2:=x2-x;
circle2cm0(x1,y1,r1,x2,y2,r2,mode,y,r,fail);
end;
procedure circle2cmy(x1,y1,r1,x2,y2,r2,y:double;mode:integer;
 var x,r:double;var fail:boolean);
{Kreis tangential an 2 Kreise und Mittelpunkt auf X Achse}
begin
circle2cmx(y1,x1,r1,y2,x2,r2,y,mode,x,r,fail);
end;
end.