Kstars

matr.cpp
1 /*
2  SPDX-FileCopyrightText: 2012 Andrew Stepanenko
3 
4  Modified by Jasem Mutlaq <[email protected]> for KStars:
5  SPDX-FileCopyrightText: 2012 Jasem Mutlaq <[email protected]>
6 
7  SPDX-License-Identifier: GPL-2.0-or-later
8 */
9 
10 #include "matr.h"
11 
12 #include "vect.h"
13 
14 #include <cmath>
15 
16 namespace GuiderUtils
17 {
18 Matrix ::Matrix(double v)
19 {
20  for (int i = 0; i < 4; i++)
21  for (int j = 0; j < 4; j++)
22  x[i][j] = (i == j) ? v : 0.0;
23  x[3][3] = 1;
24 }
25 
26 Matrix ::Matrix()
27 {
28  for (auto &row : x)
29  for (double &item : row)
30  {
31  item = 0.0;
32  }
33 
34  x[3][3] = 1;
35 }
36 
37 void Matrix ::Invert()
38 {
39  Matrix Out(1);
40  for (int i = 0; i < 4; i++)
41  {
42  double d = x[i][i];
43  if (d != 1.0)
44  {
45  for (int j = 0; j < 4; j++)
46  {
47  Out.x[i][j] /= d;
48  x[i][j] /= d;
49  }
50  }
51 
52  for (int j = 0; j < 4; j++)
53  {
54  if (j != i)
55  {
56  if (x[j][i] != 0.0)
57  {
58  double mulby = x[j][i];
59  for (int k = 0; k < 4; k++)
60  {
61  x[j][k] -= mulby * x[i][k];
62  Out.x[j][k] -= mulby * Out.x[i][k];
63  }
64  }
65  }
66  }
67  }
68  *this = Out;
69 }
70 
71 void Matrix ::Transpose()
72 {
73  double t;
74  for (int i = 0; i < 4; i++)
75  for (int j = i; j < 4; j++)
76  if (i != j)
77  {
78  t = x[i][j];
79  x[i][j] = x[j][i];
80  x[j][i] = t;
81  }
82 }
83 
84 Matrix &Matrix ::operator+=(const Matrix &A)
85 {
86  if (this == &A)
87  return *this;
88 
89  for (int i = 0; i < 4; i++)
90  for (int j = 0; j < 4; j++)
91  x[i][j] += A.x[i][j];
92 
93  return *this;
94 }
95 
96 Matrix &Matrix ::operator-=(const Matrix &A)
97 {
98  if (this == &A)
99  return *this;
100 
101  for (int i = 0; i < 4; i++)
102  for (int j = 0; j < 4; j++)
103  x[i][j] -= A.x[i][j];
104 
105  return *this;
106 }
107 
108 Matrix &Matrix ::operator*=(double v)
109 {
110  for (auto &row : x)
111  for (double &item : row)
112  {
113  item *= v;
114  }
115 
116  return *this;
117 }
118 
119 Matrix &Matrix ::operator*=(const Matrix &A)
120 {
121  if (this == &A)
122  return *this;
123 
124  Matrix res = *this;
125  for (int i = 0; i < 4; i++)
126  for (int j = 0; j < 4; j++)
127  {
128  double sum = 0;
129  for (int k = 0; k < 4; k++)
130  sum += res.x[i][k] * A.x[k][j];
131 
132  x[i][j] = sum;
133  }
134  return *this;
135 }
136 
137 Matrix operator+(const Matrix &A, const Matrix &B)
138 {
139  Matrix res;
140  for (int i = 0; i < 4; i++)
141  for (int j = 0; j < 4; j++)
142  res.x[i][j] = A.x[i][j] + B.x[i][j];
143 
144  return res;
145 }
146 
147 Matrix operator-(const Matrix &A, const Matrix &B)
148 {
149  Matrix res;
150  for (int i = 0; i < 4; i++)
151  for (int j = 0; j < 4; j++)
152  res.x[i][j] = A.x[i][j] - B.x[i][j];
153 
154  return res;
155 }
156 
157 Matrix operator*(const Matrix &A, const Matrix &B)
158 {
159  Matrix res;
160  for (int i = 0; i < 4; i++)
161  for (int j = 0; j < 4; j++)
162  {
163  double sum = 0;
164  for (int k = 0; k < 4; k++)
165  sum += A.x[i][k] * B.x[k][j];
166 
167  res.x[i][j] = sum;
168  }
169 
170  return res;
171 }
172 
173 Matrix operator*(const Matrix &A, double v)
174 {
175  Matrix res;
176  for (int i = 0; i < 4; i++)
177  for (int j = 0; j < 4; j++)
178  res.x[i][j] = A.x[i][j] * v;
179 
180  return res;
181 }
182 
183 GuiderUtils::Vector operator*(const GuiderUtils::Vector &v, const Matrix &M)
184 {
185  GuiderUtils::Vector res;
186 
187  res.x = v.x * M.x[0][0] + v.y * M.x[0][1] + v.z * M.x[0][2] + M.x[0][3];
188  res.y = v.x * M.x[1][0] + v.y * M.x[1][1] + v.z * M.x[1][2] + M.x[1][3];
189  res.z = v.x * M.x[2][0] + v.y * M.x[2][1] + v.z * M.x[2][2] + M.x[2][3];
190  /*
191  res.x = v.x*M.x[0][0] + v.y*M.x[1][0] + v.z*M.x[2][0] + M.x[3][0];
192  res.y = v.x*M.x[0][1] + v.y*M.x[1][1] + v.z*M.x[2][1] + M.x[3][1];
193  res.z = v.x*M.x[0][2] + v.y*M.x[1][2] + v.z*M.x[2][2] + M.x[3][2];
194 
195  double denom = v.x*M.x[0][3] + v.y*M.x[1][3] + v.z*M.x[2][3] + M.x[3][3];
196  if( denom != 1.0 )
197  {
198  res /= denom;
199  ShowMessage("denom="+FloatToStr(denom));
200  }
201  */
202  // ShowMessage( FloatToStr(v.x)+"\n" + FloatToStr(v.x)+"\n" + FloatToStr(v.x) );
203  // ShowMessage("res.x="+FloatToStr(res.x)+"\nres.y="+FloatToStr(res.y));
204  /*
205  ShowMessage( FloatToStr(M.x[0][0])+", " + FloatToStr(M.x[1][0])+", " + FloatToStr(M.x[2][0])+", " + FloatToStr(M.x[3][0])+"\n" +
206  FloatToStr(M.x[0][1])+", " + FloatToStr(M.x[1][1])+", " + FloatToStr(M.x[2][1])+", " + FloatToStr(M.x[3][1])+"\n" +
207  FloatToStr(M.x[0][2])+", " + FloatToStr(M.x[1][2])+", " + FloatToStr(M.x[2][2])+", " + FloatToStr(M.x[3][2])+"\n" +
208  FloatToStr(M.x[0][3])+", " + FloatToStr(M.x[1][3])+", " + FloatToStr(M.x[2][3])+", " + FloatToStr(M.x[3][3])+"\n"
209  );
210  */
211  return res;
212 }
213 
214 Matrix Translate(const GuiderUtils::Vector &Loc)
215 {
216  Matrix res(1);
217  res.x[0][3] = Loc.x;
218  res.x[1][3] = Loc.y;
219  res.x[2][3] = Loc.z;
220  /*
221  res.x[3][0] = Loc.x;
222  res.x[3][1] = Loc.y;
223  res.x[3][2] = Loc.z;
224  */
225  return res;
226 }
227 
228 Matrix Scale(const GuiderUtils::Vector &v)
229 {
230  Matrix res(1);
231  res.x[0][0] = v.x;
232  res.x[1][1] = v.y;
233  res.x[2][2] = v.z;
234 
235  return res;
236 }
237 
238 Matrix RotateX(double Angle)
239 {
240  Matrix res(1);
241  double Cosine = cos(Angle);
242  double Sine = sin(Angle);
243 
244  res.x[1][1] = Cosine;
245  res.x[2][1] = Sine;
246  res.x[1][2] = -Sine;
247  res.x[2][2] = Cosine;
248  /*
249  res.x[1][1] = Cosine;
250  res.x[2][1] = -Sine;
251  res.x[1][2] = Sine;
252  res.x[2][2] = Cosine;
253  */
254  return res;
255 }
256 
257 Matrix RotateY(double Angle)
258 {
259  Matrix res(1);
260  double Cosine = cos(Angle);
261  double Sine = sin(Angle);
262 
263  res.x[0][0] = Cosine;
264  res.x[2][0] = Sine; //-
265  res.x[0][2] = -Sine; //+
266  res.x[2][2] = Cosine;
267 
268  return res;
269 }
270 
271 Matrix RotateZ(double Angle)
272 {
273  Matrix res(1);
274  double Cosine = cos(Angle);
275  double Sine = sin(Angle);
276 
277  // ShowMessage( "sin="+FloatToStr(Sine)+"\ncos="+FloatToStr(Cosine) );
278 
279  res.x[0][0] = Cosine;
280  res.x[1][0] = Sine; //-
281  res.x[0][1] = -Sine; //+
282  res.x[1][1] = Cosine;
283 
284  return res;
285 }
286 
287 Matrix Rotate(const GuiderUtils::Vector &axis, double angle)
288 {
289  Matrix res(1);
290  double Cosine = cos(angle);
291  double Sine = sin(angle);
292 
293  res.x[0][0] = axis.x * axis.x + (1 - axis.x * axis.x) * Cosine;
294  res.x[0][1] = axis.x * axis.y * (1 - Cosine) + axis.z * Sine;
295  res.x[0][2] = axis.x * axis.z * (1 - Cosine) - axis.y * Sine;
296  res.x[0][3] = 0;
297 
298  res.x[1][0] = axis.x * axis.y * (1 - Cosine) - axis.z * Sine;
299  res.x[1][1] = axis.y * axis.y + (1 - axis.y * axis.y) * Cosine;
300  res.x[1][2] = axis.y * axis.z * (1 - Cosine) + axis.x * Sine;
301  res.x[1][3] = 0;
302 
303  res.x[2][0] = axis.x * axis.z * (1 - Cosine) + axis.y * Sine;
304  res.x[2][1] = axis.y * axis.z * (1 - Cosine) - axis.x * Sine;
305  res.x[2][2] = axis.z * axis.z + (1 - axis.z * axis.z) * Cosine;
306  res.x[2][3] = 0;
307 
308  res.x[3][0] = 0;
309  res.x[3][1] = 0;
310  res.x[3][2] = 0;
311  res.x[3][3] = 1;
312 
313  return res;
314 }
315 // Transforms V into coord sys. v1, v2, v3
316 Matrix Transform(const GuiderUtils::Vector &v1, const GuiderUtils::Vector &v2, const GuiderUtils::Vector &v3)
317 {
318  Matrix res(1);
319 
320  // Flipped columns & rows vs prev version
321  res.x[0][0] = v1.x;
322  res.x[0][1] = v1.y;
323  res.x[0][2] = v1.z;
324  res.x[0][3] = 0;
325 
326  res.x[1][0] = v2.x;
327  res.x[1][1] = v2.y;
328  res.x[1][2] = v2.z;
329  res.x[1][3] = 0;
330 
331  res.x[2][0] = v3.x;
332  res.x[2][1] = v3.y;
333  res.x[2][2] = v3.z;
334  res.x[2][3] = 0;
335 
336  res.x[3][0] = 0;
337  res.x[3][1] = 0;
338  res.x[3][2] = 0;
339  res.x[3][3] = 1;
340 
341  return res;
342 }
343 
344 Matrix MirrorX()
345 {
346  Matrix res(1);
347  res.x[0][0] = -1;
348  return res;
349 }
350 
351 Matrix MirrorY()
352 {
353  Matrix res(1);
354  res.x[1][1] = -1;
355  return res;
356 }
357 
358 Matrix MirrorZ()
359 {
360  Matrix res(1);
361  res.x[2][2] = -1;
362  return res;
363 }
364 } // namespace GuiderUtils
const QCA_EXPORT SecureArray operator+(const SecureArray &a, const SecureArray &b)
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Sat Aug 13 2022 04:01:55 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.