package math3d;

public class Eigensystem3x3Float {
   float[][] m;
   float[] eigenVectors;
   float[] eigenValues;

   public Eigensystem3x3Float(float[][] symmetricMatrix) {
      this.m = symmetricMatrix;
      if (this.m[0][1] != this.m[1][0] || this.m[0][2] != this.m[2][0] || this.m[1][2] != this.m[2][1]) {
         throw new RuntimeException("Eigensystem3x3Float only works with symmetric matrices");
      }
   }

   public void getEvalues(float[] eigenValues) {
      eigenValues[0] = this.eigenValues[0];
      eigenValues[1] = this.eigenValues[1];
      eigenValues[2] = this.eigenValues[2];
   }

   public float[] getEvaluesCopy() {
      return (float[])this.eigenValues.clone();
   }

   public float[] getEvalues() {
      return this.eigenValues;
   }

   public boolean findEvalues() {
      this.eigenValues = new float[3];
      double A = (double)this.m[0][0];
      double B = (double)this.m[0][1];
      double C = (double)this.m[0][2];
      double D = (double)this.m[1][1];
      double E = (double)this.m[1][2];
      double F = (double)this.m[2][2];
      double a = -1.0;
      double b = A + D + F;
      double c = B * B + C * C + E * E - A * D - A * F - D * F;
      double d = A * D * F - A * E * E - B * B * F + 2.0 * B * C * E - C * C * D;
      double third = 0.3333333333333333;
      double q = (3.0 * a * c - b * b) / (9.0 * a * a);
      double r = (9.0 * a * b * c - 27.0 * a * a * d - 2.0 * b * b * b) / (54.0 * a * a * a);
      double discriminant = q * q * q + r * r;
      if (discriminant > 0.0) {
         return false;
      } else if (discriminant < 0.0) {
         double rootThree = 1.7320508075688772;
         double innerSize = Math.sqrt(r * r - discriminant);
         double innerAngle;
         if (r > 0.0) {
            innerAngle = Math.atan(Math.sqrt(-discriminant) / r);
         } else {
            innerAngle = Math.PI - Math.atan(Math.sqrt(-discriminant) / -r);
         }

         double stSize = Math.pow(innerSize, 0.3333333333333333);
         double sAngle = innerAngle / 3.0;
         double tAngle = -innerAngle / 3.0;
         double sPlusT = 2.0 * stSize * Math.cos(sAngle);
         this.eigenValues[0] = (float)(sPlusT - b / (3.0 * a));
         double firstPart = -(sPlusT / 2.0) - b / 3.0 * a;
         double lastPart = -rootThree * stSize * Math.sin(sAngle);
         this.eigenValues[1] = (float)(firstPart + lastPart);
         this.eigenValues[2] = (float)(firstPart - lastPart);
         return true;
      } else {
         double sPlusT;
         if (r >= 0.0) {
            sPlusT = 2.0 * Math.pow(r, 0.3333333333333333);
         } else {
            sPlusT = -2.0 * Math.pow(-r, 0.3333333333333333);
         }

         double bOver3A = b / (3.0 * a);
         this.eigenValues[0] = (float)(sPlusT - bOver3A);
         this.eigenValues[1] = (float)(-sPlusT / 2.0 - bOver3A);
         this.eigenValues[2] = this.eigenValues[1];
         return true;
      }
   }
}