Calibration

ION MagNav Workshop, PLANS 2023, Monterey, CA

Author

Aaron Nielsen

The views expressed in this article are those of the author and do not necessarily reflect the official policy or position of the United States Government, Department of Defense, United States Air Force or Air University.
Distribution A: Authorized for public release. Distribution is unlimited. Case No. 2023-0427.

How-to

There’s a docker container with a jupyter notebook that contains this exercise:

docker run  -p 8888:8888 anielsen/magnav-plans-workshop:latest

Then select work/calibration.ipynb

Tolles-Lawson

Set up code for python

Code

import sys
sys.path.append('../src')
import numpy as np
import quaternion
from scipy import signal as sps

import vector as vec
from utils import db10
import multipole as mp

from matplotlib import pylab as plt

```{python}
import vector as vec
import numpy as np
from matplotlib import pylab as plt
import quaternion
```

Simulated magnetic measurement

Create an external Earth magnetic field vector in a North, East, Down coordinate system.

Code

b_earth_magnitude = 50_000.0 # nanoTesla
b_earth_inclination = 0.0
b_earth_declination = 0.0

# sensor flat in horizontal plane, rotating one full rotation
headings = np.linspace(0,360,100) # headings in degrees for sensor/platform orientation

Code

abg = np.radians(np.array([b_earth_declination, -b_earth_inclination, 0.0]))
qrot = quaternion.from_euler_angles(abg)
b_earth_ned = quaternion.rotate_vectors(qrot,np.array([[b_earth_magnitude,0,0]]))

# what does the Earth's core field look like to the sensor if the sensor spins in space?

# rotating the aircraft to a specific heading in the Earth reference frame is the
# same as rotating the Earth magnetic field to -heading in the aircraft reference frame.
abg = np.radians(np.array([-headings,
                           np.zeros_like(headings),
                           np.zeros_like(headings)])).T

qrot = quaternion.from_euler_angles(abg)

#b_earth_at_sensor = quaternion.rotate_vectors(qrot,b_earth_ned)
b_earth_at_sensor = quaternion.as_vector_part(
    qrot * quaternion.from_vector_part(b_earth_ned) * qrot.conjugate()
    )


# plot the measured magnetic field at the sensor
fig, ax = plt.subplots(3,1)
fig.suptitle('Earth B-field at sensor')
axlabels = [ 'x', 'y', 'z']
for sax in [0,1,2]:
    ax[sax].plot(headings, b_earth_at_sensor[:,sax], label=axlabels[sax])
    ax[sax].grid()
    ax[sax].set_xlabel('heading')
    ax[sax].legend()
    ax[sax].set_xticks(np.arange(0, 361, 30))

Add permanent moments to aircraft

Code

# permanent moments on the aircraft remain fixed with the platform and
# do not rotate/move with respect to the sensor.
b_perm = np.array([5_000.0, 0.0, 0.0] )

Solve this problem in different ways:

True direction cosines
Direction cosines corrupted by aircraft moment
Tolles-Lawson model that is “correct” - permanent only terms
Tolles-Lawson with induced moments

True direction cosines + Correct Tolles-Lawson model

Direction cosines

Compute the true direction cosines using just the Earth field:

Code

# compute direction cosines from b_earth_at_sensor
dc_earth_at_sensor = vec.unit(b_earth_at_sensor)

A-matrix

This simulation has only a permanent moment, so only include a permanent moment in the A-matrix:

Code

def AMatrix(dc,btot=1,perm_only=False):

    if perm_only:
        return dc #* vec.mag(btot)
    else:
        Aperm = dc

        Aindu = np.hstack([ dc**2 ,
                            dc[:,0,np.newaxis]*dc[:,1:],
                            dc[:,1,np.newaxis]*dc[:,2,np.newaxis]   ])
        Aindu *= vec.mag(btot)

        return np.hstack([ Aperm, Aindu] ) #* vec.mag(btot)

perm_only=True
amat_earth_at_sensor = AMatrix(
    dc_earth_at_sensor,
    b_earth_at_sensor,
    perm_only=perm_only)

Code

plot_AMatrix(headings,amat_earth_at_sensor)

(<Figure size 816x672 with 1 Axes>,
 <Axes: title={'center': 'Permanent terms'}>)

Solve with true direction cosines

We can solve using either the true direction cosines or the DC with the moment included.

Start with true external field, not corrupted by permament moment

Text(0.5, 0.98, 'Computed magnetic correction')

Compensating field

<matplotlib.legend.Legend at 0x7f33b367d6d0>

Compensated field

The compensated field has a residual of 250 nT amplitude at double the frequency of the rotation rate.

<matplotlib.legend.Legend at 0x7f33b35f4100>

True direction cosines + Tolles-Lawson with permanent and induced moments

Code

perm_only=False
amat_earth_at_sensor = AMatrix(
    dc_earth_at_sensor,
    b_earth_at_sensor,
    perm_only=perm_only)

Solve with true direction cosines

We can solve using either the true direction cosines or the DC with the moment included.

Start with true external field, not corrupted by permament moment

Text(0.5, 0.98, 'Computed magnetic correction')

Compensating field

<matplotlib.legend.Legend at 0x7f33b344f490>

Compensated field

The compensated field has a residual of 250 nT amplitude at double the frequency of the rotation rate.

<matplotlib.legend.Legend at 0x7f33b33d5040>

Code

plot_AMatrix(headings,amat_earth_at_sensor)

(<Figure size 816x672 with 2 Axes>,
 array([<Axes: title={'center': 'Permanent terms'}>,
        <Axes: title={'center': 'Induced terms'}>], dtype=object))

Corrupted direction cosines + Tolles-Lawson with permanent moments

Compute the corrupted direction cosines using just the Earth field:

Code

# compute direction cosines from b_earth_at_sensor
dc_sensor = vec.unit(b_sensor)

Code

perm_only=True
amat_sensor = AMatrix(
    dc_sensor,
    b_sensor,
    perm_only=perm_only)

Solve with true direction cosines

We can solve using either the true direction cosines or the DC with the moment included.

Start with true external field, not corrupted by permament moment

Text(0.5, 0.98, 'Computed magnetic correction')

Compensating field

<matplotlib.legend.Legend at 0x7f33b314f2e0>

Compensated field

The compensated field has a residual of 250 nT amplitude at double the frequency of the rotation rate.

<matplotlib.legend.Legend at 0x7f33b31a66d0>

Code

plot_AMatrix(headings,amat_sensor)

(<Figure size 816x672 with 1 Axes>,
 <Axes: title={'center': 'Permanent terms'}>)

Corrupted direction cosines + Tolles-Lawson with permanent and induced moments

Compute the corrupted direction cosines using just the Earth field:

Code

# compute direction cosines from b_earth_at_sensor
dc_sensor = vec.unit(b_sensor)

Code

perm_only=False
amat_sensor = AMatrix(
    dc_sensor,
    b_sensor,
    perm_only=perm_only)

Solve with true direction cosines

We can solve using either the true direction cosines or the DC with the moment included.

Start with true external field, not corrupted by permament moment

Text(0.5, 0.98, 'Computed magnetic correction')

Compensating field

<matplotlib.legend.Legend at 0x7f33b2e94520>

Compensated field

The compensated field has a residual of 250 nT amplitude at double the frequency of the rotation rate.

<matplotlib.legend.Legend at 0x7f33b2dec670>

Code

plot_AMatrix(headings,amat_sensor)

(<Figure size 816x672 with 2 Axes>,
 array([<Axes: title={'center': 'Permanent terms'}>,
        <Axes: title={'center': 'Induced terms'}>], dtype=object))