1 | |
---|
2 | #include <stdio.h> |
---|
3 | #include <math.h> |
---|
4 | #include "matrix.h" |
---|
5 | |
---|
6 | |
---|
7 | class Vector |
---|
8 | { |
---|
9 | public: |
---|
10 | Vector (float x, float y, float z) { this->x=x; this->y = y; this->z = z; }; |
---|
11 | float x, y, z; |
---|
12 | inline Vector cross (const Vector& v) const { return Vector(y * v.z - z * v.y, z * v.x - x * v.z, x * v.y - y * v.x ); } |
---|
13 | }; |
---|
14 | |
---|
15 | void Matrix::eigVl(const Matrix& mat) |
---|
16 | { |
---|
17 | |
---|
18 | float eigValue[3]; |
---|
19 | float eigVc[9]; |
---|
20 | |
---|
21 | float a = 0; |
---|
22 | float b = 0; |
---|
23 | |
---|
24 | float c[3]; |
---|
25 | |
---|
26 | // c[0] is the determinante of mat |
---|
27 | c[0] = this->m11 * this->m22 * this->m33 + |
---|
28 | 2* this->m12 * this->m13 * this->m23 - |
---|
29 | this->m11 * this->m23 * this->m23 - |
---|
30 | this->m22 * this->m13 * this->m13 - |
---|
31 | this->m33 * this->m12 * this->m12; |
---|
32 | |
---|
33 | // c[1] is the trace of a |
---|
34 | c[1] = this->m11 * this->m22 - |
---|
35 | this->m12 * this->m12 + |
---|
36 | this->m11 * this->m33 - |
---|
37 | this->m13 * this->m13 + |
---|
38 | this->m22 * this->m33 - |
---|
39 | this->m23 * this->m23; |
---|
40 | |
---|
41 | // c[2] is the sum of the diagonal elements |
---|
42 | c[2] = this->m11 + |
---|
43 | this->m22 + |
---|
44 | this->m33; |
---|
45 | |
---|
46 | |
---|
47 | // Computing the roots: |
---|
48 | a = (3.0*c[1] - c[2]*c[2]) / 3.0; |
---|
49 | b = (-2.0*c[2]*c[2]*c[2] + 9.0*c[1]*c[2] - 27.0*c[0]) / 27.0; |
---|
50 | |
---|
51 | float Q = b*b/4.0 + a*a*a/27.0; |
---|
52 | |
---|
53 | // 3 distinct Roots |
---|
54 | if (Q < 0) |
---|
55 | { |
---|
56 | printf("good\n"); |
---|
57 | float psi = atan2(sqrt(-Q), -b/2.0); |
---|
58 | float p = sqrt((b/2.0)*(b/2.0) - Q); |
---|
59 | |
---|
60 | eigValue[0] = c[2]/3.0 + 2 * pow(p, 1/3.0) * cos(psi/3.0); |
---|
61 | eigValue[1] = c[2]/3.0 - pow(p, 1/3.0) * (cos(psi/3.0) |
---|
62 | + sqrt(3.0) * sin(psi/3.0)); |
---|
63 | eigValue[2] = c[2]/3.0 - pow(p, 1/3.0) * (cos(psi/3.0) |
---|
64 | - sqrt(3.0) * sin(psi/3.0)); |
---|
65 | |
---|
66 | } |
---|
67 | // 2 Distinct Roots |
---|
68 | else if (Q == 0) |
---|
69 | { |
---|
70 | eigValue[0] = c[2]/3.0 + pow(b/2.0, 1.0/3.0); |
---|
71 | eigValue[1] = c[2]/3.0 + pow(b/2.0, 1.0/3.0); |
---|
72 | eigValue[2] = c[2]/3.0 + 2* pow(b/2.0, 1.0/3.0); |
---|
73 | } |
---|
74 | // 1 Root (not calculating anything.) |
---|
75 | else if (Q > 0) |
---|
76 | { |
---|
77 | printf("A is multiple of Identity matrix (lambda * I3))\n"); |
---|
78 | eigValue[0] = eigValue[1] = eigValue[2] = 1; |
---|
79 | } |
---|
80 | |
---|
81 | Matrix M = *this - Matrix::identity() * eigValue[0]; |
---|
82 | |
---|
83 | float u11, u12, u13, u22, u23, u33; |
---|
84 | |
---|
85 | this->debug(); |
---|
86 | |
---|
87 | printf("%f %f %f\n", eigValue[0], eigValue[1], eigValue[2]); |
---|
88 | |
---|
89 | } |
---|
90 | |
---|
91 | |
---|
92 | void Matrix::debug() const |
---|
93 | { |
---|
94 | printf("input: | %f | %f | %f |\n", this->m11, this->m12, this->m13 ); |
---|
95 | printf(" | %f | %f | %f |\n", this->m21, this->m22, this->m23 ); |
---|
96 | printf(" | %f | %f | %f |\n", this->m31, this->m32, this->m33 ); |
---|
97 | |
---|
98 | } |
---|
99 | |
---|