1 /*
2 SPDX-FileCopyrightText: 2020 Chris Rowland <chris.rowland@cherryfield.me.uk>
3
4 SPDX-License-Identifier: GPL-2.0-or-later
5 */
6
7 #include "poleaxis.h"
8
9 #include <cmath>
10
11
12 /******************************************************************
13 poleAxis(sp1, sp2, sp3) finds the mount's RA axis of rotation determined
14 by three points sampled by fixing the mount's DEC and sampling a point
15 a three different RA position.
16
17 For each SkyPoint sp, it finds its corresponding x,y,z coordinates,
18 which are points on a unit sphere. Those 3 coordinates define a plane.
19 That plane intesects the sphere, and the intersection of a plane and a
20 sphere is a circle. The center of that circle would be the axis of rotation
21 defined by the origial 3 points. So, finding the normal to the plane,
22 and pointing in that (normal) direction from the center of the sphere
23 (the origin) gives the axis of rotation of the mount.
24
25 poleAxis returns the normal. The way to interpret this vector
26 is a unit vector (direction) in the x,y,z space. To convert this back to an
27 angle useful for polar alignment, the ra and dec angles can be
28 retrieved with the primary(V) and secondary(V) functions.
29 These can then be turned into an altitude and azimuth angles using methods
30 in the SkyPoint class.
31
32 Further detail
33
34 Definitions: SkyPoint sp contains an hour angle ha and declication dec,
35 with the usual convention that dec=0 is the equator, dec=90 is the north pole
36 about which the axis spins, and ha is the spin angle
37 (positive is counter clockwise when looking north from the center of the sphere)
38 where 0-degrees would be pointing "up", e.g. toward the zenith when
39 looking at the north pole. The sphere has radius 1.0.
40
41 dirCos(sp) will return a 3D vector V whose components x,y,z
42 are 3D cartesean coordinate positions given these ha & dec angles
43 for points on the surface of a unit sphere.
44 The z direction points towards the north pole.
45 The y direction points left when looking at the north pole
46 The x direction points "up".
47
48 primary(V) where V contains the above x,y,z coordinates,
49 will return the ha angle pointing to V.
50
51 secondary(V) where V contains the above x,y,z coordinates,
52 will return the dec angle pointing to V.
53 ******************************************************************/
54
dirCos(const dms primary,const dms secondary)55 Rotations::V3 PoleAxis::dirCos(const dms primary, const dms secondary)
56 {
57 return Rotations::V3(
58 static_cast<float>(secondary.cos() * primary.cos()),
59 static_cast<float>(secondary.cos() * primary.sin()),
60 static_cast<float>(secondary.sin()));
61 }
62
dirCos(const SkyPoint sp)63 Rotations::V3 PoleAxis::dirCos(const SkyPoint sp)
64 {
65 return dirCos(sp.ra(), sp.dec());
66 }
67
primary(Rotations::V3 dirCos)68 dms PoleAxis::primary(Rotations::V3 dirCos)
69 {
70 dms p;
71 p.setRadians(static_cast<double>(std::atan2(dirCos.y(), dirCos.x())));
72 return p;
73 }
74
secondary(Rotations::V3 dirCos)75 dms PoleAxis::secondary(Rotations::V3 dirCos)
76 {
77 dms p;
78 p.setRadians(static_cast<double>(std::asin(dirCos.z())));
79 return p;
80 }
81
skyPoint(Rotations::V3 dc)82 SkyPoint PoleAxis::skyPoint(Rotations::V3 dc)
83 {
84 return SkyPoint(primary(dc), secondary(dc));
85 }
86
poleAxis(SkyPoint p1,SkyPoint p2,SkyPoint p3)87 Rotations::V3 PoleAxis::poleAxis(SkyPoint p1, SkyPoint p2, SkyPoint p3)
88 {
89 // convert the three positions to vectors, these define the plane
90 // of the Ha axis rotation
91 Rotations::V3 v1 = PoleAxis::dirCos(p1);
92 Rotations::V3 v2 = PoleAxis::dirCos(p2);
93 Rotations::V3 v3 = PoleAxis::dirCos(p3);
94
95 // the Ha axis direction is the normal to the plane
96 Rotations::V3 p = Rotations::V3::normal(v1, v2, v3);
97
98 // p points to the north or south pole depending on the rotation of the points
99 // the other pole position can be determined by reversing the sign of the Dec and
100 // adding 12hrs to the Ha value.
101 // if we want only the north then this would do it
102 //if (p.z() < 0)
103 // p = -p;
104 return p;
105 }
106