[5661] | 1 | |
---|
| 2 | #include <stdio.h> |
---|
| 3 | #include <math.h> |
---|
[5662] | 4 | #include "matrix.h" |
---|
[5661] | 5 | |
---|
[5662] | 6 | |
---|
[5665] | 7 | Vector Matrix::eigenValues() const |
---|
[5662] | 8 | { |
---|
[5665] | 9 | Vector eigVl; |
---|
[5662] | 10 | |
---|
[5661] | 11 | float a = 0; |
---|
| 12 | float b = 0; |
---|
| 13 | |
---|
| 14 | float c[3]; |
---|
| 15 | |
---|
| 16 | // c[0] is the determinante of mat |
---|
[5662] | 17 | c[0] = this->m11 * this->m22 * this->m33 + |
---|
| 18 | 2* this->m12 * this->m13 * this->m23 - |
---|
| 19 | this->m11 * this->m23 * this->m23 - |
---|
| 20 | this->m22 * this->m13 * this->m13 - |
---|
| 21 | this->m33 * this->m12 * this->m12; |
---|
[5661] | 22 | |
---|
| 23 | // c[1] is the trace of a |
---|
[5662] | 24 | c[1] = this->m11 * this->m22 - |
---|
| 25 | this->m12 * this->m12 + |
---|
| 26 | this->m11 * this->m33 - |
---|
| 27 | this->m13 * this->m13 + |
---|
| 28 | this->m22 * this->m33 - |
---|
| 29 | this->m23 * this->m23; |
---|
[5661] | 30 | |
---|
| 31 | // c[2] is the sum of the diagonal elements |
---|
[5662] | 32 | c[2] = this->m11 + |
---|
| 33 | this->m22 + |
---|
| 34 | this->m33; |
---|
[5661] | 35 | |
---|
| 36 | |
---|
| 37 | // Computing the roots: |
---|
| 38 | a = (3.0*c[1] - c[2]*c[2]) / 3.0; |
---|
| 39 | b = (-2.0*c[2]*c[2]*c[2] + 9.0*c[1]*c[2] - 27.0*c[0]) / 27.0; |
---|
| 40 | |
---|
| 41 | float Q = b*b/4.0 + a*a*a/27.0; |
---|
| 42 | |
---|
[5662] | 43 | // 3 distinct Roots |
---|
[5661] | 44 | if (Q < 0) |
---|
| 45 | { |
---|
| 46 | printf("good\n"); |
---|
| 47 | float psi = atan2(sqrt(-Q), -b/2.0); |
---|
| 48 | float p = sqrt((b/2.0)*(b/2.0) - Q); |
---|
| 49 | |
---|
[5665] | 50 | eigVl.x = c[2]/3.0 + 2 * pow(p, 1/3.0) * cos(psi/3.0); |
---|
| 51 | eigVl.y = c[2]/3.0 - pow(p, 1/3.0) * (cos(psi/3.0) |
---|
[5661] | 52 | + sqrt(3.0) * sin(psi/3.0)); |
---|
[5665] | 53 | eigVl.z = c[2]/3.0 - pow(p, 1/3.0) * (cos(psi/3.0) |
---|
[5661] | 54 | - sqrt(3.0) * sin(psi/3.0)); |
---|
| 55 | |
---|
| 56 | } |
---|
[5662] | 57 | // 2 Distinct Roots |
---|
[5661] | 58 | else if (Q == 0) |
---|
| 59 | { |
---|
[5665] | 60 | eigVl.x = c[2]/3.0 + pow(b/2.0, 1.0/3.0); |
---|
| 61 | eigVl.y = c[2]/3.0 + pow(b/2.0, 1.0/3.0); |
---|
| 62 | eigVl.z = c[2]/3.0 + 2* pow(b/2.0, 1.0/3.0); |
---|
[5661] | 63 | } |
---|
[5662] | 64 | // 1 Root (not calculating anything.) |
---|
[5661] | 65 | else if (Q > 0) |
---|
| 66 | { |
---|
| 67 | printf("A is multiple of Identity matrix (lambda * I3))\n"); |
---|
[5665] | 68 | eigVl.x = eigVl.y = eigVl.z = 1; |
---|
[5661] | 69 | } |
---|
[5665] | 70 | printf("%f %f %f\n", eigVl.x, eigVl.y, eigVl.z); |
---|
| 71 | return eigVl; |
---|
[5661] | 72 | |
---|
[5665] | 73 | } |
---|
[5661] | 74 | |
---|
[5665] | 75 | void Matrix::eigenVectors(Vector& a, Vector& b, Vector& c) const |
---|
| 76 | { |
---|
| 77 | Vector eigVl = this->eigenValues(); |
---|
| 78 | |
---|
| 79 | float eigVc[9]; |
---|
[5664] | 80 | // EigenVectors |
---|
| 81 | for (int i = 0; i < 3; ++i) |
---|
| 82 | { |
---|
| 83 | printf (":: i = %d\n", i); |
---|
[5665] | 84 | Matrix M = *this - Matrix::identity() * eigVl.x; |
---|
[5664] | 85 | Vector m1, m2, m3; |
---|
[5662] | 86 | |
---|
[5664] | 87 | M.getTransposed().toVectors(m1, m2, m3); |
---|
| 88 | |
---|
| 89 | Vector u1, u2, u3; |
---|
| 90 | |
---|
| 91 | u1 = m2.cross(m3); u1 /= u1.len(); |
---|
| 92 | u2 = m3.cross(m1); u2 /= u2.len(); |
---|
| 93 | u3 = m1.cross(m2); u3 /= u3.len(); |
---|
| 94 | |
---|
| 95 | printf("%f, %f, %f\n", u1.x, u1.y, u1.z); |
---|
| 96 | printf("%f, %f, %f\n", u2.x, u2.y, u2.z); |
---|
| 97 | printf("%f, %f, %f\n", u3.x, u3.y, u3.z); |
---|
| 98 | |
---|
| 99 | // u1 = M*u1; |
---|
| 100 | // u2 = M*u2; |
---|
| 101 | // u3 = M*u3; |
---|
[5665] | 102 | // |
---|
[5664] | 103 | // printf("%f, %f, %f\n", u1.x, u1.y, u1.z); |
---|
| 104 | // printf("%f, %f, %f\n", u2.x, u2.y, u2.z); |
---|
| 105 | // printf("%f, %f, %f\n", u3.x, u3.y, u3.z); |
---|
| 106 | // printf("\n\n"); |
---|
| 107 | } |
---|
| 108 | |
---|
[5667] | 109 | //Vector eigVc[3]; |
---|
[5666] | 110 | /* eigenvec test */ |
---|
[5667] | 111 | /* |
---|
[5666] | 112 | for(int i = 0; i < 3; i++) |
---|
| 113 | { |
---|
[5667] | 114 | eigVc[i].x =-1/this->m13*(this->m33 - eigValue[i]) + (this->m32*(-this->m31*this->m32 + this->m12*this->m33 - this->m12*eigVl[i])) / |
---|
| 115 | this->m13*(-this->m13*this->m22 - this->m12*this->m23 + this->m13*eigVl[i]); |
---|
| 116 | eigVc[i].y = -( -this->m13*this->m23 + this->m12*this->m33 - this->m12*eigVl[i]) / |
---|
| 117 | (-this->m31*this->m22 + this->m12*this->m23 + this->m13*eigVl[i]); |
---|
[5666] | 118 | eigVc[i].z = 1.0f; |
---|
| 119 | |
---|
| 120 | |
---|
| 121 | printf("home brewn: %f, %f, %f\n", eigVc[i].x, eigVc[i].y, eigVc[i].z); |
---|
| 122 | |
---|
[5667] | 123 | }*/ |
---|
[5666] | 124 | |
---|
| 125 | |
---|
| 126 | |
---|
[5662] | 127 | this->debug(); |
---|
| 128 | } |
---|
[5661] | 129 | |
---|
| 130 | |
---|
[5662] | 131 | void Matrix::debug() const |
---|
| 132 | { |
---|
| 133 | printf("input: | %f | %f | %f |\n", this->m11, this->m12, this->m13 ); |
---|
| 134 | printf(" | %f | %f | %f |\n", this->m21, this->m22, this->m23 ); |
---|
| 135 | printf(" | %f | %f | %f |\n", this->m31, this->m32, this->m33 ); |
---|
[5661] | 136 | |
---|
| 137 | } |
---|
[5662] | 138 | |
---|