Kstars

rotations.cpp
1 /*
2  SPDX-FileCopyrightText: 2021 Hy Murveit <[email protected]>
3 
4  SPDX-License-Identifier: GPL-2.0-or-later
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 
19 namespace Rotations
20 {
21 
22 double d2r(double degrees)
23 {
24  return 2 * M_PI * degrees / 360.0;
25 }
26 
27 double r2d(double radians)
28 {
29  return 360.0 * radians / (2 * M_PI);
30 }
31 
32 V3 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 
49 double V3::length()
50 {
51  return sqrt(X * X + Y * Y + Z * Z);
52 }
53 
54 V3 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 
70 QPointF 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 
82  return QPointF(r2d(azRadians), r2d(altRadians));
83 }
84 
85 V3 haDec2xyz(const QPointF &haDec, double latitude)
86 {
87  const double haRadians = d2r(haDec.x());
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 
97 QPointF 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 
108 V3 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
115 double 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
124 // https://sites.google.com/site/glennmurray/Home/rotation-matrices-and-formulas/rotation-about-an-arbitrary-axis-in-3-dimensions
125 // This rotates point around axis (which should be a unit vector) by degrees.
126 V3 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)
146 V3 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 
156 V3 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.
169 QPointF 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
QTextStream & dec(QTextStream &stream)
qreal x() const const
qreal y() const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2022 The KDE developers.
Generated on Fri Aug 12 2022 04:00:57 by doxygen 1.8.17 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.