CaRMtl/rene/zirkel/objects/QuadricObject.java

830 lines
27 KiB
Java
Raw Permalink Normal View History

/*
Copyright 2006 Rene Grothmann, modified by Eric Hakenholz
This file is part of C.a.R. software.
C.a.R. is a free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, version 3 of the License.
C.a.R. is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
package rene.zirkel.objects;
// file: QuadricObject.java
import java.util.ArrayList;
import java.util.Enumeration;
import rene.util.xml.XmlWriter;
import rene.gui.Global;
import rene.zirkel.ZirkelCanvas;
import rene.zirkel.construction.Construction;
import rene.zirkel.construction.Count;
import rene.zirkel.expression.ExpressionColor;
import rene.zirkel.expression.Quartic;
import rene.zirkel.graphics.MyGraphics;
import rene.zirkel.graphics.PolygonDrawer;
import rene.zirkel.structures.Complex;
import rene.zirkel.structures.CoordinatesXY;
public class QuadricObject extends ConstructionObject implements PointonObject,
MoveableObject {
public PointObject P[];
static Count N=new Count();
public double X[];
public QuadricObject(final Construction c, final PointObject p[]) {
super(c);
P=p;
validate();
updateText();
}
@Override
public void setDefaults() {
setShowName(Global.getParameter("options.quadric.shownames", false));
setShowValue(Global.getParameter("options.quadric.showvalues", false));
setColor(Global.getParameter("options.quadric.color", 0), Global.getParameter("options.circle.pcolor", (ExpressionColor) null, this));
setColorType(Global.getParameter("options.quadric.colortype", 0));
setHidden(Cn.Hidden);
//setObtuse(Cn.Obtuse);
//setSolid(Cn.Solid);
setLarge(Cn.LargeFont);
setBold(Cn.BoldFont);
//setPartial(Cn.Partial);
}
@Override
public void setTargetDefaults() {
setShowName(Global.getParameter("options.quadric.shownames", false));
setShowValue(Global.getParameter("options.quadric.showvalues", false));
setColor(Global.getParameter("options.quadric.color", 0), Global.getParameter("options.circle.pcolor", (ExpressionColor) null, this));
setColorType(Global.getParameter("options.quadric.colortype", 0));
}
@Override
public String getTag() {
return "Quadric";
}
@Override
public int getN() {
return N.next();
}
@Override
public void updateText() {
try {
final String Names[]=new String[P.length];
for (int i=0; i<P.length; i++) {
Names[i]=P[i].getName();
}
if (P[0].is3D()&&P[1].is3D()&&P[2].is3D()&&P[3].is3D()&&P[4].is3D()) setText(Global.name("text.circle3D"));
else setText(textAny(Global.name("text.quadric"), Names));
} catch (final Exception e) {
}
}
private static final double cub2=Math.pow(2, 1.0/3.0);
public static ArrayList<CoordinatesXY> intersect(final QuadricObject quad1,
final QuadricObject quad2) {
ArrayList<CoordinatesXY> points=new ArrayList<CoordinatesXY>();
// coeffs de x^2, xy, y^2, x, y et Cte :
double a=quad1.X[0], b=quad1.X[4], c=quad1.X[1], d=quad1.X[2], e=quad1.X[3], f=quad1.X[5];
double aa=quad2.X[0], bb=quad2.X[4], cc=quad2.X[1], dd=quad2.X[2], ee=quad2.X[3], ff=quad2.X[5];
// System.out.println("*************");
// System.out.println("a="+a);
// System.out.println("b="+b);
// System.out.println("c="+c);
// System.out.println("d="+d);
// System.out.println("e="+e);
// System.out.println("f="+f);
// System.out.println("aa="+aa);
// System.out.println("bb="+bb);
// System.out.println("cc="+cc);
// System.out.println("dd="+dd);
// System.out.println("ee="+ee);
// System.out.println("ff="+ff);
double k1=-aa*b*bb*c+a*bb*bb*c+aa*aa*c*c+aa*b*b*cc-a*b*bb*cc-2*a*aa*c*cc+a*a*cc*cc;
double k2=bb*bb*c*d-b*bb*cc*d-2*aa*c*cc*d+2*a*cc*cc*d-b*bb*c*dd+2*aa*c*c*dd+b*b*cc*dd-2*a*c*cc*dd-aa*bb*c*e+2*aa*b*cc*e-a*bb*cc*e-aa*b*c*ee+2*a*bb*c*ee-a*b*cc*ee;
double k3=cc*cc*d*d-2*c*cc*d*dd+c*c*dd*dd-bb*cc*d*e-bb*c*dd*e+2*b*cc*dd*e+aa*cc*e*e+2*bb*c*d*ee-b*cc*d*ee-b*c*dd*ee-aa*c*e*ee-a*cc*e*ee+a*c*ee*ee+bb*bb*c*f-b*bb*cc*f-2*aa*c*cc*f+2*a*cc*cc*f-b*bb*c*ff+2*aa*c*c*ff+b*b*cc*ff-2*a*c*cc*ff;
double k4=cc*dd*e*e-cc*d*e*ee-c*dd*e*ee+c*d*ee*ee+2*cc*cc*d*f-2*c*cc*dd*f-bb*cc*e*f+2*bb*c*ee*f-b*cc*ee*f-2*c*cc*d*ff+2*c*c*dd*ff-bb*c*e*ff+2*b*cc*e*ff-b*c*ee*ff;
double k5=-cc*e*ee*f+c*ee*ee*f+cc*cc*f*f+cc*e*e*ff-c*e*ee*ff-2*c*cc*f*ff+c*c*ff*ff;
double u1=k2/(4*k1);
double u2=k2*k2/(4*k1*k1)-2*k3/(3*k1);
double u3=k2*k2/(2*k1*k1)-4*k3/(3*k1);
double u4=-k2*k2*k2/(k1*k1*k1)+4*k2*k3/(k1*k1)-8*k4/k1;
double p1=k3*k3-3*k2*k4+12*k1*k5;
double p2=2*k3*k3*k3-9*k2*k3*k4+27*k1*k4*k4+27*k2*k2*k5-72*k1*k3*k5;
// System.out.println("*************");
// System.out.println("k1="+k1);
// System.out.println("k2="+k2);
// System.out.println("k3="+k3);
// System.out.println("k4="+k4);
// System.out.println("k5="+k5);
Complex q1=new Complex(-4*p1*p1*p1+p2*p2).sqrt();
q1=Complex.plus(q1, new Complex(p2, 0));
q1=q1.sqrt3();
// System.out.println("q1="+q1.real()+" + "+q1.img()+" I");
Complex r1=Complex.div(new Complex(cub2*p1), Complex.mult(new Complex(3*k1), q1));
r1=Complex.plus(r1, Complex.div(q1, new Complex(3*cub2*k1)));
Complex sa=Complex.plus(new Complex(u2), r1);
sa=Complex.div(sa.sqrt(), 2);
Complex sb=Complex.minus(new Complex(u3), r1);
sb=Complex.minus(sb, Complex.div(new Complex(u4), Complex.mult(8, sa)));
sb=Complex.div(sb.sqrt(), 2);
Complex sc=Complex.minus(new Complex(u3), r1);
sc=Complex.plus(sc, Complex.div(new Complex(u4), Complex.mult(8, sa)));
sc=Complex.div(sc.sqrt(), 2);
Complex[] X=new Complex[4];
Complex cu1=new Complex(-u1);
X[0]=Complex.minus(cu1, sa);
X[0]=Complex.minus(X[0], sb);
X[1]=Complex.minus(cu1, sa);
X[1]=Complex.plus(X[1], sb);
X[2]=Complex.plus(cu1, sa);
X[2]=Complex.minus(X[2], sc);
X[3]=Complex.plus(cu1, sa);
X[3]=Complex.plus(X[3], sc);
// System.out.println("*************");
// System.out.println("x1 = "+X[0].real()+" + "+X[0].img()+" I");
// System.out.println("x2 = "+X[1].real()+" + "+X[1].img()+" I");
// System.out.println("x3 = "+X[2].real()+" + "+X[2].img()+" I");
// System.out.println("x4 = "+X[3].real()+" + "+X[3].img()+" I");
// Recherche des ordonnées des points d'intersection :
double A=c, B, C, AA=cc, BB, CC;
for (int i=0; i<4; i++) {
if (Math.abs(X[i].img())>1E-10) {
points.add(new CoordinatesXY());
} else {
B=b*X[i].real()+e;
C=a*X[i].real()*X[i].real()+d*X[i].real()+f;
BB=bb*X[i].real()+ee;
CC=aa*X[i].real()*X[i].real()+dd*X[i].real()+ff;
double denom=A*BB-B*AA;
if (Math.abs(denom)<1E-20) {
points.add(new CoordinatesXY());
} else {
double y=(C*AA-A*CC)/denom;//formula by Dominique Tournès
points.add(new CoordinatesXY(X[i].real(), y));
}
// System.out.println("*************");
// System.out.println("x["+i+"]="+X[i].real());
// System.out.println("y["+i+"]="+y);
// System.out.println("A*BB-B*AA="+(A*BB-B*AA));
}
}
// System.out.println("*************");
// for (int i=0;i<points.size();i++){
// CoordinatesXY coords=points.get(i);
// System.out.println("point "+i+" : "+coords.X+" "+coords.Y);
// }
return points;
}
@Override
public void validate() {
for (int i=0; i<P.length; i++) {
if (!P[i].valid()) {
Valid=false;
return;
}
}
Valid=true;
// Baue Koeffizientenmatrix auf (x^2,y^2,x,y,xy,1):
final double A[][]=new double[5][6];
for (int i=0; i<5; i++) {
final double x=P[i].getX(), y=P[i].getY();
A[i][0]=x*x;
A[i][1]=y*y;
A[i][2]=x;
A[i][3]=y;
A[i][4]=x*y;
A[i][5]=1;
double sum=0;
for (int j=0; j<6; j++) {
sum+=A[i][j]*A[i][j];
}
sum=Math.sqrt(sum);
for (int j=0; j<6; j++) {
A[i][j]/=sum;
}
}
// Gauflverfahren, um auf untere Dreiecksmatrix zu kommen
int r=0;
final int colindex[]=new int[6]; // Index der Stufe oder -1 (keine
// Stufe)
// Iteration ¸ber alle Spalten:
for (int c=0; c<6; c++) {
if (r>=5) // Schema schon fertig
{
colindex[c]=-1;
continue;
}
// Berechne Pivotelement mit spaltenweiser Maximumssuche
double max=Math.abs(A[r][c]);
int imax=r;
for (int i=r+1; i<5; i++) {
final double h=Math.abs(A[i][c]);
if (h>max) {
max=h;
imax=i;
}
}
if (max>1e-13) { // Vertausche Zeilen:
if (imax!=r) {
final double h[]=A[imax];
A[imax]=A[r];
A[r]=h;
}
// Mache restliche Spalte zu 0:
for (int i=r+1; i<5; i++) {
final double lambda=A[i][c]/A[r][c];
for (int j=c+1; j<6; j++) {
A[i][j]-=lambda*A[r][j];
}
}
colindex[c]=r;
r++;
} else {
colindex[c]=-1;
}
}
// Berechne die x-Werte:
X=new double[6];
for (int j=5; j>=0; j--) {
if (colindex[j]<0) {
X[j]=1;
} else {
double h=0;
final int i=colindex[j];
for (int k=j+1; k<6; k++) {
h+=A[i][k]*X[k];
}
X[j]=-h/A[i][j];
}
}
// Normalisiere
double sum=0;
for (int i=0; i<=5; i++) {
sum+=Math.abs(X[i]);
}
if (sum<1e-10) {
Valid=false;
}
for (int i=0; i<=5; i++) {
X[i]/=sum;
// Ce qui suit ressemble à un gag, pourtant il semble que l'epsilon au lieu de 0 en coeffs permet
// de surmonter les effets de bord dans des cas particuliers (ex. hyperbole equilatère/parabole)
// sans pour autant porter atteinte à la précision des coordonnées des points d'intersections
// qui restent fiables à 1e-12, soit la précision maximale affichée du logiciel :
X[i]=n(X[i]);
}
}
private static double n(double x) {
if (Math.abs(x)<1E-16) {
double sign=(Math.signum(x)<0)?-1:1;
return 1.0E-16*sign;
}
return x;
}
@Override
public void paint(final MyGraphics g, final ZirkelCanvas zc) {
if (!Valid||mustHide(zc)) {
return;
}
g.setColor(this);
// Draw the lower part of the quadrik (minus the root):
final double start=zc.minX();
double x=start;
final double end=zc.maxX();
final double h=zc.dx(zc.getOne());
boolean valid=false, ptext=false;
double c0=0, r0=0;
double ctext=20, rtext=20;
final PolygonDrawer pd=new PolygonDrawer(true,g, this);
// Draw the lower part of the quadric (plus the root):
while (x<=end) {
try {
final double y=computeLower(x);
final double c=zc.col(x), r=zc.row(y);
if (valid) {//en ajoutant r>0, la ligne verticale disparaissait; elle est réapparue!
pd.drawTo(c, r);
if (!ptext&&r0-r>c-c0&&zc.isInside(x, y)) {
ctext=c;
rtext=r;
ptext=true;
}
} else {
pd.startPolygon(c, r);
}
c0=c;
r0=r;
valid=true;
} catch (final RuntimeException e) {
valid=false;
}
x+=h;
}
pd.finishPolygon();
// Draw the upper part of the quadric (plus the root):
x=start-2*h;
valid=false;
while (x<=end+2*h) {
try {
final double y=computeUpper(x);
final double c=zc.col(x), r=zc.row(y);
if (valid) {
pd.drawTo(c, r);
// Try to find a position for the label:
if (!ptext&&r0-r>c-c0&&zc.isInside(x, y)) {
ctext=c;
rtext=r;
ptext=true;
}
} else // left edge of quadric, connect with lower part
{
try {
final double y1=computeLower(x);
if (x>=start-h&&x<=end+h) {
g.drawLine(c, zc.row(y1), c, r, this);
}
} catch (final RuntimeException e) {
}
pd.startPolygon(c, r);
}
c0=c;
r0=r;
valid=true;
} catch (final RuntimeException e) // no points in that range
{
if (valid) // we just left the right edge of the quadric
{
try {
final double y1=computeLower(x-h);
if (x-h>=start-h&&x-h<=end+h) {
g.drawLine(c0, zc.row(y1), c0, r0, this);
}
} catch (final RuntimeException ex) {
}
}
valid=false;
}
x+=h;
}
pd.finishPolygon();
final String s=getDisplayText();
if (!s.equals("")) {
g.setLabelColor(this);
setFont(g);
DisplaysText=true;
TX1=ctext+zc.col(XcOffset)-zc.col(0);
TY1=rtext+zc.row(YcOffset)-zc.row(0);
drawLabel(g, s);
}
}
static public final String Tags[]={"x^2", "y^2", "x", "y", "xy"};
// public String getDisplayValue() {
// String s="";
// for (int i=0; i<5; i++) {
// s=s+helpDisplayValue(i==0, -X[i], Tags[i]);
// }
// return s+"="+roundDisplay(X[5]);
// }
@Override
public String getDisplayValue() {
String s="";
s+="("+Global.getLocaleNumber(-X[0], "length");
s+="*";
s+=Tags[0]+")";
for (int i=1; i<5; i++) {
s+="+";
s+="("+Global.getLocaleNumber(-X[i], "length");
s+="*";
s+=Tags[i]+")";
// s=s+helpDisplayValue(i==0, -X[i], Tags[i]);
}
return s+"="+Global.getLocaleNumber(X[5], "length");
// return s+"="+roundDisplay(X[5]);
}
@Override
public String getEquation() {
return getDisplayValue();
}
@Override
public boolean nearto(final int cc, final int rr, final ZirkelCanvas zc) {
if (!displays(zc)) {
return false;
}
final int size=(int) zc.selectionSize();
final double start=zc.minX();
double x=start;
final double end=zc.maxX();
final double h=zc.dx(zc.getOne());
while (x<=end) {
try {
final double y=computeUpper(x);
final double c=zc.col(x), r=zc.row(y);
if (Math.abs(cc-c)<=size*3/2
&&Math.abs(rr-r)<=size*3/2) {
return true;
}
} catch (final Exception e) {
}
try {
final double y=computeLower(x);
final double c=zc.col(x), r=zc.row(y);
if (Math.abs(cc-c)<=size*3/2
&&Math.abs(rr-r)<=size*3/2) {
return true;
}
} catch (final Exception e) {
}
x+=h;
}
return false;
}
public double computeUpper(final double x) {
if (Math.abs(X[1])>1e-13) {
final double p=(X[3]+x*X[4])/X[1], q=(X[0]*x*x+X[2]
*x+X[5])
/X[1];
final double h=p*p/4-q;
if (h<0) {
throw new RuntimeException("");
}
return -p/2+Math.sqrt(h);
} else {
return -(X[0]*x*x+X[2]*x+X[5])/(X[3]+X[4]*x);
}
}
public double computeLower(final double x) {
if (Math.abs(X[1])>1e-13) {
final double p=(X[3]+x*X[4])/X[1], q=(X[0]*x*x+X[2]
*x+X[5])
/X[1];
final double h=p*p/4-q;
if (h<0) {
throw new RuntimeException("");
}
return -p/2-Math.sqrt(h);
} else {
throw new RuntimeException("");
}
}
@Override
public void printArgs(final XmlWriter xml) {
for (int i=0; i<P.length; i++) {
xml.printArg("point"+(i+1), P[i].getName());
}
}
@Override
public Enumeration secondaryParams() {
DL.reset();
for (final PointObject element : P) {
DL.add(element);
}
return DL.elements();
}
@Override
public Enumeration depending() {
DL.reset();
for (final PointObject element : P) {
// System.out.println(element.getName());
DL.add(element);
}
return DL.elements();
}
@Override
public void translate() {
for (int i=0; i<P.length; i++) {
P[i]=(PointObject) P[i].getTranslation();
}
}
@Override
public ConstructionObject copy(final double x, final double y) {
try {
final QuadricObject o=(QuadricObject) clone();
setTranslation(o);
o.P=new PointObject[P.length];
for (int i=0; i<P.length; i++) {
o.P[i]=P[i];
}
o.translateConditionals();
o.translate();
o.setName();
o.updateText();
o.setBreak(false);
// o.setTarget(false); Dibs
return o;
} catch (final Exception e) {
return null;
}
}
@Override
public boolean onlynearto(final int x, final int y, final ZirkelCanvas zc) {
return false;
}
@Override
public boolean equals(final ConstructionObject o) {
if (!(o instanceof QuadricObject)||!o.valid()) {
return false;
}
try {
for (int i=0; i<6; i++) {
if (!equals(X[i], ((QuadricObject) o).X[i])) {
return false;
}
}
} catch (final RuntimeException e) {
return false;
}
return true;
}
@Override
public boolean hasUnit() {
return false;
}
public void keepBaricentricCoords(final PointObject P) {
if (!P.isPointOn()) {
return;
}
if (P.BarycentricCoordsInitialzed) {
final PointObject AA=this.P[0];
final PointObject BB=this.P[1];
final PointObject CC=this.P[2];
final double xa=AA.getX(), ya=AA.getY();
final double xb=BB.getX(), yb=BB.getY();
final double xc=CC.getX(), yc=CC.getY();
final double xm=xa+P.Gx*(xb-xa)+P.Gy*(xc-xa);
final double ym=ya+P.Gx*(yb-ya)+P.Gy*(yc-ya);
P.move(xm, ym);
} else {
P.computeBarycentricCoords();
}
}
@Override
public int getDistance(final PointObject P) {
final double a=X[0], b=X[1], c=X[2], d=X[3], e=X[4], r=X[5];
final double xc=P.getX(), yc=P.getY();
if (Math.abs(a*xc*xc+b*yc*yc+c*xc+d*yc+e*xc*yc
+r)<1e-13) // close enough
{
return 0;
}
final double t[]=new double[5], s[]=new double[5];
// Coefficients for fourth order polynomial for lambda (Lagrange factor)
// Minimize (x-xc)^2+(y-yc)^2 with a*x^2+b*y^2+c*x+d*y+e*x*y+r=0
// Computed with Maple
t[0]=a*e*e*d*d-4*a*b*b*c*c+4*a*e*d*b
*c-4*b*a*a*d*d+b*c*c*e*e-c
*Math.pow(e, 3)*d+r*Math.pow(e, 4)-8*r*e*e*b
*a+16*r*b*b*a*a;
t[1]=8*b*b*c*c+8*a*a*d*d-8*e*d*b*c-8
*a*d*c*e+8*r*e*e*b+8*a*b*c*c+8*b
*a*d*d+8*r*e*e*a-32*r*b*b*a-32*r
*b*a*a;
t[2]=12*e*d*c+16*r*b*b-4*b*d*d-8*r*e*e
+4*e*e*d*yc+16*b*b*xc*c-16*b*c*c
-16*a*d*d-4*a*c*c+16*r*a*a+16*a*a
*d*yc+4*xc*e*e*c-8*e*d*b*xc-8*e*yc
*b*c-8*a*d*xc*e-8*a*yc*c*e+16*a*b
*b*xc*xc-4*a*e*e*yc*yc+16*b*a*a*yc
*yc-4*b*xc*xc*e*e+4*Math.pow(e, 3)*yc*xc
+64*r*b*a-16*a*b*xc*e*yc;
t[3]=-32*r*b+8*d*d+8*c*c+16*e*d*xc+8*e
*e*yc*yc+8*xc*xc*e*e-32*r*a-32*b*xc
*c+16*e*yc*c-32*a*d*yc-32*a*b*xc*xc
-32*b*a*yc*yc;
t[4]=16*b*yc*yc+16*d*yc+16*c*xc+16*xc*e*yc
+16*r+16*a*xc*xc;
final int k=Quartic.solve(t, s);
// System.out.println(k+"Solutions found.");
double dmin=1e30, xmin=xc, ymin=yc;
for (int i=0; i<k; i++) // Choose closest solution of Lagrange
// equation
{
final double l=s[i];
// Solve for x,y when lambda is known.
// Computed with Maple
final double px=-(-e*d+4*b*l*xc-2*e*l*yc-4
*l*l*xc+2*b*c-2*l*c)
/(-e*e+4*b*a-4*b*l-4*l*a+4*l*l);
final double py=-(2*a*d+4*a*l*yc-2*l*d-4*l
*l*yc-2*l*xc*e-c*e)
/(-e*e+4*b*a-4*b*l-4*l*a+4*l*l);
final double dist=(px-xc)*(px-xc)+(py-yc)*(py-yc);
if (dist<dmin) {
dmin=dist;
xmin=px;
ymin=py;
}
}
final double dd=Math.sqrt((P.getX()-xmin)*(P.getX()-xmin)
+(P.getY()-ymin)*(P.getY()-ymin));
return (int) Math.round(dd*Cn.getPixel());
}
public void project(final PointObject P) {
keepBaricentricCoords(P);
final double a=X[0], b=X[1], c=X[2], d=X[3], e=X[4], r=X[5];
final double xc=P.getX(), yc=P.getY();
if (Math.abs(a*xc*xc+b*yc*yc+c*xc+d*yc+e*xc*yc
+r)<1e-13) // close enough
{
return;
}
final double t[]=new double[5], s[]=new double[5];
// Coefficients for fourth order polynomial for lambda (Lagrange factor)
// Minimize (x-xc)^2+(y-yc)^2 with a*x^2+b*y^2+c*x+d*y+e*x*y+r=0
// Computed with Maple
t[0]=a*e*e*d*d-4*a*b*b*c*c+4*a*e*d*b
*c-4*b*a*a*d*d+b*c*c*e*e-c
*Math.pow(e, 3)*d+r*Math.pow(e, 4)-8*r*e*e*b
*a+16*r*b*b*a*a;
t[1]=8*b*b*c*c+8*a*a*d*d-8*e*d*b*c-8
*a*d*c*e+8*r*e*e*b+8*a*b*c*c+8*b
*a*d*d+8*r*e*e*a-32*r*b*b*a-32*r
*b*a*a;
t[2]=12*e*d*c+16*r*b*b-4*b*d*d-8*r*e*e
+4*e*e*d*yc+16*b*b*xc*c-16*b*c*c
-16*a*d*d-4*a*c*c+16*r*a*a+16*a*a
*d*yc+4*xc*e*e*c-8*e*d*b*xc-8*e*yc
*b*c-8*a*d*xc*e-8*a*yc*c*e+16*a*b
*b*xc*xc-4*a*e*e*yc*yc+16*b*a*a*yc
*yc-4*b*xc*xc*e*e+4*Math.pow(e, 3)*yc*xc
+64*r*b*a-16*a*b*xc*e*yc;
t[3]=-32*r*b+8*d*d+8*c*c+16*e*d*xc+8*e
*e*yc*yc+8*xc*xc*e*e-32*r*a-32*b*xc
*c+16*e*yc*c-32*a*d*yc-32*a*b*xc*xc
-32*b*a*yc*yc;
t[4]=16*b*yc*yc+16*d*yc+16*c*xc+16*xc*e*yc
+16*r+16*a*xc*xc;
final int k=Quartic.solve(t, s);
// System.out.println(k+"Solutions found.");
double dmin=1e30, xmin=xc, ymin=yc;
for (int i=0; i<k; i++) // Choose closest solution of Lagrange
// equation
{
final double l=s[i];
// Solve for x,y when lambda is known.
// Computed with Maple
final double px=-(-e*d+4*b*l*xc-2*e*l*yc-4
*l*l*xc+2*b*c-2*l*c)
/(-e*e+4*b*a-4*b*l-4*l*a+4*l*l);
final double py=-(2*a*d+4*a*l*yc-2*l*d-4*l
*l*yc-2*l*xc*e-c*e)
/(-e*e+4*b*a-4*b*l-4*l*a+4*l*l);
final double dist=(px-xc)*(px-xc)+(py-yc)*(py-yc);
if (dist<dmin) {
dmin=dist;
xmin=px;
ymin=py;
}
}
P.move(xmin, ymin);
}
public void project(final PointObject P, final double alpha) {
project(P);
}
public void dragTo(final double x, final double y) {
for (int i=0; i<5; i++) {
P[i].move(xd[i]+(x-x1), yd[i]+(y-y1));
}
}
@Override
public void move(final double x, final double y) {
}
public boolean moveable() {
for (int i=0; i<5; i++) {
if (!P[i].moveable()) {
return false;
}
}
return true;
}
double xd[], yd[], x1, y1;
public void startDrag(final double x, final double y) {
if (xd==null) {
xd=new double[5];
yd=new double[5];
}
for (int i=0; i<5; i++) {
xd[i]=P[i].getX();
yd[i]=P[i].getY();
}
x1=x;
y1=y;
}
public double getOldX() {
return 0;
}
public double getOldY() {
return 0;
}
@Override
public void snap(final ZirkelCanvas zc) {
if (moveable()) {
for (int i=0; i<5; i++) {
P[i].snap(zc);
}
}
}
public boolean canInteresectWith(final ConstructionObject o) {
return true;
}
public void repulse(final PointObject P) {
project(P);
}
public PointObject[] getP() {
return P;
}
}