Kstars

matr.cpp
1/*
2 SPDX-FileCopyrightText: 2012 Andrew Stepanenko
3
4 Modified by Jasem Mutlaq <mutlaqja@ikarustech.com> for KStars:
5 SPDX-FileCopyrightText: 2012 Jasem Mutlaq <mutlaqja@ikarustech.com>
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
16namespace GuiderUtils
17{
18Matrix ::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
26Matrix ::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
37void 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
71void 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
84Matrix &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
96Matrix &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
108Matrix &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
119Matrix &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
137Matrix 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
147Matrix 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
157Matrix 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
173Matrix 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
183GuiderUtils::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
214Matrix 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
228Matrix 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
238Matrix 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
257Matrix 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
271Matrix 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
287Matrix 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
316Matrix 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
344Matrix MirrorX()
345{
346 Matrix res(1);
347 res.x[0][0] = -1;
348 return res;
349}
350
351Matrix MirrorY()
352{
353 Matrix res(1);
354 res.x[1][1] = -1;
355 return res;
356}
357
358Matrix MirrorZ()
359{
360 Matrix res(1);
361 res.x[2][2] = -1;
362 return res;
363}
364} // namespace GuiderUtils
QCA::operator+
QCA_EXPORT const SecureArray operator+(const SecureArray &a, const SecureArray &b)
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Mon Nov 4 2024 16:38:43 by doxygen 1.12.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.