# Kstars

rotations.cpp
1/*
2 SPDX-FileCopyrightText: 2021 Hy Murveit <hy@murveit.com>
3
5*/
6
7#include "rotations.h"
8
9#include <cmath>
10
11// In the coordinate system used in this file:
12// x points forward
13// y points to the left
14// z points up.
15// In this system, theta is the angle between x & y and is related
16// to our azimuth: theta = 90-azimuth.
17// Phi is the angle away from z, and is related to our altitude as 90-altitude.
18
19namespace Rotations
20{
21
22double d2r(double degrees)
23{
24 return 2 * M_PI * degrees / 360.0;
25}
26
28{
29 return 360.0 * radians / (2 * M_PI);
30}
31
32V3 V3::normal(const V3 &v1, const V3 &v2, const V3 &v3)
33{
34 // First subtract d21 = V2-V1; d32 = V3-V2
35 const V3 d21 = V3(v2.x() - v1.x(), v2.y() - v1.y(), v2.z() - v1.z());
36 const V3 d32 = V3(v3.x() - v2.x(), v3.y() - v2.y(), v3.z() - v2.z());
37 // Now take the cross-product of d21 and d31
38 const V3 cross = V3(d21.y() * d32.z() - d21.z() * d32.y(),
39 d21.z() * d32.x() - d21.x() * d32.z(),
40 d21.x() * d32.y() - d21.y() * d32.x());
41 // Finally normalize cross so that it is a unit vector.
42 const double lenSq = cross.x() * cross.x() + cross.y() * cross.y() + cross.z() * cross.z();
43 if (lenSq == 0.0) return V3();
44 const double len = sqrt(lenSq);
45 // Should we also fail if len < e.g. 5e-8 ??
46 return V3(cross.x() / len, cross.y() / len, cross.z() / len);
47}
48
49double V3::length()
50{
51 return sqrt(X * X + Y * Y + Z * Z);
52}
53
54V3 azAlt2xyz(const QPointF &azAlt)
55{
56 // Convert the new point to xyz
57 // See https://mathworld.wolfram.com/SphericalCoordinates.html
58 const double azRadians = d2r(azAlt.x());
59 const double altRadians = d2r(azAlt.y());
60
61 const double theta = -azRadians;
62 const double phi = (M_PI / 2.0) - altRadians;
63 const double x = cos(theta) * sin(phi);
64 const double y = sin(theta) * sin (phi);
65 const double z = cos(phi);
66 return V3(x, y, z);
67
68}
69
70QPointF xyz2azAlt(const V3 &xyz)
71{
72 // Deal with degenerate values for the atan along the meridian (y == 0).
73 if (xyz.y() == 0.0 && xyz.x() == 0.0)
74 {
75 // Straight overhead
76 return QPointF(0.0, 90.0);
77 }
78
79 const double azRadians = (xyz.y() == 0) ? 0.0 : -atan2(xyz.y(), xyz.x());
80 const double altRadians = (M_PI / 2.0) - acos(xyz.z());
81
83}
84
86{
88 // Since dec=90 points at the pole.
89 const double decRadians = d2r(haDec.y() - 90);
90
91 const double x = cos(decRadians);
92 const double y = sin(haRadians) * sin(decRadians);
93 const double z = cos(haRadians) * sin(decRadians);
94 return rotateAroundY(V3(x, y, z), -fabs(latitude));
95}
96
97QPointF xyz2haDec(const V3 &xyz, double latitude)
98{
99
100 const V3 pt = rotateAroundY(xyz, latitude);
101 const double ha = r2d(atan2(pt.y(), pt.z()));
102
103 V3 pole = Rotations::azAlt2xyz(QPointF(0, fabs(latitude)));
104 const double dec = 90 - getAngle(xyz, pole);
105 return QPointF(ha, dec);
106}
107
108V3 getAxis(const V3 &p1, const V3 &p2, const V3 &p3)
109{
110 return V3::normal(p1, p2, p3);
111}
112
113// Returns the angle between two points on a sphere using the spherical law of
114// cosines: https://en.wikipedia.org/wiki/Great-circle_distance
115double getAngle(const V3 &p1, const V3 &p2)
116{
117 QPointF a1 = xyz2azAlt(p1);
118 QPointF a2 = xyz2azAlt(p2);
119 return r2d(acos(sin(d2r(a1.y())) * sin(d2r(a2.y())) +
120 cos(d2r(a1.y())) * cos(d2r(a2.y())) * cos(d2r(a1.x()) - d2r(a2.x()))));
121}
122
123// Using the equations in
125// This rotates point around axis (which should be a unit vector) by degrees.
126V3 rotateAroundAxis(const V3 &point, const V3 &axis, double degrees)
127{
128 const double cosAngle = cos(d2r(degrees));
129 const double sinAngle = sin(d2r(degrees));
130 const double pointDotAxis = (point.x() * axis.x() + point.y() * axis.y() + point.z() * axis.z());
131 const double x = axis.x() * pointDotAxis * (1.0 - cosAngle) + point.x() * cosAngle + (-axis.z() * point.y() + axis.y() *
132 point.z()) * sinAngle;
133 const double y = axis.y() * pointDotAxis * (1.0 - cosAngle) + point.y() * cosAngle + (axis.z() * point.x() + -axis.x() *
134 point.z()) * sinAngle;
135 const double z = axis.z() * pointDotAxis * (1.0 - cosAngle) + point.z() * cosAngle + (-axis.y() * point.x() + axis.x() *
136 point.y()) * sinAngle;
137 return V3(x, y, z);
138}
139
140// Simpler version of above for altitude rotations, our main case.
141// Multiply [x,y,z] by the rotate-Y by "angle" rotation matrix
142// as in https://en.wikipedia.org/wiki/Rotation_matrix
143// cos(angle) 0 sin(angle)
144// 0 1 0
145// -sin(angle) 0 cos(angle)
146V3 rotateAroundY(const V3 &point, double degrees)
147{
148 const double radians = d2r(degrees);
149 const double cosAngle = cos(radians);
150 const double sinAngle = sin(radians);
151 return V3( point.x() * cosAngle + point.z() * sinAngle,
152 point.y(),
153 -point.x() * sinAngle + point.z() * cosAngle);
154}
155
156V3 rotateAroundZ(const V3 &point, double degrees)
157{
158 const double radians = d2r(degrees);
159 const double cosAngle = cos(radians);
160 const double sinAngle = sin(radians);
161 return V3( point.x() * cosAngle - point.y() * sinAngle,
162 point.x() * sinAngle + point.y() * cosAngle,
163 point.z());
164}
165
166// Rotates in altitude then azimuth, as is done to correct for polar alignment.
167// Note, NOT a single rotation along a great circle, but rather two separate
168// rotations.
169QPointF rotateRaAxis(const QPointF &azAltPoint, const QPointF &azAltRotation)
170{
171 const V3 point = azAlt2xyz(azAltPoint);
172 const V3 altRotatedPoint = rotateAroundY(point, -azAltRotation.y());
173
174 // Az rotation is simply adding in the az angle.
175 const QPointF altAz = xyz2azAlt(altRotatedPoint);
176 return QPointF(altAz.x() + azAltRotation.x(), altAz.y());
177}
178
179} // namespace rotations
qreal x() const const
qreal y() const const
QTextStream & dec(QTextStream &stream)
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Fri May 24 2024 11:49:21 by doxygen 1.10.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.