# Finding rotation for magnetometer data to astronomical coordinates

Nghia Ho  provides a python program that uses SVD to find a rotation matrix from one set of 3D points to another set. I used the provided python code to construct a rotation matrix that transforms magnetometer coordinates (that point to the magnetic north) to spherical coordinates of stars the sky. Three magnetometer observations when pointed at known stars in the sky are compared with the stars actual coordinates. Then find a rotation matrix R so R.magnetometer_coordinates = sphere coordinates.

I performed a test by holding my development Raspberry Breadboard with compass together with my phone running GoSkyWatch, pointed the phone to three stars, and noted both the magnetometer readouts (blue) and the star’s real coordinates (green). The results are plotted on the sphere below.

• Blue are magnetometer coordinates with the compass oriented to a star.
• Green are the real coordinates of the stars.
• Yellow are the magnetometer coordinates rotated with the derived matrix. Some important things to note are:

There are several conventions for spherical coordinates and also alt-azimuth angle order. I chose to adopt the convention used by Matlab’s cart2sph function, that uses the signatures (azimuth, altitude, radius) and (N, E, U) for spherical and cartersian coordinates, respectively. For example, North is sperical (0, 0, 1) and cartesian (1,0,0). West is sperical (270,0,1) and cartesian (0,-1,0).

Though there is a multi degree difference between each green and yellow star pair, this was expected since the app on my phone has a similar error when looking at the sky. What the test does show is that magnetometer readouts can be used to transform to star spherical coordinates. For the next test I have to mount the magnetometer to the telescope.

Matlab code to plot the sphere and coordinates:

```hold off
r = 0.95
[x, y, z] = sphere (24);
h = surf(r*x, r*y, r*z);
set(h, 'FaceAlpha', 0.5)
hold on

oz=[0.51 -0.86 -0.84]
oy=[0.53 0.39 0.43]
ox=[-0.67 -0.33 0.31]
scatter3(ox,oy,oz,10, 'filled')

ez=[0.29 0.83 0.86]
ey=[-0.82 -0.08 0.51]
ex=[0.50 -0.55 -0.08]
scatter3(ex,ey,ez,10, 'green','filled')

cz=[0.19 0.88 0.81]
cy=[-0.86 -0.04 0.54]
cx=[0.46 -0.48 -0.21]
scatter3(cx,cy,cz,10, 'yellow','filled')
hold off
```

Main program that collects star info from both cartesian as spherical coordinates, calculates the rotation matrix R, and uses R to calculate star positions (B2) from the magnetometer observed coordinates (A).

```import numpy as np
import altazimuth as az
import math
from transform_3d import rigid_transform_3D_0

class Star(object):
def __init__(self, *args, **kwargs):
self._azi = None
self._alt = None
self.name = kwargs.get('name')
self.azi = az.deg_to_rad(kwargs.get('azi')) if 'azi' in kwargs else None
self.alt = az.deg_to_rad(kwargs.get('alt')) if 'alt' in kwargs else None
if 'x' in kwargs and 'y' in kwargs and 'z' in kwargs:
self.azi = az.deg_to_rad(az.azimuth(kwargs.get('x'), kwargs.get('y'), kwargs.get('z')))
self.alt = az.deg_to_rad(az.altitude(kwargs.get('x'), kwargs.get('y'), kwargs.get('z')))
self.radius = np.linalg.norm([kwargs.get('x'), kwargs.get('y'), kwargs.get('z')])

@property
def azi(self):

@azi.setter
def azi(self, value):
self._azi = value

@property
def alt(self):

@alt.setter
def alt(self, value):
self._alt = value

@property

@property
def cart(self):
return az.azi_alt_to_vec(self._azi, self._alt)

def chord(self, b):
return math.sqrt((self.cart-b.cart)**2 + (self.cart-b.cart)**2 + (self.cart-b.cart)**2)

def __str__(self):
return &amp;quot;%s: x = %.2f, y = %.2f, z=%.2f&amp;quot; % (self.name, self.cart, self.cart, self.cart)

# north = Star(azi=0, alt=90)
# east = Star(azi=90, alt=0)
# south = Star(azi=180, alt=0)
# west = Star(azi=270, alt=0)

castor_obs = Star(name=&amp;quot;castor_obs&amp;quot;, x=-556.99, y=442.38, z=434.35)
arcturus_obs = Star(name=&amp;quot;arcturus_obs&amp;quot;, x=-235.79, y=279.59, z=-609.55)
vega_obs = Star(name=&amp;quot;vega_obs&amp;quot;, x=228.49, y=319, z=-620)
print castor_obs
print arcturus_obs
print vega_obs

print castor_obs.chord(arcturus_obs)
print castor_obs.chord(vega_obs)
print arcturus_obs.chord(vega_obs)

castor = Star(name=&amp;quot;castor&amp;quot;, azi=301.3, alt=16.8)
arcturus = Star(name=&amp;quot;arcturus&amp;quot;, azi=188.5, alt=56.5)
vega = Star(name=&amp;quot;vega&amp;quot;, azi=99.0, alt=59.1)
print castor
print arcturus
print vega

print castor.chord(arcturus)
print castor.chord(vega)
print arcturus.chord(vega)

# Test
A = np.matrix(np.array((castor_obs.cart, arcturus_obs.cart, vega_obs.cart)))
B = np.matrix(np.array((castor.cart, arcturus.cart, vega.cart)))

# R, t = rigid_transform_3D(A, B)
R = rigid_transform_3D_0(A, B)
# n = 2 # number of points
# B2 = (R*A.T) + np.tile(t, (1, n))
B2 = (R*A.T)
B2 = B2.T
print B
print B2
```

Rotation matrix code, from Finding optimal rotation and translation between corresponding 3D points, adapted to not perform translations (all coordinates have origin as center).

```from numpy import *
from math import sqrt

# Input: expects Nx3 matrix of points
# Returns R,t
# R = 3x3 rotation matrix
# t = 3x1 column vector

def rigid_transform_3D_0(A, B):
assert len(A) == len(B)

N = A.shape; # total points

# dot is matrix multiplication for array
H = transpose(A) * B

U, S, Vt = linalg.svd(H)

R = Vt.T * U.T

# special reflection case
if linalg.det(R) &amp;lt; 0:
print &amp;quot;Reflection detected&amp;quot;
Vt[2,:] *= -1
R = Vt.T * U.T

return R

```

Alt Azimuth to spherical coordinate conversions helper code.

```import math
import requests
import time

PIURL='http://pi:5000/'

def getAxes():
response = requests.get(PIURL)
data = response.json()
t = (data.get('x'), data.get('y'), data.get('z'))
return t

def azimuth(x,y,z):
return (180/math.pi * math.atan2(y, x)) % 360

def altitude(x,y,z):
h = math.hypot(y, x)
return (180/math.pi * math.atan2(z, h)) % 360

if isinstance(l, list):
return [e*math.pi/180 for e in l]
else:
return l*math.pi/180

if isinstance(l, list):
return [e*180/math.pi for e in l]
else:
return l*180/math.pi

def azi_alt_to_vec(azi, alt):
z  = math.sin(alt)
hyp = math.cos(alt)
x = hyp*math.cos(azi)
y = hyp*math.sin(azi)
return (x, y, z)

def printAxes(x,y,z):
print &amp;quot;x = %.2f, y = %.2f, z=%.2f&amp;quot; % (x,y,z)

def printAziAlt(azi, alt):
print &amp;quot;azi = %.2f, alt = %.2f&amp;quot; % (azi, alt)

def main():
while True:
(x,y,z) = getAxes()
printAxes(x,y,z)
printAziAlt(azimuth(x,y,z), altitude(x,y,z))
time.sleep(0.5)

if __name__ == '__main__':
main()
```