Magnetometers are key instruments on board spacecraft that probe
the plasma environments of planets and other solar system bodies. The linear
conversion of raw magnetometer outputs to fully calibrated magnetic field
measurements requires the accurate knowledge of 12 calibration parameters:
six angles, three gain factors, and three offset values. The in-flight determination of
8 of those 12 parameters is enormously supported if the spacecraft is spin-stabilized, as an incorrect choice of those parameters will lead to
systematic spin harmonic disturbances in the calibrated data. We show that
published equations and algorithms for the determination of the
eight spin-related parameters are far from optimal, as they do not take into
account the physical behavior of science-grade magnetometers and the
influence of a varying spacecraft attitude on the in-flight calibration
process. Here, we address these issues. Based on decade-long developments
and experience in calibration activities at the Braunschweig University of
Technology, we introduce advanced calibration equations, parameters, and
algorithms. With their help, it is possible to decouple different effects on
the calibration parameters, originating from the spacecraft or the
magnetometer itself. A key point of the algorithms is the bulk determination
of parameters and associated uncertainties. The lowest uncertainties are expected
under parameter-specific conditions. By application to THEMIS-C (Time History of
Events and Macroscale Interactions during Substorms) magnetometer
measurements, we show where these conditions are fulfilled along a highly
elliptical orbit around Earth.

The investigation of the plasma environment in the heliosphere, around
planets, moons, comets, or other solar system bodies, requires accurate
in situ observations of the magnetic field. Magnetometers on board spacecraft
can provide these key measurements if accurately calibrated on the ground and in
flight. The calibration process delivers the parameters needed to convert raw
magnetometer measurements into magnetic field observations B=(Bx,
By, Bz)T in physically meaningful coordinate systems and units
(usually nanotesla: nT). Commonly, a linear calibration equation is applied
for this conversion (e.g., Fornaçon et al., 1999; Balogh et al., 2001b; Auster et al., 2008):

(1)B=C⋅BS-OS.

Here BS=(BS1, BS2,
BS3)T is the raw magnetometer output in non-orthogonal
sensor coordinates, OS corrects for non-vanishing
magnetometers outputs in zero ambient fields (so-called offsets, which
include spacecraft-generated magnetic fields at the sensor position), and
C is the 3×3 coupling matrix. This matrix may have the
following form (e.g., Kepko et al., 1996):

The coupling matrix C depends on three scaling factors
(GS1, GS2, and GS3, also called the
gains) and six angles (θ1, θ2, θ3, and ϕ1,
ϕ2, ϕ3) which define the directions of the three sensor axes in
the orthogonal coordinate system to which B pertains. Calibrating a
magnetometer means finding the three gains, six angles, and three offset
components (i.e., in total 12 parameters) so that B can accurately be
obtained from BS.

Ground calibration of magnetometers is facilitated by rotating them in
Earth's magnetic field (Green, 1990). Similarly, operating a
magnetometer on a spinning spacecraft, instead of on a three-axis stabilized
spacecraft, enormously supports the in-flight determination of 8 of the
12 calibration parameters. These eight spin-related parameters are the two spin
plane offset components, five of the six sensor direction angles (all but one
defining the rotation about the spin axis), and the ratio of the spin plane
gains. The reason is that an incorrect choice in any of those eight spin-related
parameters leads to the appearance of clear, systematic signals at the spin
frequency (also called the first harmonic) and/or at twice the spin frequency
(second harmonic) in the de-spun magnetic field measurements. Hence,
minimization of these signals can be used to determine the eight calibration
parameters, as described in Farrell et al. (1995) and Kepko et al. (1996).

The other four (spin-unrelated) parameters are the absolute gains in the spin
plane and along the spin axis, the spin axis offset, and the angle of
rotation of the sensor about the spin axis. Gains and angle can be derived in
flight through comparison of magnetic field measurements with the International
Geomagnetic Reference Field (IGRF) or the Tsyganenko field models, which are
fairly accurate close to Earth (e.g., Thébault et al., 2015; Tsyganenko and Sitnov, 2007). For the determination of the spin axis offset in flight, a
list of different methods exists. Typically, the offset is obtained from
careful analysis of Alfvénic magnetic field fluctuations, present in the
pristine solar wind (e.g., Belcher, 1973; Hedgecock, 1975; Leinweber et al., 2008). If strongly compressional fluctuations are observed instead
of Alfvénic fluctuations, then the mirror mode method may be used
(Plaschke and Narita, 2016; Plaschke et al., 2017). The offset may also be obtained from
comparison with measurements from an absolute magnetometer or time-of-flight
measurements of electrons emitted and observed by an electron drift
instrument (Georgescu et al., 2006; Nakamura et al., 2014; Plaschke et al., 2014). Furthermore,
the spin axis offset may also be obtained in regions of space where the
fields are known, for instance, in diamagnetic cavities in the vicinity of
comets (Goetz et al., 2016a, b).

From the preceding paragraphs, the reader might get the impression that
in-flight calibration of magnetometers on spinning spacecraft is a solved
issue, and in theory this is the case. However, as we will show in the
following sections, the published methods for spin-aided calibration
(Farrell et al., 1995; Kepko et al., 1996) are not optimal in practice because they do
not take into account the physical behavior of the sensor package and the
influence of a varying spacecraft attitude on the in-flight calibration.

This paper aims at identifying deficiencies and suggesting improvements with
respect to the calibration equations (Eqs. 1 and 2) and the specific choice of the calibration parameters. Thereafter, we
identify optimal conditions for spin-related calibration parameter
determination. Finally, we introduce advanced algorithms for parameter
determination based on our findings that facilitate the automation and
distribution of calibration activities. A version of these algorithms is
routinely applied to calibrate magnetometer data from the Magnetospheric
Multiscale (MMS) mission (Burch et al., 2016; Torbert et al., 2016; Russell et al., 2016). The
calibration principles and algorithms described here are based on
developments at the Braunschweig University of Technology
(Fornaçon et al., 2011) that have been successfully applied for decades to
calibrate magnetometer data from, e.g., the Equator-S (Fornaçon et al., 1999),
the Cluster (Balogh et al., 2001a, b), and the Time History of
Events and Macroscale Interactions during Substorms (THEMIS) missions
(Angelopoulos, 2008; Auster et al., 2008).

Equations (1) and (2) in principle allow for any
linear conversion of BS into B. The coupling matrix
(Eq. 2) is obviously split into two components:

(3)C=Θ⋅G.

Here, the diagonal matrix G only includes the gains, and the
matrix Θ only includes the angular dependencies. Let us focus
first on the matrix Θ. The parameters θ1, θ2,
and θ3 are the angles between the three mutually non-orthogonal
sensor axes (directions S1–S3) and the spin axis in the z direction in
an orthogonal, spin-axis-aligned, and spacecraft-fixed coordinate system
(directions x, y, and z). The parameters ϕ1, ϕ2,
and ϕ3 correspond to the angles between the spacecraft-fixed x direction
in the spin plane (x–y plane, perpendicular to the spin axis) and the
projections of the sensor axes onto that plane. For simplicity, the sensor
axes S1, S2, and S3 are assumed to be approximately aligned with x,
y, and z. Note that all coordinate systems used in this paper are listed in Table 1.

Table 1List of coordinate system notations used in this paper similar to Table 1
in Kepko et al. (1996).

The individual link of the sensor axes to a spacecraft-fixed, spin-axis-aligned system is an issue here as it does not reflect the actual situation
on the spacecraft: there, the three sensor axes are typically packaged
together into one sensor system. One of the design criteria of modern
fluxgate magnetometer sensors is the temperature and long-term stability of
the sensor axis directions as defined with respect to the sensor package. The
angles between the sensor axes are usually well known from ground calibration
activities (e.g., Auster et al., 2008; Russell et al., 2016), and we can expect the
three angles between the sensor axes to be relatively stable parameters.
Consequently, in a first step, the magnetometer output in non-orthogonal
sensor coordinates should be transformed into an orthogonal sensor package
fixed coordinate system (coordinates: Px, Py, Pz; see Table 1).
The conversion matrix may have the following form:

(4)Γ=sinθS10cosθS1cosϕS12sinθS2sinϕS12sinθS2cosθS2001-1.

Here, θS1 and θS2 are the angles
between the sensor axes S1 and S2 with respect to S3=Pz, and
ϕS12 is the angle between the projections of S1 and S2 onto
a plane perpendicular to S3, the Px–Py plane. Note that
S1 lies in the Px–Pz plane (see Fig. 1a).

In the next step, the orientation of that sensor package system needs to be
defined in a spacecraft-fixed spin-axis-aligned coordinate system. This
latter transformation is expected to change every time there is a maneuver of
the spacecraft, as fuel consumption will change the tensor of inertia and,
thus, the spin axis direction in any spacecraft-fixed coordinate system. The
spin axis direction can be defined in the orthogonal sensor package system
using two parameters or angles. During maneuvers, only those two
parameters or angles should change because the geometry inside the sensor
package should not be affected. A rotation matrix Σ into a
spin-axis-aligned coordinate system dependent on the two angles σPx
and σPy can be defined as follows:

Σ=cosσPx0-sinσPx010sinσPx0cosσPx(5)⋅1000cosσPy-sinσPy0sinσPycosσPy.

Here, σPy is the angle between Pz and the
projection of the spin axis onto the Py–Pz plane,
positive towards Py; σPx is the angle between
that projection and the spin axis, positive towards Px. The angles
are illustrated in Fig. 1b. Note that the spin axis is
assumed to be approximately aligned with the Pz=S3
axis. As a result, the angles σPx and σPy
will be small and can be associated with the
Px and Py coordinates of a unit vector that points in spin axis direction.

Using the angles σPx and σPy to define
the spin axis direction is advantageous over using the angles θ3
and ϕ3, as the latter angle is badly defined if θ3 is small.
Furthermore, it should also be noted that a change in direction of the spin
axis requires an update of all angles of matrix Θ as defined
above, even though the magnetometer (sensor) itself is unaffected. Only two
parameters (σPx and σPy) need to be
changed here to adapt the matrix Σ to the new spin axis direction.

To completely orient the sensor package (system) in the spin-axis-aligned
coordinate system, a rotation about the spin axis (rotation matrix Φ)
also needs to be taken into account:

(6)Φ=cosϕa-sinϕa0sinϕacosϕa0001.

As we will show later, this rotation does not affect the spin tone content in
the de-spun magnetic field observations. The angle is affected by the
orientation of a magnetometer boom and may change due to boom bending (Farrell et al., 1995).

Altogether, we can replace the orthogonalization and reorientation
matrix Θ by Φ⋅Σ⋅Γ
in Eq. (3). Let us focus then again on the
gain matrix G in that equation. As mentioned in the Introduction,
the spacecraft spin aids the determination of the ratio g2=GS1/GS2
of the spin plane gains but not the
absolute gains in the spin plane Gp=GS1GS2 and
along the spin axis Ga=GS3.
Hence, it makes sense to use the parameters g and Gp instead
of GS1 and GS2 in the matrix G to
decouple parameters that can be frequently updated from parameters that are
only obtainable in flight from comparison to model fields or measurements of
other instruments:

(7)G=gGp000Gp/g000Ga.

Note that Kepko et al. (1996) use the difference of the inverse gains
ΔG21=1/GS2-1/GS1 instead of g. However,
later changes in the absolute gains GS1 and GS2
then necessarily require an update of ΔG21 in order to avoid
perturbations at the second harmonic in the de-spun data. The gain ratio g,
instead, is decoupled from changes in the absolute gains Gp and Ga.

The gains should be stable parameters in the absence of temperature
variations. These variations in the gains can be determined from ground
calibration, resulting in a diagonal gain correction matrix
GT(Ts, Te) that is dependent on the
magnetometer sensor (Ts) and electronics (Te)
temperatures. That matrix should be directly applied to the magnetometer
output BS, requiring the knowledge of the sensor and
electronics temperatures:

(8)BST=GTTs,Te⋅BS.

The resulting temperature-corrected output BST may then be
further converted to B via the coupling matrix C=Φ⋅Σ⋅Γ⋅G
and the offset vector O using Eq. (1), after
replacing BS with BST. This also has the
advantage that the further applied absolute gains Gp and Ga
and the gain ratio g2 should all be approximately 1 and unitless.

Altogether, we suggest using the following improved calibration equation:

(9)B=Φ⋅Σ⋅Γ⋅G⋅GTTs,Te⋅BS︸=BST-OS,

with matrices defined in Eqs. (4) to (7) instead
of the simpler Eqs. (1) and (2). The
parameters whose determination is supported by the spacecraft spin are
θS1, θS2, ϕS12, σPx, σPy, g,
OS1, and OS2.

To determine the influence of the calibration parameters on the spin tone
harmonic disturbances in the de-spun magnetic field measurements, we use a
similar mathematical approach to Kepko et al. (1996) in this section. Based on
the results, we go on to derive the optimal conditions for the determination
of each parameter in Sect. 4.

First, we compute the temperature-corrected sensor output BST
as a function of the external field B in the spinning coordinate system:

(10)BST=G-1⋅Γ-1⋅Σ-1⋅Φ-1⋅B+OS.

We linearize all the matrices, using the following simplifying assumptions.
The validity of these assumptions and the admissible deviations are discussed
in Sect. 7.

Furthermore, without loss of generality, we assume the magnetic field in the
de-spun (inertial) coordinate system (directions X, Y, and Z) to be in
the X–Z plane, and the spacecraft to spin around the Z axis, which
corresponds to the z axis in the spacecraft-fixed, spin-aligned
coordinate system (see Table 1). In that latter system, the field
rotates and has the following form:

(20)Bx=Bpcosωt=BXcosωt(21)By=Bpsinωt=BXsinωt(22)Bz=Ba=BZ.

Here, ω is the angular frequency of the spacecraft rotation, usually
determined from sun sensor or star tracker measurements, and t denotes the
time. Inserting these relations in Eq. (10) yields the
expected temperature-corrected output of the magnetometer in sensor
coordinates. By applying the de-spun rotation matrix

(23)D=cosωt-sinωt0sinωtcosωt0001

to Eq. (10) to transform BST into a
non-orthogonal, de-spun coordinate system (directions X′, Y′, and Z′;
see Table 1), after sorting by frequency and phase of the terms,
and further dropping second-order factors, we obtain the following relations.
They are structurally similar to Eqs. (11a)–(11c) in Kepko et al. (1996),
but different in detail:

These equations show how the parameters affect the signal content at the spin
tone harmonics in the de-spun measurements. The first terms in all three
Eqs. (24)–(26) are the primary
measurements terms. In the spin plane, the ambient magnetic field only has a
BX=Bp component. Consequently, the first term of BX′ is
approximately Bp, while the first term of BY′ is
approximately 0 as ϕa≈0 and δϕS12≈0.
In the spin axis, we find BZ′≈Ba with
Ga≈1 and OS3≈0. In addition,
superposed first and second harmonic signals are expected as functions of the
calibration parameters. The first harmonic signals are described by the
second and third terms in Eqs. (24)–(26). In the spin plane,
Eqs. (24) and (25), the spin tone signals are the result of spin plane
offsets OS1∕2 and projections of the spin axis field Ba
onto the spin plane. In the spin axis, Eq. (26), first
harmonic disturbances are due to the projection of spin plane fields onto the
spin axis, the reason being an incorrect description of the spin axis
direction by the angles σPx∕y. Second harmonic signals are
only expected in the de-spun spin plane components (fourth and fifth terms of
Eqs. 24 and 25). These are due to a mismatch in spin
plane gains (parameter g) or an unaccounted non-orthogonality between the
spin plane sensor axes S1 and S2 (parameter δϕS12).

4 Favorable conditions for the determination of the calibration parameters

From the factors pertaining to the first and second harmonic terms
of BX′, BY′, and BZ′ (Eqs. 24 to 26), it
is possible to derive the conditions that should be favorable for the
determination of the eight previously mentioned parameters. These factors are

As can be seen, the first factor of the latter group (Eq. 29) is
dependent on Bp, the external field in the spin plane which we
assume to be constant, on Gp, the absolute gain factor in the
spin plane which should be approximately 1, and on 1/g-g, which is 0 only
if g=1. Hence, the presence of one part of the second harmonic
disturbance, though modulated by Bp, is ultimately only dependent
on g, the ratio of spin plane gains. Consequently, this relation can be
used to determine g correctly. The signal to do that, in particular,
the signal-to-noise ratio (SNR), is larger if Bp is larger. We
capture this relation in the second line of Table 2. As
the second harmonic disturbance in the spin plane is to be minimized to get g,
the natural fluctuations around that frequency (of amplitude F2p)
should also be low in the spin plane. The uncertainty in g is then expected
to be on the order of F2p∕Bp.

The same is true for the complementary factor, on the right side of
Eq. (29): this second harmonic disturbance is also modulated
by Bp. When g is accurately determined, then the
ϕa influence vanishes, and the entire factor can only vanish by correctly
choosing δϕS12. Hence, to determine this parameter
accurately, Bp should also be large and the natural fluctuations
at the second harmonic should be of low amplitude (low F2p).
The uncertainty ΔϕS12 of δϕS12
and, ultimately, ϕS12 is expected to be on the order of
2F2p∕Bp (see line 2 in Table 2).

Let us focus on Eq. (28). The spin frequency
disturbance is clearly modulated by Bp, as Ga should
be close to 1, so Bp benefits the SNR. These disturbances vanish
if σPx and σPy become 0, i.e., if they
are precisely determined. A low amplitude in the natural fluctuations at the
spin frequency along the spin axis Fa would also support the
determination. The uncertainty in σPx
and σPy is then expected to be on the order of
ΔσPx/y≈Fa/Bp (line 1 in Table 2).

The first set of factors in Eq. (27) pertain to the spin frequency
disturbances in the spin plane components. They consist of two parts: a spin
plane offset component OS1 or OS2 and a term
that is modulated by Ba and which may vanish if the difference
(σPx−δθS1) or
(σPy−δθS2) vanishes. Obviously, if
Ba vanishes, then the spin plane spin frequency disturbances can
only come from the spin plane offset components. Hence, for their
determination, it is beneficial if the spin axis field Ba is low
and if the natural fluctuation level around the spin frequency in the spin
plane Fp is low. The uncertainty in OS1 and OS2
is then expected to be on the order of Fp+BaΔσPx/y+BaΔθS1/2
(see line 3 of Table 2).

The remaining elevation angles δθS1 and δθS2
are most difficult to determine: it is
beneficial if the spin axis field Ba is high. In addition,
however, it is necessary that the spin axis itself is determined well, as the
parameters σPx and σPy equally
influence the spin tone signal in the spin plane
as δθS1 and δθS2. Note that
σPx and σPy can be independently
determined by minimizing the spin frequency disturbances in the spin axis
component. Fp should again be low. Altogether, the uncertainty
in δθS1∕2 is on the order of
ΔθS1/2≈Fp/Ba+ΔOS1/2/Ba+ΔσPx/y
(see line 4 of Table 2).

Based on the findings from the previous section, we propose algorithms to
determine the eight spin-related parameters in an iterative manner
(Sect. 5.2 to 5.5). The algorithms are based on
computing estimates of the parameters for short intervals and evaluate the
uncertainties of those estimates based on the uncertainties indicated in
Table 2. Then, the estimates with uncertainties below a
certain acceptable threshold are chosen to form the basis of one parameter correction.

5.1 Precalibration

The temperature-dependent gains
GT(Ts, Te) determined on the ground
should be used to convert the raw magnetometer output BS to
a precalibrated, temperature-corrected intermediate product BST,
according to Eq. (8).

The offset vector OS and the calibration matrices Φ,
Σ, Γ, and G
should be initiated with the best known values at the time of calibration. At
the beginning, these will be ground-obtained values:

for θS1, θS2, ϕS12,
and OS from ground magnetometer calibration,

for ϕa from nominal spacecraft design or mirror- and laser-based alignment measurements,

for σPx and σPy from an initial
estimate of the spin axis direction (alternatively σPx=σPy=0
may be chosen),

and Gp=Ga=g=1 due to precalibration.

If in-flight calibration has already taken place, then these values will be
superseded by better in-flight determined values.

5.2 Calibration of the spin axis direction

The entire interval of magnetic field measurements should be divided into
small (overlapping) subintervals of length tint=2πn/ω,
with n∈ℕ. The factor n should not be too small;
hence, the subintervals should contain a number of spin periods, so that the
spin tone at the spin frequency and also the power around that frequency can
be accurately determined. On the other hand, subintervals should not be too
large, so that the field or environmental conditions can be assumed constant.

For each of the subintervals, the uncertainties ΔσPx/y≈Fa/Bp need to be
calculated (line 1 in Table 2). Conservatively, we
choose Bp to be the minimal modulus of the spin plane field over the
subinterval: min(Bx2+By2). Fa can be
estimated by taking the maximum of the discrete Fourier components Fa±
of the spin axis magnetic field Bz at frequencies ω± that
are slightly over and under the spin frequency:
ω±=2πn±/tint,
with n±∈N and slightly over and under n:

(30)Fa±=FBz,ω±=2N∑k=0N-1Bzt0+kδtexp-iω±kδt(31)Fa=maxFa±.

Here t0 is the start of a subinterval considered, N is the number of
magnetic field measurement samples in that subinterval, and δt is the
sampling period. The frequencies ω+ and ω− should
sufficiently differ from ω to avoid leakage from spin tone. However,
ω+ and ω− should also be close enough to ω so that
the amplitudes at those frequencies resemble the natural amplitude level at
the spin frequency. Note that the optimal choice of ω+
and ω− is subinterval-specific. In practice, however, fixed frequencies
can be used that are at some distance |ω±–ω|, if that
distance is safely larger than usual spin tone spectral peak widths.

From here on, we use ℱ(B, ω) to denote the Fourier component
of B at frequency ω. It should be noted that it may be recommended
to de-trend the B data before computing ℱ(B, ω), by simply
subtracting a linear fit. Linear trends will not occur if the external field
can be assumed to be constant. In many real applications, however, the
spacecraft will move through field gradients during subintervals considered,
and, in these cases, the linear trend in the field measurements will increase
the spectral content across the spectrum.

Parameter estimates σPx and σPy are
determined by minimization of the spin tone Sa in the spin axis
component: Sa=ℱ(Bz, ω). This minimization is
performed for each subinterval. Hence, we obtain for each subinterval one
estimate for σPx, for σPy, and for the
uncertainty ΔσPx∕y. A final parameter update
for σPx and σPy for the entire interval of
interest may be obtained by selecting the most accurate subinterval estimates
of those parameters, pertaining to minimal uncertainties ΔσPx∕y.
From those estimates, the median or average may be
computed. The selection of the best estimates can be threshold-based with
respect to ΔσPx∕y.

5.3 Calibration of gain ratio and azimuthal angle

As detailed in the previous section, Sect. 5.2, an interval of
interest is divided into short (overlapping) subintervals of length
tint=2πn/ω. For each of these subintervals, the
uncertainties Δg and ΔϕS12 are computed (see
line 2 in Table 2). Therefore, the fluctuation amplitudes
F2p=max(F(Bx2+By2, 2ω±)) need
to be computed, with 2ω±=4πn±/tint and
n±∈N slightly over and under n. Subsequently, the parameters g
and δϕS12 are determined for each subinterval by
minimization of S2p=F(Bx2+By2, 2ω).
From the set of g and δϕS12
estimates from all subintervals, those associated with the lowest uncertainties can
be chosen to yield final updates for g and δϕS12.

It should be noted that we are using here the modulus of the spin plane field
(Bx2+By2=BX2+BY2) to compute F2p
and S2p instead of any individual spin plane component (BX
or BY) in a de-spun coordinate system, as would be suggested by the
analytical treatment outlined in Sect. 3. Both approaches
(using the modulus or a de-spun component) are, however, mathematically
equivalent. To show this, we can compute BX′2+BY′2 using
Eqs. (24) and (25). From the sum BX′2+BY′2,
only those terms are large which contain the first term of
Eq. (24) because all other terms are products of multiple
small factors and, hence, can be omitted due to linearization. Taking that
into account, we obtain BX′2+BY′2≈BX′. Hence,
Bx2+By2 contains the field and variations corresponding to the
de-spun component, along which the external field is pointing. Evaluating
Bx2+By2 is hence equivalent to evaluating BX if the field
points in the X direction. This result is based on the assumption of the
spin tone and second harmonic terms being small in comparison to the constant
spin plane magnetic field, which should be fulfilled even in low field
conditions if the initial set of calibration parameters is not too
inaccurate. Note also that it is not possible to obtain additional
information with respect to the calibration parameters by evaluating the
field-perpendicular component BY because the coefficients pertaining to
the sin and cos terms of BX′ and BY′ are the same (compare
Eqs. 24 and 25).

The equivalence of the approaches (using the modulus or a de-spun component)
brings up two questions: (i) why did we not use the modulus when calculating
the influences of the spin-related parameters in Sect. 3 and
(ii) why would we prefer using the modulus over a de-spun component here and
in any practical application of the calibration algorithms outlined in this
Sect. 5? The answer to question (i) is that the mathematical
treatment of the modulus is slightly more involved than the treatment of the
individual de-spun coordinates. Furthermore, Sect. 3 follows
the approach of Kepko et al. (1996), who also use de-spun coordinates. Therefore,
our results from Sect. 3 become directly comparable to the
results of their study. In their and our analytical treatments, the
de-spinning process is exactly defined, perfectly known, and accurate. Hence,
it does not introduce additional uncertainty into the calibration process.
This latter statement is not true in any real application, which leads us to
the answer of question (ii): the modulus of the spin plane field is readily
available in any spinning coordinate system. De-spinning is not necessary for
magnetometer calibration, and it is not advised because it could introduce
additional, unnecessary uncertainty.

5.4 Calibration of spin plane offsets

The uncertainties for each subinterval are computed as suggested in line 3 of
Table 2. Therefore, the maximum of the spin axis field
(Ba=max|Bz|) over each subinterval should be
used. Fp is evaluated following Eq. (31).
Furthermore, estimates for the uncertainties ΔσPx∕y
and ΔθS1∕2 need to be obtained. These may be based on
the variability of the selected estimates of σPx∕y (see
Sect. 5.2) and δθS1∕2 (see Sect. 5.5) used to compute the final values of those parameters.

The offset estimates OS1 and OS2 are determined
for each subinterval by minimization of Sp=F(Bx2+By2, ω). From the set of
OS1 and OS2 estimates, the most accurate can be chosen to compute
final updates for the spin plane offsets. It should be noted that the offsets
are known to be the most variable parameters. Hence, it could be desirable to
compute final offset updates more often than updates of other spin-related
parameters, if possible.

5.5 Calibration of elevation angles

The uncertainties for each subinterval are computed as suggested in line 4 of
Table 2, this time using Ba=min|Bz|.
Estimates for the uncertainties ΔσPx∕y and ΔOS1∕2
need to be obtained, e.g., from the variability of
selected σPx∕y and OS1∕2 estimates.
Subsequently, the elevation angles δθS1
and δθS2 are determined for each subinterval by
minimization of Sp=F(Bx2+By2, ω).
From the set of δθS1 and δθS2
estimates, the most accurate can be chosen (lowest uncertainties) to yield
final updates of those parameters.

It should be noted that the same quantity Sp is minimized to
obtain the elevation angles δθS1
and δθS2 and the offset components OS1
and OS2. Hence, the final selection of estimates according to the
uncertainties ΔθS1∕2 and ΔOS1∕2,
which are heavily dependent on |Bz|, is very important here. In
low |Bz|, minimization of Sp yields the offset components,
whereas in high fields the offsets do not matter any more and any spin tone
may safely be attributed to an incorrect choice of the elevation angles if
the spin axis direction is precisely known.

5.6 Exclusion of data intervals

Certain intervals may be excluded from parameter determination, as some of
the underlying assumptions may not be met well. For instance, intervals
featuring large spacecraft and sensor temperature changes should be avoided, as
parameters may vary within such intervals. Hence, uncertainties in the
parameters may be significantly higher than what is reflected in the
uncertainty estimates stated in Table 2. Large temperature
variations are expected during eclipse intervals, when the spacecraft is in
shadow (e.g., of Earth), and hours after eclipse intervals as temperatures
relax to stationary values. Furthermore, magnetic field measurements at
saturation levels need to be avoided. Lastly, intervals during and after
spacecraft maneuvers may be problematic for calibration, as the spin axis
will fluctuate during maneuvers and nutation may be visible for periods of
time after maneuvers. It should be noted that all these considerations are
spacecraft- and orbit-specific.

To ascertain the accuracies that parameters may be determined with in
different regions of near-Earth space, on a highly elliptical orbit around
Earth, we apply the algorithms detailed above to two days (20 and 21 July 2007)
of THEMIS-C (Angelopoulos, 2008) fluxgate magnetometer (FGM) data
(Auster et al., 2008). The data are available at 4 Hz sampling
frequency (data product: FGL); they are already fully calibrated and the
applied calibration parameters do not change over the two days considered.
The magnitudes of the magnetic field along the spin axis |Bz| and in the
spin plane Bx2+By2 are displayed in Fig. 2a in red
and blue, respectively.

Figure 2From top to bottom: (a) magnitude of the spin axis and spin
plane magnetic fields in red and blue, (b) omnidirectional ion
spectral energy flux densities, and (c–f) uncertainties of the
estimates of the respective calibration parameters calculated in accordance
with Table 2.

The different regions that THEMIS-C went through during these two particular
days are best identified using the omnidirectional ion spectral energy flux
densities, measured by the electrostatic analyzer
(ESA; McFadden et al., 2008) and displayed in Fig. 2b. At the
beginning of 20 July 2007, THEMIC-C is located in the dayside magnetosheath.
This is clearly visible in the broad ion energy spectrum, which is
characteristic of the thermalized solar wind plasma population present
downstream of the bow shock. THEMIS-C fully transitions through the
magnetopause into the magnetosphere at about 05:06 UT, moving inbound towards
perigee at about 10:27 UT. At about 15:33 UT, THEMIS-C went back into the
magnetosheath until about 22:33 UT, when it transitioned through the bow
shock into the solar wind, characterized by a narrow energy signature,
corresponding to a cold plasma moving at solar wind speed. On 21 July,
THEMIS-C went back into the magnetosheath at about 06:07 UT, and then went
further into the magnetosphere at about 10:53 UT. The perigee pass on that
day took place at about 17:47 UT.

As can be seen in Fig. 2a, the solar wind interval is characterized by low
magnetic fields, typically below 10 nT. In the dayside
magnetosheath, the field strength is somewhat higher, on the order of a few
tens of nanotesla, and highly fluctuating. Inside the magnetosphere, the
fluctuation level is again low. The lowest field strengths of a few tens of
nanotesla are measured on the earthward side of the magnetopause, so just
inside the inner magnetosphere. The field strength continuously increases
towards Earth. On this particular THEMIS-C orbit, field strengths on the
order of 104 nT are reached along the spin axis and in the spin
plane close to perigee. As discussed above, both the fluctuation levels and
field strengths have a major influence on the expected uncertainties of the
calibration parameter estimates. We determine the parameters and the
corresponding uncertainties for overlapping subintervals of 100 spin periods
each, a spin period lasting for approximately 3 s. Hence, the
subintervals have interval lengths of approximately 5 min. Note that we
do not consider subintervals containing fields above 2×104 nT,
due to FGM instrument saturation, and also excluded
intervals in eclipse (Earth shadow) around perigee that lasted for
approximately 22 min per orbit. Over the remaining times,
temperature variations at the FGM sensor and electronics are limited to
within 3∘. In the THEMIS case, these small variations are not expected
to have a significant influence on the calibration parameters.

Subinterval lengths of 100 spin periods ensure good estimates of the power at
around (and double) the spin frequency ω=2π/(3s)≈2 rad s−1, while the calibration parameters to be determined
and the ambient magnetic field conditions may well be considered constant
over such short intervals. Estimates of the power Fa±,
Fp±, and F2p± around (double) the spin
frequency are taken at 85 % and 115 % of ω and at 185 % and
215 % of ω, respectively. Following the equations from
Table 2 and from Sect. 5.2
to 5.5 above, we determine the uncertainties for the calibration
parameter estimates for all subintervals. They are shown in Fig. 2c–f.

Figure 3Calibration parameter estimates as a function of their respective
uncertainties. Threshold levels for blue and red marked estimates are the
same as in Fig. 2. Panels (b–d), (g), and
(h) have secondary axes in orange, showing parameter values in
degrees.

Figure 2c and d show the uncertainties Δg=ΔϕS12/2
and ΔσPx∕y. For the
corresponding parameters, uncertainties on the order of
10−4 (rad) are generally acceptable. In the case of the gain
ratio parameter g, an error of 10−4 would translate into an absolute
error of 1 nT in 10 000 nT fields. With respect to the
angle ϕS12 (or δϕS12) and to the spin
axis angles σPx∕y, an error of 1×10-4 rad is equivalent to approximately 0.5 % of a degree.
Uncertainties below 10−4 (rad) are marked in blue in
Fig. 2c and d. As can be seen, estimates of g,
ϕS12, and σPx∕y with uncertainties below
this threshold can be obtained almost everywhere in the inner magnetosphere,
where fields are relatively stable, but not in the magnetosheath
(fluctuations too high) or in the solar wind (fields too low). Estimates
associated with uncertainties below 10−5 (rad) are marked in
red in Fig. 2c and d. These are only obtained in the
regions of highest ambient fields, close to perigee.

The parameter estimates themselves (g, δϕS12,
σPx, and σPy) are shown in Fig. 3a to d as a
function of their respective uncertainties.
Again, uncertainty thresholds of 10−4 and
10−5 (rad) are marked in blue and red, respectively. Taking
the averages of the estimates associated with uncertainties below
10−5 (rad), we obtain

Here, the error values are the corresponding standard deviations of the
estimates. We see that all parameters are close to 0 (or 1 in the case of g);
an update may only be advised for σPy, as its
deviation from 0 is significantly larger than the error value (see Fig. 3d).

In Fig. 3c and d, a split in values associated with low
uncertainties can be clearly seen. A closer look at this phenomenon reveals
that lower/higher σPx∕σPy values
correspond to times before/after perigee passes. Hence, the spin axis
direction in the orthogonalized sensor package coordinate system changes
during perigee. This might be related to a temperature-driven change in
spacecraft geometry, i.e., in boom alignment to the spacecraft body,
occurring in eclipse during perigee passes.

In order to calculate the uncertainties of the offset and elevation angle
estimates (ΔOS1∕2 and ΔθS1∕2; see
lines 3 and 4 in Table 2), we have to assume uncertainties
in the knowledge of the spin axis direction angles (ΔσPx∕y),
the offsets (ΔOS1∕2), and the
elevation angles (ΔθS1∕2). Based on Eq. (35),
we set ΔσPx/y=6×10-5 rad. Furthermore, as we can justify this a posteriori based on
Eqs. (37) and (39), we set
ΔOS1/2=25 pT and ΔθS1/2=7×10-4 rad.
Therefore, we obtain the uncertainty estimates per subinterval shown in Fig. 2e and f.

The offsets directly influence the absolute accuracies of the magnetic field
measurements. Typically, uncertainties on the order of or below
0.1 nT are desired in low fields. Uncertainties meeting this
threshold are marked in blue in Fig. 2e. As can be seen,
corresponding offset estimates can be routinely obtained in the solar wind,
due to the low fields, and also in the outer, low field parts of the inner
magnetosphere. Within the magnetosheath, however, many estimate uncertainties
surpass the threshold as the fluctuations levels are too high for accurate
offset determinations. Estimates with uncertainties below 10 pT
(in red) can only be obtained in the solar wind at low fields. From those
(red dots in Fig. 3e and f), we obtain average offsets of

(36)OS1=(-0.007±0.023)nT(37)OS2=(0.036±0.025)nT.

The error values here motivate the choice of ΔOS1∕2 for
the computation of the uncertainties of the elevation angles. These angles
should also be known to the order of 10−4 rad. Unfortunately,
estimates with uncertainties lower than this threshold are only obtained in
very high fields, close to perigee, as can be seen by the red dots in
Figs. 2f and 3g and h. The blue dots correspond to the
lower threshold of 10−3 rad in this case, already equivalent
to 5.7 % of a degree in angular uncertainty. From the δθS1∕2 estimates
pertaining to uncertainties lower than 10−4 rad we obtain the following averages:

(38)δθS1=(5±4)×10-4rad(39)δθS2=(3±7)×10-4rad.

Within the group of sensor orthogonalization angles (δθS1,
δθS2, and δϕS12) and spin axis angles (σPx
and σPy), these elevation angles are least accurately defined.
Apparently, it is difficult to determine them to accuracies on the order of
10−4 rad or better. When determined on the ground, in higher
fields, the parameters δθS1∕2 may however be better
determined, with lower uncertainties than 10−4 rad. Hence,
regular in-flight updating of these parameters may not be recommended, as
those updates may introduce unnecessary jitter without any benefit for the
overall accuracy of the magnetometer calibration.

The only purpose of linearizing the calibration equation in
Sect. 3 is to obtain the uncertainty relations shown in the right
column of Table 2. These uncertainties
(ΔσPx∕y, Δg, ΔϕS12,
ΔOS1∕2, and ΔθS1∕2) are used to
evaluate which calibration parameter estimates (from which subintervals) are
most suitable to compute final calibration parameter updates. The estimates
themselves are calculated using the nonlinearized calibration equation, Eq. (9). Hence, simplifications and assumptions associated with
linearization do not influence the accuracy of those estimates. They only
influence the accuracy of their uncertainties.

To obtain these uncertainties, a series of assumptions need to be made, e.g.,
in the form of Eqs. (11) to (15). This
approach raises the question of how much the calibration parameters may be
allowed to deviate from the assumed values. As shown in Sect. 6, we are interested in the orders of magnitude (factor
of 10) of the uncertainties ΔσPx∕y, Δg,
ΔϕS12, ΔOS1∕2, and ΔθS1∕2. Conservatively limiting the errors in these
uncertainties to a factor of 2, and taking into account the multiplication of
parameters due to Eq. (10), we can allocate lower limit
individual admissible error factors of 24=1.189 to gain
factors g, Gp, and Ga, as well as
angles σPx, σPy, δθS1,
δθS2, δϕS12, and ϕa.
Such error factors of 1.189 or, equivalently,
deviations by 19 % from the assumed values are extraordinarily large when
compared to the accuracy in the knowledge of the calibration parameters and
of the geometry of the sensor and spacecraft, even before performing any
in-flight calibration. For the angles, this means that deviations from 0 by
up to 11∘ are acceptable (arcsin (19 %)); note that alignment
uncertainties and deviations from sensor orthogonality should usually be
lower than 1∘. For the gains, deviations by 19 % would be
acceptable; ground calibration, however, should reduce these deviations to
less than 0.1 %. Hence, the assumptions made to linearize the calibration
equation are not restrictive at all and can easily be fulfilled in practice.

The orthogonalization angles are known to be relatively stable when compared
to the spin axis direction angles. Fortunately, as shown in
Sect. 6, the spin axis angles can be updated with high accuracy
more regularly than the sensor elevation angles δθS1
and δθS2. The parameter decoupling introduced in
Sect. 2 pays off here, as spin axis variations do not require
redetermination of the sensor elevation angles as would be the case when
using the calibration Eq. (1) with the coupling matrix
(Eq. 2) instead of Eq. (9).

It should be noted that both Eqs. (1) and (9)
assume raw magnetometer outputs to be linearly transformable into accurate
magnetic field estimates. This assumption of linearity can only be fulfilled
to a certain degree when dealing with actual magnetometer hardware.
Nonlinearities (e.g., Auster et al., 2008; Russell et al., 2016) will adversely
affect the calibration as described here if not characterized, quantified,
and corrected beforehand, as they produce spin tone and higher harmonic
signals in the magnetic field measurements. THEMIS FGM data, for instance,
suffer from slight nonlinearities in digital-to-analogue converters that are
part of the magnetometer hardware. These are known from ground
characterization of the instruments and are routinely corrected in advance of
any in-flight calibration activities and/or any conversion of magnetometer
outputs into calibrated magnetic field measurements (Auster et al., 2008).

Assuming instrument linearity, the uncertainty-based approach to determining
the spin-related calibration parameters allows for a meaningful estimation of
the error alongside any parameter updates. These errors can be compared
to the uncertainties of the already known parameters, determined either on
the
ground or in flight. Therewith, it is possible to decide whether any update
of the calibration parameters is necessary or advised or, instead, would just
introduce unnecessary variations in the calibration parameters over time.

In addition, the availability of calibration parameter estimates associated
with low uncertainties, sufficient in number and quality, determines what is
possible in terms of cadence of parameter updates. This availability depends
on the orbit of the spacecraft (the presence in regions of certain field
conditions) and also on the spin period of the spacecraft. In general, short
spin periods (high spacecraft spin frequencies) are favorable, as they
increase the number of spins that may be taken into account in subintervals
of certain length. A larger number of spin periods reduces the influence of
natural field fluctuations at (double) the spin frequency, while short
subinterval lengths ensure the constancy of the parameters and environmental
conditions. In the given THEMIS-C example, the spin plane offsets OS1∕2
may be continuously tracked while the spacecraft remained
in the solar wind and in the low field parts of the magnetosphere. The spin
axis components σPx∕y, the gain ratio g, and azimuthal
orthogonalization angle ϕS12 can easily be determined
separately before and after each perigee pass, whereas accurate determinations of
the elevation angles θS1∕2 may only be possible when
taking into account estimates from several spacecraft orbits.

Finally we would like to note that the benefits of parameter decoupling
(i.e., a sensible choice of parameters when taking into account the behavior
of the magnetometer and spacecraft hardware) and of the uncertainty-based
determination of those parameters are not tied to the exact definitions of
the calibration equation (Eq. 9) and matrices (Eqs. 4
to 7). For example, the offsets may be applied in an orthogonal,
spacecraft-fixed coordinate system instead of in the sensor coordinate
system if the main contribution to the offsets is expected from spacecraft
stray fields at the sensor position. The order of the gain,
orthogonalization, and alignment matrices (here G,
Γ, Σ, and Φ) may be changed,
and/or the 12 degrees of freedom of the calibration parameters may be
distributed over a larger number of matrices and offset vectors to account
for changes pertaining to different parts of the magnetometer–spacecraft
system in different coordinate systems (e.g., see equations
in Fornaçon et al., 1999; Auster et al., 2008). Hence, while following the principles set
out in this paper, a different set of calibration parameters and
corresponding calibration equation may be specifically selected for each
magnetometer–spacecraft combination.

FP, DF, and WM conceived the study. KHF, HUA, and
IR contributed key knowledge and experience in spacecraft magnetometer calibration.
DC and YN helped with the derivation and interpretation of equations.

We acknowledge NASA contract NAS5-02099 and Vassilis Angelopoulos for use of
data from the THEMIS mission. Specifically, we acknowledge Charles W. Carlson
and James P. McFadden for use of ESA data and Karl-Heinz Glassmeier,
Hans-Ulrich Auster, and Wolfgang Baumjohann for the use of FGM data provided
under the lead of the Technical University of Braunschweig and with financial
support from the German Ministry for Economy and Technology and the German
Center for Aviation and Space (DLR) under contract 50 OC 0302.

Raw output of spacecraft magnetometers has to be converted into meaningful units and coordinate systems before it is usable for scientific applications. This conversion is defined by 12 calibration parameters, 8 of which are more easily determined in flight if the spacecraft is spinning. We present theory and advanced algorithms to determine these eight parameters. They take into account the physical magnetometer and spacecraft behavior, making them superior to previously published algorithms.

Raw output of spacecraft magnetometers has to be converted into meaningful units and coordinate...