import java.applet.Applet;
import java.awt.event.*;
import java.awt.*;
import java.awt.geom.*;


public class Conic {

    public static double[] dualConic(double[] X) {

	Matrix M=new Matrix();
	M.a[0][0]=X[0];
	M.a[1][1]=X[2];
	M.a[2][2]=X[5];
	M.a[1][0]=.5*X[1];
	M.a[0][1]=.5*X[1];
	M.a[0][2]=.5*X[3];
	M.a[2][0]=.5*X[3];
	M.a[1][2]=.5*X[4];
	M.a[2][1]=.5*X[4];

	M=M.inverse();
	double[] Y={M.a[0][0],2*M.a[1][0],M.a[1][1],2*M.a[0][2],2*M.a[1][2],M.a[2][2]};
	return(Y);
    }


    public static double[] conicCenter(double[] X) {
	double x=2*X[3]*X[2]-X[4]*X[1];
	double y=2*X[0]*X[4]-X[3]*X[1];
	double d=X[1]*X[1]-4*X[0]*X[2];
	x=x/d;
	y=y/d;
	double[] Y={x,y};
	return(Y);
    }

    public static Complex conicCenter2(double[] X) {
                  double[] x=conicCenter(X);
                  return(new Complex(x[0],x[1]));
    }

      public static double evaluate(double[] Y,Complex z) {
                  double u=z.x;
                  double v=z.y;
    	double test=Y[0]*u*u+Y[1]*u*v+Y[2]*v*v+Y[3]*u+Y[4]*v+Y[5];
                  return(test);
      }


/*This is set up so that the evaluation on Z[0] is pos and on Z[1] is neg*/

      public static Complex[] narrow(double[] Y,Complex[] Z) {

         Complex z=Complex.average(Z[0],Z[1]);
         double test=evaluate(Y,z);   
         Complex[] W=new Complex[2];
        if(test<0) {
             W[0]=new Complex(Z[0]);
             W[1]=new Complex(z);
             return(W);
        }
        W[0]=new Complex(z);
        W[1]=new Complex(Z[1]);
        return(W);
      }


      public static Complex narrow(double[] Y,double t,int n) {
             Complex z1=conicCenter2(Y);
             double Pi=Math.PI;
             double c=Math.cos(2*Pi*t);
             double s=Math.sin(2*Pi*t);
             Complex z2=new Complex(c,s);
             Complex z3=new Complex(10000*z2.x-9999*z1.x,10000*z2.y-9999*z1.y);
             Complex[] Z={z1,z3};
             for(int i=0;i<n;++i) Z=narrow(Y,Z);
             Complex fin=Complex.average(Z[0],Z[1]);
             return(fin);
       }





    /**I adapted this code from a program written by Tom Davis.  I found it on the web.
       Given 5 vectors p0,...,p4, representing points in the plane, the routine 
       returns the valuers A,B,C,D,E,F such that  Ax^2 + B xy + C y^2 + D x + E y + F =0.*/

    public static double[] toConic(Vector p0,Vector p1,Vector p2,Vector p3,Vector p4) {

	Vector L0=Vector.cross(p0,p1);
	Vector L1=Vector.cross(p1,p2);
	Vector L2=Vector.cross(p2,p3);
	Vector L3=Vector.cross(p3,p4);
	Vector Q=Vector.cross(L0,L3);

	double A=Q.x[0];
	double B=Q.x[1];
	double C=Q.x[2];

        double a1, a2, b1, b2, c1, c2;
        double x0, x4, y0, y4, w0, w4;
        double aa, bb, cc, dd, ee, ff;
        double y4w0, w4y0, w4w0, y4y0, x4w0, w4x0, x4x0, y4x0, x4y0;
        double a1a2, a1b2, a1c2, b1a2, b1b2, b1c2, c1a2, c1b2, c1c2;


        a1 = L1.x[0]; b1 = L1.x[1]; c1 = L1.x[2];
        a2 = L2.x[0]; b2 = L2.x[1]; c2 = L2.x[2];
        x0 = p0.x[0]; y0 = p0.x[1]; w0 = p0.x[2];
        x4 = p4.x[0]; y4 = p4.x[1]; w4 = p4.x[2];
    
    y4w0 = y4*w0;
    w4y0 = w4*y0;
    w4w0 = w4*w0;
    y4y0 = y4*y0;
    x4w0 = x4*w0;
    w4x0 = w4*x0;
    x4x0 = x4*x0;
    y4x0 = y4*x0;
    x4y0 = x4*y0;
    a1a2 = a1*a2;
    a1b2 = a1*b2;
    a1c2 = a1*c2;
    b1a2 = b1*a2;
    b1b2 = b1*b2;
    b1c2 = b1*c2;
    c1a2 = c1*a2;
    c1b2 = c1*b2;
    c1c2 = c1*c2;

    aa = -A*a1a2*y4w0
	 +A*a1a2*w4y0 
	 -B*b1a2*y4w0
	 -B*c1a2*w4w0
	 +B*a1b2*w4y0
	 +B*a1c2*w4w0
	 +C*b1a2*y4y0
	 +C*c1a2*w4y0
	 -C*a1b2*y4y0
	 -C*a1c2*y4w0;

    cc =  A*c1b2*w4w0
         +A*a1b2*x4w0
	 -A*b1c2*w4w0
	 -A*b1a2*w4x0
	 +B*b1b2*x4w0
	 -B*b1b2*w4x0
	 +C*b1c2*x4w0
	 +C*b1a2*x4x0
	 -C*c1b2*w4x0
	 -C*a1b2*x4x0;

    ff =  A*c1a2*y4x0
	 +A*c1b2*y4y0
	 -A*a1c2*x4y0
	 -A*b1c2*y4y0
	 -B*c1a2*x4x0
	 -B*c1b2*x4y0
	 +B*a1c2*x4x0
	 +B*b1c2*y4x0
	 -C*c1c2*x4y0
	 +C*c1c2*y4x0;

    bb =  A*c1a2*w4w0
	 +A*a1a2*x4w0
	 -A*a1b2*y4w0
	 -A*a1c2*w4w0
	 -A*a1a2*w4x0
	 +A*b1a2*w4y0
	 +B*b1a2*x4w0
	 -B*b1b2*y4w0
	 -B*c1b2*w4w0
	 -B*a1b2*w4x0
	 +B*b1b2*w4y0
	 +B*b1c2*w4w0
	 -C*b1c2*y4w0
	 -C*b1a2*x4y0
	 -C*b1a2*y4x0
	 -C*c1a2*w4x0
	 +C*c1b2*w4y0
	 +C*a1b2*x4y0
	 +C*a1b2*y4x0
	 +C*a1c2*x4w0;

    dd = -A*c1a2*y4w0
	 +A*a1a2*y4x0
	 +A*a1b2*y4y0
	 +A*a1c2*w4y0
	 -A*a1a2*x4y0
	 -A*b1a2*y4y0
	 +B*b1a2*y4x0
	 +B*c1a2*w4x0
	 +B*c1a2*x4w0
	 +B*c1b2*w4y0
	 -B*a1b2*x4y0
	 -B*a1c2*w4x0
	 -B*a1c2*x4w0
	 -B*b1c2*y4w0
	 +C*b1c2*y4y0
	 +C*c1c2*w4y0
	 -C*c1a2*x4y0
	 -C*c1b2*y4y0
	 -C*c1c2*y4w0
	 +C*a1c2*y4x0;

    ee = -A*c1a2*w4x0
	 -A*c1b2*y4w0
	 -A*c1b2*w4y0
	 -A*a1b2*x4y0
	 +A*a1c2*x4w0
	 +A*b1c2*y4w0
	 +A*b1c2*w4y0
	 +A*b1a2*y4x0
	 -B*b1a2*x4x0
	 -B*b1b2*x4y0
	 +B*c1b2*x4w0
	 +B*a1b2*x4x0
	 +B*b1b2*y4x0
	 -B*b1c2*w4x0
	 -C*b1c2*x4y0
	 +C*c1c2*x4w0
	 +C*c1a2*x4x0
	 +C*c1b2*y4x0
	 -C*c1c2*w4x0
	 -C*a1c2*x4x0;

    aa=aa/ff;
    bb=bb/ff;
    cc=cc/ff;
    dd=dd/ff;
    ee=ee/ff;
    ff=1;

    double[] X={aa,bb,cc,dd,ee,ff};
    return(X);
    }




}