# Determine 3D orientation based on a magnetometer readout

Now I can read the x, y and z axis, my first question was: what do these values actually mean? How can I used them to determine the orientation of the magnetometer?

What do x, y and z mean: it is the direction of the magnetic north. Honeywell’s paper Compass Heading using Magnetometers from Honeywell describes the basics of the earth’s magnetic field, and how the axis readouts can be used to determine the magnetic north heading, and using a correction factor (declination angle) specific for your location to determine the true north heading (geometric north corresponding to the earth’s rotation axis).

Many example programs and projects show to use the magnetometer to build a 2D compass, such as the very interesting blog post Magnetometer II: Real Data that shows how to perform coordinate transformations with matrix operations. While the compass bearing is the Azimuth coordinate, I also need to determine the Altitude, the angle of elevation (inclination) of the telescope.

The direction is defined in the Honeywell PDF linked above. \begin{aligned} \text{let}\ \theta& = \text{atan}(x/y)\\ \text{direction}& = \begin{cases} 90 - \frac{180}{\pi} \theta,& \text{if } y > 0\\ 270 - \frac{180}{\pi} \theta,& \text{if } y < 0\\ 180,& \text{if } y = 0 \land x < 0\\ 0,& \text{if } y = 0 \land x \geq 0\\ \end{cases} \end{aligned}

The y < 0 part of this definition is wrong! As a counterexample, calculate the direction for x, y = 100, -50. All computations should use the 90-… rule for y > 0. To prevent negative values, the modulo operation can be used: $\text{direction} = (90 - \frac{180}{\pi} \theta)\ \text{mod}\ 360$

For further simplification, we can get rid of $90 - \ldots$ by swapping x and y: $\text{direction} = (\frac{180}{\pi} \text{atan}(y,x))\ \text{mod}\ 360$

Since the python math.atan2() function properly handles the corner case where the denominator (adjacent side) is zero, we do not need to check for x or y = 0 in the python code. This results in the following function to calculate the AziMuth:

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


For the inclination first calculate the length of the vector (x, y): $h = \sqrt{x^2 + y^2}$

Then the inclination is given by $90 - (\text{atan}(h/z)) * 180/\pi$

which simplifies to $\text{atan}(z/h) * 180/\pi$

In Python:

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


These values are uncorrected values: I still need to add or subtract ‘declination’ (for AziMuth) and ‘inclination’ (for Altitude) correction values for my geographical location. There are several sources on the web that provide these values, for instance http://www.magnetic-declination.com/

While the resulting AltAzimuth orientation is somewhat in the right direction, as far as I can tell based upon indoor testing with my wired Raspberry and big breadboard, I do not want to use a calculation based upon precise measurement of the magnetic north on the telescope, for the following reasons:

• Geographic location must be known to lookup Magnetic declination and inclination correction factors
• Calibration or offset errors to the magnetic North result in the same kinds of errors when pointing the telescope
• The magnetometer responds to orientation changes with sub-degree precision, which is in the order of precision that I need for the telescope. The error on the absolute measurement of the magnetic North is greater than this change, or delta precision. If I do not need the magnetic North, the absolute offset error is not a problem as long as it is constant.

Instead I am thinking of the following alternative solution: Calculate the required rotation from measured orientation to apparent orientation using the location of two known stars. The basic idea is to take two measurements with the magnetometer pointed at two different stars of known Alt-Azimuth and thus x,y,z coordinates. First calculate the rotation matrix required to move measurement 1 to star 1. Rotate measurement 2 using the same rotation matrix. Then rotate measurement 2 further by spinning on the axis defined by star 1, until the error (distance) of star 2 and rotated measurement 2 is minimal. Hopefully this error is small. The two rotations can be expressed as a single rotation, which should be useful for transformations of all magnetometer measurements.

Complete code of Alt Azimuth calculation (without magnetic correction):

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)

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

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

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

if __name__ == '__main__':
main()

1. Chris Drake zegt: