public class Quadratic extends Shape {
   double A,B,C,D,E,F,G,H,I,J, _A,_B,_C,_D,_E,_F,_G,_H,_I,_J;

   int nClipPlanes = 0;
   double cp[][] = new double[10][4], _cp[][] = new double[10][4];

   void addClipPlane(double a, double b, double c, double d) {
      cp[nClipPlanes][0] = _cp[nClipPlanes][0] = a;
      cp[nClipPlanes][1] = _cp[nClipPlanes][1] = b;
      cp[nClipPlanes][2] = _cp[nClipPlanes][2] = c;
      cp[nClipPlanes][3] = _cp[nClipPlanes][3] = d;
      nClipPlanes++;
   }

   Quadratic(double c[]) {
      A=c[0];B=c[1];C=c[2];D=c[3];E=c[4];F=c[5];G=c[6];H=c[7];I=c[8];J=c[9];
      _A=A;_B=B;_C=C;_D=D;_E=E;_F=F;_G=G;_H=H;_I=I;_J=J;
   }
   void transform() {
      double m[][] = { { cos, 0, sin, x},
	               {   0, 1,   0, y},
	               {-sin, 0, cos, z},
	               {   0, 0,   0, 1} };
      transform(m);
   }
   double m[][] = new double[4][4];
   void rotate(double theta) {
      Matrix.identity(m);
      Matrix.rotateY(m, theta);
      transform(m);
   }

   void transform(double m[][]) {
      Matrix.invert(m, matrix);

      // --------------------------------------------------------------------
      // TRANSFORM EQUATION OF SECOND ORDER POLYNOMIAL AS FOLLOWS:
      //
      // APPLY POLYNOMIAL Axx + Byy + Czz + Dyz + Ezx + Fxy + Gx + Hy + Iz + J
      //
      // TO TRANSFORMED POINT [ ax+dy+gz+j , bx+ey+hz+k , cx+fy+iz+l ] TO GET:
      //
      //    A(ax+dy+gz+j)^2 + B(bx+ey+hz+k)^2 + ... + I(cx+fy+iz+l) + J
      //
      // AND GET NEW VALUES FOR A...J BY REGROUPING COEFFICIENTS.
      // --------------------------------------------------------------------

      double a = matrix[0][0], b = matrix[1][0], c = matrix[2][0];
      double d = matrix[0][1], e = matrix[1][1], f = matrix[2][1];
      double g = matrix[0][2], h = matrix[1][2], i = matrix[2][2];
      double j = matrix[0][3], k = matrix[1][3], l = matrix[2][3];

      A = _A*a*a + _B*b*b + _C*c*c + _D*b*c + _E*c*a + _F*a*b;
      B = _A*d*d + _B*e*e + _C*f*f + _D*e*f + _E*f*d + _F*d*e;
      C = _A*g*g + _B*h*h + _C*i*i + _D*h*i + _E*i*g + _F*g*h;
      D = 2*(_A*d*g + _B*e*h + _C*f*i) + _D*(e*i+h*f) + _E*(f*g+i*d) + _F*(d*h+g*e);
      E = 2*(_A*g*a + _B*h*b + _C*i*c) + _D*(h*c+b*i) + _E*(i*a+c*g) + _F*(g*b+a*h);
      F = 2*(_A*a*d + _B*b*e + _C*c*f) + _D*(b*f+e*c) + _E*(c*d+f*a) + _F*(a*e+d*b);
      G = 2*(_A*a*j + _B*b*k + _C*c*l) + _G*a + _H*b + _I*c;
      H = 2*(_A*d*j + _B*e*k + _C*f*l) + _G*d + _H*e + _I*f;
      I = 2*(_A*g*j + _B*h*k + _C*i*l) + _G*g + _H*h + _I*i;
      J = _A*j*j + _B*k*k + _C*l*l + _D*k*l + _E*l*j + _F*j*k + _G*j + _H*k + _I*l + _J;

// TRANSFORM EQUATIONS OF CLIPPING PLANES:

      for (int I = 0 ; I < nClipPlanes ; I++) {
         cp[I][0] = _cp[I][0] * a + _cp[I][1] * b + _cp[I][2] * c;
         cp[I][1] = _cp[I][0] * d + _cp[I][1] * e + _cp[I][2] * f;
         cp[I][2] = _cp[I][0] * g + _cp[I][1] * h + _cp[I][2] * i;
         cp[I][3] = _cp[I][0] * j + _cp[I][1] * k + _cp[I][2] * l + _cp[I][3];
      }
   }
   int traceRay(double V[], double t[]) {
      int nroots = quadratic(t,                                    // SOLVE Att + Bt + C:
         A*V[0]  + B*V[1]  + C*V[2]  + D*V[3]  + E*V[4]  + F*V[5]  ,      // A

         A*V[6]  + B*V[7]  + C*V[8]  + D*V[9]  + E*V[10] + F*V[11] +
	                               G*V[12] + H*V[13] + I*V[14] ,      // B

         A*V[15] + B*V[16] + C*V[17] + D*V[18] + E*V[19] + F*V[20] +
	                               G*V[21] + H*V[22] + I*V[23] + J);  // C

      // --------------------------------------------------------------------
      // APPLY CLIP PLANE [a,b,c,d] TO POINT (v + t w) TO GET:
      //
      //    a(vx + t wx) + b(vy + t wy) + c(vz + t wz) + d >= 0
      //
      // WHICH IS THE SAME AS:
      //
      //    t * (a wx + b wy + c wz) + (a vx + b vy + c vz + d) >= 0
      // --------------------------------------------------------------------

      if (nroots > 0 && nClipPlanes > 0) {
         for (int i = 0 ; i < nClipPlanes ; i++) {
	     double A = 2*(cp[i][0]*V[12] + cp[i][1]*V[13] + cp[i][2]*V[14]);
	     double B = cp[i][0]*V[21]+cp[i][1]*V[22]+cp[i][2]*V[23]+cp[i][3];

	     if (A > 0 && t[0] * A + B < 0)   // IF NEAR CLIP PLANE IS AFTER 1ST ROOT
		t[0] = -B / A;                //    THEN REPLACE 1ST ROOT
             if (A < 0 && t[1] * A + B < 0)   // IF FAR CLIP PLANE IS BEFORE 2ND ROOT
		t[1] = -B / A;                //    THEN REPLACE 2ND ROOT
         }
         if (t[0] > t[1]) // SEE WHETHER CLIPPING HAS RESULTED IN NO ROOTS
	     nroots = 0;
      }
      return nroots;
   }
   int quadratic(double t[], double A, double B, double C) {
      double d = B * B - A * C;
      if (d < 0)
	 return 0;
      d = Math.sqrt(d);
      t[0] = (-B - d) / A;
      t[1] = (-B + d) / A;
      return 2;
   }
   void computeNormal(double p[], double n[]) {
      n[0] = 2*A*p[0] + E*p[2] + F*p[1] + G;
      n[1] = 2*B*p[1] + D*p[2] + F*p[0] + H;
      n[2] = 2*C*p[2] + D*p[1] + E*p[0] + I;
      normalize(n);
   }
   static void compute_V(double v[], double w[], double V[]) {
      double vx = v[0], vy = v[1], vz = v[2], wx = w[0], wy = w[1], wz = w[2];

      V[0] = wx*wx; V[1] = wy*wy; V[2] = wz*wz;
      V[3] = wy*wz; V[4] = wz*wx; V[5] = wx*wy;
      V[6] = vx*wx; V[7] = vy*wy; V[8] = vz*wz;
      V[ 9] = .5*(vy*wz + vz*wy);
      V[10] = .5*(vz*wx + vx*wz);
      V[11] = .5*(vx*wy + vy*wx);
      V[12] = .5*wx; V[13] = .5*wy; V[14] = .5*wz;
      V[15] = vx*vx; V[16] = vy*vy; V[17] = vz*vz;
      V[18] = vy*vz; V[19] = vz*vx; V[20] = vx*vy;
      V[21] =    vx; V[22] =    vy; V[23] =    vz;
   }
}