View Javadoc
1   /*
2    * Copyright (C) 2022 Alberto Irurueta Carro (alberto@irurueta.com)
3    *
4    * Licensed under the Apache License, Version 2.0 (the "License");
5    * you may not use this file except in compliance with the License.
6    * You may obtain a copy of the License at
7    *
8    *         http://www.apache.org/licenses/LICENSE-2.0
9    *
10   * Unless required by applicable law or agreed to in writing, software
11   * distributed under the License is distributed on an "AS IS" BASIS,
12   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13   * See the License for the specific language governing permissions and
14   * limitations under the License.
15   */
16  package com.irurueta.navigation.inertial.calibration.magnetometer;
17  
18  import com.irurueta.algebra.AlgebraException;
19  import com.irurueta.algebra.Matrix;
20  import com.irurueta.algebra.Utils;
21  import com.irurueta.algebra.WrongSizeException;
22  import com.irurueta.navigation.LockedException;
23  import com.irurueta.navigation.NotReadyException;
24  import com.irurueta.navigation.inertial.BodyKinematics;
25  import com.irurueta.navigation.inertial.BodyMagneticFluxDensity;
26  import com.irurueta.navigation.inertial.calibration.CalibrationException;
27  import com.irurueta.navigation.inertial.calibration.MagneticFluxDensityTriad;
28  import com.irurueta.navigation.inertial.calibration.StandardDeviationBodyMagneticFluxDensity;
29  import com.irurueta.numerical.EvaluationException;
30  import com.irurueta.numerical.GradientEstimator;
31  import com.irurueta.numerical.fitting.FittingException;
32  import com.irurueta.numerical.fitting.LevenbergMarquardtMultiDimensionFitter;
33  import com.irurueta.numerical.fitting.LevenbergMarquardtMultiDimensionFunctionEvaluator;
34  import com.irurueta.statistics.MaxIterationsExceededException;
35  import com.irurueta.units.MagneticFluxDensity;
36  import com.irurueta.units.MagneticFluxDensityConverter;
37  import com.irurueta.units.MagneticFluxDensityUnit;
38  
39  import java.util.Collection;
40  
41  /**
42   * Abstract class to estimate magnetometer hard-iron biases, cross couplings and scaling factors.
43   * This calibrator uses Levenberg-Marquardt to find a minimum least squared
44   * error solution.
45   * <p>
46   * To use this calibrator at least 10 measurements taken at a single unknown position and unknown orientations
47   * when common z-axis is assumed, otherwise at least 13 measurements are required.
48   * <p>
49   * Measured magnetic flux density is assumed to follow the model shown below:
50   * <pre>
51   *     mBmeas = bm + (I + Mm) * mBtrue + w
52   * </pre>
53   * Where:
54   * - mBmeas is the measured magnetic flux density. This is a 3x1 vector.
55   * - bm is magnetometer hard-iron bias. Ideally, on a perfect magnetometer,
56   * this should be a 3x1 zero vector.
57   * - I is the 3x3 identity matrix.
58   * - Mm is the 3x3 soft-iron matrix containing cross-couplings and scaling
59   * factors. Ideally, on a perfect magnetometer, this should be a 3x3 zero
60   * matrix.
61   * - mBtrue is ground-truth magnetic flux density. This is a 3x1 vector.
62   * - w is measurement noise. This is a 3x1 vector.
63   * Notice that this calibrator assumes that all measurements are taken in a short span of time,
64   * where Earth magnetic field can be assumed to be constant.
65   *
66   * @param <C> Calibrator type.
67   * @param <L> Listener type.
68   */
69  public abstract class BaseMagneticFluxDensityNormMagnetometerCalibrator<
70          C extends BaseMagneticFluxDensityNormMagnetometerCalibrator<?, ?>,
71          L extends BaseMagneticFluxDensityNormMagnetometerCalibratorListener<C>> implements
72          MagnetometerNonLinearCalibrator, UnknownHardIronNonLinearMagnetometerCalibrator,
73          UnorderedStandardDeviationBodyMagneticFluxDensityMagnetometerCalibrator {
74  
75      /**
76       * Indicates whether by default a common z-axis is assumed for the accelerometer,
77       * gyroscope and magnetometer.
78       */
79      public static final boolean DEFAULT_USE_COMMON_Z_AXIS = false;
80  
81      /**
82       * Number of unknowns when common z-axis is assumed for the accelerometer,
83       * gyroscope and magnetometer.
84       */
85      private static final int COMMON_Z_AXIS_UNKNOWNS = 9;
86  
87      /**
88       * Number of unknowns for the general case.
89       */
90      private static final int GENERAL_UNKNOWNS = 12;
91  
92      /**
93       * Required minimum number of measurements when common z-axis is assumed.
94       */
95      public static final int MINIMUM_MEASUREMENTS_COMMON_Z_AXIS = COMMON_Z_AXIS_UNKNOWNS + 1;
96  
97      /**
98       * Required minimum number of measurements for the general case.
99       */
100     public static final int MINIMUM_MEASUREMENTS_GENERAL = GENERAL_UNKNOWNS + 1;
101 
102     /**
103      * Ground truth magnetic flux density norm to be expected at location where measurements have been made,
104      * expressed in Teslas (T).
105      */
106     protected Double groundTruthMagneticFluxDensityNorm;
107 
108     /**
109      * Levenberg-Marquardt fitter to find a non-linear solution.
110      */
111     private final LevenbergMarquardtMultiDimensionFitter fitter = new LevenbergMarquardtMultiDimensionFitter();
112 
113     /**
114      * Initial x-coordinate of hard-iron bias to be used to find a solution.
115      * This is expressed in Teslas (T).
116      */
117     private double initialHardIronX;
118 
119     /**
120      * Initial y-coordinate of hard-iron bias to be used to find a solution.
121      * This is expressed in Teslas (T).
122      */
123     private double initialHardIronY;
124 
125     /**
126      * Initial z-coordinate of hard-iron bias to be used to find a solution.
127      * This is expressed in Teslas (T).
128      */
129     private double initialHardIronZ;
130 
131     /**
132      * Initial x scaling factor.
133      */
134     private double initialSx;
135 
136     /**
137      * Initial y scaling factor.
138      */
139     private double initialSy;
140 
141     /**
142      * Initial z scaling factor.
143      */
144     private double initialSz;
145 
146     /**
147      * Initial x-y cross coupling error.
148      */
149     private double initialMxy;
150 
151     /**
152      * Initial x-z cross coupling error.
153      */
154     private double initialMxz;
155 
156     /**
157      * Initial y-x cross coupling error.
158      */
159     private double initialMyx;
160 
161     /**
162      * Initial y-z cross coupling error.
163      */
164     private double initialMyz;
165 
166     /**
167      * Initial z-x cross coupling error.
168      */
169     private double initialMzx;
170 
171     /**
172      * Initial z-y cross coupling error.
173      */
174     private double initialMzy;
175 
176     /**
177      * Contains a collection of body magnetic flux density measurements taken
178      * at a given position with different unknown orientations and containing the
179      * standard deviation of magnetometer measurements.
180      */
181     private Collection<StandardDeviationBodyMagneticFluxDensity> measurements;
182 
183     /**
184      * This flag indicates whether z-axis is assumed to be common for accelerometer,
185      * gyroscope and magnetometer.
186      * When enabled, this eliminates 3 variables from Mm matrix.
187      */
188     private boolean commonAxisUsed = DEFAULT_USE_COMMON_Z_AXIS;
189 
190     /**
191      * Listener to handle events raised by this calibrator.
192      */
193     private L listener;
194 
195     /**
196      * Estimated magnetometer hard-iron biases for each magnetometer axis
197      * expressed in Teslas (T).
198      */
199     private double[] estimatedHardIron;
200 
201     /**
202      * Estimated magnetometer soft-iron matrix containing scale factors
203      * and cross coupling errors.
204      * This is the product of matrix Tm containing cross coupling errors and Km
205      * containing scaling factors.
206      * So tat:
207      * <pre>
208      *     Mm = [sx    mxy  mxz] = Tm*Km
209      *          [myx   sy   myz]
210      *          [mzx   mzy  sz ]
211      * </pre>
212      * Where:
213      * <pre>
214      *     Km = [sx 0   0 ]
215      *          [0  sy  0 ]
216      *          [0  0   sz]
217      * </pre>
218      * and
219      * <pre>
220      *     Tm = [1          -alphaXy    alphaXz ]
221      *          [alphaYx    1           -alphaYz]
222      *          [-alphaZx   alphaZy     1       ]
223      * </pre>
224      * Hence:
225      * <pre>
226      *     Mm = [sx    mxy  mxz] = Tm*Km =  [sx             -sy * alphaXy   sz * alphaXz ]
227      *          [myx   sy   myz]            [sx * alphaYx   sy              -sz * alphaYz]
228      *          [mzx   mzy  sz ]            [-sx * alphaZx  sy * alphaZy    sz           ]
229      * </pre>
230      * This instance allows any 3x3 matrix however, typically alphaYx, alphaZx and alphaZy
231      * are considered to be zero if the accelerometer z-axis is assumed to be the same
232      * as the body z-axis. When this is assumed, myx = mzx = mzy = 0 and the Mm matrix
233      * becomes upper diagonal:
234      * <pre>
235      *     Mm = [sx    mxy  mxz]
236      *          [0     sy   myz]
237      *          [0     0    sz ]
238      * </pre>
239      * Values of this matrix are unit-less.
240      */
241     private Matrix estimatedMm;
242 
243     /**
244      * Estimated covariance matrix for estimated parameters.
245      */
246     private Matrix estimatedCovariance;
247 
248     /**
249      * Estimated chi square value.
250      */
251     private double estimatedChiSq;
252 
253     /**
254      * Estimated degrees of freedom of chi square value. Degrees of freedom is equal to the number of sampled data
255      * minus the number of estimated parameters.
256      */
257     private int estimatedChiSqDegreesOfFreedom;
258 
259     /**
260      * Estimated reduced chi square value. This is equal to estimated chi square value divided by its degrees of
261      * freedom. Ideally this value should be close to 1.0.
262      */
263     private double estimatedReducedChiSq;
264 
265     /**
266      * Estimated mean square error respect to provided measurements.
267      */
268     private double estimatedMse;
269 
270     /**
271      * Estimated probability of finding a smaller chi square value expressed as a value between 0.0 and 1.0. The smaller
272      * the found chi square value is, the better the fit of the estimated parameters to the actual parameter. Thus, the
273      * smaller the chance of finding a smaller chi square value, then the better the estimated fit is.
274      */
275     private double estimatedP;
276 
277     /**
278      * Estimated measure of quality of estimated fit as a value between 0.0 and 1.0. The larger the quality value is,
279      * the better the fit that has been estimated.
280      */
281     private double estimatedQ;
282 
283     /**
284      * Indicates whether calibrator is running.
285      */
286     private boolean running;
287 
288     /**
289      * Internally holds x-coordinate of measured magnetic flux density
290      * during calibration.
291      */
292     private double bmeasX;
293 
294     /**
295      * Internally holds y-coordinate of measured magnetic flux density
296      * during calibration.
297      */
298     private double bmeasY;
299 
300     /**
301      * Internally holds z-coordinate of measured magnetic flux density
302      * during calibration.
303      */
304     private double bmeasZ;
305 
306     /**
307      * Internally holds measured magnetic flux density during calibration
308      * expressed as a column matrix.
309      */
310     private Matrix bmeas;
311 
312     /**
313      * Internally holds cross-coupling errors during calibration.
314      */
315     private Matrix m;
316 
317     /**
318      * Internally holds inverse of cross-coupling errors during calibration.
319      */
320     private Matrix invM;
321 
322     /**
323      * Internally holds biases during calibration.
324      */
325     private Matrix b;
326 
327     /**
328      * Internally holds computed true magnetic flux density during
329      * calibration.
330      */
331     private Matrix btrue;
332 
333     /**
334      * Constructor.
335      */
336     protected BaseMagneticFluxDensityNormMagnetometerCalibrator() {
337     }
338 
339     /**
340      * Constructor.
341      *
342      * @param listener listener to handle events raised by this calibrator.
343      */
344     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(final L listener) {
345         this.listener = listener;
346     }
347 
348     /**
349      * Constructor.
350      *
351      * @param measurements collection of body magnetic flux density
352      *                     measurements with standard deviation of
353      *                     magnetometer measurements taken at the same
354      *                     position with zero velocity and unknown different
355      *                     orientations.
356      */
357     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
358             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements) {
359         this.measurements = measurements;
360     }
361 
362     /**
363      * Constructor.
364      *
365      * @param commonAxisUsed indicates whether z-axis is assumed to be common
366      *                       for the accelerometer, gyroscope and magnetometer.
367      */
368     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(final boolean commonAxisUsed) {
369         this.commonAxisUsed = commonAxisUsed;
370     }
371 
372     /**
373      * Constructor.
374      *
375      * @param initialHardIron initial hard-iron to find a solution.
376      * @throws IllegalArgumentException if provided hard-iron array does
377      *                                  not have length 3.
378      */
379     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(final double[] initialHardIron) {
380         try {
381             setInitialHardIron(initialHardIron);
382         } catch (final LockedException ignore) {
383             // never happens
384         }
385     }
386 
387     /**
388      * Constructor.
389      *
390      * @param initialHardIron initial hard-iron to find a solution.
391      * @throws IllegalArgumentException if provided hard-iron matrix is not
392      *                                  3x1.
393      */
394     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(final Matrix initialHardIron) {
395         try {
396             setInitialHardIron(initialHardIron);
397         } catch (final LockedException ignore) {
398             // never happens
399         }
400     }
401 
402     /**
403      * Constructor.
404      *
405      * @param initialHardIron initial hard-iron to find a solution.
406      * @param initialMm       initial soft-iron matrix containing scale factors
407      *                        and cross coupling errors.
408      * @throws IllegalArgumentException if provided hard-iron matrix is not
409      *                                  3x1 or if soft-iron matrix is not
410      *                                  3x3.
411      */
412     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(final Matrix initialHardIron, final Matrix initialMm) {
413         this(initialHardIron);
414         try {
415             setInitialMm(initialMm);
416         } catch (final LockedException ignore) {
417             // never happens
418         }
419     }
420 
421     /**
422      * Constructor.
423      *
424      * @param measurements collection of body magnetic flux density
425      *                     measurements with standard deviation of
426      *                     magnetometer measurements taken at the same
427      *                     position with zero velocity and unknown different
428      *                     orientations.
429      * @param listener     listener to handle events raised by this calibrator.
430      */
431     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
432             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final L listener) {
433         this(measurements);
434         this.listener = listener;
435     }
436 
437     /**
438      * Constructor.
439      *
440      * @param measurements   collection of body magnetic flux density
441      *                       measurements with standard deviation of
442      *                       magnetometer measurements taken at the same
443      *                       position with zero velocity and unknown different
444      *                       orientations.
445      * @param commonAxisUsed indicates whether z-axis is assumed to be common
446      *                       for the accelerometer, gyroscope and magnetometer.
447      */
448     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
449             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed) {
450         this(measurements);
451         this.commonAxisUsed = commonAxisUsed;
452     }
453 
454     /**
455      * Constructor.
456      *
457      * @param measurements   collection of body magnetic flux density
458      *                       measurements with standard deviation of
459      *                       magnetometer measurements taken at the same
460      *                       position with zero velocity and unknown different
461      *                       orientations.
462      * @param commonAxisUsed indicates whether z-axis is assumed to be common
463      *                       for the accelerometer, gyroscope and magnetometer.
464      * @param listener       listener to handle events raised by this calibrator.
465      */
466     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
467             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed,
468             final L listener) {
469         this(measurements, commonAxisUsed);
470         this.listener = listener;
471     }
472 
473     /**
474      * Constructor.
475      *
476      * @param measurements    collection of body magnetic flux density
477      *                        measurements with standard deviation of
478      *                        magnetometer measurements taken at the same
479      *                        position with zero velocity and unknown different
480      *                        orientations.
481      * @param initialHardIron initial hard-iron to find a solution.
482      * @throws IllegalArgumentException if provided hard-iron array does
483      *                                  not have length 3.
484      */
485     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
486             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final double[] initialHardIron) {
487         this(initialHardIron);
488         this.measurements = measurements;
489     }
490 
491     /**
492      * Constructor.
493      *
494      * @param measurements    collection of body magnetic flux density
495      *                        measurements with standard deviation of
496      *                        magnetometer measurements taken at the same
497      *                        position with zero velocity and unknown different
498      *                        orientations.
499      * @param initialHardIron initial hard-iron to find a solution.
500      * @param listener        listener to handle events raised by this calibrator.
501      * @throws IllegalArgumentException if provided hard-iron array does
502      *                                  not have length 3.
503      */
504     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
505             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final double[] initialHardIron,
506             final L listener) {
507         this(measurements, initialHardIron);
508         this.listener = listener;
509     }
510 
511     /**
512      * Constructor.
513      *
514      * @param measurements    collection of body magnetic flux density
515      *                        measurements with standard deviation of
516      *                        magnetometer measurements taken at the same
517      *                        position with zero velocity and unknown different
518      *                        orientations.
519      * @param commonAxisUsed  indicates whether z-axis is assumed to be common
520      *                        for the accelerometer, gyroscope and magnetometer.
521      * @param initialHardIron initial hard-iron to find a solution.
522      * @throws IllegalArgumentException if provided hard-iron array does
523      *                                  not have length 3.
524      */
525     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
526             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed,
527             final double[] initialHardIron) {
528         this(measurements, initialHardIron);
529         this.commonAxisUsed = commonAxisUsed;
530     }
531 
532     /**
533      * Constructor.
534      *
535      * @param measurements    collection of body magnetic flux density
536      *                        measurements with standard deviation of
537      *                        magnetometer measurements taken at the same
538      *                        position with zero velocity and unknown different
539      *                        orientations.
540      * @param commonAxisUsed  indicates whether z-axis is assumed to be common
541      *                        for the accelerometer, gyroscope and magnetometer.
542      * @param initialHardIron initial hard-iron to find a solution.
543      * @param listener        listener to handle events raised by this calibrator.
544      * @throws IllegalArgumentException if provided hard-iron array does
545      *                                  not have length 3.
546      */
547     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
548             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed,
549             final double[] initialHardIron, final L listener) {
550         this(measurements, commonAxisUsed, initialHardIron);
551         this.listener = listener;
552     }
553 
554     /**
555      * Constructor.
556      *
557      * @param measurements    collection of body magnetic flux density
558      *                        measurements with standard deviation of
559      *                        magnetometer measurements taken at the same
560      *                        position with zero velocity and unknown different
561      *                        orientations.
562      * @param initialHardIron initial hard-iron to find a solution.
563      * @throws IllegalArgumentException if provided hard-iron matrix is not
564      *                                  3x1.
565      */
566     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
567             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final Matrix initialHardIron) {
568         this(initialHardIron);
569         this.measurements = measurements;
570     }
571 
572     /**
573      * Constructor.
574      *
575      * @param measurements    collection of body magnetic flux density
576      *                        measurements with standard deviation of
577      *                        magnetometer measurements taken at the same
578      *                        position with zero velocity and unknown different
579      *                        orientations.
580      * @param initialHardIron initial hard-iron to find a solution.
581      * @param listener        listener to handle events raised by this calibrator.
582      * @throws IllegalArgumentException if provided hard-iron matrix is not
583      *                                  3x1.
584      */
585     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
586             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final Matrix initialHardIron,
587             final L listener) {
588         this(measurements, initialHardIron);
589         this.listener = listener;
590     }
591 
592     /**
593      * Constructor.
594      *
595      * @param measurements    collection of body magnetic flux density
596      *                        measurements with standard deviation of
597      *                        magnetometer measurements taken at the same
598      *                        position with zero velocity and unknown different
599      *                        orientations.
600      * @param commonAxisUsed  indicates whether z-axis is assumed to be common
601      *                        for the accelerometer, gyroscope and magnetometer.
602      * @param initialHardIron initial hard-iron to find a solution.
603      * @throws IllegalArgumentException if provided hard-iron matrix is not
604      *                                  3x1.
605      */
606     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
607             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed,
608             final Matrix initialHardIron) {
609         this(measurements, initialHardIron);
610         this.commonAxisUsed = commonAxisUsed;
611     }
612 
613     /**
614      * Constructor.
615      *
616      * @param measurements    collection of body magnetic flux density
617      *                        measurements with standard deviation of
618      *                        magnetometer measurements taken at the same
619      *                        position with zero velocity and unknown different
620      *                        orientations.
621      * @param commonAxisUsed  indicates whether z-axis is assumed to be common
622      *                        for the accelerometer, gyroscope and magnetometer.
623      * @param initialHardIron initial hard-iron to find a solution.
624      * @param listener        listener to handle events raised by this calibrator.
625      * @throws IllegalArgumentException if provided hard-iron matrix is not
626      *                                  3x1.
627      */
628     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
629             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed,
630             final Matrix initialHardIron, final L listener) {
631         this(measurements, commonAxisUsed, initialHardIron);
632         this.listener = listener;
633     }
634 
635     /**
636      * Constructor.
637      *
638      * @param measurements    collection of body magnetic flux density
639      *                        measurements with standard deviation of
640      *                        magnetometer measurements taken at the same
641      *                        position with zero velocity and unknown different
642      *                        orientations.
643      * @param initialHardIron initial hard-iron to find a solution.
644      * @param initialMm       initial soft-iron matrix containing scale factors
645      *                        and cross coupling errors.
646      * @throws IllegalArgumentException if provided hard-iron matrix is not
647      *                                  3x1 or if soft-iron matrix is not
648      *                                  3x3.
649      */
650     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
651             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final Matrix initialHardIron,
652             final Matrix initialMm) {
653         this(initialHardIron, initialMm);
654         this.measurements = measurements;
655     }
656 
657     /**
658      * Constructor.
659      *
660      * @param measurements    collection of body magnetic flux density
661      *                        measurements with standard deviation of
662      *                        magnetometer measurements taken at the same
663      *                        position with zero velocity and unknown different
664      *                        orientations.
665      * @param initialHardIron initial hard-iron to find a solution.
666      * @param initialMm       initial soft-iron matrix containing scale factors
667      *                        and cross coupling errors.
668      * @param listener        listener to handle events raised by this calibrator.
669      * @throws IllegalArgumentException if provided hard-iron matrix is not
670      *                                  3x1 or if soft-iron matrix is not
671      *                                  3x3.
672      */
673     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
674             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final Matrix initialHardIron,
675             final Matrix initialMm, final L listener) {
676         this(measurements, initialHardIron, initialMm);
677         this.listener = listener;
678     }
679 
680     /**
681      * Constructor.
682      *
683      * @param measurements    collection of body magnetic flux density
684      *                        measurements with standard deviation of
685      *                        magnetometer measurements taken at the same
686      *                        position with zero velocity and unknown different
687      *                        orientations.
688      * @param commonAxisUsed  indicates whether z-axis is assumed to be common
689      *                        for the accelerometer, gyroscope and magnetometer.
690      * @param initialHardIron initial hard-iron to find a solution.
691      * @param initialMm       initial soft-iron matrix containing scale factors
692      *                        and cross coupling errors.
693      * @throws IllegalArgumentException if provided hard-iron matrix is not
694      *                                  3x1 or if soft-iron matrix is not
695      *                                  3x3.
696      */
697     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
698             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed,
699             final Matrix initialHardIron, final Matrix initialMm) {
700         this(measurements, initialHardIron, initialMm);
701         this.commonAxisUsed = commonAxisUsed;
702     }
703 
704     /**
705      * Constructor.
706      *
707      * @param measurements    collection of body magnetic flux density
708      *                        measurements with standard deviation of
709      *                        magnetometer measurements taken at the same
710      *                        position with zero velocity and unknown different
711      *                        orientations.
712      * @param commonAxisUsed  indicates whether z-axis is assumed to be common
713      *                        for the accelerometer, gyroscope and magnetometer.
714      * @param initialHardIron initial hard-iron to find a solution.
715      * @param initialMm       initial soft-iron matrix containing scale factors
716      *                        and cross coupling errors.
717      * @param listener        listener to handle events raised by this calibrator.
718      * @throws IllegalArgumentException if provided hard-iron matrix is not
719      *                                  3x1 or if soft-iron matrix is not
720      *                                  3x3.
721      */
722     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
723             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed,
724             final Matrix initialHardIron, final Matrix initialMm, final L listener) {
725         this(measurements, commonAxisUsed, initialHardIron, initialMm);
726         this.listener = listener;
727     }
728 
729     /**
730      * Constructor.
731      *
732      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
733      * @throws IllegalArgumentException if provided magnetic flux norm value is negative.
734      */
735     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(final Double groundTruthMagneticFluxDensityNorm) {
736         internalSetGroundTruthMagneticFluxDensityNorm(groundTruthMagneticFluxDensityNorm);
737     }
738 
739     /**
740      * Constructor.
741      *
742      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
743      * @param listener                           listener to handle events raised by this calibrator.
744      * @throws IllegalArgumentException if provided magnetic flux norm value is negative.
745      */
746     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
747             final Double groundTruthMagneticFluxDensityNorm, final L listener) {
748         this(groundTruthMagneticFluxDensityNorm);
749         this.listener = listener;
750     }
751 
752     /**
753      * Constructor.
754      *
755      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
756      * @param measurements                       collection of body magnetic flux density
757      *                                           measurements with standard deviation of
758      *                                           magnetometer measurements taken at the same
759      *                                           position with zero velocity and unknown different
760      *                                           orientations.
761      * @throws IllegalArgumentException if provided magnetic flux norm value is negative.
762      */
763     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
764             final Double groundTruthMagneticFluxDensityNorm,
765             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements) {
766         this(groundTruthMagneticFluxDensityNorm);
767         this.measurements = measurements;
768     }
769 
770     /**
771      * Constructor.
772      *
773      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
774      * @param commonAxisUsed                     indicates whether z-axis is assumed to be common
775      *                                           for the accelerometer, gyroscope and magnetometer.
776      * @throws IllegalArgumentException if provided magnetic flux norm value is negative.
777      */
778     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
779             final Double groundTruthMagneticFluxDensityNorm, final boolean commonAxisUsed) {
780         this(groundTruthMagneticFluxDensityNorm);
781         this.commonAxisUsed = commonAxisUsed;
782     }
783 
784     /**
785      * Constructor.
786      *
787      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
788      * @param initialHardIron                    initial hard-iron to find a solution.
789      * @throws IllegalArgumentException if provided magnetic flux norm value is
790      *                                  negative, or if provided hard-iron array does
791      *                                  not have length 3.
792      */
793     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
794             final Double groundTruthMagneticFluxDensityNorm, final double[] initialHardIron) {
795         this(groundTruthMagneticFluxDensityNorm);
796         try {
797             setInitialHardIron(initialHardIron);
798         } catch (final LockedException ignore) {
799             // never happens
800         }
801     }
802 
803     /**
804      * Constructor.
805      *
806      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
807      * @param initialHardIron                    initial hard-iron to find a solution.
808      * @throws IllegalArgumentException if provided magnetic flux norm value is
809      *                                  negative, or if provided hard-iron matrix is not
810      *                                  3x1.
811      */
812     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
813             final Double groundTruthMagneticFluxDensityNorm, final Matrix initialHardIron) {
814         this(groundTruthMagneticFluxDensityNorm);
815         try {
816             setInitialHardIron(initialHardIron);
817         } catch (final LockedException ignore) {
818             // never happens
819         }
820     }
821 
822     /**
823      * Constructor.
824      *
825      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
826      * @param initialHardIron                    initial hard-iron to find a solution.
827      * @param initialMm                          initial soft-iron matrix containing scale factors
828      *                                           and cross coupling errors.
829      * @throws IllegalArgumentException if provided magnetic flux norm value is
830      *                                  negative, or if provided hard-iron matrix is not
831      *                                  3x1 or if soft-iron matrix is not
832      *                                  3x3.
833      */
834     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
835             final Double groundTruthMagneticFluxDensityNorm, final Matrix initialHardIron, final Matrix initialMm) {
836         this(groundTruthMagneticFluxDensityNorm, initialHardIron);
837         try {
838             setInitialMm(initialMm);
839         } catch (final LockedException ignore) {
840             // never happens
841         }
842     }
843 
844     /**
845      * Constructor.
846      *
847      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
848      * @param measurements                       collection of body magnetic flux density
849      *                                           measurements with standard deviation of
850      *                                           magnetometer measurements taken at the same
851      *                                           position with zero velocity and unknown different
852      *                                           orientations.
853      * @param listener                           listener to handle events raised by this calibrator.
854      * @throws IllegalArgumentException if provided magnetic flux norm value is negative.
855      */
856     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
857             final Double groundTruthMagneticFluxDensityNorm,
858             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final L listener) {
859         this(groundTruthMagneticFluxDensityNorm, measurements);
860         this.listener = listener;
861     }
862 
863     /**
864      * Constructor.
865      *
866      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
867      * @param measurements                       collection of body magnetic flux density
868      *                                           measurements with standard deviation of
869      *                                           magnetometer measurements taken at the same
870      *                                           position with zero velocity and unknown different
871      *                                           orientations.
872      * @param commonAxisUsed                     indicates whether z-axis is assumed to be common
873      *                                           for the accelerometer, gyroscope and magnetometer.
874      * @throws IllegalArgumentException if provided magnetic flux norm value is negative.
875      */
876     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
877             final Double groundTruthMagneticFluxDensityNorm,
878             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed) {
879         this(groundTruthMagneticFluxDensityNorm, measurements);
880         this.commonAxisUsed = commonAxisUsed;
881     }
882 
883     /**
884      * Constructor.
885      *
886      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
887      * @param measurements                       collection of body magnetic flux density
888      *                                           measurements with standard deviation of
889      *                                           magnetometer measurements taken at the same
890      *                                           position with zero velocity and unknown different
891      *                                           orientations.
892      * @param commonAxisUsed                     indicates whether z-axis is assumed to be common
893      *                                           for the accelerometer, gyroscope and magnetometer.
894      * @param listener                           listener to handle events raised by this calibrator.
895      * @throws IllegalArgumentException if provided magnetic flux norm value is negative.
896      */
897     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
898             final Double groundTruthMagneticFluxDensityNorm,
899             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed,
900             final L listener) {
901         this(groundTruthMagneticFluxDensityNorm, measurements, commonAxisUsed);
902         this.listener = listener;
903     }
904 
905     /**
906      * Constructor.
907      *
908      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
909      * @param measurements                       collection of body magnetic flux density
910      *                                           measurements with standard deviation of
911      *                                           magnetometer measurements taken at the same
912      *                                           position with zero velocity and unknown different
913      *                                           orientations.
914      * @param initialHardIron                    initial hard-iron to find a solution.
915      * @throws IllegalArgumentException if provided magnetic flux norm value is negative,
916      *                                  or if provided hard-iron array does not have length 3.
917      */
918     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
919             final Double groundTruthMagneticFluxDensityNorm,
920             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final double[] initialHardIron) {
921         this(groundTruthMagneticFluxDensityNorm, initialHardIron);
922         this.measurements = measurements;
923     }
924 
925     /**
926      * Constructor.
927      *
928      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
929      * @param measurements                       collection of body magnetic flux density
930      *                                           measurements with standard deviation of
931      *                                           magnetometer measurements taken at the same
932      *                                           position with zero velocity and unknown different
933      *                                           orientations.
934      * @param initialHardIron                    initial hard-iron to find a solution.
935      * @param listener                           listener to handle events raised by this calibrator.
936      * @throws IllegalArgumentException if provided magnetic flux norm value is negative,
937      *                                  or if provided hard-iron array does not have length 3.
938      */
939     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
940             final Double groundTruthMagneticFluxDensityNorm,
941             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final double[] initialHardIron,
942             final L listener) {
943         this(groundTruthMagneticFluxDensityNorm, measurements, initialHardIron);
944         this.listener = listener;
945     }
946 
947     /**
948      * Constructor.
949      *
950      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
951      * @param measurements                       collection of body magnetic flux density
952      *                                           measurements with standard deviation of
953      *                                           magnetometer measurements taken at the same
954      *                                           position with zero velocity and unknown different
955      *                                           orientations.
956      * @param commonAxisUsed                     indicates whether z-axis is assumed to be common
957      *                                           for the accelerometer, gyroscope and magnetometer.
958      * @param initialHardIron                    initial hard-iron to find a solution.
959      * @throws IllegalArgumentException if provided magnetic flux norm value is negative,
960      *                                  or if provided hard-iron array does not have length 3.
961      */
962     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
963             final Double groundTruthMagneticFluxDensityNorm,
964             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed,
965             final double[] initialHardIron) {
966         this(groundTruthMagneticFluxDensityNorm, measurements, initialHardIron);
967         this.commonAxisUsed = commonAxisUsed;
968     }
969 
970     /**
971      * Constructor.
972      *
973      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
974      * @param measurements                       collection of body magnetic flux density
975      *                                           measurements with standard deviation of
976      *                                           magnetometer measurements taken at the same
977      *                                           position with zero velocity and unknown different
978      *                                           orientations.
979      * @param commonAxisUsed                     indicates whether z-axis is assumed to be common
980      *                                           for the accelerometer, gyroscope and magnetometer.
981      * @param initialHardIron                    initial hard-iron to find a solution.
982      * @param listener                           listener to handle events raised by this calibrator.
983      * @throws IllegalArgumentException if provided magnetic flux norm value is negative,
984      *                                  or if provided hard-iron array does not have length 3.
985      */
986     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
987             final Double groundTruthMagneticFluxDensityNorm,
988             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed,
989             final double[] initialHardIron, final L listener) {
990         this(groundTruthMagneticFluxDensityNorm, measurements, commonAxisUsed, initialHardIron);
991         this.listener = listener;
992     }
993 
994     /**
995      * Constructor.
996      *
997      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
998      * @param measurements                       collection of body magnetic flux density
999      *                                           measurements with standard deviation of
1000      *                                           magnetometer measurements taken at the same
1001      *                                           position with zero velocity and unknown different
1002      *                                           orientations.
1003      * @param initialHardIron                    initial hard-iron to find a solution.
1004      * @throws IllegalArgumentException if provided magnetic flux norm value is negative,
1005      *                                  or if provided hard-iron matrix is not 3x1.
1006      */
1007     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
1008             final Double groundTruthMagneticFluxDensityNorm,
1009             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final Matrix initialHardIron) {
1010         this(groundTruthMagneticFluxDensityNorm, initialHardIron);
1011         this.measurements = measurements;
1012     }
1013 
1014     /**
1015      * Constructor.
1016      *
1017      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
1018      * @param measurements                       collection of body magnetic flux density
1019      *                                           measurements with standard deviation of
1020      *                                           magnetometer measurements taken at the same
1021      *                                           position with zero velocity and unknown different
1022      *                                           orientations.
1023      * @param initialHardIron                    initial hard-iron to find a solution.
1024      * @param listener                           listener to handle events raised by this calibrator.
1025      * @throws IllegalArgumentException if provided magnetic flux norm value is negative,
1026      *                                  or if provided hard-iron matrix is not 3x1.
1027      */
1028     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
1029             final Double groundTruthMagneticFluxDensityNorm,
1030             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final Matrix initialHardIron,
1031             final L listener) {
1032         this(groundTruthMagneticFluxDensityNorm, measurements, initialHardIron);
1033         this.listener = listener;
1034     }
1035 
1036     /**
1037      * Constructor.
1038      *
1039      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
1040      * @param measurements                       collection of body magnetic flux density
1041      *                                           measurements with standard deviation of
1042      *                                           magnetometer measurements taken at the same
1043      *                                           position with zero velocity and unknown different
1044      *                                           orientations.
1045      * @param commonAxisUsed                     indicates whether z-axis is assumed to be common
1046      *                                           for the accelerometer, gyroscope and magnetometer.
1047      * @param initialHardIron                    initial hard-iron to find a solution.
1048      * @throws IllegalArgumentException if provided magnetic flux norm value is negative,
1049      *                                  or if provided hard-iron matrix is not 3x1.
1050      */
1051     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
1052             final Double groundTruthMagneticFluxDensityNorm,
1053             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed,
1054             final Matrix initialHardIron) {
1055         this(groundTruthMagneticFluxDensityNorm, measurements, initialHardIron);
1056         this.commonAxisUsed = commonAxisUsed;
1057     }
1058 
1059     /**
1060      * Constructor.
1061      *
1062      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
1063      * @param measurements                       collection of body magnetic flux density
1064      *                                           measurements with standard deviation of
1065      *                                           magnetometer measurements taken at the same
1066      *                                           position with zero velocity and unknown different
1067      *                                           orientations.
1068      * @param commonAxisUsed                     indicates whether z-axis is assumed to be common
1069      *                                           for the accelerometer, gyroscope and magnetometer.
1070      * @param initialHardIron                    initial hard-iron to find a solution.
1071      * @param listener                           listener to handle events raised by this calibrator.
1072      * @throws IllegalArgumentException if provided magnetic flux norm value is negative,
1073      *                                  or if provided hard-iron matrix is not 3x1.
1074      */
1075     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
1076             final Double groundTruthMagneticFluxDensityNorm,
1077             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed,
1078             final Matrix initialHardIron, final L listener) {
1079         this(groundTruthMagneticFluxDensityNorm, measurements, commonAxisUsed, initialHardIron);
1080         this.listener = listener;
1081     }
1082 
1083     /**
1084      * Constructor.
1085      *
1086      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
1087      * @param measurements                       collection of body magnetic flux density
1088      *                                           measurements with standard deviation of
1089      *                                           magnetometer measurements taken at the same
1090      *                                           position with zero velocity and unknown different
1091      *                                           orientations.
1092      * @param initialHardIron                    initial hard-iron to find a solution.
1093      * @param initialMm                          initial soft-iron matrix containing scale factors
1094      *                                           and cross coupling errors.
1095      * @throws IllegalArgumentException if provided magnetic flux norm value is negative,
1096      *                                  or if provided hard-iron matrix is not 3x1 or if
1097      *                                  soft-iron matrix is not 3x3.
1098      */
1099     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
1100             final Double groundTruthMagneticFluxDensityNorm,
1101             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final Matrix initialHardIron,
1102             final Matrix initialMm) {
1103         this(groundTruthMagneticFluxDensityNorm, initialHardIron, initialMm);
1104         this.measurements = measurements;
1105     }
1106 
1107     /**
1108      * Constructor.
1109      *
1110      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
1111      * @param measurements                       collection of body magnetic flux density
1112      *                                           measurements with standard deviation of
1113      *                                           magnetometer measurements taken at the same
1114      *                                           position with zero velocity and unknown different
1115      *                                           orientations.
1116      * @param initialHardIron                    initial hard-iron to find a solution.
1117      * @param initialMm                          initial soft-iron matrix containing scale factors
1118      *                                           and cross coupling errors.
1119      * @param listener                           listener to handle events raised by this calibrator.
1120      * @throws IllegalArgumentException if provided magnetic flux norm value is negative,
1121      *                                  or if provided hard-iron matrix is not 3x1 or if
1122      *                                  soft-iron matrix is not 3x3.
1123      */
1124     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
1125             final Double groundTruthMagneticFluxDensityNorm,
1126             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final Matrix initialHardIron,
1127             final Matrix initialMm, final L listener) {
1128         this(groundTruthMagneticFluxDensityNorm, measurements, initialHardIron, initialMm);
1129         this.listener = listener;
1130     }
1131 
1132     /**
1133      * Constructor.
1134      *
1135      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
1136      * @param measurements                       collection of body magnetic flux density
1137      *                                           measurements with standard deviation of
1138      *                                           magnetometer measurements taken at the same
1139      *                                           position with zero velocity and unknown different
1140      *                                           orientations.
1141      * @param commonAxisUsed                     indicates whether z-axis is assumed to be common
1142      *                                           for the accelerometer, gyroscope and magnetometer.
1143      * @param initialHardIron                    initial hard-iron to find a solution.
1144      * @param initialMm                          initial soft-iron matrix containing scale factors
1145      *                                           and cross coupling errors.
1146      * @throws IllegalArgumentException if provided magnetic flux norm value is negative,
1147      *                                  or if provided hard-iron matrix is not 3x1
1148      *                                  or if soft-iron matrix is not 3x3.
1149      */
1150     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
1151             final Double groundTruthMagneticFluxDensityNorm,
1152             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed,
1153             final Matrix initialHardIron, final Matrix initialMm) {
1154         this(groundTruthMagneticFluxDensityNorm, measurements, initialHardIron, initialMm);
1155         this.commonAxisUsed = commonAxisUsed;
1156     }
1157 
1158     /**
1159      * Constructor.
1160      *
1161      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm expressed in Teslas (T).
1162      * @param measurements                       collection of body magnetic flux density
1163      *                                           measurements with standard deviation of
1164      *                                           magnetometer measurements taken at the same
1165      *                                           position with zero velocity and unknown different
1166      *                                           orientations.
1167      * @param commonAxisUsed                     indicates whether z-axis is assumed to be common
1168      *                                           for the accelerometer, gyroscope and magnetometer.
1169      * @param initialHardIron                    initial hard-iron to find a solution.
1170      * @param initialMm                          initial soft-iron matrix containing scale factors
1171      *                                           and cross coupling errors.
1172      * @param listener                           listener to handle events raised by this calibrator.
1173      * @throws IllegalArgumentException if provided magnetic flux norm value is negative,
1174      *                                  or if provided hard-iron matrix is not 3x1
1175      *                                  or if soft-iron matrix is not 3x3.
1176      */
1177     protected BaseMagneticFluxDensityNormMagnetometerCalibrator(
1178             final Double groundTruthMagneticFluxDensityNorm,
1179             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements, final boolean commonAxisUsed,
1180             final Matrix initialHardIron, final Matrix initialMm, final L listener) {
1181         this(groundTruthMagneticFluxDensityNorm, measurements, commonAxisUsed, initialHardIron, initialMm);
1182         this.listener = listener;
1183     }
1184 
1185 
1186     /**
1187      * Gets ground truth magnetic flux density norm to be expected at location where measurements have been made,
1188      * expressed in Teslas (T).
1189      *
1190      * @return ground truth magnetic flux density or null.
1191      */
1192     public Double getGroundTruthMagneticFluxDensityNorm() {
1193         return groundTruthMagneticFluxDensityNorm;
1194     }
1195 
1196     /**
1197      * Gets ground truth magnetic flux density norm to be expected at location where measurements have been made.
1198      *
1199      * @return ground truth magnetic flux density or null.
1200      */
1201     public MagneticFluxDensity getGroundTruthMagneticFluxDensityNormAsMagneticFluxDensity() {
1202         return groundTruthMagneticFluxDensityNorm != null
1203                 ? new MagneticFluxDensity(groundTruthMagneticFluxDensityNorm, MagneticFluxDensityUnit.TESLA) : null;
1204     }
1205 
1206     /**
1207      * Gets ground truth magnetic flux density norm to be expected at location where measurements have been made.
1208      *
1209      * @param result instance where result will be stored.
1210      * @return true if ground truth magnetic flux density norm has been defined, false if it is not available yet.
1211      */
1212     public boolean getGroundTruthMagneticFluxDensityNormAsMagneticFluxDensity(final MagneticFluxDensity result) {
1213         if (groundTruthMagneticFluxDensityNorm != null) {
1214             result.setValue(groundTruthMagneticFluxDensityNorm);
1215             result.setUnit(MagneticFluxDensityUnit.TESLA);
1216             return true;
1217         } else {
1218             return false;
1219         }
1220     }
1221 
1222     /**
1223      * Gets initial x-coordinate of magnetometer hard-iron bias to be used
1224      * to find a solution.
1225      * This is expressed in Teslas (T).
1226      *
1227      * @return initial x-coordinate of magnetometer hard-iron bias.
1228      */
1229     @Override
1230     public double getInitialHardIronX() {
1231         return initialHardIronX;
1232     }
1233 
1234     /**
1235      * Sets initial x-coordinate of magnetometer hard-iron bias to be used
1236      * to find a solution.
1237      * This is expressed in Teslas (T).
1238      *
1239      * @param initialHardIronX initial x-coordinate of magnetometer
1240      *                         hard-iron bias.
1241      * @throws LockedException if calibrator is currently running.
1242      */
1243     @Override
1244     public void setInitialHardIronX(final double initialHardIronX) throws LockedException {
1245         if (running) {
1246             throw new LockedException();
1247         }
1248         this.initialHardIronX = initialHardIronX;
1249     }
1250 
1251     /**
1252      * Gets initial y-coordinate of magnetometer hard-iron bias to be used
1253      * to find a solution.
1254      * This is expressed in Teslas (T).
1255      *
1256      * @return initial y-coordinate of magnetometer hard-iron bias.
1257      */
1258     @Override
1259     public double getInitialHardIronY() {
1260         return initialHardIronY;
1261     }
1262 
1263     /**
1264      * Sets initial y-coordinate of magnetometer hard-iron bias to be used
1265      * to find a solution.
1266      * This is expressed in Teslas (T).
1267      *
1268      * @param initialHardIronY initial y-coordinate of magnetometer
1269      *                         hard-iron bias.
1270      * @throws LockedException if calibrator is currently running.
1271      */
1272     @Override
1273     public void setInitialHardIronY(final double initialHardIronY) throws LockedException {
1274         if (running) {
1275             throw new LockedException();
1276         }
1277         this.initialHardIronY = initialHardIronY;
1278     }
1279 
1280     /**
1281      * Gets initial z-coordinate of magnetometer hard-iron bias to be used
1282      * to find a solution.
1283      * This is expressed in Teslas (T).
1284      *
1285      * @return initial z-coordinate of magnetometer hard-iron bias.
1286      */
1287     @Override
1288     public double getInitialHardIronZ() {
1289         return initialHardIronZ;
1290     }
1291 
1292     /**
1293      * Sets initial z-coordinate of magnetometer hard-iron bias to be used
1294      * to find a solution.
1295      * This is expressed in meters Teslas (T).
1296      *
1297      * @param initialHardIronZ initial z-coordinate of magnetometer
1298      *                         hard-iron bias.
1299      * @throws LockedException if calibrator is currently running.
1300      */
1301     @Override
1302     public void setInitialHardIronZ(final double initialHardIronZ) throws LockedException {
1303         if (running) {
1304             throw new LockedException();
1305         }
1306         this.initialHardIronZ = initialHardIronZ;
1307     }
1308 
1309     /**
1310      * Gets initial x-coordinate of magnetometer hard iron bias to be used
1311      * to find a solution.
1312      *
1313      * @return initial x-coordinate of magnetometer hard-iron bias.
1314      */
1315     @Override
1316     public MagneticFluxDensity getInitialHardIronXAsMagneticFluxDensity() {
1317         return new MagneticFluxDensity(initialHardIronX, MagneticFluxDensityUnit.TESLA);
1318     }
1319 
1320     /**
1321      * Gets initial x-coordinate of magnetometer hard iron bias to be used
1322      * to find a solution.
1323      *
1324      * @param result instance where result will be stored.
1325      */
1326     @Override
1327     public void getInitialHardIronXAsMagneticFluxDensity(final MagneticFluxDensity result) {
1328         result.setValue(initialHardIronX);
1329         result.setUnit(MagneticFluxDensityUnit.TESLA);
1330     }
1331 
1332     /**
1333      * Sets initial x-coordinate of magnetometer hard iron bias to be used
1334      * to find a solution.
1335      *
1336      * @param initialHardIronX initial x-coordinate of magnetometer bias.
1337      * @throws LockedException if calibrator is currently running.
1338      */
1339     @Override
1340     public void setInitialHardIronX(final MagneticFluxDensity initialHardIronX) throws LockedException {
1341         if (running) {
1342             throw new LockedException();
1343         }
1344         this.initialHardIronX = convertMagneticFluxDensity(initialHardIronX);
1345     }
1346 
1347     /**
1348      * Gets initial y-coordinate of magnetometer hard iron bias to be used
1349      * to find a solution.
1350      *
1351      * @return initial y-coordinate of magnetometer hard-iron bias.
1352      */
1353     @Override
1354     public MagneticFluxDensity getInitialHardIronYAsMagneticFluxDensity() {
1355         return new MagneticFluxDensity(initialHardIronY, MagneticFluxDensityUnit.TESLA);
1356     }
1357 
1358     /**
1359      * Gets initial y-coordinate of magnetometer hard iron bias to be used
1360      * to find a solution.
1361      *
1362      * @param result instance where result will be stored.
1363      */
1364     @Override
1365     public void getInitialHardIronYAsMagneticFluxDensity(final MagneticFluxDensity result) {
1366         result.setValue(initialHardIronY);
1367         result.setUnit(MagneticFluxDensityUnit.TESLA);
1368     }
1369 
1370     /**
1371      * Sets initial y-coordinate of magnetometer hard iron bias to be used
1372      * to find a solution.
1373      *
1374      * @param initialHardIronY initial y-coordinate of magnetometer bias.
1375      * @throws LockedException if calibrator is currently running.
1376      */
1377     @Override
1378     public void setInitialHardIronY(final MagneticFluxDensity initialHardIronY) throws LockedException {
1379         if (running) {
1380             throw new LockedException();
1381         }
1382         this.initialHardIronY = convertMagneticFluxDensity(initialHardIronY);
1383     }
1384 
1385     /**
1386      * Gets initial z-coordinate of magnetometer hard iron bias to be used
1387      * to find a solution.
1388      *
1389      * @return initial z-coordinate of magnetometer hard-iron bias.
1390      */
1391     @Override
1392     public MagneticFluxDensity getInitialHardIronZAsMagneticFluxDensity() {
1393         return new MagneticFluxDensity(initialHardIronZ, MagneticFluxDensityUnit.TESLA);
1394     }
1395 
1396     /**
1397      * Gets initial z-coordinate of magnetometer hard iron bias to be used
1398      * to find a solution.
1399      *
1400      * @param result instance where result will be stored.
1401      */
1402     @Override
1403     public void getInitialHardIronZAsMagneticFluxDensity(final MagneticFluxDensity result) {
1404         result.setValue(initialHardIronZ);
1405         result.setUnit(MagneticFluxDensityUnit.TESLA);
1406     }
1407 
1408     /**
1409      * Sets initial z-coordinate of magnetometer hard iron bias to be used
1410      * to find a solution.
1411      *
1412      * @param initialHardIronZ initial z-coordinate of magnetometer bias.
1413      * @throws LockedException if calibrator is currently running.
1414      */
1415     @Override
1416     public void setInitialHardIronZ(final MagneticFluxDensity initialHardIronZ) throws LockedException {
1417         if (running) {
1418             throw new LockedException();
1419         }
1420         this.initialHardIronZ = convertMagneticFluxDensity(initialHardIronZ);
1421     }
1422 
1423     /**
1424      * Sets initial hard-iron bias coordinates of magnetometer used to find
1425      * a solution expressed in Teslas (T).
1426      *
1427      * @param initialHardIronX initial x-coordinate of magnetometer
1428      *                         hard-iron bias.
1429      * @param initialHardIronY initial y-coordinate of magnetometer
1430      *                         hard-iron bias.
1431      * @param initialHardIronZ initial z-coordinate of magnetometer
1432      *                         hard-iron bias.
1433      * @throws LockedException if calibrator is currently running.
1434      */
1435     @Override
1436     public void setInitialHardIron(
1437             final double initialHardIronX, final double initialHardIronY, final double initialHardIronZ)
1438             throws LockedException {
1439         if (running) {
1440             throw new LockedException();
1441         }
1442         this.initialHardIronX = initialHardIronX;
1443         this.initialHardIronY = initialHardIronY;
1444         this.initialHardIronZ = initialHardIronZ;
1445     }
1446 
1447     /**
1448      * Sets initial hard iron coordinates of magnetometer used to find a solution.
1449      *
1450      * @param initialHardIronX initial x-coordinate of magnetometer bias.
1451      * @param initialHardIronY initial y-coordinate of magnetometer bias.
1452      * @param initialHardIronZ initial z-coordinate of magnetometer bias.
1453      * @throws LockedException if calibrator is currently running.
1454      */
1455     @Override
1456     public void setInitialHardIron(
1457             final MagneticFluxDensity initialHardIronX, final MagneticFluxDensity initialHardIronY,
1458             final MagneticFluxDensity initialHardIronZ) throws LockedException {
1459         if (running) {
1460             throw new LockedException();
1461         }
1462 
1463         this.initialHardIronX = convertMagneticFluxDensity(initialHardIronX);
1464         this.initialHardIronY = convertMagneticFluxDensity(initialHardIronY);
1465         this.initialHardIronZ = convertMagneticFluxDensity(initialHardIronZ);
1466     }
1467 
1468     /**
1469      * Gets initial hard-iron used to find a solution.
1470      *
1471      * @return initial hard-iron.
1472      */
1473     @Override
1474     public MagneticFluxDensityTriad getInitialHardIronAsTriad() {
1475         return new MagneticFluxDensityTriad(MagneticFluxDensityUnit.TESLA,
1476                 initialHardIronX, initialHardIronY, initialHardIronZ);
1477     }
1478 
1479     /**
1480      * Gets initial hard-iron used to find a solution.
1481      *
1482      * @param result instance where result will be stored.
1483      */
1484     @Override
1485     public void getInitialHardIronAsTriad(final MagneticFluxDensityTriad result) {
1486         result.setValueCoordinatesAndUnit(initialHardIronX, initialHardIronY, initialHardIronZ,
1487                 MagneticFluxDensityUnit.TESLA);
1488     }
1489 
1490     /**
1491      * Sets initial hard-iron used to find a solution.
1492      *
1493      * @param initialHardIron initial hard-iron to be set.
1494      * @throws LockedException if calibrator is currently running.
1495      */
1496     @Override
1497     public void setInitialHardIron(final MagneticFluxDensityTriad initialHardIron) throws LockedException {
1498         if (running) {
1499             throw new LockedException();
1500         }
1501 
1502         initialHardIronX = convertMagneticFluxDensity(initialHardIron.getValueX(), initialHardIron.getUnit());
1503         initialHardIronY = convertMagneticFluxDensity(initialHardIron.getValueY(), initialHardIron.getUnit());
1504         initialHardIronZ = convertMagneticFluxDensity(initialHardIron.getValueZ(), initialHardIron.getUnit());
1505     }
1506 
1507     /**
1508      * Gets initial x scaling factor.
1509      *
1510      * @return initial x scaling factor.
1511      */
1512     @Override
1513     public double getInitialSx() {
1514         return initialSx;
1515     }
1516 
1517     /**
1518      * Sets initial x scaling factor.
1519      *
1520      * @param initialSx initial x scaling factor.
1521      * @throws LockedException if calibrator is currently running.
1522      */
1523     @Override
1524     public void setInitialSx(final double initialSx) throws LockedException {
1525         if (running) {
1526             throw new LockedException();
1527         }
1528         this.initialSx = initialSx;
1529     }
1530 
1531     /**
1532      * Gets initial y scaling factor.
1533      *
1534      * @return initial y scaling factor.
1535      */
1536     @Override
1537     public double getInitialSy() {
1538         return initialSy;
1539     }
1540 
1541     /**
1542      * Sets initial y scaling factor.
1543      *
1544      * @param initialSy initial y scaling factor.
1545      * @throws LockedException if calibrator is currently running.
1546      */
1547     @Override
1548     public void setInitialSy(final double initialSy) throws LockedException {
1549         if (running) {
1550             throw new LockedException();
1551         }
1552         this.initialSy = initialSy;
1553     }
1554 
1555     /**
1556      * Gets initial z scaling factor.
1557      *
1558      * @return initial z scaling factor.
1559      */
1560     @Override
1561     public double getInitialSz() {
1562         return initialSz;
1563     }
1564 
1565     /**
1566      * Sets initial z scaling factor.
1567      *
1568      * @param initialSz initial z scaling factor.
1569      * @throws LockedException if calibrator is currently running.
1570      */
1571     @Override
1572     public void setInitialSz(final double initialSz) throws LockedException {
1573         if (running) {
1574             throw new LockedException();
1575         }
1576         this.initialSz = initialSz;
1577     }
1578 
1579     /**
1580      * Gets initial x-y cross coupling error.
1581      *
1582      * @return initial x-y cross coupling error.
1583      */
1584     @Override
1585     public double getInitialMxy() {
1586         return initialMxy;
1587     }
1588 
1589     /**
1590      * Sets initial x-y cross coupling error.
1591      *
1592      * @param initialMxy initial x-y cross coupling error.
1593      * @throws LockedException if calibrator is currently running.
1594      */
1595     @Override
1596     public void setInitialMxy(final double initialMxy) throws LockedException {
1597         if (running) {
1598             throw new LockedException();
1599         }
1600         this.initialMxy = initialMxy;
1601     }
1602 
1603     /**
1604      * Gets initial x-z cross coupling error.
1605      *
1606      * @return initial x-z cross coupling error.
1607      */
1608     @Override
1609     public double getInitialMxz() {
1610         return initialMxz;
1611     }
1612 
1613     /**
1614      * Sets initial x-z cross coupling error.
1615      *
1616      * @param initialMxz initial x-z cross coupling error.
1617      * @throws LockedException if calibrator is currently running.
1618      */
1619     @Override
1620     public void setInitialMxz(final double initialMxz) throws LockedException {
1621         if (running) {
1622             throw new LockedException();
1623         }
1624         this.initialMxz = initialMxz;
1625     }
1626 
1627     /**
1628      * Gets initial y-x cross coupling error.
1629      *
1630      * @return initial y-x cross coupling error.
1631      */
1632     @Override
1633     public double getInitialMyx() {
1634         return initialMyx;
1635     }
1636 
1637     /**
1638      * Sets initial y-x cross coupling error.
1639      *
1640      * @param initialMyx initial y-x cross coupling error.
1641      * @throws LockedException if calibrator is currently running.
1642      */
1643     @Override
1644     public void setInitialMyx(final double initialMyx) throws LockedException {
1645         if (running) {
1646             throw new LockedException();
1647         }
1648         this.initialMyx = initialMyx;
1649     }
1650 
1651     /**
1652      * Gets initial y-z cross coupling error.
1653      *
1654      * @return initial y-z cross coupling error.
1655      */
1656     @Override
1657     public double getInitialMyz() {
1658         return initialMyz;
1659     }
1660 
1661     /**
1662      * Sets initial y-z cross coupling error.
1663      *
1664      * @param initialMyz initial y-z cross coupling error.
1665      * @throws LockedException if calibrator is currently running.
1666      */
1667     @Override
1668     public void setInitialMyz(final double initialMyz) throws LockedException {
1669         if (running) {
1670             throw new LockedException();
1671         }
1672         this.initialMyz = initialMyz;
1673     }
1674 
1675     /**
1676      * Gets initial z-x cross coupling error.
1677      *
1678      * @return initial z-x cross coupling error.
1679      */
1680     @Override
1681     public double getInitialMzx() {
1682         return initialMzx;
1683     }
1684 
1685     /**
1686      * Sets initial z-x cross coupling error.
1687      *
1688      * @param initialMzx initial z-x cross coupling error.
1689      * @throws LockedException if calibrator is currently running.
1690      */
1691     @Override
1692     public void setInitialMzx(final double initialMzx) throws LockedException {
1693         if (running) {
1694             throw new LockedException();
1695         }
1696         this.initialMzx = initialMzx;
1697     }
1698 
1699     /**
1700      * Gets initial z-y cross coupling error.
1701      *
1702      * @return initial z-y cross coupling error.
1703      */
1704     @Override
1705     public double getInitialMzy() {
1706         return initialMzy;
1707     }
1708 
1709     /**
1710      * Sets initial z-y cross coupling error.
1711      *
1712      * @param initialMzy initial z-y cross coupling error.
1713      * @throws LockedException if calibrator is currently running.
1714      */
1715     @Override
1716     public void setInitialMzy(final double initialMzy) throws LockedException {
1717         if (running) {
1718             throw new LockedException();
1719         }
1720         this.initialMzy = initialMzy;
1721     }
1722 
1723     /**
1724      * Sets initial scaling factors.
1725      *
1726      * @param initialSx initial x scaling factor.
1727      * @param initialSy initial y scaling factor.
1728      * @param initialSz initial z scaling factor.
1729      * @throws LockedException if calibrator is currently running.
1730      */
1731     @Override
1732     public void setInitialScalingFactors(
1733             final double initialSx, final double initialSy, final double initialSz) throws LockedException {
1734         if (running) {
1735             throw new LockedException();
1736         }
1737         this.initialSx = initialSx;
1738         this.initialSy = initialSy;
1739         this.initialSz = initialSz;
1740     }
1741 
1742     /**
1743      * Sets initial cross coupling errors.
1744      *
1745      * @param initialMxy initial x-y cross coupling error.
1746      * @param initialMxz initial x-z cross coupling error.
1747      * @param initialMyx initial y-x cross coupling error.
1748      * @param initialMyz initial y-z cross coupling error.
1749      * @param initialMzx initial z-x cross coupling error.
1750      * @param initialMzy initial z-y cross coupling error.
1751      * @throws LockedException if calibrator is currently running.
1752      */
1753     @Override
1754     public void setInitialCrossCouplingErrors(
1755             final double initialMxy, final double initialMxz, final double initialMyx,
1756             final double initialMyz, final double initialMzx, final double initialMzy) throws LockedException {
1757         if (running) {
1758             throw new LockedException();
1759         }
1760         this.initialMxy = initialMxy;
1761         this.initialMxz = initialMxz;
1762         this.initialMyx = initialMyx;
1763         this.initialMyz = initialMyz;
1764         this.initialMzx = initialMzx;
1765         this.initialMzy = initialMzy;
1766     }
1767 
1768     /**
1769      * Sets initial scaling factors and cross coupling errors.
1770      *
1771      * @param initialSx  initial x scaling factor.
1772      * @param initialSy  initial y scaling factor.
1773      * @param initialSz  initial z scaling factor.
1774      * @param initialMxy initial x-y cross coupling error.
1775      * @param initialMxz initial x-z cross coupling error.
1776      * @param initialMyx initial y-x cross coupling error.
1777      * @param initialMyz initial y-z cross coupling error.
1778      * @param initialMzx initial z-x cross coupling error.
1779      * @param initialMzy initial z-y cross coupling error.
1780      * @throws LockedException if calibrator is currently running.
1781      */
1782     @Override
1783     public void setInitialScalingFactorsAndCrossCouplingErrors(
1784             final double initialSx, final double initialSy, final double initialSz,
1785             final double initialMxy, final double initialMxz, final double initialMyx,
1786             final double initialMyz, final double initialMzx, final double initialMzy) throws LockedException {
1787         if (running) {
1788             throw new LockedException();
1789         }
1790         setInitialScalingFactors(initialSx, initialSy, initialSz);
1791         setInitialCrossCouplingErrors(initialMxy, initialMxz, initialMyx, initialMyz, initialMzx, initialMzy);
1792     }
1793 
1794     /**
1795      * Gets initial hard-iron bias to be used to find a solution as an array.
1796      * Array values are expressed in Teslas (T).
1797      *
1798      * @return array containing coordinates of initial bias.
1799      */
1800     @Override
1801     public double[] getInitialHardIron() {
1802         final var result = new double[BodyMagneticFluxDensity.COMPONENTS];
1803         getInitialHardIron(result);
1804         return result;
1805     }
1806 
1807     /**
1808      * Gets initial hard-iron  bias to be used to find a solution as an array.
1809      * Array values are expressed in Teslas (T).
1810      *
1811      * @param result instance where result data will be copied to.
1812      * @throws IllegalArgumentException if provided array does not have
1813      *                                  length 3.
1814      */
1815     @Override
1816     public void getInitialHardIron(final double[] result) {
1817         if (result.length != BodyMagneticFluxDensity.COMPONENTS) {
1818             throw new IllegalArgumentException();
1819         }
1820         result[0] = initialHardIronX;
1821         result[1] = initialHardIronY;
1822         result[2] = initialHardIronZ;
1823     }
1824 
1825     /**
1826      * Sets initial hard-iron bias to be used to find a solution as an array.
1827      * Array values are expressed in Teslas (T).
1828      *
1829      * @param initialHardIron initial hard-iron to find a solution.
1830      * @throws LockedException          if calibrator is currently running.
1831      * @throws IllegalArgumentException if provided array does not have length 3.
1832      */
1833     @Override
1834     public void setInitialHardIron(final double[] initialHardIron) throws LockedException {
1835         if (running) {
1836             throw new LockedException();
1837         }
1838 
1839         if (initialHardIron.length != BodyMagneticFluxDensity.COMPONENTS) {
1840             throw new IllegalArgumentException();
1841         }
1842         initialHardIronX = initialHardIron[0];
1843         initialHardIronY = initialHardIron[1];
1844         initialHardIronZ = initialHardIron[2];
1845     }
1846 
1847     /**
1848      * Gets initial hard-iron bias to be used to find a solution as a
1849      * column matrix.
1850      * Values are expressed in Teslas (T).
1851      *
1852      * @return initial hard-iron bias to be used to find a solution as a
1853      * column matrix.
1854      */
1855     @Override
1856     public Matrix getInitialHardIronAsMatrix() {
1857         Matrix result;
1858         try {
1859             result = new Matrix(BodyMagneticFluxDensity.COMPONENTS, 1);
1860             getInitialHardIronAsMatrix(result);
1861         } catch (final WrongSizeException ignore) {
1862             // never happens
1863             result = null;
1864         }
1865         return result;
1866     }
1867 
1868     /**
1869      * Gets initial hard-iron bias to be used to find a solution as a
1870      * column matrix.
1871      * Values are expressed in Teslas (T).
1872      *
1873      * @param result instance where result data will be copied to.
1874      * @throws IllegalArgumentException if provided matrix is not 3x1.
1875      */
1876     @Override
1877     public void getInitialHardIronAsMatrix(final Matrix result) {
1878         if (result.getRows() != BodyMagneticFluxDensity.COMPONENTS || result.getColumns() != 1) {
1879             throw new IllegalArgumentException();
1880         }
1881         result.setElementAtIndex(0, initialHardIronX);
1882         result.setElementAtIndex(1, initialHardIronY);
1883         result.setElementAtIndex(2, initialHardIronZ);
1884     }
1885 
1886     /**
1887      * Sets initial hard-iron bias to be used to find a solution as a column
1888      * matrix with values expressed in Teslas (T).
1889      *
1890      * @param initialHardIron initial hard-iron bias to find a solution.
1891      * @throws LockedException          if calibrator is currently running.
1892      * @throws IllegalArgumentException if provided matrix is not 3x1.
1893      */
1894     @Override
1895     public void setInitialHardIron(final Matrix initialHardIron) throws LockedException {
1896         if (running) {
1897             throw new LockedException();
1898         }
1899         if (initialHardIron.getRows() != BodyMagneticFluxDensity.COMPONENTS || initialHardIron.getColumns() != 1) {
1900             throw new IllegalArgumentException();
1901         }
1902 
1903         initialHardIronX = initialHardIron.getElementAtIndex(0);
1904         initialHardIronY = initialHardIron.getElementAtIndex(1);
1905         initialHardIronZ = initialHardIron.getElementAtIndex(2);
1906     }
1907 
1908     /**
1909      * Gets initial scale factors and cross coupling errors matrix.
1910      *
1911      * @return initial scale factors and cross coupling errors matrix.
1912      */
1913     @Override
1914     public Matrix getInitialMm() {
1915         Matrix result;
1916         try {
1917             result = new Matrix(BodyMagneticFluxDensity.COMPONENTS, BodyMagneticFluxDensity.COMPONENTS);
1918             getInitialMm(result);
1919         } catch (final WrongSizeException ignore) {
1920             // never happens
1921             result = null;
1922         }
1923         return result;
1924     }
1925 
1926     /**
1927      * Gets initial scale factors and cross coupling errors matrix.
1928      *
1929      * @param result instance where data will be stored.
1930      * @throws IllegalArgumentException if provided matrix is not 3x3.
1931      */
1932     @Override
1933     public void getInitialMm(final Matrix result) {
1934         if (result.getRows() != BodyKinematics.COMPONENTS || result.getColumns() != BodyKinematics.COMPONENTS) {
1935             throw new IllegalArgumentException();
1936         }
1937         result.setElementAtIndex(0, initialSx);
1938         result.setElementAtIndex(1, initialMyx);
1939         result.setElementAtIndex(2, initialMzx);
1940 
1941         result.setElementAtIndex(3, initialMxy);
1942         result.setElementAtIndex(4, initialSy);
1943         result.setElementAtIndex(5, initialMzy);
1944 
1945         result.setElementAtIndex(6, initialMxz);
1946         result.setElementAtIndex(7, initialMyz);
1947         result.setElementAtIndex(8, initialSz);
1948     }
1949 
1950     /**
1951      * Sets initial scale factors and cross coupling errors matrix.
1952      *
1953      * @param initialMm initial scale factors and cross coupling errors matrix.
1954      * @throws IllegalArgumentException if provided matrix is not 3x3.
1955      * @throws LockedException          if calibrator is currently running.
1956      */
1957     @Override
1958     public void setInitialMm(final Matrix initialMm) throws LockedException {
1959         if (running) {
1960             throw new LockedException();
1961         }
1962         if (initialMm.getRows() != BodyKinematics.COMPONENTS || initialMm.getColumns() != BodyKinematics.COMPONENTS) {
1963             throw new IllegalArgumentException();
1964         }
1965 
1966         initialSx = initialMm.getElementAtIndex(0);
1967         initialMyx = initialMm.getElementAtIndex(1);
1968         initialMzx = initialMm.getElementAtIndex(2);
1969 
1970         initialMxy = initialMm.getElementAtIndex(3);
1971         initialSy = initialMm.getElementAtIndex(4);
1972         initialMzy = initialMm.getElementAtIndex(5);
1973 
1974         initialMxz = initialMm.getElementAtIndex(6);
1975         initialMyz = initialMm.getElementAtIndex(7);
1976         initialSz = initialMm.getElementAtIndex(8);
1977     }
1978 
1979     /**
1980      * Gets collection of body magnetic flux density measurements taken
1981      * at a given position with different unknown orientations and containing the
1982      * standard deviation of magnetometer measurements.
1983      *
1984      * @return collection of body magnetic flux density measurements at
1985      * a known position and timestamp with unknown orientations.
1986      */
1987     @Override
1988     public Collection<StandardDeviationBodyMagneticFluxDensity> getMeasurements() {
1989         return measurements;
1990     }
1991 
1992     /**
1993      * Sets collection of body magnetic flux density measurements taken
1994      * at a given position with different unknown orientations and containing the
1995      * standard deviation of magnetometer measurements.
1996      *
1997      * @param measurements collection of body magnetic flux density
1998      *                     measurements at a known position and timestamp
1999      *                     with unknown orientations.
2000      * @throws LockedException if calibrator is currently running.
2001      */
2002     @Override
2003     public void setMeasurements(
2004             final Collection<StandardDeviationBodyMagneticFluxDensity> measurements) throws LockedException {
2005         if (running) {
2006             throw new LockedException();
2007         }
2008         this.measurements = measurements;
2009     }
2010 
2011     /**
2012      * Indicates the type of measurement used by this calibrator.
2013      *
2014      * @return type of measurement used by this calibrator.
2015      */
2016     @Override
2017     public MagnetometerCalibratorMeasurementType getMeasurementType() {
2018         return MagnetometerCalibratorMeasurementType.STANDARD_DEVIATION_BODY_MAGNETIC_FLUX_DENSITY;
2019     }
2020 
2021     /**
2022      * Indicates whether this calibrator requires ordered measurements in a
2023      * list or not.
2024      *
2025      * @return true if measurements must be ordered, false otherwise.
2026      */
2027     @Override
2028     public boolean isOrderedMeasurementsRequired() {
2029         return false;
2030     }
2031 
2032     /**
2033      * Indicates whether this calibrator requires quality scores for each
2034      * measurement or not.
2035      *
2036      * @return true if quality scores are required, false otherwise.
2037      */
2038     @Override
2039     public boolean isQualityScoresRequired() {
2040         return false;
2041     }
2042 
2043     /**
2044      * Indicates whether z-axis is assumed to be common for accelerometer,
2045      * gyroscope and magnetometer.
2046      * When enabled, this eliminates 3 variables from Mm (soft-iron) matrix.
2047      *
2048      * @return true if z-axis is assumed to be common for accelerometer,
2049      * gyroscope and magnetometer, false otherwise.
2050      */
2051     @Override
2052     public boolean isCommonAxisUsed() {
2053         return commonAxisUsed;
2054     }
2055 
2056     /**
2057      * Specifies whether z-axis is assumed to be common for accelerometer and
2058      * gyroscope.
2059      * When enabled, this eliminates 3 variables from Mm matrix.
2060      *
2061      * @param commonAxisUsed true if z-axis is assumed to be common for
2062      *                       accelerometer, gyroscope and magnetometer, false
2063      *                       otherwise.
2064      * @throws LockedException if estimator is currently running.
2065      */
2066     @Override
2067     public void setCommonAxisUsed(final boolean commonAxisUsed) throws LockedException {
2068         if (running) {
2069             throw new LockedException();
2070         }
2071 
2072         this.commonAxisUsed = commonAxisUsed;
2073     }
2074 
2075     /**
2076      * Gets listener to handle events raised by this calibrator.
2077      *
2078      * @return listener to handle events raised by this calibrator.
2079      */
2080     public L getListener() {
2081         return listener;
2082     }
2083 
2084     /**
2085      * Sets listener to handle events raised by this calibrator.
2086      *
2087      * @param listener listener to handle events raised by this calibrator.
2088      * @throws LockedException if calibrator is currently running.
2089      */
2090     public void setListener(final L listener) throws LockedException {
2091         if (running) {
2092             throw new LockedException();
2093         }
2094 
2095         this.listener = listener;
2096     }
2097 
2098     /**
2099      * Gets minimum number of required measurements.
2100      *
2101      * @return minimum number of required measurements.
2102      */
2103     @Override
2104     public int getMinimumRequiredMeasurements() {
2105         return commonAxisUsed ? MINIMUM_MEASUREMENTS_COMMON_Z_AXIS : MINIMUM_MEASUREMENTS_GENERAL;
2106     }
2107 
2108     /**
2109      * Indicates whether calibrator is ready to start.
2110      *
2111      * @return true if calibrator is ready, false otherwise.
2112      */
2113     @Override
2114     public boolean isReady() {
2115         return measurements != null && measurements.size() >= getMinimumRequiredMeasurements();
2116     }
2117 
2118     /**
2119      * Indicates whether calibrator is currently running or not.
2120      *
2121      * @return true if calibrator is running, false otherwise.
2122      */
2123     @Override
2124     public boolean isRunning() {
2125         return running;
2126     }
2127 
2128     /**
2129      * Estimates magnetometer calibration parameters containing scale factors
2130      * and cross-coupling errors.
2131      *
2132      * @throws LockedException      if calibrator is currently running.
2133      * @throws NotReadyException    if calibrator is not ready.
2134      * @throws CalibrationException if calibration fails for numerical reasons.
2135      */
2136     @Override
2137     public void calibrate() throws LockedException, NotReadyException, CalibrationException {
2138         if (running) {
2139             throw new LockedException();
2140         }
2141 
2142         if (!isReady()) {
2143             throw new NotReadyException();
2144         }
2145 
2146         try {
2147             running = true;
2148 
2149             onBeforeCalibrate();
2150 
2151             if (listener != null) {
2152                 //noinspection unchecked
2153                 listener.onCalibrateStart((C) this);
2154             }
2155 
2156             if (commonAxisUsed) {
2157                 calibrateCommonAxis();
2158             } else {
2159                 calibrateGeneral();
2160             }
2161 
2162             if (listener != null) {
2163                 //noinspection unchecked
2164                 listener.onCalibrateEnd((C) this);
2165             }
2166 
2167         } catch (final AlgebraException | FittingException | com.irurueta.numerical.NotReadyException e) {
2168             throw new CalibrationException(e);
2169         } finally {
2170             running = false;
2171         }
2172     }
2173 
2174     /**
2175      * Gets array containing x,y,z components of estimated magnetometer
2176      * hard-iron biases expressed in Teslas (T).
2177      *
2178      * @return array containing x,y,z components of estimated magnetometer
2179      * hard-iron biases.
2180      */
2181     @Override
2182     public double[] getEstimatedHardIron() {
2183         return estimatedHardIron;
2184     }
2185 
2186     /**
2187      * Gets array containing x,y,z components of estimated magnetometer
2188      * hard-iron biases expressed in Teslas (T).
2189      *
2190      * @param result instance where estimated magnetometer biases will be
2191      *               stored.
2192      * @return true if result instance was updated, false otherwise (when
2193      * estimation is not yet available).
2194      */
2195     @Override
2196     public boolean getEstimatedHardIron(final double[] result) {
2197         if (estimatedHardIron != null) {
2198             System.arraycopy(estimatedHardIron, 0, result, 0, estimatedHardIron.length);
2199             return true;
2200         } else {
2201             return false;
2202         }
2203     }
2204 
2205     /**
2206      * Gets column matrix containing x,y,z components of estimated
2207      * magnetometer hard-iron biases expressed in Teslas (T).
2208      *
2209      * @return column matrix containing x,y,z components of estimated
2210      * magnetometer hard-iron biases.
2211      */
2212     @Override
2213     public Matrix getEstimatedHardIronAsMatrix() {
2214         return estimatedHardIron != null ? Matrix.newFromArray(estimatedHardIron) : null;
2215     }
2216 
2217     /**
2218      * Gets column matrix containing x,y,z components of estimated
2219      * magnetometer hard-iron biases expressed in Teslas (T).
2220      *
2221      * @param result instance where result data will be stored.
2222      * @return true if result was updated, false otherwise.
2223      * @throws WrongSizeException if provided result instance has invalid size.
2224      */
2225     @Override
2226     public boolean getEstimatedHardIronAsMatrix(final Matrix result) throws WrongSizeException {
2227         if (estimatedHardIron != null) {
2228             result.fromArray(estimatedHardIron);
2229             return true;
2230         } else {
2231             return false;
2232         }
2233     }
2234 
2235     /**
2236      * Gets x coordinate of estimated magnetometer bias expressed in
2237      * Teslas (T).
2238      *
2239      * @return x coordinate of estimated magnetometer bias or null if not
2240      * available.
2241      */
2242     @Override
2243     public Double getEstimatedHardIronX() {
2244         return estimatedHardIron != null ? estimatedHardIron[0] : null;
2245     }
2246 
2247     /**
2248      * Gets y coordinate of estimated magnetometer bias expressed in
2249      * Teslas (T).
2250      *
2251      * @return y coordinate of estimated magnetometer bias or null if not
2252      * available.
2253      */
2254     @Override
2255     public Double getEstimatedHardIronY() {
2256         return estimatedHardIron != null ? estimatedHardIron[1] : null;
2257     }
2258 
2259     /**
2260      * Gets z coordinate of estimated magnetometer bias expressed in
2261      * Teslas (T).
2262      *
2263      * @return z coordinate of estimated magnetometer bias or null if not
2264      * available.
2265      */
2266     @Override
2267     public Double getEstimatedHardIronZ() {
2268         return estimatedHardIron != null ? estimatedHardIron[2] : null;
2269     }
2270 
2271     /**
2272      * Gets x coordinate of estimated magnetometer bias.
2273      *
2274      * @return x coordinate of estimated magnetometer bias.
2275      */
2276     @Override
2277     public MagneticFluxDensity getEstimatedHardIronXAsMagneticFluxDensity() {
2278         return estimatedHardIron != null
2279                 ? new MagneticFluxDensity(estimatedHardIron[0], MagneticFluxDensityUnit.TESLA) : null;
2280     }
2281 
2282     /**
2283      * Gets x coordinate of estimated magnetometer bias.
2284      *
2285      * @param result instance where result will be stored.
2286      * @return true if estimated magnetometer bias is available, false otherwise.
2287      */
2288     @Override
2289     public boolean getEstimatedHardIronXAsMagneticFluxDensity(final MagneticFluxDensity result) {
2290         if (estimatedHardIron != null) {
2291             result.setValue(estimatedHardIron[0]);
2292             result.setUnit(MagneticFluxDensityUnit.TESLA);
2293             return true;
2294         } else {
2295             return false;
2296         }
2297     }
2298 
2299     /**
2300      * Gets y coordinate of estimated magnetometer bias.
2301      *
2302      * @return y coordinate of estimated magnetometer bias.
2303      */
2304     @Override
2305     public MagneticFluxDensity getEstimatedHardIronYAsMagneticFluxDensity() {
2306         return estimatedHardIron != null
2307                 ? new MagneticFluxDensity(estimatedHardIron[1], MagneticFluxDensityUnit.TESLA) : null;
2308     }
2309 
2310     /**
2311      * Gets y coordinate of estimated magnetometer bias.
2312      *
2313      * @param result instance where result will be stored.
2314      * @return true if estimated magnetometer bias is available, false otherwise.
2315      */
2316     @Override
2317     public boolean getEstimatedHardIronYAsMagneticFluxDensity(final MagneticFluxDensity result) {
2318         if (estimatedHardIron != null) {
2319             result.setValue(estimatedHardIron[1]);
2320             result.setUnit(MagneticFluxDensityUnit.TESLA);
2321             return true;
2322         } else {
2323             return false;
2324         }
2325     }
2326 
2327     /**
2328      * Gets z coordinate of estimated magnetometer bias.
2329      *
2330      * @return z coordinate of estimated magnetometer bias.
2331      */
2332     @Override
2333     public MagneticFluxDensity getEstimatedHardIronZAsMagneticFluxDensity() {
2334         return estimatedHardIron != null
2335                 ? new MagneticFluxDensity(estimatedHardIron[2], MagneticFluxDensityUnit.TESLA) : null;
2336     }
2337 
2338     /**
2339      * Gets z coordinate of estimated magnetometer bias.
2340      *
2341      * @param result instance where result will be stored.
2342      * @return true if estimated magnetometer bias is available, false otherwise.
2343      */
2344     @Override
2345     public boolean getEstimatedHardIronZAsMagneticFluxDensity(final MagneticFluxDensity result) {
2346         if (estimatedHardIron != null) {
2347             result.setValue(estimatedHardIron[2]);
2348             result.setUnit(MagneticFluxDensityUnit.TESLA);
2349             return true;
2350         } else {
2351             return false;
2352         }
2353     }
2354 
2355     /**
2356      * Gets estimated magnetometer bias.
2357      *
2358      * @return estimated magnetometer bias or null if not available.
2359      */
2360     @Override
2361     public MagneticFluxDensityTriad getEstimatedHardIronAsTriad() {
2362         return estimatedHardIron != null
2363                 ? new MagneticFluxDensityTriad(MagneticFluxDensityUnit.TESLA,
2364                 estimatedHardIron[0], estimatedHardIron[1], estimatedHardIron[2])
2365                 : null;
2366     }
2367 
2368     /**
2369      * Gets estimated magnetometer bias.
2370      *
2371      * @param result instance where result will be stored.
2372      * @return true if estimated magnetometer bias is available and result was
2373      * modified, false otherwise.
2374      */
2375     @Override
2376     public boolean getEstimatedHardIronAsTriad(final MagneticFluxDensityTriad result) {
2377         if (estimatedHardIron != null) {
2378             result.setValueCoordinatesAndUnit(estimatedHardIron[0], estimatedHardIron[1], estimatedHardIron[2],
2379                     MagneticFluxDensityUnit.TESLA);
2380             return true;
2381         } else {
2382             return false;
2383         }
2384     }
2385 
2386     /**
2387      * Gets estimated magnetometer soft-iron matrix containing scale factors
2388      * and cross coupling errors.
2389      * This is the product of matrix Tm containing cross coupling errors and Km
2390      * containing scaling factors.
2391      * So tat:
2392      * <pre>
2393      *     Mm = [sx    mxy  mxz] = Tm*Km
2394      *          [myx   sy   myz]
2395      *          [mzx   mzy  sz ]
2396      * </pre>
2397      * Where:
2398      * <pre>
2399      *     Km = [sx 0   0 ]
2400      *          [0  sy  0 ]
2401      *          [0  0   sz]
2402      * </pre>
2403      * and
2404      * <pre>
2405      *     Tm = [1          -alphaXy    alphaXz ]
2406      *          [alphaYx    1           -alphaYz]
2407      *          [-alphaZx   alphaZy     1       ]
2408      * </pre>
2409      * Hence:
2410      * <pre>
2411      *     Mm = [sx    mxy  mxz] = Tm*Km =  [sx             -sy * alphaXy   sz * alphaXz ]
2412      *          [myx   sy   myz]            [sx * alphaYx   sy              -sz * alphaYz]
2413      *          [mzx   mzy  sz ]            [-sx * alphaZx  sy * alphaZy    sz           ]
2414      * </pre>
2415      * This instance allows any 3x3 matrix however, typically alphaYx, alphaZx and alphaZy
2416      * are considered to be zero if the accelerometer z-axis is assumed to be the same
2417      * as the body z-axis. When this is assumed, myx = mzx = mzy = 0 and the Mm matrix
2418      * becomes upper diagonal:
2419      * <pre>
2420      *     Mm = [sx    mxy  mxz]
2421      *          [0     sy   myz]
2422      *          [0     0    sz ]
2423      * </pre>
2424      * Values of this matrix are unit-less.
2425      *
2426      * @return estimated magnetometer soft-iron scale factors and cross coupling errors,
2427      * or null if not available.
2428      */
2429     @Override
2430     public Matrix getEstimatedMm() {
2431         return estimatedMm;
2432     }
2433 
2434     /**
2435      * Gets estimated x-axis scale factor.
2436      *
2437      * @return estimated x-axis scale factor or null if not available.
2438      */
2439     @Override
2440     public Double getEstimatedSx() {
2441         return estimatedMm != null ? estimatedMm.getElementAt(0, 0) : null;
2442     }
2443 
2444     /**
2445      * Gets estimated y-axis scale factor.
2446      *
2447      * @return estimated y-axis scale factor or null if not available.
2448      */
2449     @Override
2450     public Double getEstimatedSy() {
2451         return estimatedMm != null ? estimatedMm.getElementAt(1, 1) : null;
2452     }
2453 
2454     /**
2455      * Gets estimated z-axis scale factor.
2456      *
2457      * @return estimated z-axis scale factor or null if not available.
2458      */
2459     @Override
2460     public Double getEstimatedSz() {
2461         return estimatedMm != null ? estimatedMm.getElementAt(2, 2) : null;
2462     }
2463 
2464     /**
2465      * Gets estimated x-y cross-coupling error.
2466      *
2467      * @return estimated x-y cross-coupling error or null if not available.
2468      */
2469     @Override
2470     public Double getEstimatedMxy() {
2471         return estimatedMm != null ? estimatedMm.getElementAt(0, 1) : null;
2472     }
2473 
2474     /**
2475      * Gets estimated x-z cross-coupling error.
2476      *
2477      * @return estimated x-z cross-coupling error or null if not available.
2478      */
2479     @Override
2480     public Double getEstimatedMxz() {
2481         return estimatedMm != null ? estimatedMm.getElementAt(0, 2) : null;
2482     }
2483 
2484     /**
2485      * Gets estimated y-x cross-coupling error.
2486      *
2487      * @return estimated y-x cross-coupling error or null if not available.
2488      */
2489     @Override
2490     public Double getEstimatedMyx() {
2491         return estimatedMm != null ? estimatedMm.getElementAt(1, 0) : null;
2492     }
2493 
2494     /**
2495      * Gets estimated y-z cross-coupling error.
2496      *
2497      * @return estimated y-z cross-coupling error or null if not available.
2498      */
2499     @Override
2500     public Double getEstimatedMyz() {
2501         return estimatedMm != null ? estimatedMm.getElementAt(1, 2) : null;
2502     }
2503 
2504     /**
2505      * Gets estimated z-x cross-coupling error.
2506      *
2507      * @return estimated z-x cross-coupling error or null if not available.
2508      */
2509     @Override
2510     public Double getEstimatedMzx() {
2511         return estimatedMm != null ? estimatedMm.getElementAt(2, 0) : null;
2512     }
2513 
2514     /**
2515      * Gets estimated z-y cross-coupling error.
2516      *
2517      * @return estimated z-y cross-coupling error or null if not available.
2518      */
2519     @Override
2520     public Double getEstimatedMzy() {
2521         return estimatedMm != null ? estimatedMm.getElementAt(2, 1) : null;
2522     }
2523 
2524     /**
2525      * Gets estimated covariance matrix for estimated calibration parameters.
2526      * Diagonal elements of the matrix contains variance for the following
2527      * parameters (following indicated order): bx, by, bz, sx, sy, sz,
2528      * mxy, mxz, myx, myz, mzx, mzy.
2529      *
2530      * @return estimated covariance matrix for estimated calibration parameters.
2531      */
2532     @Override
2533     public Matrix getEstimatedCovariance() {
2534         return estimatedCovariance;
2535     }
2536 
2537     /**
2538      * Gets estimated chi square value.
2539      *
2540      * @return estimated chi square value.
2541      */
2542     @Override
2543     public double getEstimatedChiSq() {
2544         return estimatedChiSq;
2545     }
2546 
2547     /**
2548      * Gets estimated chi square degrees of freedom. Degrees of freedom is equal to the number of sampled data minus the
2549      * number of estimated parameters.
2550      *
2551      * @return estimated degrees of freedom of chi square value
2552      */
2553     @Override
2554     public int getEstimatedChiSqDegreesOfFreedom() {
2555         return estimatedChiSqDegreesOfFreedom;
2556     }
2557 
2558     /**
2559      * Gets estimated reduced chi square value. This is equal to estimated chi square value divided by its degrees of
2560      * freedom. Ideally this value should be close to 1.0, indicating that fit is optimal.
2561      * A value larger than 1.0 indicates that fit is not good or noise has been underestimated, and a value smaller than
2562      * 1.0 indicates that there is overfitting or noise has been overestimated.
2563      *
2564      * @return estimated reduced chi square value
2565      */
2566     @Override
2567     public double getEstimatedReducedChiSq() {
2568         return estimatedReducedChiSq;
2569     }
2570 
2571     /**
2572      * Gets estimated mean square error respect to provided measurements.
2573      *
2574      * @return estimated mean square error respect to provided measurements.
2575      */
2576     @Override
2577     public double getEstimatedMse() {
2578         return estimatedMse;
2579     }
2580 
2581     /**
2582      * Gets estimated probability of finding a smaller chi square value expressed as a value between 0.0 and 1.0. The
2583      * smaller the found chi square value is, the better the fit of the estimated parameters to the actual parameter.
2584      * Thus, the smaller the chance of finding a smaller chi square value, then the better the estimated fit is.
2585      *
2586      * @return estimated probability of finding a smaller chi square value.
2587      */
2588     @Override
2589     public double getEstimatedP() {
2590         return estimatedP;
2591     }
2592 
2593     /**
2594      * Gets estimated measure of quality of estimated fit as a value between 0.0 and 1.0. The larger the quality value
2595      * is, the better the fit that has been estimated.
2596      *
2597      * @return estimated measure of quality of estimated fit.
2598      */
2599     @Override
2600     public double getEstimatedQ() {
2601         return estimatedQ;
2602     }
2603 
2604     /**
2605      * Gets variance of estimated x coordinate of magnetometer bias expressed in
2606      * squared Teslas (T^2).
2607      *
2608      * @return variance of estimated x coordinate of magnetometer bias or null if
2609      * not available.
2610      */
2611     public Double getEstimatedHardIronXVariance() {
2612         return estimatedCovariance != null ? estimatedCovariance.getElementAt(0, 0) : null;
2613     }
2614 
2615     /**
2616      * Gets standard deviation of estimated x coordinate of magnetometer bias
2617      * expressed in Teslas (T).
2618      *
2619      * @return standard deviation of estimated x coordinate of magnetometer bias
2620      * or null if not available.
2621      */
2622     public Double getEstimatedHardIronXStandardDeviation() {
2623         final var variance = getEstimatedHardIronXVariance();
2624         return variance != null ? Math.sqrt(variance) : null;
2625     }
2626 
2627     /**
2628      * Gets standard deviation of estimated x coordinate of magnetometer bias.
2629      *
2630      * @return standard deviation of estimated x coordinate of magnetometer bias
2631      * or null if not available.
2632      */
2633     public MagneticFluxDensity getEstimatedHardIronXStandardDeviationAsMagneticFluxDensity() {
2634         return estimatedCovariance != null
2635                 ? new MagneticFluxDensity(getEstimatedHardIronXStandardDeviation(), MagneticFluxDensityUnit.TESLA)
2636                 : null;
2637     }
2638 
2639     /**
2640      * Gets standard deviation of estimated x coordinate of magnetometer bias.
2641      *
2642      * @param result instance where result will be stored.
2643      * @return true if standard deviation of estimated x coordinate of
2644      * magnetometer bias is available, false otherwise.
2645      */
2646     public boolean getEstimatedHardIronXStandardDeviationAsMagneticFluxDensity(
2647             final MagneticFluxDensity result) {
2648         if (estimatedCovariance != null) {
2649             result.setValue(getEstimatedHardIronXStandardDeviation());
2650             result.setUnit(MagneticFluxDensityUnit.TESLA);
2651             return true;
2652         } else {
2653             return false;
2654         }
2655     }
2656 
2657     /**
2658      * Gets variance of estimated y coordinate of magnetometer bias expressed in
2659      * squared Teslas (T^2).
2660      *
2661      * @return variance of estimated y coordinate of magnetometer bias or null if
2662      * not available.
2663      */
2664     public Double getEstimatedHardIronYVariance() {
2665         return estimatedCovariance != null ? estimatedCovariance.getElementAt(1, 1) : null;
2666     }
2667 
2668     /**
2669      * Gets standard deviation of estimated y coordinate of magnetometer bias
2670      * expressed in Teslas (T).
2671      *
2672      * @return standard deviation of estimated y coordinate of magnetometer bias
2673      * or null if not available.
2674      */
2675     public Double getEstimatedHardIronYStandardDeviation() {
2676         final var variance = getEstimatedHardIronYVariance();
2677         return variance != null ? Math.sqrt(variance) : null;
2678     }
2679 
2680     /**
2681      * Gets standard deviation of estimated y coordinate of magnetometer bias.
2682      *
2683      * @return standard deviation of estimated y coordinate of magnetometer bias
2684      * or null if not available.
2685      */
2686     public MagneticFluxDensity getEstimatedHardIronYStandardDeviationAsMagneticFluxDensity() {
2687         return estimatedCovariance != null
2688                 ? new MagneticFluxDensity(getEstimatedHardIronYStandardDeviation(), MagneticFluxDensityUnit.TESLA)
2689                 : null;
2690     }
2691 
2692     /**
2693      * Gets standard deviation of estimated y coordinate of magnetometer bias.
2694      *
2695      * @param result instance where result will be stored.
2696      * @return true if standard deviation of estimated y coordinate of
2697      * magnetometer bias is available, false otherwise.
2698      */
2699     public boolean getEstimatedHardIronYStandardDeviationAsMagneticFluxDensity(
2700             final MagneticFluxDensity result) {
2701         if (estimatedCovariance != null) {
2702             result.setValue(getEstimatedHardIronYStandardDeviation());
2703             result.setUnit(MagneticFluxDensityUnit.TESLA);
2704             return true;
2705         } else {
2706             return false;
2707         }
2708     }
2709 
2710     /**
2711      * Gets variance of estimated z coordinate of magnetometer bias expressed in
2712      * squared Teslas (T^2).
2713      *
2714      * @return variance of estimated z coordinate of magnetometer bias or null if
2715      * not available.
2716      */
2717     public Double getEstimatedHardIronZVariance() {
2718         return estimatedCovariance != null ? estimatedCovariance.getElementAt(2, 2) : null;
2719     }
2720 
2721     /**
2722      * Gets standard deviation of estimated z coordinate of magnetometer bias
2723      * expressed in Teslas (T).
2724      *
2725      * @return standard deviation of estimated z coordinate of magnetometer bias
2726      * or null if not available.
2727      */
2728     public Double getEstimatedHardIronZStandardDeviation() {
2729         final var variance = getEstimatedHardIronZVariance();
2730         return variance != null ? Math.sqrt(variance) : null;
2731     }
2732 
2733     /**
2734      * Gets standard deviation of estimated z coordinate of magnetometer bias.
2735      *
2736      * @return standard deviation of estimated z coordinate of magnetometer bias
2737      * or null if not available.
2738      */
2739     public MagneticFluxDensity getEstimatedHardIronZStandardDeviationAsMagneticFluxDensity() {
2740         return estimatedCovariance != null
2741                 ? new MagneticFluxDensity(getEstimatedHardIronZStandardDeviation(), MagneticFluxDensityUnit.TESLA)
2742                 : null;
2743     }
2744 
2745     /**
2746      * Gets standard deviation of estimated z coordinate of magnetometer bias.
2747      *
2748      * @param result instance where result will be stored.
2749      * @return true if standard deviation of estimated z coordinate of
2750      * magnetometer bias is available, false otherwise.
2751      */
2752     public boolean getEstimatedHardIronZStandardDeviationAsMagneticFluxDensity(final MagneticFluxDensity result) {
2753         if (estimatedCovariance != null) {
2754             result.setValue(getEstimatedHardIronZStandardDeviation());
2755             result.setUnit(MagneticFluxDensityUnit.TESLA);
2756             return true;
2757         } else {
2758             return false;
2759         }
2760     }
2761 
2762     /**
2763      * Gets standard deviation of estimated magnetometer bias coordinates.
2764      *
2765      * @return standard deviation of estimated magnetometer bias coordinates.
2766      */
2767     public MagneticFluxDensityTriad getEstimatedHardIronStandardDeviation() {
2768         return estimatedCovariance != null
2769                 ? new MagneticFluxDensityTriad(MagneticFluxDensityUnit.TESLA,
2770                 getEstimatedHardIronXStandardDeviation(),
2771                 getEstimatedHardIronYStandardDeviation(),
2772                 getEstimatedHardIronZStandardDeviation())
2773                 : null;
2774     }
2775 
2776     /**
2777      * Gets standard deviation of estimated magnetometer bias coordinates.
2778      *
2779      * @param result instance where result will be stored.
2780      * @return true if standard deviation of magnetometer bias was available,
2781      * false otherwise.
2782      */
2783     public boolean getEstimatedHardIronStandardDeviation(final MagneticFluxDensityTriad result) {
2784         if (estimatedCovariance != null) {
2785             result.setValueCoordinatesAndUnit(
2786                     getEstimatedHardIronXStandardDeviation(),
2787                     getEstimatedHardIronYStandardDeviation(),
2788                     getEstimatedHardIronZStandardDeviation(),
2789                     MagneticFluxDensityUnit.TESLA);
2790             return true;
2791         } else {
2792             return false;
2793         }
2794     }
2795 
2796     /**
2797      * Gets average of estimated standard deviation of magnetometer bias coordinates
2798      * expressed in Teslas (T).
2799      *
2800      * @return average of estimated standard deviation of magnetometer bias coordinates,
2801      * or null if not available.
2802      */
2803     public Double getEstimatedHardIronStandardDeviationAverage() {
2804         return estimatedCovariance != null
2805                 ? (getEstimatedHardIronXStandardDeviation() + getEstimatedHardIronYStandardDeviation()
2806                 + getEstimatedHardIronZStandardDeviation()) / 3.0
2807                 : null;
2808     }
2809 
2810     /**
2811      * Gets average of estimated standard deviation of magnetometer bias coordinates.
2812      *
2813      * @return average of estimated standard deviation of magnetometer bias coordinates,
2814      * or null if not available.
2815      */
2816     public MagneticFluxDensity getEstimatedHardIronStandardDeviationAverageAsMagneticFluxDensity() {
2817         return estimatedCovariance != null
2818                 ? new MagneticFluxDensity(getEstimatedHardIronStandardDeviationAverage(), MagneticFluxDensityUnit.TESLA)
2819                 : null;
2820     }
2821 
2822     /**
2823      * Gets average of estimated standard deviation of magnetometer bias coordinates.
2824      *
2825      * @param result instance where result will be stored.
2826      * @return true if average of estimated standard deviation of magnetometer bias is available,
2827      * false otherwise.
2828      */
2829     public boolean getEstimatedHardIronStandardDeviationAverageAsMagneticFluxDensity(final MagneticFluxDensity result) {
2830         if (estimatedCovariance != null) {
2831             result.setValue(getEstimatedHardIronStandardDeviationAverage());
2832             result.setUnit(MagneticFluxDensityUnit.TESLA);
2833             return true;
2834         } else {
2835             return false;
2836         }
2837     }
2838 
2839     /**
2840      * Gets norm of estimated standard deviation of magnetometer bias expressed in
2841      * Teslas (T).
2842      *
2843      * @return norm of estimated standard deviation of magnetometer bias or null
2844      * if not available.
2845      */
2846     public Double getEstimatedHardIronStandardDeviationNorm() {
2847         return estimatedCovariance != null
2848                 ? Math.sqrt(getEstimatedHardIronXVariance() + getEstimatedHardIronYVariance()
2849                 + getEstimatedHardIronZVariance()) : null;
2850     }
2851 
2852     /**
2853      * Gets norm of estimated standard deviation of magnetometer bias.
2854      *
2855      * @return norm of estimated standard deviation of magnetometer bias or null
2856      * if not available.
2857      */
2858     public MagneticFluxDensity getEstimatedHardIronStandardDeviationNormAsMagneticFluxDensity() {
2859         return estimatedCovariance != null
2860                 ? new MagneticFluxDensity(getEstimatedHardIronStandardDeviationNorm(), MagneticFluxDensityUnit.TESLA)
2861                 : null;
2862     }
2863 
2864     /**
2865      * Gets norm of estimated standard deviation of magnetometer bias.
2866      *
2867      * @param result instance where result will be stored.
2868      * @return true if norm of estimated standard deviation of magnetometer bias
2869      * is available, false otherwise.
2870      */
2871     public boolean getEstimatedHardIronStandardDeviationNormAsMagneticFluxDensity(
2872             final MagneticFluxDensity result) {
2873         if (estimatedCovariance != null) {
2874             result.setValue(getEstimatedHardIronStandardDeviationNorm());
2875             result.setUnit(MagneticFluxDensityUnit.TESLA);
2876             return true;
2877         } else {
2878             return false;
2879         }
2880     }
2881 
2882     /**
2883      * Called before calibration occurs.
2884      * This can be overridden by subclasses.
2885      *
2886      * @throws CalibrationException if anything fails.
2887      */
2888     protected void onBeforeCalibrate() throws CalibrationException {
2889     }
2890 
2891     /**
2892      * Internally sets ground truth magnetic flux density norm to be expected at location where
2893      * measurements have been made, expressed in Teslas (T).
2894      *
2895      * @param groundTruthMagneticFluxDensityNorm ground truth magnetic flux density norm or null if undefined.
2896      * @throws IllegalArgumentException if provided value is negative.
2897      */
2898     protected void internalSetGroundTruthMagneticFluxDensityNorm(final Double groundTruthMagneticFluxDensityNorm) {
2899         if (groundTruthMagneticFluxDensityNorm != null && groundTruthMagneticFluxDensityNorm < 0.0) {
2900             throw new IllegalArgumentException();
2901         }
2902         this.groundTruthMagneticFluxDensityNorm = groundTruthMagneticFluxDensityNorm;
2903     }
2904 
2905     /**
2906      * Sets input data into Levenberg-Marquardt fitter.
2907      *
2908      * @throws WrongSizeException never happens.
2909      */
2910     protected void setInputData() throws WrongSizeException {
2911         final var gtb = groundTruthMagneticFluxDensityNorm;
2912         final var gtb2 = gtb * gtb;
2913 
2914         final var numMeasurements = measurements.size();
2915         final var x = new Matrix(numMeasurements, BodyMagneticFluxDensity.COMPONENTS);
2916         final var y = new double[numMeasurements];
2917         final var specificForceStandardDeviations = new double[numMeasurements];
2918         var i = 0;
2919         for (final var measurement : measurements) {
2920             final var measuredMagneticFluxDensity = measurement.getMagneticFluxDensity();
2921 
2922             final var bmeasuredX = measuredMagneticFluxDensity.getBx();
2923             final var bmeasuredY = measuredMagneticFluxDensity.getBy();
2924             final var bmeasuredZ = measuredMagneticFluxDensity.getBz();
2925 
2926             x.setElementAt(i, 0, bmeasuredX);
2927             x.setElementAt(i, 1, bmeasuredY);
2928             x.setElementAt(i, 2, bmeasuredZ);
2929 
2930             y[i] = gtb2;
2931 
2932             specificForceStandardDeviations[i] = measurement.getMagneticFluxDensityStandardDeviation();
2933 
2934             i++;
2935         }
2936 
2937         fitter.setInputData(x, y, specificForceStandardDeviations);
2938     }
2939 
2940     /**
2941      * Internal method to perform general calibration.
2942      *
2943      * @throws FittingException                         if Levenberg-Marquardt fails for numerical reasons.
2944      * @throws AlgebraException                         if there are numerical instabilities that prevent
2945      *                                                  matrix inversion.
2946      * @throws com.irurueta.numerical.NotReadyException never happens.
2947      */
2948     private void calibrateGeneral() throws AlgebraException, FittingException,
2949             com.irurueta.numerical.NotReadyException {
2950         // The magnetometer model is:
2951         // bmeas = ba + (I + Mm) * btrue + w
2952 
2953         // Ideally a least squares solution tries to minimize noise component, so:
2954         // bmeas = ba + (I + Mm) * btrue
2955 
2956         // For convergence purposes of the Levenberg-Marquardt algorithm, the
2957         // magnetometer model can be better expressed as:
2958         // bmeas = T*K*(btrue + b)
2959         // bmeas = M*(btrue + b)
2960         // bmeas = M*btrue + M*b
2961 
2962         // where:
2963         // M = I + Mm
2964         // bm = M*b = (I + Mm)*b --> b = M^-1*bm
2965 
2966         // We know that the norm of the true body magnetic flux density
2967         // is equal to the amount of Earth magnetic flux density at provided
2968         // position and timestamp
2969         // ||btrue|| = ||bEarth|| --> from 30 µT to 60 µT
2970 
2971         // Hence:
2972         // bmeas - M*b = M*btrue
2973 
2974         // M^-1 * (bmeas - M*b) = btrue
2975 
2976         // ||bEarth||^2 = ||btrue||^2 = (M^-1 * (bmeas - M*b))^T * (M^-1 * (bmeas - M*b))
2977         // ||bEarth||^2 = (bmeas - M*b)^T*(M^-1)^T * M^-1 * (bmeas - M*b)
2978         // ||bEarth||^2 = (bmeas - M * b)^T * ||M^-1||^2 * (bmeas - M * b)
2979         // ||bEarth||^2 = ||bmeas - M * b||^2 * ||M^-1||^2
2980 
2981         // Where:
2982 
2983         // b = [bx]
2984         //     [by]
2985         //     [bz]
2986 
2987         // M = [m11 	m12 	m13]
2988         //     [m21 	m22 	m23]
2989         //     [m31 	m32 	m33]
2990 
2991         final var gradientEstimator = new GradientEstimator(this::evaluateGeneral);
2992 
2993         final var initialM = Matrix.identity(BodyMagneticFluxDensity.COMPONENTS, BodyMagneticFluxDensity.COMPONENTS);
2994         initialM.add(getInitialMm());
2995 
2996         final var invInitialM = Utils.inverse(initialM);
2997         final var initialBm = getInitialHardIronAsMatrix();
2998         final var initialB = invInitialM.multiplyAndReturnNew(initialBm);
2999 
3000         fitter.setFunctionEvaluator(new LevenbergMarquardtMultiDimensionFunctionEvaluator() {
3001             @Override
3002             public int getNumberOfDimensions() {
3003                 // Input points are measured magnetic flux density coordinates
3004                 return BodyMagneticFluxDensity.COMPONENTS;
3005             }
3006 
3007             @Override
3008             public double[] createInitialParametersArray() {
3009                 final var initial = new double[GENERAL_UNKNOWNS];
3010 
3011                 // biases b
3012                 for (var i = 0; i < BodyMagneticFluxDensity.COMPONENTS; i++) {
3013                     initial[i] = initialB.getElementAtIndex(i);
3014                 }
3015 
3016                 // cross coupling errors M
3017                 final var num = BodyMagneticFluxDensity.COMPONENTS * BodyMagneticFluxDensity.COMPONENTS;
3018                 for (int i = 0, j = BodyMagneticFluxDensity.COMPONENTS; i < num; i++, j++) {
3019                     initial[j] = initialM.getElementAtIndex(i);
3020                 }
3021 
3022                 return initial;
3023             }
3024 
3025             @Override
3026             public double evaluate(
3027                     final int i, final double[] point, final double[] params, final double[] derivatives)
3028                     throws EvaluationException {
3029 
3030                 bmeasX = point[0];
3031                 bmeasY = point[1];
3032                 bmeasZ = point[2];
3033 
3034                 gradientEstimator.gradient(params, derivatives);
3035 
3036                 return evaluateGeneral(params);
3037             }
3038         });
3039 
3040         setInputData();
3041 
3042         fitter.fit();
3043 
3044         final var result = fitter.getA();
3045 
3046         final var bx = result[0];
3047         final var by = result[1];
3048         final var bz = result[2];
3049 
3050         final var m11 = result[3];
3051         final var m21 = result[4];
3052         final var m31 = result[5];
3053 
3054         final var m12 = result[6];
3055         final var m22 = result[7];
3056         final var m32 = result[8];
3057 
3058         final var m13 = result[9];
3059         final var m23 = result[10];
3060         final var m33 = result[11];
3061 
3062         final var mb = new Matrix(BodyMagneticFluxDensity.COMPONENTS, 1);
3063         mb.setElementAtIndex(0, bx);
3064         mb.setElementAtIndex(1, by);
3065         mb.setElementAtIndex(2, bz);
3066 
3067         final var mm = new Matrix(BodyMagneticFluxDensity.COMPONENTS, BodyMagneticFluxDensity.COMPONENTS);
3068         mm.setElementAtIndex(0, m11);
3069         mm.setElementAtIndex(1, m21);
3070         mm.setElementAtIndex(2, m31);
3071 
3072         mm.setElementAtIndex(3, m12);
3073         mm.setElementAtIndex(4, m22);
3074         mm.setElementAtIndex(5, m32);
3075 
3076         mm.setElementAtIndex(6, m13);
3077         mm.setElementAtIndex(7, m23);
3078         mm.setElementAtIndex(8, m33);
3079 
3080         setResult(mm, mb);
3081 
3082         // at this point covariance is expressed in terms of b and M, and must
3083         // be expressed in terms of ba and Ma.
3084         // We know that:
3085 
3086         // b = [bx]
3087         //     [by]
3088         //     [bz]
3089 
3090         // M = [m11 	m12 	m13]
3091         //     [m21 	m22 	m23]
3092         //     [m31 	m32 	m33]
3093 
3094         // and that ba and Ma are expressed as:
3095         // Mm = M - I
3096         // bm = M * b
3097 
3098         // Mm = [m11 - 1    m12         m13    ] =  [sx     mxy     mxz]
3099         //      [m21        m22 - 1     m23    ]    [myx    sy      myz]
3100         //      [m31        m32         m33 - 1]    [mzx    mzy     sz ]
3101 
3102         // bm = [m11 * bx + m12 * by + m13 * bz] = 	[bmx]
3103         //      [m21 * bx + m22 * by + m23 * bz]	[bmy]
3104         //      [m31 * bx + m32 * by + m33 * bz]	[bmz]
3105 
3106         // Defining the linear application:
3107         // F(b, M) = F(bx, by, bz, m11, m21, m31, m12, m22, m32, m13, m23, m33)
3108         // as:
3109         // [bmx] = [m11 * bx + m12 * by + m13 * bz]
3110         // [bmy]   [m21 * bx + m22 * by + m23 * bz]
3111         // [bmz]   [m31 * bx + m32 * by + m33 * bz]
3112         // [sx]	   [m11 - 1]
3113         // [sy]	   [m22 - 1]
3114         // [sz]	   [m33 -1]
3115         // [mxy]   [m12]
3116         // [mxz]   [m13]
3117         // [myx]   [m21]
3118         // [myz]   [m23]
3119         // [mzx]   [m31]
3120         // [mzy]   [m32]
3121 
3122         // Then the Jacobian of F(b, M) is:
3123         // J = 	[m11  m12  m13  bx  0   0   by  0   0   bz  0   0 ]
3124         //	    [m21  m22  m23  0   bx  0   0   by  0   0   bz  0 ]
3125         //	    [m31  m32  m33  0   0   bx  0   0   by  0   0   bz]
3126         //	    [0    0    0    1   0   0   0   0   0   0   0   0 ]
3127         //	    [0    0    0    0   0   0   0   1   0   0   0   0 ]
3128         //	    [0    0    0    0   0   0   0   0   0   0   0   1 ]
3129         //	    [0    0    0    0   0   0   1   0   0   0   0   0 ]
3130         //	    [0    0    0    0   0   0   0   0   0   1   0   0 ]
3131         //	    [0    0    0    0   1   0   0   0   0   0   0   0 ]
3132         //	    [0    0    0    0   0   0   0   0   0   0   1   0 ]
3133         //	    [0    0    0    0   0   1   0   0   0   0   0   0 ]
3134         //	    [0    0    0    0   0   0   0   0   1   0   0   0 ]
3135 
3136         // We know that the propagated covariance is J * Cov * J', hence:
3137         final var jacobian = new Matrix(GENERAL_UNKNOWNS, GENERAL_UNKNOWNS);
3138 
3139         jacobian.setElementAt(0, 0, m11);
3140         jacobian.setElementAt(1, 0, m21);
3141         jacobian.setElementAt(2, 0, m31);
3142 
3143         jacobian.setElementAt(0, 1, m12);
3144         jacobian.setElementAt(1, 1, m22);
3145         jacobian.setElementAt(2, 1, m32);
3146 
3147         jacobian.setElementAt(0, 2, m13);
3148         jacobian.setElementAt(1, 2, m23);
3149         jacobian.setElementAt(2, 2, m33);
3150 
3151         jacobian.setElementAt(0, 3, bx);
3152         jacobian.setElementAt(3, 3, 1.0);
3153 
3154         jacobian.setElementAt(1, 4, bx);
3155         jacobian.setElementAt(8, 4, 1.0);
3156 
3157         jacobian.setElementAt(2, 5, bx);
3158         jacobian.setElementAt(10, 5, 1.0);
3159 
3160         jacobian.setElementAt(0, 6, by);
3161         jacobian.setElementAt(6, 6, 1.0);
3162 
3163         jacobian.setElementAt(1, 7, by);
3164         jacobian.setElementAt(4, 7, 1.0);
3165 
3166         jacobian.setElementAt(2, 8, by);
3167         jacobian.setElementAt(11, 8, 1.0);
3168 
3169         jacobian.setElementAt(0, 9, bz);
3170         jacobian.setElementAt(7, 9, 1.0);
3171 
3172         jacobian.setElementAt(1, 10, bz);
3173         jacobian.setElementAt(9, 10, 1.0);
3174 
3175         jacobian.setElementAt(2, 11, bz);
3176         jacobian.setElementAt(5, 11, 1.0);
3177 
3178         final var jacobianTrans = jacobian.transposeAndReturnNew();
3179         jacobian.multiply(estimatedCovariance);
3180         jacobian.multiply(jacobianTrans);
3181         estimatedCovariance = jacobian;
3182     }
3183 
3184     /**
3185      * Internal method to perform calibration when common z-axis is assumed for both
3186      * the accelerometer and gyroscope.
3187      *
3188      * @throws FittingException                         if Levenberg-Marquardt fails for numerical reasons.
3189      * @throws AlgebraException                         if there are numerical instabilities that prevent
3190      *                                                  matrix inversion.
3191      * @throws com.irurueta.numerical.NotReadyException never happens.
3192      */
3193     private void calibrateCommonAxis() throws AlgebraException, FittingException,
3194             com.irurueta.numerical.NotReadyException {
3195         // The magnetometer model is:
3196         // bmeas = bm + (I + Mm) * btrue + w
3197 
3198         // Ideally a least squares solution tries to minimize noise component, so:
3199         // bmeas = bm + (I + Mm) * btrue
3200 
3201         // For convergence purposes of the Levenberg-Marquardt algorithm, the
3202         // magnetometer model can be better expressed as:
3203         // bmeas = T*K*(btrue + b)
3204         // bmeas = M*(btrue + b)
3205         // bmeas = M*btrue + M*b
3206 
3207         //where:
3208         // M = I + Mm
3209         // bm = M*b = (I + Mm)*b --> b = M^-1*bm
3210 
3211         // We know that the norm of the true body magnetic flux density
3212         // is equal to the amount of Earth magnetic flux density at provided
3213         // position and timestamp
3214         // ||btrue|| = ||bEarth|| --> from 30 µT to 60 µT
3215 
3216         // Hence:
3217         // bmeas - M*b = M*btrue
3218 
3219         // M^-1 * (bmeas - M*b) = btrue
3220 
3221         // ||bEarth||^2 = ||btrue||^2 = (M^-1 * (bmeas - M*b))^T * (M^-1 * (bmeas - M*b))
3222         // ||bEarth||^2 = (bmeas - M*b)^T*(M^-1)^T * M^-1 * (bmeas - M*b)
3223         // ||bEarth||^2 = (bmeas - M * b)^T * ||M^-1||^2 * (bmeas - M * b)
3224         // ||bEarth||^2 = ||bmeas - M * b||^2 * ||M^-1||^2
3225 
3226         // Where:
3227 
3228         // b = [bx]
3229         //     [by]
3230         //     [bz]
3231 
3232         // M = [m11 	m12 	m13]
3233         //     [0 		m22 	m23]
3234         //     [0 	 	0 		m33]
3235 
3236 
3237         final var gradientEstimator = new GradientEstimator(this::evaluateCommonAxis);
3238 
3239         final var initialM = Matrix.identity(BodyMagneticFluxDensity.COMPONENTS, BodyMagneticFluxDensity.COMPONENTS);
3240         initialM.add(getInitialMm());
3241 
3242         // Force initial M to be upper diagonal
3243         initialM.setElementAt(1, 0, 0.0);
3244         initialM.setElementAt(2, 0, 0.0);
3245         initialM.setElementAt(2, 1, 0.0);
3246 
3247         final var invInitialM = Utils.inverse(initialM);
3248         final var initialBm = getInitialHardIronAsMatrix();
3249         final var initialB = invInitialM.multiplyAndReturnNew(initialBm);
3250 
3251         fitter.setFunctionEvaluator(new LevenbergMarquardtMultiDimensionFunctionEvaluator() {
3252             @Override
3253             public int getNumberOfDimensions() {
3254                 // Input points are measured magnetic flux density coordinates
3255                 return BodyKinematics.COMPONENTS;
3256             }
3257 
3258             @Override
3259             public double[] createInitialParametersArray() {
3260                 final var initial = new double[COMMON_Z_AXIS_UNKNOWNS];
3261 
3262                 // biases b
3263                 for (var i = 0; i < BodyMagneticFluxDensity.COMPONENTS; i++) {
3264                     initial[i] = initialB.getElementAtIndex(i);
3265                 }
3266 
3267                 // upper diagonal cross coupling errors M
3268                 var k = BodyMagneticFluxDensity.COMPONENTS;
3269                 for (var j = 0; j < BodyMagneticFluxDensity.COMPONENTS; j++) {
3270                     for (var i = 0; i < BodyMagneticFluxDensity.COMPONENTS; i++) {
3271                         if (i <= j) {
3272                             initial[k] = initialM.getElementAt(i, j);
3273                             k++;
3274                         }
3275                     }
3276                 }
3277 
3278                 return initial;
3279             }
3280 
3281             @Override
3282             public double evaluate(
3283                     final int i, final double[] point, final double[] params, final double[] derivatives)
3284                     throws EvaluationException {
3285 
3286                 bmeasX = point[0];
3287                 bmeasY = point[1];
3288                 bmeasZ = point[2];
3289 
3290                 gradientEstimator.gradient(params, derivatives);
3291 
3292                 return evaluateCommonAxis(params);
3293             }
3294         });
3295 
3296         setInputData();
3297 
3298         fitter.fit();
3299 
3300         final var result = fitter.getA();
3301 
3302         final var bx = result[0];
3303         final var by = result[1];
3304         final var bz = result[2];
3305 
3306         final var m11 = result[3];
3307 
3308         final var m12 = result[4];
3309         final var m22 = result[5];
3310 
3311         final var m13 = result[6];
3312         final var m23 = result[7];
3313         final var m33 = result[8];
3314 
3315         final var mb = new Matrix(BodyMagneticFluxDensity.COMPONENTS, 1);
3316         mb.setElementAtIndex(0, bx);
3317         mb.setElementAtIndex(1, by);
3318         mb.setElementAtIndex(2, bz);
3319 
3320         final var mm = new Matrix(BodyMagneticFluxDensity.COMPONENTS, BodyMagneticFluxDensity.COMPONENTS);
3321         mm.setElementAtIndex(0, m11);
3322         mm.setElementAtIndex(1, 0.0);
3323         mm.setElementAtIndex(2, 0.0);
3324 
3325         mm.setElementAtIndex(3, m12);
3326         mm.setElementAtIndex(4, m22);
3327         mm.setElementAtIndex(5, 0.0);
3328 
3329         mm.setElementAtIndex(6, m13);
3330         mm.setElementAtIndex(7, m23);
3331         mm.setElementAtIndex(8, m33);
3332 
3333         setResult(mm, mb);
3334 
3335         // at this point covariance is expressed in terms of b and M, and must
3336         // be expressed in terms of ba and Ma.
3337         // We know that:
3338 
3339         // b = [bx]
3340         //     [by]
3341         //     [bz]
3342 
3343         // M = [m11 	m12 	m13]
3344         //     [0    	m22 	m23]
3345         //     [0    	0    	m33]
3346 
3347         // m21 = m31 = m32 = 0
3348 
3349         // and that ba and Ma are expressed as:
3350         // Ma = M - I
3351         // ba = M * b
3352 
3353         // Ma = [m11 - 1    m12         m13    ] =  [sx     mxy     mxz]
3354         //      [0          m22 - 1     m23    ]    [0      sy      myz]
3355         //      [0          0           m33 - 1]    [0      0       sz ]
3356 
3357         // ba = [m11 * bx + m12 * by + m13 * bz] = 	[bax]
3358         //      [           m22 * by + m23 * bz]	[bay]
3359         //      [                      m33 * bz]	[baz]
3360 
3361         // Defining the linear application:
3362         // F(b, M) = F(bx, by, bz, m11, m12, m22, m13, m23, m33)
3363         // as:
3364         // [bax] = [m11 * bx + m12 * by + m13 * bz]
3365         // [bay]   [m22 * by + m23 * bz]
3366         // [baz]   [m33 * bz]
3367         // [sx]	   [m11 - 1]
3368         // [sy]	   [m22 - 1]
3369         // [sz]	   [m33 -1]
3370         // [mxy]   [m12]
3371         // [mxz]   [m13]
3372         // [myx]   [0]
3373         // [myz]   [m23]
3374         // [mzx]   [0]
3375         // [mzy]   [0]
3376 
3377         // Then the Jacobian of F(b, M) is:
3378         // J = 	[m11  m12  m13  bx  by  0   bz  0   0 ]
3379         //	    [0    m22  m23  0   0   by  0   bz  0 ]
3380         //	    [0    0    m33  0   0   0   0   0   bz]
3381         //	    [0    0    0    1   0   0   0   0   0 ]
3382         //	    [0    0    0    0   0   1   0   0   0 ]
3383         //	    [0    0    0    0   0   0   0   0   1 ]
3384         //	    [0    0    0    0   1   0   0   0   0 ]
3385         //	    [0    0    0    0   0   0   1   0   0 ]
3386         //	    [0    0    0    0   0   0   0   0   0 ]
3387         //	    [0    0    0    0   0   0   0   1   0 ]
3388         //	    [0    0    0    0   0   0   0   0   0 ]
3389         //	    [0    0    0    0   0   0   0   0   0 ]
3390 
3391         // We know that the propagated covariance is J * Cov * J', hence:
3392         final var jacobian = new Matrix(GENERAL_UNKNOWNS, COMMON_Z_AXIS_UNKNOWNS);
3393 
3394         jacobian.setElementAt(0, 0, m11);
3395 
3396         jacobian.setElementAt(0, 1, m12);
3397         jacobian.setElementAt(1, 1, m22);
3398 
3399         jacobian.setElementAt(0, 2, m13);
3400         jacobian.setElementAt(1, 2, m23);
3401         jacobian.setElementAt(2, 2, m33);
3402 
3403         jacobian.setElementAt(0, 3, bx);
3404         jacobian.setElementAt(3, 3, 1.0);
3405 
3406         jacobian.setElementAt(0, 4, by);
3407         jacobian.setElementAt(6, 4, 1.0);
3408 
3409         jacobian.setElementAt(1, 5, by);
3410         jacobian.setElementAt(4, 5, 1.0);
3411 
3412         jacobian.setElementAt(0, 6, bz);
3413         jacobian.setElementAt(7, 6, 1.0);
3414 
3415         jacobian.setElementAt(1, 7, bz);
3416         jacobian.setElementAt(9, 7, 1.0);
3417 
3418         jacobian.setElementAt(2, 8, bz);
3419         jacobian.setElementAt(5, 8, 1.0);
3420 
3421         final var jacobianTrans = jacobian.transposeAndReturnNew();
3422         jacobian.multiply(estimatedCovariance);
3423         jacobian.multiply(jacobianTrans);
3424         estimatedCovariance = jacobian;
3425     }
3426 
3427     /**
3428      * Makes proper conversion of internal cross-coupling and bias matrices.
3429      *
3430      * @param m internal cross-coupling matrix.
3431      * @param b internal bias matrix.
3432      * @throws AlgebraException if a numerical instability occurs.
3433      */
3434     private void setResult(final Matrix m, final Matrix b) throws AlgebraException {
3435         // Because:
3436         // M = I + Mm
3437         // b = M^-1*bm
3438 
3439         // Then:
3440         // Mm = M - I
3441         // bm = M*b
3442 
3443         if (estimatedHardIron == null) {
3444             estimatedHardIron = new double[BodyMagneticFluxDensity.COMPONENTS];
3445         }
3446 
3447         final var bm = m.multiplyAndReturnNew(b);
3448         bm.toArray(estimatedHardIron);
3449 
3450         if (estimatedMm == null) {
3451             estimatedMm = m;
3452         } else {
3453             estimatedMm.copyFrom(m);
3454         }
3455 
3456         for (var i = 0; i < BodyMagneticFluxDensity.COMPONENTS; i++) {
3457             estimatedMm.setElementAt(i, i, estimatedMm.getElementAt(i, i) - 1.0);
3458         }
3459 
3460         estimatedCovariance = fitter.getCovar();
3461         estimatedChiSq = fitter.getChisq();
3462         estimatedChiSqDegreesOfFreedom = fitter.getChisqDegreesOfFreedom();
3463         estimatedReducedChiSq = fitter.getReducedChisq();
3464         estimatedMse = fitter.getMse();
3465         try {
3466             estimatedP = fitter.getP();
3467             estimatedQ = fitter.getQ();
3468         } catch (final MaxIterationsExceededException ignore) {
3469             // if numerical instabilities arise, we assume worst case (no fit at all)
3470             // probability of finding a smaller chi square value is 1.0
3471             // quality of fit is 0.0
3472             estimatedP = 1.0;
3473             estimatedQ = 0.0;
3474         }
3475     }
3476 
3477     /**
3478      * Computes estimated true magnetic flux density squared norm using current measured
3479      * body magnetic flux density and provided parameters for the general case.
3480      * This method is internally executed during gradient estimation and
3481      * Levenberg-Marquardt fitting needed for calibration computation.
3482      *
3483      * @param params array containing current parameters for the general purpose case.
3484      *               Must have length 12.
3485      * @return estimated true specific force squared norm.
3486      * @throws EvaluationException if there are numerical instabilities.
3487      */
3488     private double evaluateGeneral(final double[] params) throws EvaluationException {
3489         final var bx = params[0];
3490         final var by = params[1];
3491         final var bz = params[2];
3492 
3493         final var m11 = params[3];
3494         final var m21 = params[4];
3495         final var m31 = params[5];
3496 
3497         final var m12 = params[6];
3498         final var m22 = params[7];
3499         final var m32 = params[8];
3500 
3501         final var m13 = params[9];
3502         final var m23 = params[10];
3503         final var m33 = params[11];
3504 
3505         return evaluate(bx, by, bz, m11, m21, m31, m12, m22, m32, m13, m23, m33);
3506     }
3507 
3508     /**
3509      * Computes estimated true magnetic flux density squared norm using current measured
3510      * body magnetic flux density and provided parameters when common z-axis is assumed.
3511      * This method is internally executed during gradient estimation and
3512      * Levenberg-Marquardt fitting needed for calibration computation.
3513      *
3514      * @param params array containing current parameters for the common z-axis case.
3515      *               Must have length 9.
3516      * @return estimated true specific force squared norm.
3517      * @throws EvaluationException if there are numerical instabilities.
3518      */
3519     private double evaluateCommonAxis(final double[] params) throws EvaluationException {
3520         final var bx = params[0];
3521         final var by = params[1];
3522         final var bz = params[2];
3523 
3524         final var m11 = params[3];
3525 
3526         final var m12 = params[4];
3527         final var m22 = params[5];
3528 
3529         final var m13 = params[6];
3530         final var m23 = params[7];
3531         final var m33 = params[8];
3532 
3533         return evaluate(bx, by, bz, m11, 0.0, 0.0, m12, m22, 0.0, m13, m23, m33);
3534     }
3535 
3536     /**
3537      * Computes estimated true magnetic flux density squared norm using current measured
3538      * body magnetic flux density and provided parameters.
3539      * This method is internally executed during gradient estimation and
3540      * Levenberg-Marquardt fitting needed for calibration computation.
3541      *
3542      * @param bx  x-coordinate of bias.
3543      * @param by  y-coordinate of bias.
3544      * @param bz  z-coordinate of bias.
3545      * @param m11 element 1,1 of cross-coupling error matrix.
3546      * @param m21 element 2,1 of cross-coupling error matrix.
3547      * @param m31 element 3,1 of cross-coupling error matrix.
3548      * @param m12 element 1,2 of cross-coupling error matrix.
3549      * @param m22 element 2,2 of cross-coupling error matrix.
3550      * @param m32 element 3,2 of cross-coupling error matrix.
3551      * @param m13 element 1,3 of cross-coupling error matrix.
3552      * @param m23 element 2,3 of cross-coupling error matrix.
3553      * @param m33 element 3,3 of cross-coupling error matrix.
3554      * @return estimated true specific force squared norm.
3555      * @throws EvaluationException if there are numerical instabilities.
3556      */
3557     private double evaluate(final double bx, final double by, final double bz,
3558                             final double m11, final double m21, final double m31,
3559                             final double m12, final double m22, final double m32,
3560                             final double m13, final double m23, final double m33) throws EvaluationException {
3561 
3562         // bmeas = M*(btrue + b)
3563 
3564         // btrue = M^-1*bmeas - b
3565 
3566         try {
3567             if (bmeas == null) {
3568                 bmeas = new Matrix(BodyMagneticFluxDensity.COMPONENTS, 1);
3569             }
3570             if (m == null) {
3571                 m = new Matrix(BodyMagneticFluxDensity.COMPONENTS, BodyMagneticFluxDensity.COMPONENTS);
3572             }
3573             if (invM == null) {
3574                 invM = new Matrix(BodyMagneticFluxDensity.COMPONENTS, BodyMagneticFluxDensity.COMPONENTS);
3575             }
3576             if (b == null) {
3577                 b = new Matrix(BodyMagneticFluxDensity.COMPONENTS, 1);
3578             }
3579             if (btrue == null) {
3580                 btrue = new Matrix(BodyMagneticFluxDensity.COMPONENTS, 1);
3581             }
3582 
3583             bmeas.setElementAtIndex(0, bmeasX);
3584             bmeas.setElementAtIndex(1, bmeasY);
3585             bmeas.setElementAtIndex(2, bmeasZ);
3586 
3587             m.setElementAt(0, 0, m11);
3588             m.setElementAt(1, 0, m21);
3589             m.setElementAt(2, 0, m31);
3590 
3591             m.setElementAt(0, 1, m12);
3592             m.setElementAt(1, 1, m22);
3593             m.setElementAt(2, 1, m32);
3594 
3595             m.setElementAt(0, 2, m13);
3596             m.setElementAt(1, 2, m23);
3597             m.setElementAt(2, 2, m33);
3598 
3599             Utils.inverse(m, invM);
3600 
3601             b.setElementAtIndex(0, bx);
3602             b.setElementAtIndex(1, by);
3603             b.setElementAtIndex(2, bz);
3604 
3605             invM.multiply(bmeas, btrue);
3606             btrue.subtract(b);
3607 
3608             final var norm = Utils.normF(btrue);
3609             return norm * norm;
3610 
3611         } catch (final AlgebraException e) {
3612             throw new EvaluationException(e);
3613         }
3614     }
3615 
3616     /**
3617      * Converts magnetic flux density value and unit to Teslas.
3618      *
3619      * @param value magnetic flux density value.
3620      * @param unit  unit of magnetic flux density value.
3621      * @return converted value.
3622      */
3623     private static double convertMagneticFluxDensity(final double value, final MagneticFluxDensityUnit unit) {
3624         return MagneticFluxDensityConverter.convert(value, unit, MagneticFluxDensityUnit.TESLA);
3625     }
3626 
3627     /**
3628      * Converts magnetic flux density instance to Teslas.
3629      *
3630      * @param magneticFluxDensity magnetic flux density instance to be converted.
3631      * @return converted value.
3632      */
3633     private static double convertMagneticFluxDensity(final MagneticFluxDensity magneticFluxDensity) {
3634         return convertMagneticFluxDensity(magneticFluxDensity.getValue().doubleValue(), magneticFluxDensity.getUnit());
3635     }
3636 }