GOAT (Geometrical optics application tool) 0.1
Loading...
Searching...
No Matches
matrix.h
Go to the documentation of this file.
1/*****************************************************************/
8#pragma once
9#include "vector.h"
10#include <iostream>
11
12namespace GOAT
13{
14 namespace maths
15 {
22 template<class T> class Matrix
23 {
24 public:
26 {
27 for (int i = 0; i < 3; i++)
28 for (int j = 0; j < 3; j++)
29 M[i][j] = 0.0;
30 }
31
32 void binWrite(std::ofstream& os)
33 {
34 for (int i = 0; i < 3; i++)
35 for (int j = 0; j < 3; j++)
36 os.write((char*)&M[i][j], sizeof(M[i][j]));
37 }
38
39 void binRead(std::ifstream& is)
40 {
41 for (int i = 0; i < 3; i++)
42 for (int j = 0; j < 3; j++)
43 is.read((char*)&M[i][j], sizeof(M[i][j]));
44 }
45
46
55 {
56 for (int i = 0; i < 3; i++)
57 M[i][0] = a[i];
58 for (int i = 0; i < 3; i++)
59 M[i][1] = b[i];
60 for (int i = 0; i < 3; i++)
61 M[i][2] = c[i];
62 }
63
65 {
66 if (this == &A) return *this;
67 for (int i = 0; i < 3; i++)
68 for (int j = 0; j < 3; j++)
69 M[i][j] = A.M[i][j];
70 return *this;
71 }
72
73 bool operator == (const Matrix& A)
74 {
75 bool Erg = true;
76 for (int i = 0; i < 3; i++)
77 for (int j = 0; j < 3; j++)
78 Erg = Erg && (M[i][j] == A.M[i][j]);
79 return Erg;
80 }
81
82 bool operator != (const Matrix& M)
83 {
84 bool Erg = false;
85 for (int i = 0; i < 3; i++)
86 for (int j = 0; j < 3; j++)
87 Erg = Erg || (M[i][j] != M.M[i][j]);
88 return Erg;
89 }
90
96 T& operator () (int i, int j) { return M[i][j]; }
97
103 {
104 Matrix<T> Erg;
105 for (int i = 0; i < 3; i++)
106 for (int j = 0; j < 3; j++)
107 Erg.M[i][j] = M[i][j] + A.M[i][j];
108 return Erg;
109 }
110
112 {
113 for (int i = 0; i < 3; i++)
114 for (int j = 0; j < 3; j++)
115 M[i][j] += A.M[i][j];
116 return *this;
117 }
118
120 {
121 Matrix Erg;
122 for (int i = 0; i < 3; i++)
123 for (int j = 0; j < 3; j++)
124 Erg.M[i][j] = -M[i][j];
125 return Erg;
126 }
127
129 {
130 for (int i = 0; i < 3; i++)
131 for (int j = 0; j < 3; j++)
132 M[i][j] -= A.M[i][j];
133 return *this;
134 }
135
137 {
138 T m[3][3];
139 for (int i = 0; i < 3; i++)
140 for (int j = 0; j < 3; j++)
141 {
142 m[i][j] = 0;
143 for (int l = 0; l < 3; l++)
144 m[i][j] += M[i][l] * A.M[l][j];
145 }
146 for (int i = 0; i < 3; i++)
147 for (int j = 0; j < 3; j++)
148 M[i][j] = m[i][j];
149 return *this;
150 }
151
152 Matrix& operator /=(const T& x)
153 {
154 for (int i = 0; i < 3; i++)
155 for (int j = 0; j < 3; j++)
156 M[i][j] /= x;
157 return *this;
158 }
159
160 Matrix& operator *=(const T& x)
161 {
162 for (int i = 0; i < 3; i++)
163 for (int j = 0; j < 3; j++)
164 M[i][j] *= x;
165 return *this;
166 }
167
168
169 friend Matrix operator - (const Matrix& A, const Matrix& B)
170 {
171 Matrix Erg;
172 for (int i = 0; i < 3; i++)
173 for (int j = 0; j < 3; j++)
174 Erg.M[i][j] = A.M[i][j] - B.M[i][j];
175 return Erg;
176 }
177
178 friend Matrix operator * (const Matrix& A, const Matrix& B)
179 {
180 Matrix Erg;
181 for (int i = 0; i < 3; i++)
182 for (int j = 0; j < 3; j++)
183 for (int l = 0; l < 3; l++)
184 Erg.M[i][j] += A.M[i][l] * B.M[l][j];
185 return Erg;
186 }
187
188 friend Matrix operator * (const Matrix& A, const T& x)
189 {
190 Matrix Erg;
191 for (int i = 0; i < 3; i++)
192 for (int j = 0; j < 3; j++)
193 Erg.M[i][j] = A.M[i][j] * x;
194 return Erg;
195 }
196
197 friend Matrix operator * (const T& x, const Matrix& A)
198 {
199 Matrix Erg;
200 for (int i = 0; i < 3; i++)
201 for (int j = 0; j < 3; j++)
202 Erg.M[i][j] = x * A.M[i][j];
203 return Erg;
204 }
205
206 friend Vector<T> operator * (const Matrix& A, const Vector<T>& r)
207 {
208 Vector<T> Erg;
209 for (int i = 0; i < 3; i++)
210 for (int j = 0; j < 3; j++)
211 Erg[i] += A.M[i][j] * r[j];
212 return Erg;
213 }
214
215 friend Vector<T> operator * (const Vector<T>& r, const Matrix& A)
216 {
217 Vector<T> Erg;
218 for (int i = 0; i < 3; i++)
219 for (int j = 0; j < 3; j++)
220 Erg[i] += A.M[i][j] * r[j];
221 return Erg;
222 }
223
224 friend Matrix operator / (const Matrix& A, const T& x)
225 {
226 Matrix Erg;
227 for (int i = 0; i < 3; i++)
228 for (int j = 0; j < 3; j++)
229 Erg.M[i][j] = A.M[i][j] / x;
230 return Erg;
231 }
232
233
244 friend std::ostream& operator << (std::ostream& os, const Matrix& A)
245 {
246 for (int i = 0; i < 3; i++)
247 os << A.M[i][0] << " " << A.M[i][1] << " " << A.M[i][2] << std::endl;
248 return os;
249 }
250
258 friend std::istream& operator >> (std::istream& is, Matrix& A)
259 {
260 for (int i = 0; i < 3; i++)
261 for (int j = 0; j < 3; j++)
262 is >> A.M[i][j];
263 return is;
264 }
265
266
273 friend Matrix transpose(const Matrix& A)
274 {
275 Matrix Erg;
276 for (int i = 0; i < 3; i++)
277 for (int j = 0; j < 3; j++)
278 Erg.M[i][j] = A.M[j][i];
279 return Erg;
280 }
281
288 friend T trace(const Matrix& A)
289 {
290 T Erg = 0;
291 for (int i = 0; i < 3; i++)
292 Erg += A.M[i][i];
293 return Erg;
294 }
295
302 friend T det(const Matrix& A)
303 {
304 T Erg;
305 Erg = A.M[0][0] * (A.M[1][1] * A.M[2][2] - A.M[2][1] * A.M[1][2])
306 - A.M[0][1] * (A.M[1][0] * A.M[2][2] - A.M[2][0] * A.M[1][2])
307 + A.M[0][2] * (A.M[1][0] * A.M[2][1] - A.M[2][0] * A.M[1][1]);
308 return Erg;
309 }
310
319 friend Matrix cutMatrix(const Matrix& M, int i, int j)
320 {
321 Matrix Erg;
322 for (int k = 0; k < 3; k++)
323 for (int l = 0; l < 3; l++)
324 {
325 if ((k != i) && (l != j)) Erg.M[k][l] = M.M[k][l];
326 if ((l == j) && (k == i)) Erg.M[k][l] = 1;
327 else
328 {
329 if (l == j) Erg.M[k][l] = 0;
330 if (k == i) Erg.M[k][l] = 0;
331 }
332 }
333 return Erg;
334 }
335
342 friend Matrix invert(const Matrix<T>& A, bool& invertable)
343 {
344 Matrix Erg;
345 T C;
346 C = det(A);
347 if (C == 0.0) invertable = false;
348 else
349 {
350 Erg.M[0][0] = A.M[1][1] * A.M[2][2] - A.M[2][1] * A.M[1][2];
351 Erg.M[1][0] = A.M[2][0] * A.M[1][2] - A.M[1][0] * A.M[2][2];
352 Erg.M[2][0] = A.M[1][0] * A.M[2][1] - A.M[2][0] * A.M[1][1];
353
354 Erg.M[0][1] = A.M[2][1] * A.M[0][2] - A.M[0][1] * A.M[2][2];
355 Erg.M[1][1] = A.M[0][0] * A.M[2][2] - A.M[2][0] * A.M[0][2];
356 Erg.M[2][1] = A.M[2][0] * A.M[0][1] - A.M[0][0] * A.M[2][1];
357
358 Erg.M[0][2] = A.M[0][1] * A.M[1][2] - A.M[1][1] * A.M[0][2];
359 Erg.M[1][2] = A.M[1][0] * A.M[0][2] - A.M[0][0] * A.M[1][2];
360 Erg.M[2][2] = A.M[0][0] * A.M[1][1] - A.M[0][1] * A.M[1][0];
361 invertable = true;
362 }
363 return Erg / C;
364 }
365
371
372 friend Matrix invert(const Matrix<T>& A)
373 {
374 bool dummy;
375 return invert(A, dummy);
376 }
377
378 void clear()
379 {
380 for (int i = 0; i < 3; i++)
381 for (int j = 0; j < 3; j++)
382 M[i][j] = 0;
383 }
384
385 T M[3][3];
386 };
387
392 // Matrix-Matrix Operatoren mit gemischten Typen (std::complex<double> /double)
393 Matrix<std::complex<double> > operator * (const Matrix<double>&, const Matrix<std::complex<double> >&);
394 Matrix<std::complex<double> > operator * (const Matrix<std::complex<double> >&, const Matrix<double>&);
395 Matrix<std::complex<double> > operator + (const Matrix<double>&, const Matrix<std::complex<double> >&);
396 Matrix<std::complex<double> > operator + (const Matrix<std::complex<double> >&, const Matrix<double>&);
397 Matrix<std::complex<double> > operator - (const Matrix<double>&, const Matrix<std::complex<double> >&);
398 Matrix<std::complex<double> > operator - (const Matrix<std::complex<double> >&, const Matrix<double>&);
399
400 // Matrix-Vektor Operatoren mit gemischten Typen (std::complex<double> /double)
401 Vector<std::complex<double> > operator * (const Matrix<double>&, const Vector<std::complex<double> >&);
402 Vector<std::complex<double> > operator * (const Matrix<std::complex<double> >&, const Vector<double>&);
403
404 // Matrix-Skalar Operatoren mit gemischten Typen (std::complex<double> /double)
405 Matrix<std::complex<double> > operator * (const Matrix<double>&, const std::complex<double>&);
406 Matrix<std::complex<double> > operator * (const Matrix<std::complex<double> >&, const double&);
407 Matrix<std::complex<double> > operator * (const double&, const Matrix<std::complex<double> >&);
408 Matrix<std::complex<double> > operator * (const std::complex<double>&, const Matrix<double>&);
409 Matrix<std::complex<double> > operator / (const Matrix<double>&, const std::complex<double>&);
410 Matrix<std::complex<double> > operator / (const Matrix<std::complex<double> >&, const double&);
412
413 // Drehmatrizen um die 3 Raumachsen
423 Matrix<double> Dx(double phi);
430 Matrix<double> Dy(double phi);
437 Matrix<double> Dz(double phi);
438 Matrix<double> rotMatrix(const Vector<double> a, double gamma);
446 Matrix<double> rotMatrix(double alpha, double beta, double gamma);
447
456 Matrix<double> rotMatrix(Vector<double> P, double dtheta, double dphi);
469 void trafo(const Vector<double>& e0,
470 const Vector<double>& e1,
471 const Vector<double>& e2,
473 Matrix<double>& R);
474 const Matrix<double> UNITY = Matrix<double>(Vector<double>(1.0, 0.0, 0.0), Vector<double>(0.0, 1.0, 0.0), Vector<double>(0.0, 0.0, 1.0));
476
477 }
478}
479
This class represents a threedimensional (numeric) Matrix as a template.
Definition matrix.h:23
Matrix operator+(const Matrix &A)
Definition matrix.h:102
bool operator!=(const Matrix &M)
Definition matrix.h:82
Matrix & operator/=(const T &x)
Definition matrix.h:152
Matrix & operator+=(const Matrix &A)
Definition matrix.h:111
bool operator==(const Matrix &A)
Definition matrix.h:73
T & operator()(int i, int j)
This operator returns the value of the j-th element in the i-th column.
Definition matrix.h:96
void binWrite(std::ofstream &os)
Definition matrix.h:32
Matrix & operator*=(const Matrix &A)
Definition matrix.h:136
void binRead(std::ifstream &is)
Definition matrix.h:39
friend Matrix operator*(const Matrix &A, const Matrix &B)
Definition matrix.h:178
friend Matrix invert(const Matrix< T > &A)
Calculates the inverse of a matrix, if possible.
Definition matrix.h:372
Matrix(Vector< T > a, Vector< T > b, Vector< T > c)
Main constructor for the template class Matrix. This constructor uses three vectors as parameters,...
Definition matrix.h:54
friend T det(const Matrix &A)
Calculates the determinant of the matrix A.
Definition matrix.h:302
friend std::istream & operator>>(std::istream &is, Matrix &A)
Input of the matrix from an istream.
Definition matrix.h:258
friend Matrix cutMatrix(const Matrix &M, int i, int j)
This function returns a matrix where all elements of the i-th column and the j-th row set to zero.
Definition matrix.h:319
friend Matrix transpose(const Matrix &A)
Calculates the transpose of the matrix A.
Definition matrix.h:273
friend Matrix invert(const Matrix< T > &A, bool &invertable)
Calculates the inverse of a matrix, if possible.
Definition matrix.h:342
Matrix & operator-=(const Matrix &A)
Definition matrix.h:128
Matrix operator-()
Definition matrix.h:119
friend std::ostream & operator<<(std::ostream &os, const Matrix &A)
Output of the matrix into an ostream Output in the form: M00 M01 M02 M10 M11 M12 M20 M21 M22.
Definition matrix.h:244
friend Matrix operator/(const Matrix &A, const T &x)
Definition matrix.h:224
friend T trace(const Matrix &A)
Calculates the trace of matrix A, i.e. the sum over all diagonal elements.
Definition matrix.h:288
Matrix & operator=(const Matrix &A)
Definition matrix.h:64
Template class for threedimensional vectors.
Definition vector.h:57
Matrix< double > null()
Null-matrix (all components are zero)
const Matrix< double > UNITY
Unity matrix (double-valued)
Definition matrix.h:474
Matrix< double > rotMatrixA(Vector< double > n, Vector< double > k, double gamma)
Calculates rotation matrix for a rotation around an axis passing through a given point and pointing i...
void trafo(const Vector< double > &e0, const Vector< double > &e1, const Vector< double > &e2, Matrix< double > &H, Matrix< double > &R)
Calculates the transformation matrices from the laboratory coordinate system into a local system and ...
Matrix< std::complex< double > > cunity()
unity matrix (complex)
Matrix< double > Dx(double phi)
Rotation matrix around x axis This function returns the matrix for a rotation around the x-axis:
Matrix< double > Dy(double phi)
Rotation matrix around y axis This function returns the matrix for a rotation around the y-axis:
Matrix< std::complex< double > > operator-(const Matrix< double > &, const Matrix< std::complex< double > > &)
Matrix< std::complex< double > > operator/(const Matrix< double > &, const std::complex< double > &)
Matrix< std::complex< double > > operator*(const Matrix< double > &, const Matrix< std::complex< double > > &)
Matrix< double > Dz(double phi)
Rotation matrix around z axis This function returns the matrix for a rotation around the z-axis:
Matrix< std::complex< double > > operator+(const Matrix< double > &, const Matrix< std::complex< double > > &)
Matrix< double > rotMatrix(const Vector< double > a, double gamma)
Rotation matrix for rotation around the axis a by the angle gamma.
const Matrix< std::complex< double > > CUNITY
Unity matrix (complex-valued)
Definition matrix.h:475
Matrix< double > unity()
unity matrix (double precision)
This class is used for the iray class. This class is intended for internal use only....
Definition fresnel.h:7
This file contains the Vector template class and some useful functions around this class.