Finding rotation for magnetometer data to astronomical coordinates


  1. Finding optimal rotation and translation between corresponding 3D points which references A Comparison of Four Algorithms for Estimating 3-D Rigid Transformations
  2. Magnetometer II: Real Data
  3. Spherical coordinate system
  4. MathWorks cart2sph function
  5. Celestial Coordinates

Nghia Ho [1] 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.

Screen Shot 2016-09-04 at 20.27.01

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 = 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')])

    def azi(self):
        return az.rad_to_deg(self._azi)

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

    def alt(self):
        return az.rad_to_deg(self._alt)

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

    def radius(self):
        return self._radius

    def radius(self, value):
        self._radius = value

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

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

    def __str__(self):
        return "%s: x = %.2f, y = %.2f, z=%.2f" % (, self.cart[0], self.cart[1], self.cart[2])

# 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="castor_obs", x=-556.99, y=442.38, z=434.35)
arcturus_obs = Star(name="arcturus_obs", x=-235.79, y=279.59, z=-609.55)
vega_obs = Star(name="vega_obs", 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="castor", azi=301.3, alt=16.8)
arcturus = Star(name="arcturus", azi=188.5, alt=56.5)
vega = Star(name="vega", 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[0]; # 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) < 0:
       print "Reflection detected"
       Vt[2,:] *= -1
       R = Vt.T * U.T

    return R

Alt Azimuth to spherical coordinate conversions helper code.

import math
import requests
import time


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

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

def rad_to_deg(l):
    if isinstance(l, list):
        return [e*180/math.pi for e in l]
        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 "x = %.2f, y = %.2f, z=%.2f" % (x,y,z)

def printAziAlt(azi, alt):
    print "azi = %.2f, alt = %.2f" % (azi, alt)

def main():
    while True:
        (x,y,z) = getAxes()
        printAziAlt(azimuth(x,y,z), altitude(x,y,z))

if __name__ == '__main__':

Geef een reactie

Vul je gegevens in of klik op een icoon om in te loggen. logo

Je reageert onder je account. Log uit / Bijwerken )


Je reageert onder je Twitter account. Log uit / Bijwerken )

Facebook foto

Je reageert onder je Facebook account. Log uit / Bijwerken )

Google+ photo

Je reageert onder je Google+ account. Log uit / Bijwerken )

Verbinden met %s