View Javadoc
1   /*
2    * Copyright (C) 2020 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.gyroscope;
17  
18  import com.irurueta.algebra.AlgebraException;
19  import com.irurueta.algebra.ArrayUtils;
20  import com.irurueta.algebra.Matrix;
21  import com.irurueta.algebra.Utils;
22  import com.irurueta.algebra.WrongSizeException;
23  import com.irurueta.geometry.InhomogeneousPoint3D;
24  import com.irurueta.geometry.Quaternion;
25  import com.irurueta.geometry.RotationException;
26  import com.irurueta.navigation.LockedException;
27  import com.irurueta.navigation.NotReadyException;
28  import com.irurueta.navigation.inertial.BodyKinematics;
29  import com.irurueta.navigation.inertial.INSLooselyCoupledKalmanInitializerConfig;
30  import com.irurueta.navigation.inertial.INSTightlyCoupledKalmanInitializerConfig;
31  import com.irurueta.navigation.inertial.calibration.AccelerationFixer;
32  import com.irurueta.navigation.inertial.calibration.AngularSpeedTriad;
33  import com.irurueta.navigation.inertial.calibration.BodyKinematicsSequence;
34  import com.irurueta.navigation.inertial.calibration.CalibrationException;
35  import com.irurueta.navigation.inertial.calibration.GyroscopeBiasUncertaintySource;
36  import com.irurueta.navigation.inertial.calibration.GyroscopeCalibrationSource;
37  import com.irurueta.navigation.inertial.calibration.StandardDeviationTimedBodyKinematics;
38  import com.irurueta.numerical.EvaluationException;
39  import com.irurueta.numerical.GradientEstimator;
40  import com.irurueta.numerical.fitting.FittingException;
41  import com.irurueta.numerical.fitting.LevenbergMarquardtMultiDimensionFitter;
42  import com.irurueta.numerical.fitting.LevenbergMarquardtMultiDimensionFunctionEvaluator;
43  import com.irurueta.statistics.MaxIterationsExceededException;
44  import com.irurueta.units.Acceleration;
45  import com.irurueta.units.AccelerationConverter;
46  import com.irurueta.units.AccelerationUnit;
47  import com.irurueta.units.AngularSpeed;
48  import com.irurueta.units.AngularSpeedConverter;
49  import com.irurueta.units.AngularSpeedUnit;
50  
51  import java.util.ArrayList;
52  import java.util.List;
53  
54  /**
55   * Estimates gyroscope biases, cross couplings, scaling factors and
56   * G-dependent cross biases introduced on the gyroscope by the specific
57   * forces sensed by the accelerometer.
58   * <p>
59   * This calibrator assumes that the IMU is at a more or less fixed location on
60   * Earth, and evaluates sequences of measured body kinematics to perform
61   * calibration for unknown orientations on those provided sequences.
62   * Each provided sequence will be preceded by a static period where mean
63   * specific force will be measured to determine gravity (and hence partial
64   * body attitude).
65   * <p>
66   * To use this calibrator at least 10 sequences are needed (each containing
67   * at least 3 body kinematics measurements) when common z-axis is assumed and
68   * G-dependant cross biases are ignored, otherwise at least 13 sequences are
69   * required (each containing at least 3 body kinematics measurements) when
70   * common z-axis is not assumed.
71   * If G-dependent cross biases are being estimated, then at least 19
72   * measurements are needed when common z-axis is assumed, otherwise at least
73   * 22 sequences are required (each containing at least 3 body kinematics
74   * measurements) when common z-axis is not assumed.
75   * <p>
76   * Measured gyroscope angular rates is assumed to follow the model shown below:
77   * <pre>
78   *     Ωmeas = bg + (I + Mg) * Ωtrue + Gg * ftrue + w
79   * </pre>
80   * Where:
81   * - Ωmeas is the measured gyroscope angular rates. This is a 3x1 vector.
82   * - bg is the gyroscope bias. Ideally, on a perfect gyroscope, this should be a
83   * 3x1 zero vector.
84   * - I is the 3x3 identity matrix.
85   * - Mg is the 3x3 matrix containing cross-couplings and scaling factors. Ideally, on
86   * a perfect gyroscope, this should be a 3x3 zero matrix.
87   * - Ωtrue is ground-truth gyroscope angular rates.
88   * - Gg is the G-dependent cross biases introduced by the specific forces sensed
89   * by the accelerometer. Ideally, on a perfect gyroscope, this should be a 3x3
90   * zero matrix.
91   * - ftrue is ground-truth specific force. This is a 3x1 vector.
92   * - w is measurement noise. This is a 3x1 vector.
93   * <p>
94   * This calibrator is based on the ideas of
95   * David Tedaldi, Alberto Pretto, Emmanuelle Menegatti. A Robust and Easy to
96   * Implement Method for IMU Calibration without External Equipments.
97   */
98  public class EasyGyroscopeCalibrator implements
99          GyroscopeNonLinearCalibrator, UnknownBiasGyroscopeCalibrator, GyroscopeCalibrationSource,
100         GyroscopeBiasUncertaintySource, OrderedBodyKinematicsSequenceGyroscopeCalibrator,
101         AccelerometerDependentGyroscopeCalibrator {
102 
103     /**
104      * Indicates whether by default a common z-axis is assumed for both the accelerometer
105      * and gyroscope.
106      */
107     public static final boolean DEFAULT_USE_COMMON_Z_AXIS = true;
108 
109     /**
110      * Indicates that by default G-dependent cross biases introduced
111      * by the accelerometer on the gyroscope are estimated.
112      */
113     public static final boolean DEFAULT_ESTIMATE_G_DEPENDENT_CROSS_BIASES = true;
114 
115     /**
116      * Number of unknowns when common z-axis is assumed for both the accelerometer
117      * and gyroscope when G-dependent cross biases are being estimated.
118      */
119     public static final int COMMON_Z_AXIS_UNKNOWNS_AND_CROSS_BIASES = 18;
120 
121     /**
122      * Number of unknowns for the general case when G-dependent cross
123      * biases are being estimated.
124      */
125     public static final int GENERAL_UNKNOWNS_AND_CROSS_BIASES = 21;
126 
127     /**
128      * Number of unknowns when common z-axis is assumed for both
129      * the accelerometer and gyroscope when G-dependent cross biases
130      * are not being estimated.
131      */
132     public static final int COMMON_Z_AXIS_UNKNOWNS = 9;
133 
134     /**
135      * Number of unknowns for the general case when G-dependent cross
136      * biases are not being estimated.
137      */
138     public static final int GENERAL_UNKNOWNS = 12;
139 
140     /**
141      * Required minimum number of sequences when common z-axis is assumed
142      * and G-dependent cross biases are being estimated.
143      */
144     public static final int MINIMUM_SEQUENCES_COMMON_Z_AXIS_AND_CROSS_BIASES =
145             COMMON_Z_AXIS_UNKNOWNS_AND_CROSS_BIASES + 1;
146 
147     /**
148      * Required minimum number of sequences for the general case and
149      * G-dependent cross biases are being estimated.
150      */
151     public static final int MINIMUM_SEQUENCES_GENERAL_AND_CROSS_BIASES =
152             GENERAL_UNKNOWNS_AND_CROSS_BIASES + 1;
153 
154     /**
155      * Required minimum number of sequences when common z-axis is assumed
156      * and G-dependent cross biases are being ignored.
157      */
158     public static final int MINIMUM_SEQUENCES_COMMON_Z_AXIS = COMMON_Z_AXIS_UNKNOWNS + 1;
159 
160     /**
161      * Required minimum number of sequences for the general case and
162      * G-dependent cross biases are being ignored.
163      */
164     public static final int MINIMUM_SEQUENCES_GENERAL = GENERAL_UNKNOWNS + 1;
165 
166     /**
167      * Levenberg-Marquardt fitter to find a non-linear solution.
168      */
169     private final LevenbergMarquardtMultiDimensionFitter fitter = new LevenbergMarquardtMultiDimensionFitter();
170 
171     /**
172      * Known x-coordinate of accelerometer bias to be used to fix measured
173      * specific force and find cross biases introduced by the accelerometer.
174      * This is expressed in meters per squared second (m/s^2).
175      */
176     private double accelerometerBiasX;
177 
178     /**
179      * Known y-coordinate of accelerometer bias to be used to fix measured
180      * specific force and find cross biases introduced by the accelerometer.
181      * This is expressed in meters per squared second (m/s^2).
182      */
183     private double accelerometerBiasY;
184 
185     /**
186      * Known z-coordinate of accelerometer bias to be used to fix measured
187      * specific force and find cross biases introduced by the accelerometer.
188      * This is expressed in meters per squared second (m/s^2).
189      */
190     private double accelerometerBiasZ;
191 
192     /**
193      * Known accelerometer x scaling factor to be used to fix measured
194      * specific force and find cross biases introduced by the accelerometer.
195      */
196     private double accelerometerSx;
197 
198     /**
199      * Known accelerometer y scaling factor to be used to fix measured
200      * specific force and find cross biases introduced by the accelerometer.
201      */
202     private double accelerometerSy;
203 
204     /**
205      * Known accelerometer z scaling factor to be used to fix measured
206      * specific force and find cross biases introduced by the accelerometer.
207      */
208     private double accelerometerSz;
209 
210     /**
211      * Known accelerometer x-y cross coupling error to be used to fix measured
212      * specific force and find cross biases introduced by the accelerometer.
213      */
214     private double accelerometerMxy;
215 
216     /**
217      * Know accelerometer x-z cross coupling error to be used to fix measured
218      * specific force and find cross biases introduced by the accelerometer.
219      */
220     private double accelerometerMxz;
221 
222     /**
223      * Known accelerometer y-x cross coupling error to be used to fix measured
224      * specific force and find cross biases introduced by the accelerometer.
225      */
226     private double accelerometerMyx;
227 
228     /**
229      * Known accelerometer y-z cross coupling error to be used to fix measured
230      * specific force and find cross biases introduced by the accelerometer.
231      */
232     private double accelerometerMyz;
233 
234     /**
235      * Known accelerometer z-x cross coupling error to be used to fix measured
236      * specific force and find cross biases introduced by the accelerometer.
237      */
238     private double accelerometerMzx;
239 
240     /**
241      * Known accelerometer z-y cross coupling error to be used to fix measured
242      * specific force and find cross biases introduced by the accelerometer.
243      */
244     private double accelerometerMzy;
245 
246     /**
247      * Initial x-coordinate of gyroscope bias to be used to find a solution.
248      * This is expressed in radians per second (rad/s).
249      */
250     private double initialBiasX;
251 
252     /**
253      * Initial y-coordinate of gyroscope bias to be used to find a solution.
254      * This is expressed in radians per second (rad/s).
255      */
256     private double initialBiasY;
257 
258     /**
259      * Initial z-coordinate of gyroscope bias to be used to find a solution.
260      * This is expressed in radians per second (rad/s).
261      */
262     private double initialBiasZ;
263 
264     /**
265      * Initial gyroscope x scaling factor.
266      */
267     private double initialSx;
268 
269     /**
270      * Initial gyroscope y scaling factor.
271      */
272     private double initialSy;
273 
274     /**
275      * Initial gyroscope z scaling factor.
276      */
277     private double initialSz;
278 
279     /**
280      * Initial gyroscope x-y cross coupling error.
281      */
282     private double initialMxy;
283 
284     /**
285      * Initial gyroscope x-z cross coupling error.
286      */
287     private double initialMxz;
288 
289     /**
290      * Initial gyroscope y-x cross coupling error.
291      */
292     private double initialMyx;
293 
294     /**
295      * Initial gyroscope y-z cross coupling error.
296      */
297     private double initialMyz;
298 
299     /**
300      * Initial gyroscope z-x cross coupling error.
301      */
302     private double initialMzx;
303 
304     /**
305      * Initial gyroscope z-y cross coupling error.
306      */
307     private double initialMzy;
308 
309     /**
310      * Initial G-dependent cross biases introduced on the gyroscope by the
311      * specific forces sensed by the accelerometer.
312      */
313     private Matrix initialGg;
314 
315     /**
316      * Contains a collection of sequences of timestamped body kinematics
317      * measurements taken at a given position where the device moves freely
318      * with different orientations.
319      */
320     private List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences;
321 
322     /**
323      * This flag indicates whether z-axis is assumed to be common for accelerometer
324      * and gyroscope.
325      * When enabled, this eliminates 3 variables from Mg matrix.
326      */
327     private boolean commonAxisUsed = DEFAULT_USE_COMMON_Z_AXIS;
328 
329     /**
330      * This flag indicates whether G-dependent cross biases are being
331      * estimated or not.
332      * When enabled, this adds 9 variables from Gg matrix.
333      */
334     private boolean estimateGDependentCrossBiases = DEFAULT_ESTIMATE_G_DEPENDENT_CROSS_BIASES;
335 
336     /**
337      * Listener to handle events raised by this calibrator.
338      */
339     private EasyGyroscopeCalibratorListener listener;
340 
341     /**
342      * Estimated angular rate biases for each IMU axis expressed in radians per
343      * second (rad/s).
344      */
345     private double[] estimatedBiases;
346 
347     /**
348      * Estimated gyroscope scale factors and cross coupling errors.
349      * This is the product of matrix Tg containing cross coupling errors and Kg
350      * containing scaling factors.
351      * So that:
352      * <pre>
353      *     Mg = [sx    mxy  mxz] = Tg*Kg
354      *          [myx   sy   myz]
355      *          [mzx   mzy  sz ]
356      * </pre>
357      * Where:
358      * <pre>
359      *     Kg = [sx 0   0 ]
360      *          [0  sy  0 ]
361      *          [0  0   sz]
362      * </pre>
363      * and
364      * <pre>
365      *     Tg = [1          -alphaXy    alphaXz ]
366      *          [alphaYx    1           -alphaYz]
367      *          [-alphaZx   alphaZy     1       ]
368      * </pre>
369      * Hence:
370      * <pre>
371      *     Mg = [sx    mxy  mxz] = Tg*Kg =  [sx             -sy * alphaXy   sz * alphaXz ]
372      *          [myx   sy   myz]            [sx * alphaYx   sy              -sz * alphaYz]
373      *          [mzx   mzy  sz ]            [-sx * alphaZx  sy * alphaZy    sz           ]
374      * </pre>
375      * This instance allows any 3x3 matrix however, typically alphaYx, alphaZx and alphaZy
376      * are considered to be zero if the gyroscope z-axis is assumed to be the same
377      * as the body z-axis. When this is assumed, myx = mzx = mzy = 0 and the Mg matrix
378      * becomes upper diagonal:
379      * <pre>
380      *     Mg = [sx    mxy  mxz]
381      *          [0     sy   myz]
382      *          [0     0    sz ]
383      * </pre>
384      * Values of this matrix are unit-less.
385      */
386     private Matrix estimatedMg;
387 
388     /**
389      * Estimated G-dependent cross biases introduced on the gyroscope by the
390      * specific forces sensed by the accelerometer.
391      * This instance allows any 3x3 matrix.
392      */
393     private Matrix estimatedGg;
394 
395     /**
396      * Estimated covariance matrix for estimated parameters.
397      */
398     private Matrix estimatedCovariance;
399 
400     /**
401      * Estimated chi square value.
402      */
403     private double estimatedChiSq;
404 
405     /**
406      * Estimated degrees of freedom of chi square value. Degrees of freedom is equal to the number of sampled data
407      * minus the number of estimated parameters.
408      */
409     private int estimatedChiSqDegreesOfFreedom;
410 
411     /**
412      * Estimated reduced chi square value. This is equal to estimated chi square value divided by its degrees of
413      * freedom. Ideally this value should be close to 1.0.
414      */
415     private double estimatedReducedChiSq;
416 
417     /**
418      * Estimated mean square error respect to provided measurements.
419      */
420     private double estimatedMse;
421 
422     /**
423      * Estimated probability of finding a smaller chi square value expressed as a value between 0.0 and 1.0. The smaller
424      * the found chi square value is, the better the fit of the estimated parameters to the actual parameter. Thus, the
425      * smaller the chance of finding a smaller chi square value, then the better the estimated fit is.
426      */
427     private double estimatedP;
428 
429     /**
430      * Estimated measure of quality of estimated fit as a value between 0.0 and 1.0. The larger the quality value is,
431      * the better the fit that has been estimated.
432      */
433     private double estimatedQ;
434 
435     /**
436      * Indicates whether calibrator is running.
437      */
438     private boolean running;
439 
440     /**
441      * Acceleration fixer.
442      */
443     private final AccelerationFixer accelerationFixer = new AccelerationFixer();
444 
445     /**
446      * Index of current point being evaluated.
447      */
448     private int i;
449 
450     /**
451      * Holds integrated rotation of a sequence.
452      */
453     private final Quaternion q = new Quaternion();
454 
455     /**
456      * Contains a copy of input sequences where fixed body kinematics will be updated.
457      */
458     private List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> fixedSequences;
459 
460     /**
461      * Contains measured specific force on a sample within a sequence
462      * expressed as a 3x1 matrix.
463      */
464     private Matrix measuredSpecificForce;
465 
466     /**
467      * Contains fixed specific force for a sample within a sequence
468      * using provided accelerometer parameters and expressed as a 3x1
469      * matrix.
470      */
471     private Matrix trueSpecificForce;
472 
473     /**
474      * Contains measured angular rate on a sample within a sequence
475      * expressed as a 3x1 matrix.
476      */
477     private Matrix measuredAngularRate;
478 
479     /**
480      * Contains fixed specific force for a sample within a sequence
481      * using current parameters being estimated.
482      */
483     private Matrix trueAngularRate;
484 
485     /**
486      * Internally holds cross-coupling errors during calibration.
487      */
488     private Matrix m;
489 
490     /**
491      * Internally holds inverse of cross-coupling errors during calibration.
492      */
493     private Matrix invM;
494 
495     /**
496      * Internally holds biases during calibration.
497      */
498     private Matrix b;
499 
500     /**
501      * Internally hold g-dependent cross biases during calibration.
502      */
503     private Matrix g;
504 
505     /**
506      * Internally holds angular rate bias due to g-dependent cross biases
507      */
508     private Matrix tmp;
509 
510     /**
511      * Holds normalized gravity versor at the start of a sequence.
512      * This is used to compute rotations of gravity versor during
513      * calibration.
514      */
515     private final InhomogeneousPoint3D startPoint = new InhomogeneousPoint3D();
516 
517     /**
518      * Holds normalized gravity versor at the end of a sequence.
519      * This is used to compute rotations of gravity versor during
520      * calibration.
521      */
522     private final InhomogeneousPoint3D endPoint = new InhomogeneousPoint3D();
523 
524     /**
525      * Contains expected gravity versor obtained from measurements fixed
526      * using known accelerometer parameters.
527      * This is reused for memory efficiency purposes during calibration.
528      */
529     private final InhomogeneousPoint3D expectedEndPoint = new InhomogeneousPoint3D();
530 
531     /**
532      * Array containing normalized gravity versor before (former 3 values)
533      * and after (latter 3 values) a given sequence.
534      * This is used during calibration.
535      */
536     private double[] point;
537 
538     /**
539      * Constructor.
540      */
541     public EasyGyroscopeCalibrator() {
542         try {
543             initialGg = new Matrix(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
544         } catch (final WrongSizeException ignore) {
545             // never happens
546         }
547     }
548 
549     /**
550      * Constructor.
551      *
552      * @param sequences   collection of sequences containing timestamped body
553      *                    kinematics measurements.
554      * @param initialBias initial gyroscope bias to be used to find a solution.
555      *                    This must be 3x1 and is expressed in radians per
556      *                    second (rad/s).
557      * @param initialMg   initial gyroscope scale factors and cross coupling
558      *                    errors matrix. Must be 3x3.
559      * @param initialGg   initial gyroscope G-dependent cross biases
560      *                    introduced on the gyroscope by the specific forces
561      *                    sensed by the accelerometer. Must be 3x3.
562      * @throws IllegalArgumentException if any of the provided values does
563      *                                  not have proper size.
564      */
565     public EasyGyroscopeCalibrator(
566             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
567             final Matrix initialBias, final Matrix initialMg, final Matrix initialGg) {
568         this();
569         this.sequences = sequences;
570         try {
571             setInitialBias(initialBias);
572             setInitialMg(initialMg);
573             setInitialGg(initialGg);
574         } catch (final LockedException ignore) {
575             // never happens
576         }
577     }
578 
579     /**
580      * Constructor.
581      *
582      * @param sequences   collection of sequences containing timestamped body
583      *                    kinematics measurements.
584      * @param initialBias initial gyroscope bias to be used to find a solution.
585      *                    This must be 3x1 and is expressed in radians per
586      *                    second (rad/s).
587      * @param initialMg   initial gyroscope scale factors and cross coupling
588      *                    errors matrix. Must be 3x3.
589      * @param initialGg   initial gyroscope G-dependent cross biases
590      *                    introduced on the gyroscope by the specific forces
591      *                    sensed by the accelerometer. Must be 3x3.
592      * @param listener    listener to handle events raised by this
593      *                    calibrator.
594      * @throws IllegalArgumentException if any of the provided values does
595      *                                  not have proper size.
596      */
597     public EasyGyroscopeCalibrator(
598             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
599             final Matrix initialBias, final Matrix initialMg, final Matrix initialGg,
600             final EasyGyroscopeCalibratorListener listener) {
601         this(sequences, initialBias, initialMg, initialGg);
602         this.listener = listener;
603     }
604 
605     /**
606      * Constructor.
607      *
608      * @param sequences   collection of sequences containing timestamped body
609      *                    kinematics measurements.
610      * @param initialBias initial gyroscope bias to be used to find a
611      *                    solution. This must have length 3 and is expressed
612      *                    in radians per second (rad/s).
613      * @param initialMg   initial gyroscope scale factors and cross coupling
614      *                    errors matrix. Must be 3x3.
615      * @param initialGg   initial gyroscope G-dependent cross biases
616      *                    introduced on the gyroscope by the specific forces
617      *                    sensed by the accelerometer. Must be 3x3.
618      * @throws IllegalArgumentException if any of the provided values does
619      *                                  not have proper size.
620      */
621     public EasyGyroscopeCalibrator(
622             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
623             final double[] initialBias, final Matrix initialMg, final Matrix initialGg) {
624         this();
625         this.sequences = sequences;
626         try {
627             setInitialBias(initialBias);
628             setInitialMg(initialMg);
629             setInitialGg(initialGg);
630         } catch (final LockedException ignore) {
631             // never happens
632         }
633     }
634 
635     /**
636      * Constructor.
637      *
638      * @param sequences   collection of sequences containing timestamped body
639      *                    kinematics measurements.
640      * @param initialBias initial gyroscope bias to be used to find a
641      *                    solution. This must have length 3 and is expressed
642      *                    in radians per second (rad/s).
643      * @param initialMg   initial gyroscope scale factors and cross coupling
644      *                    errors matrix. Must be 3x3.
645      * @param initialGg   initial gyroscope G-dependent cross biases
646      *                    introduced on the gyroscope by the specific forces
647      *                    sensed by the accelerometer. Must be 3x3.
648      * @param listener    listener to handle events raised by this
649      *                    calibrator.
650      * @throws IllegalArgumentException if any of the provided values does
651      *                                  not have proper size.
652      */
653     public EasyGyroscopeCalibrator(
654             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
655             final double[] initialBias, final Matrix initialMg, final Matrix initialGg,
656             final EasyGyroscopeCalibratorListener listener) {
657         this(sequences, initialBias, initialMg, initialGg);
658         this.listener = listener;
659     }
660 
661     /**
662      * Constructor.
663      *
664      * @param sequences         collection of sequences containing timestamped body
665      *                          kinematics measurements.
666      * @param initialBias       initial gyroscope bias to be used to find a
667      *                          solution. This must have length 3 and is expressed
668      *                          in radians per second (rad/s).
669      * @param initialMg         initial gyroscope scale factors and cross coupling
670      *                          errors matrix. Must be 3x3.
671      * @param initialGg         initial gyroscope G-dependent cross biases
672      *                          introduced on the gyroscope by the specific forces
673      *                          sensed by the accelerometer. Must be 3x3.
674      * @param accelerometerBias known accelerometer bias. This must
675      *                          have length 3 and is expressed in
676      *                          meters per squared second
677      *                          (m/s^2).
678      * @param accelerometerMa   known accelerometer scale factors and
679      *                          cross coupling matrix. Must be 3x3.
680      * @throws IllegalArgumentException if any of the provided values does
681      *                                  not have proper size.
682      */
683     public EasyGyroscopeCalibrator(
684             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
685             final double[] initialBias, final Matrix initialMg, final Matrix initialGg,
686             final double[] accelerometerBias, final Matrix accelerometerMa) {
687         this(sequences, initialBias, initialMg, initialGg);
688         try {
689             setAccelerometerBias(accelerometerBias);
690             setAccelerometerMa(accelerometerMa);
691         } catch (final LockedException ignore) {
692             // never happens
693         }
694     }
695 
696     /**
697      * Constructor.
698      *
699      * @param sequences         collection of sequences containing timestamped body
700      *                          kinematics measurements.
701      * @param initialBias       initial gyroscope bias to be used to find a
702      *                          solution. This must have length 3 and is expressed
703      *                          in radians per second (rad/s).
704      * @param initialMg         initial gyroscope scale factors and cross coupling
705      *                          errors matrix. Must be 3x3.
706      * @param initialGg         initial gyroscope G-dependent cross biases
707      *                          introduced on the gyroscope by the specific forces
708      *                          sensed by the accelerometer. Must be 3x3.
709      * @param accelerometerBias known accelerometer bias. This must
710      *                          have length 3 and is expressed in
711      *                          meters per squared second
712      *                          (m/s^2).
713      * @param accelerometerMa   known accelerometer scale factors and
714      *                          cross coupling matrix. Must be 3x3.
715      * @param listener          listener to handle events raised by this
716      *                          calibrator.
717      * @throws IllegalArgumentException if any of the provided values does
718      *                                  not have proper size.
719      */
720     public EasyGyroscopeCalibrator(
721             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
722             final double[] initialBias, final Matrix initialMg, final Matrix initialGg,
723             final double[] accelerometerBias, final Matrix accelerometerMa,
724             final EasyGyroscopeCalibratorListener listener) {
725         this(sequences, initialBias, initialMg, initialGg, accelerometerBias, accelerometerMa);
726         this.listener = listener;
727     }
728 
729     /**
730      * Constructor.
731      *
732      * @param sequences         collection of sequences containing timestamped body
733      *                          kinematics measurements.
734      * @param initialBias       initial gyroscope bias to be used to find a
735      *                          solution. This must be 3x1 and is expressed
736      *                          in radians per second (rad/s).
737      * @param initialMg         initial gyroscope scale factors and cross coupling
738      *                          errors matrix. Must be 3x3.
739      * @param initialGg         initial gyroscope G-dependent cross biases
740      *                          introduced on the gyroscope by the specific forces
741      *                          sensed by the accelerometer. Must be 3x3.
742      * @param accelerometerBias known accelerometer bias. This must be 3x1
743      *                          and is expressed in meters per squared
744      *                          second (m/s^2).
745      * @param accelerometerMa   known accelerometer scale factors and
746      *                          cross coupling matrix. Must be 3x3.
747      * @throws IllegalArgumentException if any of the provided values does
748      *                                  not have proper size.
749      */
750     public EasyGyroscopeCalibrator(
751             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
752             final Matrix initialBias, final Matrix initialMg, final Matrix initialGg,
753             final Matrix accelerometerBias, final Matrix accelerometerMa) {
754         this(sequences, initialBias, initialMg, initialGg);
755         try {
756             setAccelerometerBias(accelerometerBias);
757             setAccelerometerMa(accelerometerMa);
758         } catch (final LockedException ignore) {
759             // never happens
760         }
761     }
762 
763     /**
764      * Constructor.
765      *
766      * @param sequences         collection of sequences containing timestamped body
767      *                          kinematics measurements.
768      * @param initialBias       initial gyroscope bias to be used to find a
769      *                          solution. This must be 3x1 and is expressed
770      *                          in radians per second (rad/s).
771      * @param initialMg         initial gyroscope scale factors and cross coupling
772      *                          errors matrix. Must be 3x3.
773      * @param initialGg         initial gyroscope G-dependent cross biases
774      *                          introduced on the gyroscope by the specific forces
775      *                          sensed by the accelerometer. Must be 3x3.
776      * @param accelerometerBias known accelerometer bias. This must be 3x1
777      *                          and is expressed in meters per squared
778      *                          second (m/s^2).
779      * @param accelerometerMa   known accelerometer scale factors and
780      *                          cross coupling matrix. Must be 3x3.
781      * @param listener          listener to handle events raised by this
782      *                          calibrator.
783      * @throws IllegalArgumentException if any of the provided values does
784      *                                  not have proper size.
785      */
786     public EasyGyroscopeCalibrator(
787             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
788             final Matrix initialBias, final Matrix initialMg, final Matrix initialGg, final Matrix accelerometerBias,
789             final Matrix accelerometerMa, final EasyGyroscopeCalibratorListener listener) {
790         this(sequences, initialBias, initialMg, initialGg, accelerometerBias, accelerometerMa);
791         this.listener = listener;
792     }
793 
794     /**
795      * Constructor.
796      *
797      * @param sequences                     collection of sequences containing timestamped body
798      *                                      kinematics measurements.
799      * @param commonAxisUsed                indicates whether z-axis is
800      *                                      assumed to be common for
801      *                                      accelerometer and gyroscope.
802      * @param estimateGDependentCrossBiases true if G-dependent cross biases
803      *                                      will be estimated, false
804      *                                      otherwise.
805      * @param initialBias                   initial gyroscope bias to be used to find a
806      *                                      solution. This must be 3x1 and is expressed
807      *                                      in radians per second (rad/s).
808      * @param initialMg                     initial gyroscope scale factors and cross coupling
809      *                                      errors matrix. Must be 3x3.
810      * @param initialGg                     initial gyroscope G-dependent cross biases
811      *                                      introduced on the gyroscope by the specific forces
812      *                                      sensed by the accelerometer. Must be 3x3.
813      * @throws IllegalArgumentException if any of the provided values does
814      *                                  not have proper size.
815      */
816     public EasyGyroscopeCalibrator(
817             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
818             final boolean commonAxisUsed, final boolean estimateGDependentCrossBiases,
819             final Matrix initialBias, final Matrix initialMg, final Matrix initialGg) {
820         this(sequences, initialBias, initialMg, initialGg);
821         this.commonAxisUsed = commonAxisUsed;
822         this.estimateGDependentCrossBiases = estimateGDependentCrossBiases;
823     }
824 
825     /**
826      * Constructor.
827      *
828      * @param sequences                     collection of sequences containing timestamped body
829      *                                      kinematics measurements.
830      * @param commonAxisUsed                indicates whether z-axis is
831      *                                      assumed to be common for
832      *                                      accelerometer and gyroscope.
833      * @param estimateGDependentCrossBiases true if G-dependent cross biases
834      *                                      will be estimated, false
835      *                                      otherwise.
836      * @param initialBias                   initial gyroscope bias to be used to find a
837      *                                      solution. This must be 3x1 and is expressed
838      *                                      in radians per second (rad/s).
839      * @param initialMg                     initial gyroscope scale factors and cross coupling
840      *                                      errors matrix. Must be 3x3.
841      * @param initialGg                     initial gyroscope G-dependent cross biases
842      *                                      introduced on the gyroscope by the specific forces
843      *                                      sensed by the accelerometer. Must be 3x3.
844      * @param listener                      listener to handle events raised by this
845      *                                      calibrator.
846      * @throws IllegalArgumentException if any of the provided values does
847      *                                  not have proper size.
848      */
849     public EasyGyroscopeCalibrator(
850             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
851             final boolean commonAxisUsed, final boolean estimateGDependentCrossBiases, final Matrix initialBias,
852             final Matrix initialMg, final Matrix initialGg, final EasyGyroscopeCalibratorListener listener) {
853         this(sequences, commonAxisUsed, estimateGDependentCrossBiases, initialBias, initialMg, initialGg);
854         this.listener = listener;
855     }
856 
857     /**
858      * Constructor.
859      *
860      * @param sequences                     collection of sequences containing timestamped body
861      *                                      kinematics measurements.
862      * @param commonAxisUsed                indicates whether z-axis is
863      *                                      assumed to be common for
864      *                                      accelerometer and gyroscope.
865      * @param estimateGDependentCrossBiases true if G-dependent cross biases
866      *                                      will be estimated, false
867      *                                      otherwise.
868      * @param initialBias                   initial gyroscope bias to be used to find a
869      *                                      solution. This must have length 3 and is expressed
870      *                                      in radians per second (rad/s).
871      * @param initialMg                     initial gyroscope scale factors and cross coupling
872      *                                      errors matrix. Must be 3x3.
873      * @param initialGg                     initial gyroscope G-dependent cross biases
874      *                                      introduced on the gyroscope by the specific forces
875      *                                      sensed by the accelerometer. Must be 3x3.
876      * @throws IllegalArgumentException if any of the provided values does
877      *                                  not have proper size.
878      */
879     public EasyGyroscopeCalibrator(
880             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
881             final boolean commonAxisUsed, final boolean estimateGDependentCrossBiases, final double[] initialBias,
882             final Matrix initialMg, final Matrix initialGg) {
883         this(sequences, initialBias, initialMg, initialGg);
884         this.commonAxisUsed = commonAxisUsed;
885         this.estimateGDependentCrossBiases = estimateGDependentCrossBiases;
886     }
887 
888     /**
889      * Constructor.
890      *
891      * @param sequences                     collection of sequences containing timestamped body
892      *                                      kinematics measurements.
893      * @param commonAxisUsed                indicates whether z-axis is
894      *                                      assumed to be common for
895      *                                      accelerometer and gyroscope.
896      * @param estimateGDependentCrossBiases true if G-dependent cross biases
897      *                                      will be estimated, false
898      *                                      otherwise.
899      * @param initialBias                   initial gyroscope bias to be used to find a
900      *                                      solution. This must have length 3 and is expressed
901      *                                      in radians per second (rad/s).
902      * @param initialMg                     initial gyroscope scale factors and cross coupling
903      *                                      errors matrix. Must be 3x3.
904      * @param initialGg                     initial gyroscope G-dependent cross biases
905      *                                      introduced on the gyroscope by the specific forces
906      *                                      sensed by the accelerometer. Must be 3x3.
907      * @param listener                      listener to handle events raised by this
908      *                                      calibrator.
909      * @throws IllegalArgumentException if any of the provided values does
910      *                                  not have proper size.
911      */
912     public EasyGyroscopeCalibrator(
913             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
914             final boolean commonAxisUsed, final boolean estimateGDependentCrossBiases, final double[] initialBias,
915             final Matrix initialMg, final Matrix initialGg, final EasyGyroscopeCalibratorListener listener) {
916         this(sequences, commonAxisUsed, estimateGDependentCrossBiases, initialBias, initialMg, initialGg);
917         this.listener = listener;
918     }
919 
920     /**
921      * Constructor.
922      *
923      * @param sequences                     collection of sequences containing timestamped body
924      *                                      kinematics measurements.
925      * @param commonAxisUsed                indicates whether z-axis is
926      *                                      assumed to be common for
927      *                                      accelerometer and gyroscope.
928      * @param estimateGDependentCrossBiases true if G-dependent cross biases
929      *                                      will be estimated, false
930      *                                      otherwise.
931      * @param initialBias                   initial gyroscope bias to be used to find a
932      *                                      solution. This must have length 3 and is expressed
933      *                                      in radians per second (rad/s).
934      * @param initialMg                     initial gyroscope scale factors and cross coupling
935      *                                      errors matrix. Must be 3x3.
936      * @param initialGg                     initial gyroscope G-dependent cross biases
937      *                                      introduced on the gyroscope by the specific forces
938      *                                      sensed by the accelerometer. Must be 3x3.
939      * @param accelerometerBias             known accelerometer bias. This
940      *                                      must have length 3 and is
941      *                                      expressed in meters per squared
942      *                                      second (m/s^2).
943      * @param accelerometerMa               known accelerometer scale factors
944      *                                      and cross coupling matrix. Must
945      *                                      be 3x3.
946      * @throws IllegalArgumentException if any of the provided values does
947      *                                  not have proper size.
948      */
949     public EasyGyroscopeCalibrator(
950             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
951             final boolean commonAxisUsed, final boolean estimateGDependentCrossBiases, final double[] initialBias,
952             final Matrix initialMg, final Matrix initialGg, final double[] accelerometerBias,
953             final Matrix accelerometerMa) {
954         this(sequences, initialBias, initialMg, initialGg, accelerometerBias, accelerometerMa);
955         this.commonAxisUsed = commonAxisUsed;
956         this.estimateGDependentCrossBiases = estimateGDependentCrossBiases;
957     }
958 
959     /**
960      * Constructor.
961      *
962      * @param sequences                     collection of sequences containing timestamped body
963      *                                      kinematics measurements.
964      * @param commonAxisUsed                indicates whether z-axis is
965      *                                      assumed to be common for
966      *                                      accelerometer and gyroscope.
967      * @param estimateGDependentCrossBiases true if G-dependent cross biases
968      *                                      will be estimated, false
969      *                                      otherwise.
970      * @param initialBias                   initial gyroscope bias to be used to find a
971      *                                      solution. This must have length 3 and is expressed
972      *                                      in radians per second (rad/s).
973      * @param initialMg                     initial gyroscope scale factors and cross coupling
974      *                                      errors matrix. Must be 3x3.
975      * @param initialGg                     initial gyroscope G-dependent cross biases
976      *                                      introduced on the gyroscope by the specific forces
977      *                                      sensed by the accelerometer. Must be 3x3.
978      * @param accelerometerBias             known accelerometer bias. This
979      *                                      must have length 3 and is
980      *                                      expressed in meters per squared
981      *                                      second (m/s^2).
982      * @param accelerometerMa               known accelerometer scale factors
983      *                                      and cross coupling matrix. Must
984      *                                      be 3x3.
985      * @param listener                      listener to handle events raised by this
986      *                                      calibrator.
987      * @throws IllegalArgumentException if any of the provided values does
988      *                                  not have proper size.
989      */
990     public EasyGyroscopeCalibrator(
991             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
992             final boolean commonAxisUsed, final boolean estimateGDependentCrossBiases, final double[] initialBias,
993             final Matrix initialMg, final Matrix initialGg, final double[] accelerometerBias,
994             final Matrix accelerometerMa, final EasyGyroscopeCalibratorListener listener) {
995         this(sequences, commonAxisUsed, estimateGDependentCrossBiases, initialBias, initialMg, initialGg,
996                 accelerometerBias, accelerometerMa);
997         this.listener = listener;
998     }
999 
1000     /**
1001      * Constructor.
1002      *
1003      * @param sequences                     collection of sequences containing timestamped body
1004      *                                      kinematics measurements.
1005      * @param commonAxisUsed                indicates whether z-axis is
1006      *                                      assumed to be common for
1007      *                                      accelerometer and gyroscope.
1008      * @param estimateGDependentCrossBiases true if G-dependent cross biases
1009      *                                      will be estimated, false
1010      *                                      otherwise.
1011      * @param initialBias                   initial gyroscope bias to be used to find a
1012      *                                      solution. This must be 3x1 and is expressed
1013      *                                      in radians per second (rad/s).
1014      * @param initialMg                     initial gyroscope scale factors and cross coupling
1015      *                                      errors matrix. Must be 3x3.
1016      * @param initialGg                     initial gyroscope G-dependent cross biases
1017      *                                      introduced on the gyroscope by the specific forces
1018      *                                      sensed by the accelerometer. Must be 3x3.
1019      * @param accelerometerBias             known accelerometer bias. This
1020      *                                      must have length 3 and is
1021      *                                      expressed in meters per squared
1022      *                                      second (m/s^2).
1023      * @param accelerometerMa               known accelerometer scale factors
1024      *                                      and cross coupling matrix. Must
1025      *                                      be 3x3.
1026      * @throws IllegalArgumentException if any of the provided values does
1027      *                                  not have proper size.
1028      */
1029     public EasyGyroscopeCalibrator(
1030             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
1031             final boolean commonAxisUsed, final boolean estimateGDependentCrossBiases, final Matrix initialBias,
1032             final Matrix initialMg, final Matrix initialGg, final Matrix accelerometerBias,
1033             final Matrix accelerometerMa) {
1034         this(sequences, initialBias, initialMg, initialGg, accelerometerBias, accelerometerMa);
1035         this.commonAxisUsed = commonAxisUsed;
1036         this.estimateGDependentCrossBiases = estimateGDependentCrossBiases;
1037     }
1038 
1039     /**
1040      * Constructor.
1041      *
1042      * @param sequences                     collection of sequences containing timestamped body
1043      *                                      kinematics measurements.
1044      * @param commonAxisUsed                indicates whether z-axis is
1045      *                                      assumed to be common for
1046      *                                      accelerometer and gyroscope.
1047      * @param estimateGDependentCrossBiases true if G-dependent cross biases
1048      *                                      will be estimated, false
1049      *                                      otherwise.
1050      * @param initialBias                   initial gyroscope bias to be used to find a
1051      *                                      solution. This must be 3x1 and is expressed
1052      *                                      in radians per second (rad/s).
1053      * @param initialMg                     initial gyroscope scale factors and cross coupling
1054      *                                      errors matrix. Must be 3x3.
1055      * @param initialGg                     initial gyroscope G-dependent cross biases
1056      *                                      introduced on the gyroscope by the specific forces
1057      *                                      sensed by the accelerometer. Must be 3x3.
1058      * @param accelerometerBias             known accelerometer bias. This
1059      *                                      must have length 3 and is
1060      *                                      expressed in meters per squared
1061      *                                      second (m/s^2).
1062      * @param accelerometerMa               known accelerometer scale factors
1063      *                                      and cross coupling matrix. Must
1064      *                                      be 3x3.
1065      * @param listener                      listener to handle events raised by this
1066      *                                      calibrator.
1067      * @throws IllegalArgumentException if any of the provided values does
1068      *                                  not have proper size.
1069      */
1070     public EasyGyroscopeCalibrator(
1071             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences,
1072             final boolean commonAxisUsed, final boolean estimateGDependentCrossBiases, final Matrix initialBias,
1073             final Matrix initialMg, final Matrix initialGg, final Matrix accelerometerBias,
1074             final Matrix accelerometerMa, final EasyGyroscopeCalibratorListener listener) {
1075         this(sequences, commonAxisUsed, estimateGDependentCrossBiases, initialBias, initialMg, initialGg,
1076                 accelerometerBias, accelerometerMa);
1077         this.listener = listener;
1078     }
1079 
1080     /**
1081      * Gets known x-coordinate of accelerometer bias to be used to fix
1082      * measured specific force and find cross biases introduced by the
1083      * accelerometer.
1084      * This is expressed in meters per squared second (m/s^2).
1085      *
1086      * @return known x-coordinate of accelerometer bias.
1087      */
1088     @Override
1089     public double getAccelerometerBiasX() {
1090         return accelerometerBiasX;
1091     }
1092 
1093     /**
1094      * Sets known x-coordinate of accelerometer bias to be used to fix
1095      * measured specific force and find cross biases introduced by the
1096      * accelerometer.
1097      * This is expressed in meters per squared second (m/s^2).
1098      *
1099      * @param accelerometerBiasX known x-coordinate of accelerometer bias.
1100      * @throws LockedException if calibrator is currently running.
1101      */
1102     @Override
1103     public void setAccelerometerBiasX(final double accelerometerBiasX) throws LockedException {
1104         if (running) {
1105             throw new LockedException();
1106         }
1107         this.accelerometerBiasX = accelerometerBiasX;
1108     }
1109 
1110     /**
1111      * Gets known y-coordinate of accelerometer bias to be used to fix
1112      * measured specific force and find cross biases introduced by the
1113      * accelerometer.
1114      * This is expressed in meters per squared second (m/s^2).
1115      *
1116      * @return known y-coordinate of accelerometer bias.
1117      */
1118     @Override
1119     public double getAccelerometerBiasY() {
1120         return accelerometerBiasY;
1121     }
1122 
1123     /**
1124      * Sets known y-coordinate of accelerometer bias to be used to fix
1125      * measured specific force and find cross biases introduced by the
1126      * accelerometer.
1127      * This is expressed in meters per squared second (m/s^2).
1128      *
1129      * @param accelerometerBiasY known y-coordinate of accelerometer bias.
1130      * @throws LockedException if calibrator is currently running.
1131      */
1132     @Override
1133     public void setAccelerometerBiasY(final double accelerometerBiasY) throws LockedException {
1134         if (running) {
1135             throw new LockedException();
1136         }
1137         this.accelerometerBiasY = accelerometerBiasY;
1138     }
1139 
1140     /**
1141      * Gets known z-coordinate of accelerometer bias to be used to fix
1142      * measured specific force and find cross biases introduced by the
1143      * accelerometer.
1144      * This is expressed in meters per squared second (m/s^2).
1145      *
1146      * @return known z-coordinate of accelerometer bias.
1147      */
1148     @Override
1149     public double getAccelerometerBiasZ() {
1150         return accelerometerBiasZ;
1151     }
1152 
1153     /**
1154      * Sets known z-coordinate of accelerometer bias to be used to fix
1155      * measured specific force and find cross biases introduced by the
1156      * accelerometer.
1157      * This is expressed in meters per squared second (m/s^2).
1158      *
1159      * @param accelerometerBiasZ known z-coordinate of accelerometer bias.
1160      * @throws LockedException if calibrator is currently running.
1161      */
1162     @Override
1163     public void setAccelerometerBiasZ(final double accelerometerBiasZ) throws LockedException {
1164         if (running) {
1165             throw new LockedException();
1166         }
1167         this.accelerometerBiasZ = accelerometerBiasZ;
1168     }
1169 
1170     /**
1171      * Gets known x-coordinate of accelerometer bias to be used to fix
1172      * measured specific force and find cross biases introduced by the
1173      * accelerometer.
1174      *
1175      * @return known x-coordinate of accelerometer bias.
1176      */
1177     @Override
1178     public Acceleration getAccelerometerBiasXAsAcceleration() {
1179         return new Acceleration(accelerometerBiasX, AccelerationUnit.METERS_PER_SQUARED_SECOND);
1180     }
1181 
1182     /**
1183      * Gets known x-coordinate of accelerometer bias to be used to fix
1184      * measured specific force and find cross biases introduced by the
1185      * accelerometer.
1186      *
1187      * @param result instance where result data will be stored.
1188      */
1189     @Override
1190     public void getAccelerometerBiasXAsAcceleration(final Acceleration result) {
1191         result.setValue(accelerometerBiasX);
1192         result.setUnit(AccelerationUnit.METERS_PER_SQUARED_SECOND);
1193     }
1194 
1195     /**
1196      * Sets known x-coordinate of accelerometer bias to be used to fix
1197      * measured specific force and find cross biases introduced by the
1198      * accelerometer.
1199      *
1200      * @param accelerometerBiasX x-coordinate of accelerometer bias.
1201      * @throws LockedException if calibrator is currently running.
1202      */
1203     @Override
1204     public void setAccelerometerBiasX(final Acceleration accelerometerBiasX) throws LockedException {
1205         if (running) {
1206             throw new LockedException();
1207         }
1208         this.accelerometerBiasX = convertAcceleration(accelerometerBiasX);
1209     }
1210 
1211     /**
1212      * Gets known y-coordinate of accelerometer bias to be used to fix
1213      * measured specific force and find cross biases introduced by the
1214      * accelerometer.
1215      *
1216      * @return known y-coordinate of accelerometer bias.
1217      */
1218     @Override
1219     public Acceleration getAccelerometerBiasYAsAcceleration() {
1220         return new Acceleration(accelerometerBiasY, AccelerationUnit.METERS_PER_SQUARED_SECOND);
1221     }
1222 
1223     /**
1224      * Gets known y-coordinate of accelerometer bias to be used to fix
1225      * measured specific force and find cross biases introduced by the
1226      * accelerometer.
1227      *
1228      * @param result instance where result data will be stored.
1229      */
1230     @Override
1231     public void getAccelerometerBiasYAsAcceleration(final Acceleration result) {
1232         result.setValue(accelerometerBiasY);
1233         result.setUnit(AccelerationUnit.METERS_PER_SQUARED_SECOND);
1234     }
1235 
1236     /**
1237      * Sets known y-coordinate of accelerometer bias to be used to fix
1238      * measured specific force and find cross biases introduced by the
1239      * accelerometer.
1240      *
1241      * @param accelerometerBiasY y-coordinate of accelerometer bias.
1242      * @throws LockedException if calibrator is currently running.
1243      */
1244     @Override
1245     public void setAccelerometerBiasY(final Acceleration accelerometerBiasY) throws LockedException {
1246         if (running) {
1247             throw new LockedException();
1248         }
1249         this.accelerometerBiasY = convertAcceleration(accelerometerBiasY);
1250     }
1251 
1252     /**
1253      * Gets known z-coordinate of accelerometer bias to be used to fix
1254      * measured specific force and find cross biases introduced by the
1255      * accelerometer.
1256      *
1257      * @return known z-coordinate of accelerometer bias.
1258      */
1259     @Override
1260     public Acceleration getAccelerometerBiasZAsAcceleration() {
1261         return new Acceleration(accelerometerBiasZ, AccelerationUnit.METERS_PER_SQUARED_SECOND);
1262     }
1263 
1264     /**
1265      * Gets known z-coordinate of accelerometer bias to be used to fix
1266      * measured specific force and find cross biases introduced by the
1267      * accelerometer.
1268      *
1269      * @param result instance where result data will be stored.
1270      */
1271     @Override
1272     public void getAccelerometerBiasZAsAcceleration(final Acceleration result) {
1273         result.setValue(accelerometerBiasZ);
1274         result.setUnit(AccelerationUnit.METERS_PER_SQUARED_SECOND);
1275     }
1276 
1277     /**
1278      * Sets known z-coordinate of accelerometer bias to be used to fix
1279      * measured specific force and find cross biases introduced by the
1280      * accelerometer.
1281      *
1282      * @param accelerometerBiasZ z-coordinate of accelerometer bias.
1283      * @throws LockedException if calibrator is currently running.
1284      */
1285     @Override
1286     public void setAccelerometerBiasZ(final Acceleration accelerometerBiasZ) throws LockedException {
1287         if (running) {
1288             throw new LockedException();
1289         }
1290         this.accelerometerBiasZ = convertAcceleration(accelerometerBiasZ);
1291     }
1292 
1293     /**
1294      * Sets known accelerometer bias to be used to fix measured specific
1295      * force and find cross biases introduced by the accelerometer.
1296      * This is expressed in meters per squared second (m/s^2).
1297      *
1298      * @param accelerometerBiasX x-coordinate of accelerometer bias.
1299      * @param accelerometerBiasY y-coordinate of accelerometer bias.
1300      * @param accelerometerBiasZ z-coordinate of accelerometer bias.
1301      * @throws LockedException if calibrator is currently running.
1302      */
1303     @Override
1304     public void setAccelerometerBias(
1305             final double accelerometerBiasX, final double accelerometerBiasY, final double accelerometerBiasZ)
1306             throws LockedException {
1307         if (running) {
1308             throw new LockedException();
1309         }
1310 
1311         this.accelerometerBiasX = accelerometerBiasX;
1312         this.accelerometerBiasY = accelerometerBiasY;
1313         this.accelerometerBiasZ = accelerometerBiasZ;
1314     }
1315 
1316     /**
1317      * Sets known accelerometer bias to be used to fix measured specific
1318      * force and find cross biases introduced by the accelerometer.
1319      *
1320      * @param accelerometerBiasX x-coordinate of accelerometer bias.
1321      * @param accelerometerBiasY y-coordinate of accelerometer bias.
1322      * @param accelerometerBiasZ z-coordinate of accelerometer bias.
1323      * @throws LockedException if calibrator is currently running.
1324      */
1325     @Override
1326     public void setAccelerometerBias(
1327             final Acceleration accelerometerBiasX,
1328             final Acceleration accelerometerBiasY,
1329             final Acceleration accelerometerBiasZ) throws LockedException {
1330         if (running) {
1331             throw new LockedException();
1332         }
1333 
1334         this.accelerometerBiasX = convertAcceleration(accelerometerBiasX);
1335         this.accelerometerBiasY = convertAcceleration(accelerometerBiasY);
1336         this.accelerometerBiasZ = convertAcceleration(accelerometerBiasZ);
1337     }
1338 
1339     /**
1340      * Gets known accelerometer bias to be used to fix measured specific
1341      * force and find cross biases introduced by the accelerometer.
1342      * This is expressed in meters per squared second (m/s^2).
1343      *
1344      * @return known accelerometer bias.
1345      */
1346     @Override
1347     public double[] getAccelerometerBias() {
1348         final var result = new double[BodyKinematics.COMPONENTS];
1349         getAccelerometerBias(result);
1350         return result;
1351     }
1352 
1353     /**
1354      * Gets known accelerometer bias to be used to fix measured specific
1355      * force and find cross biases introduced by the accelerometer.
1356      * This is expressed in meters per squared second (m/s^2).
1357      *
1358      * @param result instance where result data will be copied to.
1359      * @throws IllegalArgumentException if provided array does not have
1360      *                                  length 3.
1361      */
1362     @Override
1363     public void getAccelerometerBias(final double[] result) {
1364         if (result.length != BodyKinematics.COMPONENTS) {
1365             throw new IllegalArgumentException();
1366         }
1367 
1368         result[0] = accelerometerBiasX;
1369         result[1] = accelerometerBiasY;
1370         result[2] = accelerometerBiasZ;
1371     }
1372 
1373     /**
1374      * Sets known accelerometer bias to be used to fix measured specific
1375      * force and find cross biases introduced by the accelerometer.
1376      * This is expressed in meters per squared second (m/s^2).
1377      *
1378      * @param accelerometerBias known accelerometer bias.
1379      * @throws LockedException          if calibrator is currently running.
1380      * @throws IllegalArgumentException if provided array does not have
1381      *                                  length 3.
1382      */
1383     @Override
1384     public void setAccelerometerBias(final double[] accelerometerBias) throws LockedException {
1385         if (running) {
1386             throw new LockedException();
1387         }
1388 
1389         if (accelerometerBias.length != BodyKinematics.COMPONENTS) {
1390             throw new IllegalArgumentException();
1391         }
1392 
1393         accelerometerBiasX = accelerometerBias[0];
1394         accelerometerBiasY = accelerometerBias[1];
1395         accelerometerBiasZ = accelerometerBias[2];
1396     }
1397 
1398     /**
1399      * Gets known accelerometer bias to be used to fix measured specific
1400      * force and find cross biases introduced by the accelerometer.
1401      * This is expressed in meters per squared second (m/s^2).
1402      *
1403      * @return known accelerometer bias.
1404      */
1405     @Override
1406     public Matrix getAccelerometerBiasAsMatrix() {
1407         Matrix result;
1408         try {
1409             result = new Matrix(BodyKinematics.COMPONENTS, 1);
1410             getAccelerometerBiasAsMatrix(result);
1411         } catch (final WrongSizeException ignore) {
1412             // never happens
1413             result = null;
1414         }
1415         return result;
1416     }
1417 
1418     /**
1419      * Gets known accelerometer bias to be used to fix measured specific
1420      * force and find cross biases introduced by the accelerometer.
1421      * This is expressed in meters per squared second (m/s^2).
1422      *
1423      * @param result instance where result data will be copied to.
1424      * @throws IllegalArgumentException if provided matrix is not 3x1.
1425      */
1426     @Override
1427     public void getAccelerometerBiasAsMatrix(final Matrix result) {
1428         if (result.getRows() != BodyKinematics.COMPONENTS || result.getColumns() != 1) {
1429             throw new IllegalArgumentException();
1430         }
1431         result.setElementAtIndex(0, accelerometerBiasX);
1432         result.setElementAtIndex(1, accelerometerBiasY);
1433         result.setElementAtIndex(2, accelerometerBiasZ);
1434     }
1435 
1436     /**
1437      * Sets known accelerometer bias to be used to fix measured specific
1438      * force and find cross biases introduced by the accelerometer.
1439      * This is expressed in meters per squared second (m/s^2).
1440      *
1441      * @param accelerometerBias known accelerometer bias. Must be 3x1.
1442      * @throws LockedException          if calibrator is currently running.
1443      * @throws IllegalArgumentException if provided matrix is not 3x1.
1444      */
1445     @Override
1446     public void setAccelerometerBias(final Matrix accelerometerBias) throws LockedException {
1447         if (running) {
1448             throw new LockedException();
1449         }
1450         if (accelerometerBias.getRows() != BodyKinematics.COMPONENTS || accelerometerBias.getColumns() != 1) {
1451             throw new IllegalArgumentException();
1452         }
1453 
1454         accelerometerBiasX = accelerometerBias.getElementAtIndex(0);
1455         accelerometerBiasY = accelerometerBias.getElementAtIndex(1);
1456         accelerometerBiasZ = accelerometerBias.getElementAtIndex(2);
1457     }
1458 
1459     /**
1460      * Gets known accelerometer x scaling factor to be used to fix measured
1461      * specific force and find cross biases introduced by the accelerometer.
1462      *
1463      * @return known accelerometer x scaling factor.
1464      */
1465     @Override
1466     public double getAccelerometerSx() {
1467         return accelerometerSx;
1468     }
1469 
1470     /**
1471      * Sets known accelerometer x scaling factor to be used to fix measured
1472      * specific force and find cross biases introduced by the accelerometer.
1473      *
1474      * @param accelerometerSx known accelerometer x scaling factor.
1475      * @throws LockedException if calibrator is currently running.
1476      */
1477     @Override
1478     public void setAccelerometerSx(final double accelerometerSx) throws LockedException {
1479         if (running) {
1480             throw new LockedException();
1481         }
1482         this.accelerometerSx = accelerometerSx;
1483     }
1484 
1485     /**
1486      * Gets known accelerometer y scaling factor to be used to fix measured
1487      * specific force and find cross biases introduced by the accelerometer.
1488      *
1489      * @return known accelerometer y scaling factor.
1490      */
1491     @Override
1492     public double getAccelerometerSy() {
1493         return accelerometerSy;
1494     }
1495 
1496     /**
1497      * Sets known accelerometer y scaling factor to be used to fix measured
1498      * specific force and find cross biases introduced by the accelerometer.
1499      *
1500      * @param accelerometerSy known accelerometer y scaling factor.
1501      * @throws LockedException if calibrator is currently running.
1502      */
1503     @Override
1504     public void setAccelerometerSy(final double accelerometerSy) throws LockedException {
1505         if (running) {
1506             throw new LockedException();
1507         }
1508         this.accelerometerSy = accelerometerSy;
1509     }
1510 
1511     /**
1512      * Gets known accelerometer z scaling factor to be used to fix measured
1513      * specific force and find cross biases introduced by the accelerometer.
1514      *
1515      * @return known accelerometer z scaling factor.
1516      */
1517     @Override
1518     public double getAccelerometerSz() {
1519         return accelerometerSz;
1520     }
1521 
1522     /**
1523      * Sets known accelerometer z scaling factor to be used to fix measured
1524      * specific force and find cross biases introduced by the accelerometer.
1525      *
1526      * @param accelerometerSz known accelerometer z scaling factor.
1527      * @throws LockedException if calibrator is currently running.
1528      */
1529     @Override
1530     public void setAccelerometerSz(final double accelerometerSz) throws LockedException {
1531         if (running) {
1532             throw new LockedException();
1533         }
1534         this.accelerometerSz = accelerometerSz;
1535     }
1536 
1537     /**
1538      * Gets known accelerometer x-y cross coupling error to be used to fix
1539      * measured specific force and find cross biases introduced by the
1540      * accelerometer.
1541      *
1542      * @return known accelerometer x-y cross coupling error.
1543      */
1544     @Override
1545     public double getAccelerometerMxy() {
1546         return accelerometerMxy;
1547     }
1548 
1549     /**
1550      * Sets known accelerometer x-y cross coupling error to be used to fix
1551      * measured specific force and find cross biases introduced by the
1552      * accelerometer.
1553      *
1554      * @param accelerometerMxy known accelerometer x-y cross coupling error.
1555      * @throws LockedException if calibrator is currently running.
1556      */
1557     @Override
1558     public void setAccelerometerMxy(final double accelerometerMxy) throws LockedException {
1559         if (running) {
1560             throw new LockedException();
1561         }
1562         this.accelerometerMxy = accelerometerMxy;
1563     }
1564 
1565     /**
1566      * Gets known accelerometer x-z cross coupling error to be used to fix
1567      * measured specific force and find cross biases introduced by the
1568      * accelerometer.
1569      *
1570      * @return known accelerometer x-z cross coupling error.
1571      */
1572     @Override
1573     public double getAccelerometerMxz() {
1574         return accelerometerMxz;
1575     }
1576 
1577     /**
1578      * Sets known accelerometer x-z cross coupling error to be used to fix
1579      * measured specific force and find cross biases introduced by the
1580      * accelerometer.
1581      *
1582      * @param accelerometerMxz known accelerometer x-z cross coupling error.
1583      * @throws LockedException if calibrator is currently running.
1584      */
1585     @Override
1586     public void setAccelerometerMxz(final double accelerometerMxz) throws LockedException {
1587         if (running) {
1588             throw new LockedException();
1589         }
1590         this.accelerometerMxz = accelerometerMxz;
1591     }
1592 
1593     /**
1594      * Gets known accelerometer y-x cross coupling error to be used to fix
1595      * measured specific force and find cross biases introduced by the
1596      * accelerometer.
1597      *
1598      * @return known accelerometer y-x cross coupling error.
1599      */
1600     @Override
1601     public double getAccelerometerMyx() {
1602         return accelerometerMyx;
1603     }
1604 
1605     /**
1606      * Sets known accelerometer y-x cross coupling error to be used to fix
1607      * measured specific force and find cross biases introduced by the
1608      * accelerometer.
1609      *
1610      * @param accelerometerMyx known accelerometer y-x cross coupling
1611      *                         error.
1612      * @throws LockedException if calibrator is currently running.
1613      */
1614     @Override
1615     public void setAccelerometerMyx(final double accelerometerMyx) throws LockedException {
1616         if (running) {
1617             throw new LockedException();
1618         }
1619         this.accelerometerMyx = accelerometerMyx;
1620     }
1621 
1622     /**
1623      * Gets known accelerometer y-z cross coupling error to be used to fix
1624      * measured specific force and find cross biases introduced by the
1625      * accelerometer.
1626      *
1627      * @return known accelerometer y-z cross coupling error.
1628      */
1629     @Override
1630     public double getAccelerometerMyz() {
1631         return accelerometerMyz;
1632     }
1633 
1634     /**
1635      * Sets known accelerometer y-z cross coupling error to be used to fix
1636      * measured specific force and find cross biases introduced by the
1637      * accelerometer.
1638      *
1639      * @param accelerometerMyz known accelerometer y-z cross coupling
1640      *                         error.
1641      * @throws LockedException if calibrator is currently running.
1642      */
1643     @Override
1644     public void setAccelerometerMyz(final double accelerometerMyz) throws LockedException {
1645         if (running) {
1646             throw new LockedException();
1647         }
1648         this.accelerometerMyz = accelerometerMyz;
1649     }
1650 
1651     /**
1652      * Gets known accelerometer z-x cross coupling error to be used to fix
1653      * measured specific force and find cross biases introduced by the
1654      * accelerometer.
1655      *
1656      * @return known accelerometer z-x cross coupling error.
1657      */
1658     @Override
1659     public double getAccelerometerMzx() {
1660         return accelerometerMzx;
1661     }
1662 
1663     /**
1664      * Sets known accelerometer z-x cross coupling error to be used to fix
1665      * measured specific force and find cross biases introduced by the
1666      * accelerometer.
1667      *
1668      * @param accelerometerMzx known accelerometer z-x cross coupling
1669      *                         error.
1670      * @throws LockedException if calibrator is currently running.
1671      */
1672     @Override
1673     public void setAccelerometerMzx(final double accelerometerMzx) throws LockedException {
1674         if (running) {
1675             throw new LockedException();
1676         }
1677         this.accelerometerMzx = accelerometerMzx;
1678     }
1679 
1680     /**
1681      * Gets known accelerometer z-y cross coupling error to be used to fix
1682      * measured specific force and find cross biases introduced by the
1683      * accelerometer.
1684      *
1685      * @return known accelerometer z-y cross coupling error.
1686      */
1687     @Override
1688     public double getAccelerometerMzy() {
1689         return accelerometerMzy;
1690     }
1691 
1692     /**
1693      * Sets known accelerometer z-y cross coupling error to be used to fix
1694      * measured specific force and find cross biases introduced by the
1695      * accelerometer.
1696      *
1697      * @param accelerometerMzy known accelerometer z-y cross coupling
1698      *                         error.
1699      * @throws LockedException if calibrator is currently running.
1700      */
1701     @Override
1702     public void setAccelerometerMzy(final double accelerometerMzy) throws LockedException {
1703         if (running) {
1704             throw new LockedException();
1705         }
1706         this.accelerometerMzy = accelerometerMzy;
1707     }
1708 
1709     /**
1710      * Sets known accelerometer scaling factors to be used to fix measured
1711      * specific force and find cross biases introduced by the
1712      * accelerometer.
1713      *
1714      * @param accelerometerSx known accelerometer x scaling factor.
1715      * @param accelerometerSy known accelerometer y scaling factor.
1716      * @param accelerometerSz known accelerometer z scaling factor.
1717      * @throws LockedException if calibrator is currently running.
1718      */
1719     @Override
1720     public void setAccelerometerScalingFactors(
1721             final double accelerometerSx, final double accelerometerSy, final double accelerometerSz)
1722             throws LockedException {
1723         if (running) {
1724             throw new LockedException();
1725         }
1726         this.accelerometerSx = accelerometerSx;
1727         this.accelerometerSy = accelerometerSy;
1728         this.accelerometerSz = accelerometerSz;
1729     }
1730 
1731     /**
1732      * Sets known accelerometer cross coupling errors to be used to fix
1733      * measured specific force and find cross biases introduced by the
1734      * accelerometer.
1735      *
1736      * @param accelerometerMxy known accelerometer x-y cross coupling
1737      *                         error.
1738      * @param accelerometerMxz known accelerometer x-z cross coupling
1739      *                         error.
1740      * @param accelerometerMyx known accelerometer y-x cross coupling
1741      *                         error.
1742      * @param accelerometerMyz known accelerometer y-z cross coupling
1743      *                         error.
1744      * @param accelerometerMzx known accelerometer z-x cross coupling
1745      *                         error.
1746      * @param accelerometerMzy known accelerometer z-y cross coupling
1747      *                         error.
1748      * @throws LockedException if calibrator is currently running.
1749      */
1750     @Override
1751     public void setAccelerometerCrossCouplingErrors(
1752             final double accelerometerMxy, final double accelerometerMxz,
1753             final double accelerometerMyx, final double accelerometerMyz,
1754             final double accelerometerMzx, final double accelerometerMzy) throws LockedException {
1755         if (running) {
1756             throw new LockedException();
1757         }
1758         this.accelerometerMxy = accelerometerMxy;
1759         this.accelerometerMxz = accelerometerMxz;
1760         this.accelerometerMyx = accelerometerMyx;
1761         this.accelerometerMyz = accelerometerMyz;
1762         this.accelerometerMzx = accelerometerMzx;
1763         this.accelerometerMzy = accelerometerMzy;
1764     }
1765 
1766     /**
1767      * Sets known accelerometer scaling factors and cross coupling errors
1768      * to be used to fix measured specific force and find cross biases
1769      * introduced by the accelerometer.
1770      *
1771      * @param accelerometerSx  known accelerometer x scaling factor.
1772      * @param accelerometerSy  known accelerometer y scaling factor.
1773      * @param accelerometerSz  known accelerometer z scaling factor.
1774      * @param accelerometerMxy known accelerometer x-y cross coupling
1775      *                         error.
1776      * @param accelerometerMxz known accelerometer x-z cross coupling
1777      *                         error.
1778      * @param accelerometerMyx known accelerometer y-x cross coupling
1779      *                         error.
1780      * @param accelerometerMyz known accelerometer y-z cross coupling
1781      *                         error.
1782      * @param accelerometerMzx known accelerometer z-x cross coupling
1783      *                         error.
1784      * @param accelerometerMzy known accelerometer z-y cross coupling
1785      *                         error.
1786      * @throws LockedException if calibrator is currently running.
1787      */
1788     @Override
1789     public void setAccelerometerScalingFactorsAndCrossCouplingErrors(
1790             final double accelerometerSx, final double accelerometerSy,
1791             final double accelerometerSz, final double accelerometerMxy,
1792             final double accelerometerMxz, final double accelerometerMyx,
1793             final double accelerometerMyz, final double accelerometerMzx,
1794             final double accelerometerMzy) throws LockedException {
1795         if (running) {
1796             throw new LockedException();
1797         }
1798         setAccelerometerScalingFactors(accelerometerSx, accelerometerSy, accelerometerSz);
1799         setAccelerometerCrossCouplingErrors(accelerometerMxy, accelerometerMxz, accelerometerMyx,
1800                 accelerometerMyz, accelerometerMzx, accelerometerMzy);
1801     }
1802 
1803     /**
1804      * Gets known accelerometer scale factors and cross coupling
1805      * errors matrix.
1806      *
1807      * @return known accelerometer scale factors and cross coupling
1808      * errors matrix.
1809      */
1810     @Override
1811     public Matrix getAccelerometerMa() {
1812         Matrix result;
1813         try {
1814             result = new Matrix(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
1815             getAccelerometerMa(result);
1816         } catch (final WrongSizeException ignore) {
1817             // never happens
1818             result = null;
1819         }
1820         return result;
1821     }
1822 
1823     /**
1824      * Gets known accelerometer scale factors and cross coupling
1825      * errors matrix.
1826      *
1827      * @param result instance where data will be stored.
1828      * @throws IllegalArgumentException if provided matrix is not 3x3.
1829      */
1830     @Override
1831     public void getAccelerometerMa(final Matrix result) {
1832         if (result.getRows() != BodyKinematics.COMPONENTS || result.getColumns() != BodyKinematics.COMPONENTS) {
1833             throw new IllegalArgumentException();
1834         }
1835         result.setElementAtIndex(0, accelerometerSx);
1836         result.setElementAtIndex(1, accelerometerMyx);
1837         result.setElementAtIndex(2, accelerometerMzx);
1838 
1839         result.setElementAtIndex(3, accelerometerMxy);
1840         result.setElementAtIndex(4, accelerometerSy);
1841         result.setElementAtIndex(5, accelerometerMzy);
1842 
1843         result.setElementAtIndex(6, accelerometerMxz);
1844         result.setElementAtIndex(7, accelerometerMyz);
1845         result.setElementAtIndex(8, accelerometerSz);
1846     }
1847 
1848     /**
1849      * Sets known accelerometer scale factors and cross coupling
1850      * errors matrix.
1851      *
1852      * @param accelerometerMa known accelerometer scale factors and
1853      *                        cross coupling errors matrix. Must be 3x3.
1854      * @throws LockedException          if calibrator is currently running.
1855      * @throws IllegalArgumentException if provided matrix is not 3x3.
1856      */
1857     @Override
1858     public void setAccelerometerMa(final Matrix accelerometerMa) throws LockedException {
1859         if (running) {
1860             throw new LockedException();
1861         }
1862         if (accelerometerMa.getRows() != BodyKinematics.COMPONENTS
1863                 || accelerometerMa.getColumns() != BodyKinematics.COMPONENTS) {
1864             throw new IllegalArgumentException();
1865         }
1866 
1867         accelerometerSx = accelerometerMa.getElementAtIndex(0);
1868         accelerometerMyx = accelerometerMa.getElementAtIndex(1);
1869         accelerometerMzx = accelerometerMa.getElementAtIndex(2);
1870 
1871         accelerometerMxy = accelerometerMa.getElementAtIndex(3);
1872         accelerometerSy = accelerometerMa.getElementAtIndex(4);
1873         accelerometerMzy = accelerometerMa.getElementAtIndex(5);
1874 
1875         accelerometerMxz = accelerometerMa.getElementAtIndex(6);
1876         accelerometerMyz = accelerometerMa.getElementAtIndex(7);
1877         accelerometerSz = accelerometerMa.getElementAtIndex(8);
1878     }
1879 
1880     /**
1881      * Gets initial x-coordinate of gyroscope bias to be used to find
1882      * a solution.
1883      * This is expressed in radians per second (rad/s).
1884      *
1885      * @return initial x-coordinate of gyroscope bias.
1886      */
1887     public double getInitialBiasX() {
1888         return initialBiasX;
1889     }
1890 
1891     /**
1892      * Sets initial x-coordinate of gyroscope bias to be used to find
1893      * a solution.
1894      * This is expressed in radians per second (rad/s).
1895      *
1896      * @param initialBiasX initial x-coordinate of gyroscope bias.
1897      * @throws LockedException if calibrator is currently running.
1898      */
1899     public void setInitialBiasX(final double initialBiasX) throws LockedException {
1900         if (running) {
1901             throw new LockedException();
1902         }
1903         this.initialBiasX = initialBiasX;
1904     }
1905 
1906     /**
1907      * Gets initial y-coordinate of gyroscope bias to be used to find
1908      * a solution.
1909      * This is expressed in radians per second (rad/s).
1910      *
1911      * @return initial y-coordinate of gyroscope bias.
1912      */
1913     public double getInitialBiasY() {
1914         return initialBiasY;
1915     }
1916 
1917     /**
1918      * Sets initial y-coordinate of gyroscope bias to be used to find
1919      * a solution.
1920      * This is expressed in radians per second (rad/s).
1921      *
1922      * @param initialBiasY initial y-coordinate of gyroscope bias.
1923      * @throws LockedException if calibrator is currently running.
1924      */
1925     public void setInitialBiasY(final double initialBiasY) throws LockedException {
1926         if (running) {
1927             throw new LockedException();
1928         }
1929         this.initialBiasY = initialBiasY;
1930     }
1931 
1932     /**
1933      * Gets initial z-coordinate of gyroscope bias ot be used to find
1934      * a solution.
1935      * This is expressed in radians per second (rad/s).
1936      *
1937      * @return initial z-coordinate of gyroscope bias.
1938      */
1939     public double getInitialBiasZ() {
1940         return initialBiasZ;
1941     }
1942 
1943     /**
1944      * Sets initial z-coordinate of gyroscope bias to be used to find
1945      * a solution.
1946      * This is expressed in radians per second (rad/s).
1947      *
1948      * @param initialBiasZ initial z-coordinate of gyroscope bias.
1949      * @throws LockedException if calibrator is currently running.
1950      */
1951     public void setInitialBiasZ(final double initialBiasZ) throws LockedException {
1952         if (running) {
1953             throw new LockedException();
1954         }
1955         this.initialBiasZ = initialBiasZ;
1956     }
1957 
1958     /**
1959      * Gets initial x-coordinate of gyroscope bias to be used to find a
1960      * solution.
1961      *
1962      * @return initial x-coordinate of gyroscope bias.
1963      */
1964     public AngularSpeed getInitialBiasAngularSpeedX() {
1965         return new AngularSpeed(initialBiasX, AngularSpeedUnit.RADIANS_PER_SECOND);
1966     }
1967 
1968     /**
1969      * Gets initial x-coordinate of gyroscope bias to be used to find a
1970      * solution.
1971      *
1972      * @param result instance where result data will be stored.
1973      */
1974     public void getInitialBiasAngularSpeedX(final AngularSpeed result) {
1975         result.setValue(initialBiasX);
1976         result.setUnit(AngularSpeedUnit.RADIANS_PER_SECOND);
1977     }
1978 
1979     /**
1980      * Sets initial x-coordinate of gyroscope bias to be used to find a
1981      * solution.
1982      *
1983      * @param initialBiasX initial x-coordinate of gyroscope bias.
1984      * @throws LockedException if calibrator is currently running.
1985      */
1986     public void setInitialBiasX(final AngularSpeed initialBiasX) throws LockedException {
1987         if (running) {
1988             throw new LockedException();
1989         }
1990         this.initialBiasX = convertAngularSpeed(initialBiasX);
1991     }
1992 
1993     /**
1994      * Gets initial y-coordinate of gyroscope bias to be used to find a
1995      * solution.
1996      *
1997      * @return initial y-coordinate of gyroscope bias.
1998      */
1999     public AngularSpeed getInitialBiasAngularSpeedY() {
2000         return new AngularSpeed(initialBiasY, AngularSpeedUnit.RADIANS_PER_SECOND);
2001     }
2002 
2003     /**
2004      * Gets initial y-coordinate of gyroscope bias to be used to find a
2005      * solution.
2006      *
2007      * @param result instance where result data will be stored.
2008      */
2009     public void getInitialBiasAngularSpeedY(final AngularSpeed result) {
2010         result.setValue(initialBiasY);
2011         result.setUnit(AngularSpeedUnit.RADIANS_PER_SECOND);
2012     }
2013 
2014     /**
2015      * Sets initial y-coordinate of gyroscope bias to be used to find a
2016      * solution.
2017      *
2018      * @param initialBiasY initial y-coordinate of gyroscope bias.
2019      * @throws LockedException if calibrator is currently running.
2020      */
2021     public void setInitialBiasY(final AngularSpeed initialBiasY) throws LockedException {
2022         if (running) {
2023             throw new LockedException();
2024         }
2025         this.initialBiasY = convertAngularSpeed(initialBiasY);
2026     }
2027 
2028     /**
2029      * Gets initial z-coordinate of gyroscope bias to be used to find a
2030      * solution.
2031      *
2032      * @return initial z-coordinate of gyroscope bias.
2033      */
2034     public AngularSpeed getInitialBiasAngularSpeedZ() {
2035         return new AngularSpeed(initialBiasZ, AngularSpeedUnit.RADIANS_PER_SECOND);
2036     }
2037 
2038     /**
2039      * Gets initial z-coordinate of gyroscope bias to be used to find a
2040      * solution.
2041      *
2042      * @param result instance where result data will be stored.
2043      */
2044     public void getInitialBiasAngularSpeedZ(final AngularSpeed result) {
2045         result.setValue(initialBiasZ);
2046         result.setUnit(AngularSpeedUnit.RADIANS_PER_SECOND);
2047     }
2048 
2049     /**
2050      * Sets initial z-coordinate of gyroscope bias to be used to find a
2051      * solution.
2052      *
2053      * @param initialBiasZ initial z-coordinate of gyroscope bias.
2054      * @throws LockedException if calibrator is currently running.
2055      */
2056     public void setInitialBiasZ(final AngularSpeed initialBiasZ) throws LockedException {
2057         if (running) {
2058             throw new LockedException();
2059         }
2060         this.initialBiasZ = convertAngularSpeed(initialBiasZ);
2061     }
2062 
2063     /**
2064      * Sets initial bias coordinates of gyroscope used to find a solution
2065      * expressed in radians per second (rad/s).
2066      *
2067      * @param initialBiasX initial x-coordinate of gyroscope bias.
2068      * @param initialBiasY initial y-coordinate of gyroscope bias.
2069      * @param initialBiasZ initial z-coordinate of gyroscope bias.
2070      * @throws LockedException if calibrator is currently running.
2071      */
2072     public void setInitialBias(
2073             final double initialBiasX, final double initialBiasY, final double initialBiasZ) throws LockedException {
2074         if (running) {
2075             throw new LockedException();
2076         }
2077         this.initialBiasX = initialBiasX;
2078         this.initialBiasY = initialBiasY;
2079         this.initialBiasZ = initialBiasZ;
2080     }
2081 
2082     /**
2083      * Sets initial bias coordinates of gyroscope used to find a solution.
2084      *
2085      * @param initialBiasX initial x-coordinate of gyroscope bias.
2086      * @param initialBiasY initial y-coordinate of gyroscope bias.
2087      * @param initialBiasZ initial z-coordinate of gyroscope bias.
2088      * @throws LockedException if calibrator is currently running.
2089      */
2090     public void setInitialBias(
2091             final AngularSpeed initialBiasX, final AngularSpeed initialBiasY, final AngularSpeed initialBiasZ)
2092             throws LockedException {
2093         if (running) {
2094             throw new LockedException();
2095         }
2096         this.initialBiasX = convertAngularSpeed(initialBiasX);
2097         this.initialBiasY = convertAngularSpeed(initialBiasY);
2098         this.initialBiasZ = convertAngularSpeed(initialBiasZ);
2099     }
2100 
2101     /**
2102      * Gets initial x scaling factor of gyroscope.
2103      *
2104      * @return initial x scaling factor of gyroscope.
2105      */
2106     @Override
2107     public double getInitialSx() {
2108         return initialSx;
2109     }
2110 
2111     /**
2112      * Sets initial x scaling factor of gyroscope.
2113      *
2114      * @param initialSx initial x scaling factor of gyroscope.
2115      * @throws LockedException if calibrator is currently running.
2116      */
2117     @Override
2118     public void setInitialSx(final double initialSx) throws LockedException {
2119         if (running) {
2120             throw new LockedException();
2121         }
2122         this.initialSx = initialSx;
2123     }
2124 
2125     /**
2126      * Gets initial y scaling factor of gyroscope.
2127      *
2128      * @return initial y scaling factor of gyroscope.
2129      */
2130     @Override
2131     public double getInitialSy() {
2132         return initialSy;
2133     }
2134 
2135     /**
2136      * Sets initial y scaling factor of gyroscope.
2137      *
2138      * @param initialSy initial y scaling factor of gyroscope.
2139      * @throws LockedException if calibrator is currently running.
2140      */
2141     @Override
2142     public void setInitialSy(final double initialSy) throws LockedException {
2143         if (running) {
2144             throw new LockedException();
2145         }
2146         this.initialSy = initialSy;
2147     }
2148 
2149     /**
2150      * Gets initial z scaling factor of gyroscope.
2151      *
2152      * @return initial z scaling factor of gyroscope.
2153      */
2154     @Override
2155     public double getInitialSz() {
2156         return initialSz;
2157     }
2158 
2159     /**
2160      * Sets initial z scaling factor of gyroscope.
2161      *
2162      * @param initialSz initial z scaling factor of gyroscope.
2163      * @throws LockedException if calibrator is currently running.
2164      */
2165     @Override
2166     public void setInitialSz(final double initialSz) throws LockedException {
2167         if (running) {
2168             throw new LockedException();
2169         }
2170         this.initialSz = initialSz;
2171     }
2172 
2173     /**
2174      * Gets initial x-y cross coupling error of gyroscope.
2175      *
2176      * @return initial x-y cross coupling error of gyroscope.
2177      */
2178     @Override
2179     public double getInitialMxy() {
2180         return initialMxy;
2181     }
2182 
2183     /**
2184      * Sets initial x-y cross coupling error of gyroscope.
2185      *
2186      * @param initialMxy initial x-y cross coupling error of gyroscope.
2187      * @throws LockedException if calibrator is currently running.
2188      */
2189     @Override
2190     public void setInitialMxy(final double initialMxy) throws LockedException {
2191         if (running) {
2192             throw new LockedException();
2193         }
2194         this.initialMxy = initialMxy;
2195     }
2196 
2197     /**
2198      * Gets initial x-z cross coupling error of gyroscope.
2199      *
2200      * @return initial x-z cross coupling error of gyroscope.
2201      */
2202     @Override
2203     public double getInitialMxz() {
2204         return initialMxz;
2205     }
2206 
2207     /**
2208      * Sets initial x-z cross coupling error of gyroscope.
2209      *
2210      * @param initialMxz initial x-z cross coupling error of gyroscope.
2211      * @throws LockedException if calibrator is currently running.
2212      */
2213     @Override
2214     public void setInitialMxz(final double initialMxz) throws LockedException {
2215         if (running) {
2216             throw new LockedException();
2217         }
2218         this.initialMxz = initialMxz;
2219     }
2220 
2221     /**
2222      * Gets initial y-x cross coupling error of gyroscope.
2223      *
2224      * @return initial y-x cross coupling error of gyroscope.
2225      */
2226     @Override
2227     public double getInitialMyx() {
2228         return initialMyx;
2229     }
2230 
2231     /**
2232      * Sets initial y-x cross coupling error of gyroscope.
2233      *
2234      * @param initialMyx initial y-x cross coupling error of gyroscope.
2235      * @throws LockedException if calibrator is currently running.
2236      */
2237     @Override
2238     public void setInitialMyx(final double initialMyx) throws LockedException {
2239         if (running) {
2240             throw new LockedException();
2241         }
2242         this.initialMyx = initialMyx;
2243     }
2244 
2245     /**
2246      * Gets initial y-z cross coupling error of gyroscope.
2247      *
2248      * @return initial y-z cross coupling error of gyroscope.
2249      */
2250     @Override
2251     public double getInitialMyz() {
2252         return initialMyz;
2253     }
2254 
2255     /**
2256      * Sets initial y-z cross coupling error of gyroscope.
2257      *
2258      * @param initialMyz initial y-z cross coupling error of gyroscope.
2259      * @throws LockedException if calibrator is currently running.
2260      */
2261     @Override
2262     public void setInitialMyz(final double initialMyz) throws LockedException {
2263         if (running) {
2264             throw new LockedException();
2265         }
2266         this.initialMyz = initialMyz;
2267     }
2268 
2269     /**
2270      * Gets initial z-x cross coupling error of gyroscope.
2271      *
2272      * @return initial z-x cross coupling error of gyroscope.
2273      */
2274     @Override
2275     public double getInitialMzx() {
2276         return initialMzx;
2277     }
2278 
2279     /**
2280      * Sets initial z-x cross coupling error of gyroscope.
2281      *
2282      * @param initialMzx initial z-x cross coupling error of gyroscope.
2283      * @throws LockedException if calibrator is currently running.
2284      */
2285     @Override
2286     public void setInitialMzx(final double initialMzx) throws LockedException {
2287         if (running) {
2288             throw new LockedException();
2289         }
2290         this.initialMzx = initialMzx;
2291     }
2292 
2293     /**
2294      * Gets initial z-y cross coupling error of gyroscope.
2295      *
2296      * @return initial z-y cross coupling error of gyroscope.
2297      */
2298     @Override
2299     public double getInitialMzy() {
2300         return initialMzy;
2301     }
2302 
2303     /**
2304      * Sets initial z-y cross coupling error of gyroscope.
2305      *
2306      * @param initialMzy initial z-y cross coupling error of gyroscope.
2307      * @throws LockedException if calibrator is currently running.
2308      */
2309     @Override
2310     public void setInitialMzy(final double initialMzy) throws LockedException {
2311         if (running) {
2312             throw new LockedException();
2313         }
2314         this.initialMzy = initialMzy;
2315     }
2316 
2317     /**
2318      * Sets initial scaling factors of gyroscope.
2319      *
2320      * @param initialSx initial x scaling factor of gyroscope.
2321      * @param initialSy initial y scaling factor of gyroscope.
2322      * @param initialSz initial z scaling factor of gyroscope.
2323      * @throws LockedException if calibrator is currently running.
2324      */
2325     @Override
2326     public void setInitialScalingFactors(
2327             final double initialSx, final double initialSy, final double initialSz) throws LockedException {
2328         if (running) {
2329             throw new LockedException();
2330         }
2331         this.initialSx = initialSx;
2332         this.initialSy = initialSy;
2333         this.initialSz = initialSz;
2334     }
2335 
2336     /**
2337      * Sets initial cross coupling errors of gyroscope.
2338      *
2339      * @param initialMxy initial x-y cross coupling error of gyroscope.
2340      * @param initialMxz initial x-z cross coupling error of gyroscope.
2341      * @param initialMyx initial y-x cross coupling error of gyroscope.
2342      * @param initialMyz initial y-z cross coupling error of gyroscope.
2343      * @param initialMzx initial z-x cross coupling error of gyroscope.
2344      * @param initialMzy initial z-y cross coupling error of gyroscope.
2345      * @throws LockedException if calibrator is currently running.
2346      */
2347     @Override
2348     public void setInitialCrossCouplingErrors(
2349             final double initialMxy, final double initialMxz, final double initialMyx,
2350             final double initialMyz, final double initialMzx, final double initialMzy) throws LockedException {
2351         if (running) {
2352             throw new LockedException();
2353         }
2354         this.initialMxy = initialMxy;
2355         this.initialMxz = initialMxz;
2356         this.initialMyx = initialMyx;
2357         this.initialMyz = initialMyz;
2358         this.initialMzx = initialMzx;
2359         this.initialMzy = initialMzy;
2360     }
2361 
2362     /**
2363      * Sets initial scaling factors and cross coupling errors of
2364      * gyroscope.
2365      *
2366      * @param initialSx  initial x scaling factor of gyroscope.
2367      * @param initialSy  initial y scaling factor of gyroscope.
2368      * @param initialSz  initial z scaling factor of gyroscope.
2369      * @param initialMxy initial x-y cross coupling error of gyroscope.
2370      * @param initialMxz initial x-z cross coupling error of gyroscope.
2371      * @param initialMyx initial y-x cross coupling error of gyroscope.
2372      * @param initialMyz initial y-z cross coupling error of gyroscope.
2373      * @param initialMzx initial z-x cross coupling error of gyroscope.
2374      * @param initialMzy initial z-y cross coupling error of gyroscope.
2375      * @throws LockedException if calibrator is currently running.
2376      */
2377     @Override
2378     public void setInitialScalingFactorsAndCrossCouplingErrors(
2379             final double initialSx, final double initialSy, final double initialSz,
2380             final double initialMxy, final double initialMxz, final double initialMyx,
2381             final double initialMyz, final double initialMzx, final double initialMzy) throws LockedException {
2382         if (running) {
2383             throw new LockedException();
2384         }
2385         setInitialScalingFactors(initialSx, initialSy, initialSz);
2386         setInitialCrossCouplingErrors(initialMxy, initialMxz, initialMyx, initialMyz, initialMzx, initialMzy);
2387     }
2388 
2389     /**
2390      * Gets initial gyroscope bias to be used to find a solution as
2391      * an array.
2392      * Array values are expressed in radians per second (rad/s).
2393      *
2394      * @return array containing coordinates of initial gyroscope bias.
2395      */
2396     public double[] getInitialBias() {
2397         final var result = new double[BodyKinematics.COMPONENTS];
2398         getInitialBias(result);
2399         return result;
2400     }
2401 
2402     /**
2403      * Gets initial gyroscope bias to be used to find a solution as
2404      * an array.
2405      * Array values are expressed in radians per second (rad/s).
2406      *
2407      * @param result instance where result data will be copied to.
2408      * @throws IllegalArgumentException if provided array does not have length 3.
2409      */
2410     public void getInitialBias(final double[] result) {
2411         if (result.length != BodyKinematics.COMPONENTS) {
2412             throw new IllegalArgumentException();
2413         }
2414         result[0] = initialBiasX;
2415         result[1] = initialBiasY;
2416         result[2] = initialBiasZ;
2417     }
2418 
2419     /**
2420      * Sets initial gyroscope bias to be used to find a solution as
2421      * an array.
2422      * Array values are expressed in radians per second (rad/s).
2423      *
2424      * @param initialBias initial bias to find a solution.
2425      * @throws LockedException          if calibrator is currently running.
2426      * @throws IllegalArgumentException if provided array does not have length 3.
2427      */
2428     public void setInitialBias(final double[] initialBias) throws LockedException {
2429         if (running) {
2430             throw new LockedException();
2431         }
2432 
2433         if (initialBias.length != BodyKinematics.COMPONENTS) {
2434             throw new IllegalArgumentException();
2435         }
2436         initialBiasX = initialBias[0];
2437         initialBiasY = initialBias[1];
2438         initialBiasZ = initialBias[2];
2439     }
2440 
2441     /**
2442      * Gets initial gyroscope bias to be used to find a solution as a
2443      * column matrix.
2444      * Values are expressed in radians per second (rad/s).
2445      *
2446      * @return initial gyroscope bias to be used to find a solution as a
2447      * column matrix.
2448      */
2449     public Matrix getInitialBiasAsMatrix() {
2450         Matrix result;
2451         try {
2452             result = new Matrix(BodyKinematics.COMPONENTS, 1);
2453             getInitialBiasAsMatrix(result);
2454         } catch (final WrongSizeException ignore) {
2455             // never happens
2456             result = null;
2457         }
2458         return result;
2459     }
2460 
2461     /**
2462      * Gets initial gyroscope bias to be used to find a solution as a
2463      * column matrix.
2464      * Values are expressed in radians per second (rad/s).
2465      *
2466      * @param result instance where result data will be copied to.
2467      * @throws IllegalArgumentException if provided matrix is not 3x1.
2468      */
2469     public void getInitialBiasAsMatrix(final Matrix result) {
2470         if (result.getRows() != BodyKinematics.COMPONENTS || result.getColumns() != 1) {
2471             throw new IllegalArgumentException();
2472         }
2473         result.setElementAtIndex(0, initialBiasX);
2474         result.setElementAtIndex(1, initialBiasY);
2475         result.setElementAtIndex(2, initialBiasZ);
2476     }
2477 
2478     /**
2479      * Sets initial gyroscope bias to be used to find a solution as
2480      * a column matrix.
2481      * Values are expressed in radians per second (rad/s).
2482      *
2483      * @param initialBias initial gyroscope bias to find a solution.
2484      * @throws LockedException          if calibrator is currently running.
2485      * @throws IllegalArgumentException if provided matrix is not 3x1.
2486      */
2487     public void setInitialBias(final Matrix initialBias) throws LockedException {
2488         if (running) {
2489             throw new LockedException();
2490         }
2491         if (initialBias.getRows() != BodyKinematics.COMPONENTS || initialBias.getColumns() != 1) {
2492             throw new IllegalArgumentException();
2493         }
2494 
2495         initialBiasX = initialBias.getElementAtIndex(0);
2496         initialBiasY = initialBias.getElementAtIndex(1);
2497         initialBiasZ = initialBias.getElementAtIndex(2);
2498     }
2499 
2500     /**
2501      * Gets initial bias coordinates of gyroscope used to find a solution.
2502      *
2503      * @return initial bias coordinates.
2504      */
2505     public AngularSpeedTriad getInitialBiasAsTriad() {
2506         return new AngularSpeedTriad(AngularSpeedUnit.RADIANS_PER_SECOND, initialBiasX, initialBiasY, initialBiasZ);
2507     }
2508 
2509     /**
2510      * Gets initial bias coordinates of gyroscope used to find a solution.
2511      *
2512      * @param result instance where result will be stored.
2513      */
2514     public void getInitialBiasAsTriad(AngularSpeedTriad result) {
2515         result.setValueCoordinatesAndUnit(initialBiasX, initialBiasY, initialBiasZ,
2516                 AngularSpeedUnit.RADIANS_PER_SECOND);
2517     }
2518 
2519     /**
2520      * Sets initial bias coordinates of gyroscope used to find a solution.
2521      *
2522      * @param initialBias initial bias coordinates to be set.
2523      * @throws LockedException if calibrator is currently running.
2524      */
2525     public void setInitialBias(AngularSpeedTriad initialBias) throws LockedException {
2526         if (running) {
2527             throw new LockedException();
2528         }
2529 
2530         initialBiasX = convertAngularSpeed(initialBias.getValueX(), initialBias.getUnit());
2531         initialBiasY = convertAngularSpeed(initialBias.getValueY(), initialBias.getUnit());
2532         initialBiasZ = convertAngularSpeed(initialBias.getValueZ(), initialBias.getUnit());
2533     }
2534 
2535     /**
2536      * Gets initial gyroscope scale factors and cross coupling errors
2537      * matrix.
2538      *
2539      * @return initial gyroscope scale factors and cross coupling errors
2540      * matrix.
2541      */
2542     @Override
2543     public Matrix getInitialMg() {
2544         Matrix result;
2545         try {
2546             result = new Matrix(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
2547             getInitialMg(result);
2548         } catch (final WrongSizeException ignore) {
2549             // never happens
2550             result = null;
2551         }
2552         return result;
2553     }
2554 
2555     /**
2556      * Gets initial gyroscope scale factors and cross coupling errors
2557      * matrix.
2558      *
2559      * @param result instance where data will be stored.
2560      * @throws IllegalArgumentException if provided matrix is not 3x3.
2561      */
2562     @Override
2563     public void getInitialMg(final Matrix result) {
2564         if (result.getRows() != BodyKinematics.COMPONENTS || result.getColumns() != BodyKinematics.COMPONENTS) {
2565             throw new IllegalArgumentException();
2566         }
2567         result.setElementAtIndex(0, initialSx);
2568         result.setElementAtIndex(1, initialMyx);
2569         result.setElementAtIndex(2, initialMzx);
2570 
2571         result.setElementAtIndex(3, initialMxy);
2572         result.setElementAtIndex(4, initialSy);
2573         result.setElementAtIndex(5, initialMzy);
2574 
2575         result.setElementAtIndex(6, initialMxz);
2576         result.setElementAtIndex(7, initialMyz);
2577         result.setElementAtIndex(8, initialSz);
2578     }
2579 
2580     /**
2581      * Sets initial gyroscope scale factors and cross coupling errors matrix.
2582      *
2583      * @param initialMg initial scale factors and cross coupling errors matrix.
2584      * @throws IllegalArgumentException if provided matrix is not 3x3.
2585      * @throws LockedException          if calibrator is currently running.
2586      */
2587     @Override
2588     public void setInitialMg(final Matrix initialMg) throws LockedException {
2589         if (running) {
2590             throw new LockedException();
2591         }
2592         if (initialMg.getRows() != BodyKinematics.COMPONENTS || initialMg.getColumns() != BodyKinematics.COMPONENTS) {
2593             throw new IllegalArgumentException();
2594         }
2595 
2596         initialSx = initialMg.getElementAtIndex(0);
2597         initialMyx = initialMg.getElementAtIndex(1);
2598         initialMzx = initialMg.getElementAtIndex(2);
2599 
2600         initialMxy = initialMg.getElementAtIndex(3);
2601         initialSy = initialMg.getElementAtIndex(4);
2602         initialMzy = initialMg.getElementAtIndex(5);
2603 
2604         initialMxz = initialMg.getElementAtIndex(6);
2605         initialMyz = initialMg.getElementAtIndex(7);
2606         initialSz = initialMg.getElementAtIndex(8);
2607     }
2608 
2609     /**
2610      * Gets initial G-dependent cross biases introduced on the gyroscope by the
2611      * specific forces sensed by the accelerometer.
2612      *
2613      * @return a 3x3 matrix containing initial g-dependent cross biases.
2614      */
2615     @Override
2616     public Matrix getInitialGg() {
2617         return new Matrix(initialGg);
2618     }
2619 
2620     /**
2621      * Gets initial G-dependent cross biases introduced on the gyroscope by the
2622      * specific forces sensed by the accelerometer.
2623      *
2624      * @param result instance where data will be stored.
2625      * @throws IllegalArgumentException if provided matrix is not 3x3.
2626      */
2627     @Override
2628     public void getInitialGg(final Matrix result) {
2629         if (result.getRows() != BodyKinematics.COMPONENTS || result.getColumns() != BodyKinematics.COMPONENTS) {
2630             throw new IllegalArgumentException();
2631         }
2632 
2633         result.copyFrom(initialGg);
2634     }
2635 
2636     /**
2637      * Sets initial G-dependent cross biases introduced on the gyroscope by the
2638      * specific forces sensed by the accelerometer.
2639      *
2640      * @param initialGg g-dependent cross biases.
2641      * @throws LockedException          if calibrator is currently running.
2642      * @throws IllegalArgumentException if provided matrix is not 3x3.
2643      */
2644     @Override
2645     public void setInitialGg(final Matrix initialGg) throws LockedException {
2646         if (running) {
2647             throw new LockedException();
2648         }
2649 
2650         if (initialGg.getRows() != BodyKinematics.COMPONENTS || initialGg.getColumns() != BodyKinematics.COMPONENTS) {
2651             throw new IllegalArgumentException();
2652         }
2653 
2654         initialGg.copyTo(this.initialGg);
2655     }
2656 
2657     /**
2658      * Gets collection of sequences of timestamped body kinematics
2659      * measurements taken at a given position where the device moves freely
2660      * with different orientations.
2661      *
2662      * @return collection of sequences of timestamped body kinematics
2663      * measurements.
2664      */
2665     @Override
2666     public List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> getSequences() {
2667         return sequences;
2668     }
2669 
2670     /**
2671      * Sets collection of sequences of timestamped body kinematics
2672      * measurements taken at a given position where the device moves freely
2673      * with different orientations.
2674      *
2675      * @param sequences collection of sequences of timestamped body
2676      *                  kinematics measurements.
2677      * @throws LockedException if calibrator is currently running.
2678      */
2679     @Override
2680     public void setSequences(
2681             final List<BodyKinematicsSequence<StandardDeviationTimedBodyKinematics>> sequences) throws LockedException {
2682         if (running) {
2683             throw new LockedException();
2684         }
2685         this.sequences = sequences;
2686     }
2687 
2688     /**
2689      * Indicates the type of measurement or sequence used by this calibrator.
2690      *
2691      * @return type of measurement or sequence used by this calibrator.
2692      */
2693     @Override
2694     public GyroscopeCalibratorMeasurementOrSequenceType getMeasurementOrSequenceType() {
2695         return GyroscopeCalibratorMeasurementOrSequenceType.BODY_KINEMATICS_SEQUENCE;
2696     }
2697 
2698     /**
2699      * Indicates whether this calibrator requires ordered measurements or sequences
2700      * in a list or not.
2701      *
2702      * @return true if measurements or sequences must be ordered, false otherwise.
2703      */
2704     @Override
2705     public boolean isOrderedMeasurementsOrSequencesRequired() {
2706         return true;
2707     }
2708 
2709     /**
2710      * Indicates whether this calibrator requires quality scores for each
2711      * measurement/sequence or not.
2712      *
2713      * @return true if quality scores are required, false otherwise.
2714      */
2715     @Override
2716     public boolean isQualityScoresRequired() {
2717         return false;
2718     }
2719 
2720     /**
2721      * Indicates whether z-axis is assumed to be common for accelerometer and
2722      * gyroscope.
2723      * When enabled, this eliminates 3 variables from Ma matrix.
2724      *
2725      * @return true if z-axis is assumed to be common for accelerometer and gyroscope,
2726      * false otherwise.
2727      */
2728     @Override
2729     public boolean isCommonAxisUsed() {
2730         return commonAxisUsed;
2731     }
2732 
2733     /**
2734      * Specifies whether z-axis is assumed to be common for accelerometer and
2735      * gyroscope.
2736      * When enabled, this eliminates 3 variables from Ma matrix.
2737      *
2738      * @param commonAxisUsed true if z-axis is assumed to be common for accelerometer
2739      *                       and gyroscope, false otherwise.
2740      * @throws LockedException if calibrator is currently running.
2741      */
2742     @Override
2743     public void setCommonAxisUsed(final boolean commonAxisUsed) throws LockedException {
2744         if (running) {
2745             throw new LockedException();
2746         }
2747 
2748         this.commonAxisUsed = commonAxisUsed;
2749     }
2750 
2751     /**
2752      * Indicates whether G-dependent cross biases are being estimated
2753      * or not.
2754      * When enabled, this adds 9 variables from Gg matrix.
2755      *
2756      * @return true if G-dependent cross biases will be estimated,
2757      * false otherwise.
2758      */
2759     public boolean isGDependentCrossBiasesEstimated() {
2760         return estimateGDependentCrossBiases;
2761     }
2762 
2763     /**
2764      * Specifies whether G-dependent cross biases are being estimated
2765      * or not.
2766      * When enabled, this adds 9 variables from Gg matrix.
2767      *
2768      * @param estimateGDependentCrossBiases true if G-dependent cross
2769      *                                      biases will be estimated,
2770      *                                      false otherwise.
2771      * @throws LockedException if calibrator is currently running.
2772      */
2773     public void setGDependentCrossBiasesEstimated(final boolean estimateGDependentCrossBiases) throws LockedException {
2774         if (running) {
2775             throw new LockedException();
2776         }
2777 
2778         this.estimateGDependentCrossBiases = estimateGDependentCrossBiases;
2779     }
2780 
2781     /**
2782      * Gets listener to handle events raised by this estimator.
2783      *
2784      * @return listener to handle events raised by this estimator.
2785      */
2786     public EasyGyroscopeCalibratorListener getListener() {
2787         return listener;
2788     }
2789 
2790     /**
2791      * Sets listener to handle events raised by this estimator.
2792      *
2793      * @param listener listener to handle events raised by this estimator.
2794      * @throws LockedException if calibrator is currently running.
2795      */
2796     public void setListener(final EasyGyroscopeCalibratorListener listener) throws LockedException {
2797         if (running) {
2798             throw new LockedException();
2799         }
2800 
2801         this.listener = listener;
2802     }
2803 
2804     /**
2805      * Gets minimum number of required sequences.
2806      *
2807      * @return minimum number of required sequences.
2808      */
2809     @Override
2810     public int getMinimumRequiredMeasurementsOrSequences() {
2811         if (commonAxisUsed) {
2812             if (estimateGDependentCrossBiases) {
2813                 return MINIMUM_SEQUENCES_COMMON_Z_AXIS_AND_CROSS_BIASES;
2814             } else {
2815                 return MINIMUM_SEQUENCES_COMMON_Z_AXIS;
2816             }
2817         } else {
2818             if (estimateGDependentCrossBiases) {
2819                 return MINIMUM_SEQUENCES_GENERAL_AND_CROSS_BIASES;
2820             } else {
2821                 return MINIMUM_SEQUENCES_GENERAL;
2822             }
2823         }
2824     }
2825 
2826     /**
2827      * Indicates whether calibrator is ready to start.
2828      *
2829      * @return true if calibrator is ready, false otherwise.
2830      */
2831     @Override
2832     public boolean isReady() {
2833         return sequences != null && sequences.size() >= getMinimumRequiredMeasurementsOrSequences();
2834     }
2835 
2836     /**
2837      * Indicates whether calibrator is currently running or not.
2838      *
2839      * @return true if calibrator is running, false otherwise.
2840      */
2841     @Override
2842     public boolean isRunning() {
2843         return running;
2844     }
2845 
2846     /**
2847      * Estimates gyroscope calibration parameters containing bias, scale factors,
2848      * cross-coupling errors and G-dependent coupling.
2849      *
2850      * @throws LockedException      if calibrator is currently running.
2851      * @throws NotReadyException    if calibrator is not ready.
2852      * @throws CalibrationException if estimation fails for numerical reasons.
2853      */
2854     @Override
2855     public void calibrate() throws LockedException, NotReadyException, CalibrationException {
2856         if (running) {
2857             throw new LockedException();
2858         }
2859 
2860         if (!isReady()) {
2861             throw new NotReadyException();
2862         }
2863 
2864         try {
2865             running = true;
2866 
2867             reset();
2868 
2869             if (listener != null) {
2870                 listener.onCalibrateStart(this);
2871             }
2872 
2873             if (commonAxisUsed) {
2874                 if (estimateGDependentCrossBiases) {
2875                     calibrateCommonAxisAndGDependentCrossBiases();
2876                 } else {
2877                     calibrateCommonAxis();
2878                 }
2879             } else {
2880                 if (estimateGDependentCrossBiases) {
2881                     calibrateGeneralAndGDependentCrossBiases();
2882                 } else {
2883                     calibrateGeneral();
2884                 }
2885             }
2886 
2887             if (listener != null) {
2888                 listener.onCalibrateEnd(this);
2889             }
2890 
2891         } catch (final FittingException | AlgebraException | com.irurueta.numerical.NotReadyException e) {
2892             throw new CalibrationException(e);
2893         } finally {
2894             running = false;
2895         }
2896     }
2897 
2898     /**
2899      * Gets array containing x,y,z components of estimated gyroscope biases
2900      * expressed in radians per second (rad/s).
2901      *
2902      * @return array containing x,y,z components of estimated gyroscope biases.
2903      */
2904     @Override
2905     public double[] getEstimatedBiases() {
2906         return estimatedBiases;
2907     }
2908 
2909     /**
2910      * Gets array containing x,y,z components of estimated gyroscope biases
2911      * expressed in radians per second (rad/s).
2912      *
2913      * @param result instance where estimated gyroscope biases will be stored.
2914      * @return true if result instance was updated, false otherwise (when estimation
2915      * is not yet available).
2916      */
2917     @Override
2918     public boolean getEstimatedBiases(final double[] result) {
2919         if (estimatedBiases != null) {
2920             System.arraycopy(estimatedBiases, 0, result, 0, estimatedBiases.length);
2921             return true;
2922         } else {
2923             return false;
2924         }
2925     }
2926 
2927     /**
2928      * Gets column matrix containing x,y,z components of estimated gyroscope biases
2929      * expressed in radians per second (rad/s).
2930      *
2931      * @return column matrix containing x,y,z components of estimated gyroscope
2932      * biases.
2933      */
2934     @Override
2935     public Matrix getEstimatedBiasesAsMatrix() {
2936         return estimatedBiases != null ? Matrix.newFromArray(estimatedBiases) : null;
2937     }
2938 
2939     /**
2940      * Gets column matrix containing x,y,z components of estimated gyroscope biases
2941      * expressed in radians per second (rad/s).
2942      *
2943      * @param result instance where result data will be stored.
2944      * @return true if result was updated, false otherwise.
2945      * @throws WrongSizeException if provided result instance has invalid size.
2946      */
2947     @Override
2948     public boolean getEstimatedBiasesAsMatrix(final Matrix result) throws WrongSizeException {
2949         if (estimatedBiases != null) {
2950             result.fromArray(estimatedBiases);
2951             return true;
2952         } else {
2953             return false;
2954         }
2955     }
2956 
2957     /**
2958      * Gets x coordinate of estimated gyroscope bias expressed in radians per
2959      * second (rad/s).
2960      *
2961      * @return x coordinate of estimated gyroscope bias or null if not available.
2962      */
2963     @Override
2964     public Double getEstimatedBiasX() {
2965         return estimatedBiases != null ? estimatedBiases[0] : null;
2966     }
2967 
2968     /**
2969      * Gets y coordinate of estimated gyroscope bias expressed in radians per
2970      * second (rad/s).
2971      *
2972      * @return y coordinate of estimated gyroscope bias or null if not available.
2973      */
2974     @Override
2975     public Double getEstimatedBiasY() {
2976         return estimatedBiases != null ? estimatedBiases[1] : null;
2977     }
2978 
2979     /**
2980      * Gets z coordinate of estimated gyroscope bias expressed in radians per
2981      * second (rad/s).
2982      *
2983      * @return z coordinate of estimated gyroscope bias or null if not available.
2984      */
2985     @Override
2986     public Double getEstimatedBiasZ() {
2987         return estimatedBiases != null ? estimatedBiases[2] : null;
2988     }
2989 
2990     /**
2991      * Gets x coordinate of estimated gyroscope bias.
2992      *
2993      * @return x coordinate of estimated gyroscope bias or null if not available.
2994      */
2995     @Override
2996     public AngularSpeed getEstimatedBiasAngularSpeedX() {
2997         return estimatedBiases != null
2998                 ? new AngularSpeed(estimatedBiases[0], AngularSpeedUnit.RADIANS_PER_SECOND) : null;
2999     }
3000 
3001     /**
3002      * Gets x coordinate of estimated gyroscope bias.
3003      *
3004      * @param result instance where result will be stored.
3005      * @return true if result was updated, false if estimation is not available.
3006      */
3007     @Override
3008     public boolean getEstimatedBiasAngularSpeedX(final AngularSpeed result) {
3009         if (estimatedBiases != null) {
3010             result.setValue(estimatedBiases[0]);
3011             result.setUnit(AngularSpeedUnit.RADIANS_PER_SECOND);
3012             return true;
3013         } else {
3014             return false;
3015         }
3016     }
3017 
3018     /**
3019      * Gets y coordinate of estimated gyroscope bias.
3020      *
3021      * @return y coordinate of estimated gyroscope bias or null if not available.
3022      */
3023     @Override
3024     public AngularSpeed getEstimatedBiasAngularSpeedY() {
3025         return estimatedBiases != null
3026                 ? new AngularSpeed(estimatedBiases[1], AngularSpeedUnit.RADIANS_PER_SECOND) : null;
3027     }
3028 
3029     /**
3030      * Gets y coordinate of estimated gyroscope bias.
3031      *
3032      * @param result instance where result will be stored.
3033      * @return true if result was updated, false if estimation is not available.
3034      */
3035     @Override
3036     public boolean getEstimatedBiasAngularSpeedY(final AngularSpeed result) {
3037         if (estimatedBiases != null) {
3038             result.setValue(estimatedBiases[1]);
3039             result.setUnit(AngularSpeedUnit.RADIANS_PER_SECOND);
3040             return true;
3041         } else {
3042             return false;
3043         }
3044     }
3045 
3046     /**
3047      * Gets z coordinate of estimated gyroscope bias.
3048      *
3049      * @return z coordinate of estimated gyroscope bias or null if not available.
3050      */
3051     @Override
3052     public AngularSpeed getEstimatedBiasAngularSpeedZ() {
3053         return estimatedBiases != null
3054                 ? new AngularSpeed(estimatedBiases[2], AngularSpeedUnit.RADIANS_PER_SECOND) : null;
3055     }
3056 
3057     /**
3058      * Gets z coordinate of estimated gyroscope bias.
3059      *
3060      * @param result instance where result will be stored.
3061      * @return true if result was updated, false if estimation is not available.
3062      */
3063     @Override
3064     public boolean getEstimatedBiasAngularSpeedZ(final AngularSpeed result) {
3065         if (estimatedBiases != null) {
3066             result.setValue(estimatedBiases[2]);
3067             result.setUnit(AngularSpeedUnit.RADIANS_PER_SECOND);
3068             return true;
3069         } else {
3070             return false;
3071         }
3072     }
3073 
3074     /**
3075      * Gets estimated gyroscope bias.
3076      *
3077      * @return estimated gyroscope bias or null if not available.
3078      */
3079     @Override
3080     public AngularSpeedTriad getEstimatedBiasAsTriad() {
3081         return estimatedBiases != null
3082                 ? new AngularSpeedTriad(AngularSpeedUnit.RADIANS_PER_SECOND,
3083                 estimatedBiases[0], estimatedBiases[1], estimatedBiases[2])
3084                 : null;
3085     }
3086 
3087     /**
3088      * Gets estimated gyroscope bias.
3089      *
3090      * @param result instance where result will be stored.
3091      * @return true if estimated gyroscope bias is available and result was
3092      * modified, false otherwise.
3093      */
3094     @Override
3095     public boolean getEstimatedBiasAsTriad(final AngularSpeedTriad result) {
3096         if (estimatedBiases != null) {
3097             result.setValueCoordinatesAndUnit(
3098                     estimatedBiases[0], estimatedBiases[1], estimatedBiases[2], AngularSpeedUnit.RADIANS_PER_SECOND);
3099             return true;
3100         } else {
3101             return false;
3102         }
3103     }
3104 
3105     /**
3106      * Gets estimated gyroscope scale factors and cross coupling errors.
3107      * This is the product of matrix Tg containing cross coupling errors and Kg
3108      * containing scaling factors.
3109      * So that:
3110      * <pre>
3111      *     Mg = [sx    mxy  mxz] = Tg*Kg
3112      *          [myx   sy   myz]
3113      *          [mzx   mzy  sz ]
3114      * </pre>
3115      * Where:
3116      * <pre>
3117      *     Kg = [sx 0   0 ]
3118      *          [0  sy  0 ]
3119      *          [0  0   sz]
3120      * </pre>
3121      * and
3122      * <pre>
3123      *     Tg = [1          -alphaXy    alphaXz ]
3124      *          [alphaYx    1           -alphaYz]
3125      *          [-alphaZx   alphaZy     1       ]
3126      * </pre>
3127      * Hence:
3128      * <pre>
3129      *     Mg = [sx    mxy  mxz] = Tg*Kg =  [sx             -sy * alphaXy   sz * alphaXz ]
3130      *          [myx   sy   myz]            [sx * alphaYx   sy              -sz * alphaYz]
3131      *          [mzx   mzy  sz ]            [-sx * alphaZx  sy * alphaZy    sz           ]
3132      * </pre>
3133      * This instance allows any 3x3 matrix however, typically alphaYx, alphaZx and alphaZy
3134      * are considered to be zero if the gyroscope z-axis is assumed to be the same
3135      * as the body z-axis. When this is assumed, myx = mzx = mzy = 0 and the Mg matrix
3136      * becomes upper diagonal:
3137      * <pre>
3138      *     Mg = [sx    mxy  mxz]
3139      *          [0     sy   myz]
3140      *          [0     0    sz ]
3141      * </pre>
3142      * Values of this matrix are unit-less.
3143      *
3144      * @return estimated gyroscope scale factors and cross coupling errors, or null
3145      * if not available.
3146      */
3147     @Override
3148     public Matrix getEstimatedMg() {
3149         return estimatedMg;
3150     }
3151 
3152     /**
3153      * Gets estimated gyroscope x-axis scale factor.
3154      *
3155      * @return estimated gyroscope x-axis scale factor or null
3156      * if not available.
3157      */
3158     @Override
3159     public Double getEstimatedSx() {
3160         return estimatedMg != null ? estimatedMg.getElementAt(0, 0) : null;
3161     }
3162 
3163     /**
3164      * Gets estimated gyroscope y-axis scale factor.
3165      *
3166      * @return estimated gyroscope y-axis scale factor or null
3167      * if not available.
3168      */
3169     @Override
3170     public Double getEstimatedSy() {
3171         return estimatedMg != null ? estimatedMg.getElementAt(1, 1) : null;
3172     }
3173 
3174     /**
3175      * Gets estimated gyroscope z-axis scale factor.
3176      *
3177      * @return estimated gyroscope z-axis scale factor or null
3178      * if not available.
3179      */
3180     @Override
3181     public Double getEstimatedSz() {
3182         return estimatedMg != null ? estimatedMg.getElementAt(2, 2) : null;
3183     }
3184 
3185     /**
3186      * Gets estimated gyroscope x-y cross-coupling error.
3187      *
3188      * @return estimated gyroscope x-y cross-coupling error or null
3189      * if not available.
3190      */
3191     @Override
3192     public Double getEstimatedMxy() {
3193         return estimatedMg != null ? estimatedMg.getElementAt(0, 1) : null;
3194     }
3195 
3196     /**
3197      * Gets estimated gyroscope x-z cross-coupling error.
3198      *
3199      * @return estimated gyroscope x-z cross-coupling error or null
3200      * if not available.
3201      */
3202     @Override
3203     public Double getEstimatedMxz() {
3204         return estimatedMg != null ? estimatedMg.getElementAt(0, 2) : null;
3205     }
3206 
3207     /**
3208      * Gets estimated gyroscope y-x cross-coupling error.
3209      *
3210      * @return estimated gyroscope y-x cross-coupling error or null
3211      * if not available.
3212      */
3213     @Override
3214     public Double getEstimatedMyx() {
3215         return estimatedMg != null ? estimatedMg.getElementAt(1, 0) : null;
3216     }
3217 
3218     /**
3219      * Gets estimated gyroscope y-z cross-coupling error.
3220      *
3221      * @return estimated gyroscope y-z cross-coupling error or null
3222      * if not available.
3223      */
3224     @Override
3225     public Double getEstimatedMyz() {
3226         return estimatedMg != null ? estimatedMg.getElementAt(1, 2) : null;
3227     }
3228 
3229     /**
3230      * Gets estimated gyroscope z-x cross-coupling error.
3231      *
3232      * @return estimated gyroscope z-x cross-coupling error or null
3233      * if not available.
3234      */
3235     @Override
3236     public Double getEstimatedMzx() {
3237         return estimatedMg != null ? estimatedMg.getElementAt(2, 0) : null;
3238     }
3239 
3240     /**
3241      * Gets estimated gyroscope z-y cross-coupling error.
3242      *
3243      * @return estimated gyroscope z-y cross-coupling error or null
3244      * if not available.
3245      */
3246     @Override
3247     public Double getEstimatedMzy() {
3248         return estimatedMg != null ? estimatedMg.getElementAt(2, 1) : null;
3249     }
3250 
3251     /**
3252      * Gets estimated G-dependent cross biases introduced on the gyroscope by the
3253      * specific forces sensed by the accelerometer.
3254      * This instance allows any 3x3 matrix.
3255      *
3256      * @return estimated G-dependent cross biases.
3257      */
3258     @Override
3259     public Matrix getEstimatedGg() {
3260         return estimatedGg;
3261     }
3262 
3263     /**
3264      * Gets estimated covariance matrix for estimated parameters.
3265      * Diagonal elements of the matrix contains variance for the following
3266      * parameters (following indicated order): bgx, bgy, bgz, sx, sy, sz,
3267      * mxy, mxz, myx, myz, mzx, mzy, gg11, gg21, gg31, gg12, gg22, gg32,
3268      * gg13, gg23, gg33.
3269      *
3270      * @return estimated covariance matrix for estimated parameters.
3271      */
3272     @Override
3273     public Matrix getEstimatedCovariance() {
3274         return estimatedCovariance;
3275     }
3276 
3277     /**
3278      * Gets estimated chi square value.
3279      *
3280      * @return estimated chi square value.
3281      */
3282     @Override
3283     public double getEstimatedChiSq() {
3284         return estimatedChiSq;
3285     }
3286 
3287     /**
3288      * Gets estimated chi square degrees of freedom. Degrees of freedom is equal to the number of sampled data minus the
3289      * number of estimated parameters.
3290      *
3291      * @return estimated degrees of freedom of chi square value
3292      */
3293     @Override
3294     public int getEstimatedChiSqDegreesOfFreedom() {
3295         return estimatedChiSqDegreesOfFreedom;
3296     }
3297 
3298     /**
3299      * Gets estimated reduced chi square value. This is equal to estimated chi square value divided by its degrees of
3300      * freedom. Ideally this value should be close to 1.0, indicating that fit is optimal.
3301      * A value larger than 1.0 indicates that fit is not good or noise has been underestimated, and a value smaller than
3302      * 1.0 indicates that there is overfitting or noise has been overestimated.
3303      *
3304      * @return estimated reduced chi square value
3305      */
3306     @Override
3307     public double getEstimatedReducedChiSq() {
3308         return estimatedReducedChiSq;
3309     }
3310 
3311     /**
3312      * Gets estimated probability of finding a smaller chi square value expressed as a value between 0.0 and 1.0. The
3313      * smaller the found chi square value is, the better the fit of the estimated parameters to the actual parameter.
3314      * Thus, the smaller the chance of finding a smaller chi square value, then the better the estimated fit is.
3315      *
3316      * @return estimated probability of finding a smaller chi square value.
3317      */
3318     @Override
3319     public double getEstimatedP() {
3320         return estimatedP;
3321     }
3322 
3323     /**
3324      * Gets estimated measure of quality of estimated fit as a value between 0.0 and 1.0. The larger the quality value
3325      * is, the better the fit that has been estimated.
3326      *
3327      * @return estimated measure of quality of estimated fit.
3328      */
3329     @Override
3330     public double getEstimatedQ() {
3331         return estimatedQ;
3332     }
3333 
3334     /**
3335      * Gets estimated mean square error respect to provided measurements.
3336      *
3337      * @return estimated mean square error respect to provided measurements.
3338      */
3339     @Override
3340     public double getEstimatedMse() {
3341         return estimatedMse;
3342     }
3343 
3344     /**
3345      * Gets variance of estimated x coordinate of gyroscope bias expressed in (rad^2/s^2).
3346      *
3347      * @return variance of estimated x coordinate of gyroscope bias or null if not available.
3348      */
3349     public Double getEstimatedBiasXVariance() {
3350         return estimatedCovariance != null ? estimatedCovariance.getElementAt(0, 0) : null;
3351     }
3352 
3353     /**
3354      * Gets standard deviation of estimated x coordinate of gyroscope bias expressed in
3355      * radians per second (rad/s).
3356      *
3357      * @return standard deviation of estimated x coordinate of gyroscope bias or null if not
3358      * available.
3359      */
3360     public Double getEstimatedBiasXStandardDeviation() {
3361         final var variance = getEstimatedBiasXVariance();
3362         return variance != null ? Math.sqrt(variance) : null;
3363     }
3364 
3365     /**
3366      * Gets standard deviation of estimated x coordinate of gyroscope bias.
3367      *
3368      * @return standard deviation of estimated x coordinate of gyroscope bias or null if not
3369      * available.
3370      */
3371     public AngularSpeed getEstimatedBiasXStandardDeviationAsAngularSpeed() {
3372         return estimatedCovariance != null
3373                 ? new AngularSpeed(getEstimatedBiasXStandardDeviation(), AngularSpeedUnit.RADIANS_PER_SECOND)
3374                 : null;
3375     }
3376 
3377     /**
3378      * Gets standard deviation of estimated x coordinate of gyroscope bias.
3379      *
3380      * @param result instance where result will be stored.
3381      * @return true if standard deviation of estimated x coordinate of gyroscope bias is available,
3382      * false otherwise.
3383      */
3384     public boolean getEstimatedBiasXStandardDeviationAsAngularSpeed(final AngularSpeed result) {
3385         if (estimatedCovariance != null) {
3386             result.setValue(getEstimatedBiasXStandardDeviation());
3387             result.setUnit(AngularSpeedUnit.RADIANS_PER_SECOND);
3388             return true;
3389         } else {
3390             return false;
3391         }
3392     }
3393 
3394     /**
3395      * Gets variance of estimated y coordinate of gyroscope bias expressed in (rad^2/s^2).
3396      *
3397      * @return variance of estimated y coordinate of gyroscope bias or null if not available.
3398      */
3399     public Double getEstimatedBiasYVariance() {
3400         return estimatedCovariance != null ? estimatedCovariance.getElementAt(1, 1) : null;
3401     }
3402 
3403     /**
3404      * Gets standard deviation of estimated y coordinate of gyroscope bias expressed in
3405      * radians per second (rad/s).
3406      *
3407      * @return standard deviation of estimated y coordinate of gyroscope bias or null if not
3408      * available.
3409      */
3410     public Double getEstimatedBiasYStandardDeviation() {
3411         final var variance = getEstimatedBiasYVariance();
3412         return variance != null ? Math.sqrt(variance) : null;
3413     }
3414 
3415     /**
3416      * Gets standard deviation of estimated y coordinate of gyroscope bias.
3417      *
3418      * @return standard deviation of estimated y coordinate of gyroscope bias or null if not
3419      * available.
3420      */
3421     public AngularSpeed getEstimatedBiasYStandardDeviationAsAngularSpeed() {
3422         return estimatedCovariance != null
3423                 ? new AngularSpeed(getEstimatedBiasYStandardDeviation(), AngularSpeedUnit.RADIANS_PER_SECOND)
3424                 : null;
3425     }
3426 
3427     /**
3428      * Gets standard deviation of estimated y coordinate of gyroscope bias.
3429      *
3430      * @param result instance where result will be stored.
3431      * @return true if standard deviation of estimated y coordinate of gyroscope bias is available,
3432      * false otherwise.
3433      */
3434     public boolean getEstimatedBiasYStandardDeviationAsAngularSpeed(final AngularSpeed result) {
3435         if (estimatedCovariance != null) {
3436             result.setValue(getEstimatedBiasYStandardDeviation());
3437             result.setUnit(AngularSpeedUnit.RADIANS_PER_SECOND);
3438             return true;
3439         } else {
3440             return false;
3441         }
3442     }
3443 
3444     /**
3445      * Gets variance of estimated z coordinate of gyroscope bias expressed in (rad^2/s^2).
3446      *
3447      * @return variance of estimated z coordinate of gyroscope bias or null if not available.
3448      */
3449     public Double getEstimatedBiasZVariance() {
3450         return estimatedCovariance != null ? estimatedCovariance.getElementAt(2, 2) : null;
3451     }
3452 
3453     /**
3454      * Gets standard deviation of estimated z coordinate of gyroscope bias expressed in
3455      * radians per second (rad/s).
3456      *
3457      * @return standard deviation of estimated z coordinate of gyroscope bias or null if not
3458      * available.
3459      */
3460     public Double getEstimatedBiasZStandardDeviation() {
3461         final var variance = getEstimatedBiasZVariance();
3462         return variance != null ? Math.sqrt(variance) : null;
3463     }
3464 
3465     /**
3466      * Gets standard deviation of estimated z coordinate of gyroscope bias.
3467      *
3468      * @return standard deviation of estimated z coordinate of gyroscope bias or null if not
3469      * available.
3470      */
3471     public AngularSpeed getEstimatedBiasZStandardDeviationAsAngularSpeed() {
3472         return estimatedCovariance != null
3473                 ? new AngularSpeed(getEstimatedBiasZStandardDeviation(), AngularSpeedUnit.RADIANS_PER_SECOND)
3474                 : null;
3475     }
3476 
3477     /**
3478      * Gets standard deviation of estimated z coordinate of gyroscope bias.
3479      *
3480      * @param result instance where result will be stored.
3481      * @return true if standard deviation of estimated z coordinate of gyroscope bias is available,
3482      * false otherwise.
3483      */
3484     public boolean getEstimatedBiasZStandardDeviationAsAngularSpeed(final AngularSpeed result) {
3485         if (estimatedCovariance != null) {
3486             result.setValue(getEstimatedBiasZStandardDeviation());
3487             result.setUnit(AngularSpeedUnit.RADIANS_PER_SECOND);
3488             return true;
3489         } else {
3490             return false;
3491         }
3492     }
3493 
3494     /**
3495      * Gets standard deviation of estimated gyroscope bias coordinates.
3496      *
3497      * @return standard deviation of estimated gyroscope bias coordinates.
3498      */
3499     public AngularSpeedTriad getEstimatedBiasStandardDeviation() {
3500         return estimatedCovariance != null
3501                 ? new AngularSpeedTriad(AngularSpeedUnit.RADIANS_PER_SECOND,
3502                 getEstimatedBiasXStandardDeviation(),
3503                 getEstimatedBiasYStandardDeviation(),
3504                 getEstimatedBiasZStandardDeviation())
3505                 : null;
3506     }
3507 
3508     /**
3509      * Gets standard deviation of estimated gyroscope bias coordinates.
3510      *
3511      * @param result instance where result will be stored.
3512      * @return true if standard deviation of gyroscope bias was available, false
3513      * otherwise.
3514      */
3515     public boolean getEstimatedBiasStandardDeviation(final AngularSpeedTriad result) {
3516         if (estimatedCovariance != null) {
3517             result.setValueCoordinatesAndUnit(
3518                     getEstimatedBiasXStandardDeviation(),
3519                     getEstimatedBiasYStandardDeviation(),
3520                     getEstimatedBiasZStandardDeviation(),
3521                     AngularSpeedUnit.RADIANS_PER_SECOND);
3522             return true;
3523         } else {
3524             return false;
3525         }
3526     }
3527 
3528     /**
3529      * Gets average of estimated standard deviation of gyroscope bias coordinates expressed
3530      * in radians per second (rad/s).
3531      *
3532      * @return average of estimated standard deviation of gyroscope bias coordinates or null
3533      * if not available.
3534      */
3535     public Double getEstimatedBiasStandardDeviationAverage() {
3536         return estimatedCovariance != null
3537                 ? (getEstimatedBiasXStandardDeviation() + getEstimatedBiasYStandardDeviation()
3538                 + getEstimatedBiasZStandardDeviation()) / 3.0 : null;
3539     }
3540 
3541     /**
3542      * Gets average of estimated standard deviation of gyroscope bias coordinates.
3543      *
3544      * @return average of estimated standard deviation of gyroscope bias coordinates or null.
3545      */
3546     public AngularSpeed getEstimatedBiasStandardDeviationAverageAsAngularSpeed() {
3547         return estimatedCovariance != null
3548                 ? new AngularSpeed(getEstimatedBiasStandardDeviationAverage(), AngularSpeedUnit.RADIANS_PER_SECOND)
3549                 : null;
3550     }
3551 
3552     /**
3553      * Gets average of estimated standard deviation of gyroscope bias coordinates.
3554      *
3555      * @param result instance where result will be stored.
3556      * @return true if average of estimated standard deviation of gyroscope bias is available,
3557      * false otherwise.
3558      */
3559     public boolean getEstimatedBiasStandardDeviationAverageAsAngularSpeed(final AngularSpeed result) {
3560         if (estimatedCovariance != null) {
3561             result.setValue(getEstimatedBiasStandardDeviationAverage());
3562             result.setUnit(AngularSpeedUnit.RADIANS_PER_SECOND);
3563             return true;
3564         } else {
3565             return false;
3566         }
3567     }
3568 
3569     /**
3570      * Gets norm of estimated standard deviation of gyroscope bias expressed in
3571      * radians per second (rad/s).
3572      * This can be used as the initial gyroscope bias uncertainty for
3573      * {@link INSLooselyCoupledKalmanInitializerConfig} or {@link INSTightlyCoupledKalmanInitializerConfig}.
3574      *
3575      * @return norm of estimated standard deviation of gyroscope bias or null
3576      * if not available.
3577      */
3578     @Override
3579     public Double getEstimatedBiasStandardDeviationNorm() {
3580         return estimatedCovariance != null
3581                 ? Math.sqrt(getEstimatedBiasXVariance() + getEstimatedBiasYVariance() + getEstimatedBiasZVariance())
3582                 : null;
3583     }
3584 
3585     /**
3586      * Gets norm of estimated standard deviation of gyroscope bias.
3587      * This can be used as the initial gyroscope bias uncertainty for
3588      * {@link INSLooselyCoupledKalmanInitializerConfig} or {@link INSTightlyCoupledKalmanInitializerConfig}.
3589      *
3590      * @return norm of estimated standard deviation of gyroscope bias or null
3591      * if not available.
3592      */
3593     public AngularSpeed getEstimatedBiasStandardDeviationNormAsAngularSpeed() {
3594         return estimatedCovariance != null
3595                 ? new AngularSpeed(getEstimatedBiasStandardDeviationNorm(), AngularSpeedUnit.RADIANS_PER_SECOND)
3596                 : null;
3597     }
3598 
3599     /**
3600      * Gets norm of estimated standard deviation of gyroscope bias coordinates.
3601      * This can be used as the initial gyroscope bias uncertainty for
3602      * {@link INSLooselyCoupledKalmanInitializerConfig} or {@link INSTightlyCoupledKalmanInitializerConfig}.
3603      *
3604      * @param result instance where result will be stored.
3605      * @return true if norm of estimated standard deviation of gyroscope bias is
3606      * available, false otherwise.
3607      */
3608     public boolean getEstimatedBiasStandardDeviationNormAsAngularSpeed(final AngularSpeed result) {
3609         if (estimatedCovariance != null) {
3610             result.setValue(getEstimatedBiasStandardDeviationNorm());
3611             result.setUnit(AngularSpeedUnit.RADIANS_PER_SECOND);
3612             return true;
3613         } else {
3614             return false;
3615         }
3616     }
3617 
3618     /**
3619      * Internal method to perform calibration when common z-axis is assumed
3620      * for both the accelerometer and gyroscope and when G-dependent cross
3621      * biases are being estimated.
3622      *
3623      * @throws AlgebraException                         if accelerometer parameters prevent fixing
3624      *                                                  measured accelerometer values due to numerical instabilities.
3625      * @throws FittingException                         if no convergence to solution is found.
3626      * @throws com.irurueta.numerical.NotReadyException if fitter is not ready.
3627      */
3628     private void calibrateCommonAxisAndGDependentCrossBiases() throws AlgebraException, FittingException,
3629             com.irurueta.numerical.NotReadyException {
3630         // The gyroscope model is
3631         // Ωmeas = bg + (I + Mg) * Ωtrue + Gg * ftrue + w
3632 
3633         // Ideally a least squares solution tries to minimize noise component, so:
3634         // Ωmeas = bg + (I + Mg) * Ωtrue + Gg * ftrue
3635 
3636         // Where Mg is upper triangular
3637 
3638         // For convergence purposes of the Levenberg-Marquardt algorithm, we
3639         // take common factor M = I + Mg
3640 
3641         // and the gyroscope model can be better expressed as:
3642 
3643         // Ωmeas = M*(Ωtrue + b + G * ftrue)
3644 
3645         // where:
3646         // bg = M*b --> b = M^-1*bg
3647         // Gg = M*G --> G = M^-1*Gg
3648 
3649         // Hence
3650         // Ωmeas - M*b - M*G*ftrue = M*Ωtrue
3651         // M^-1 * (Ωmeas - M*b - M*G*ftrue) = Ωtrue
3652 
3653         // Notice that M is upper diagonal because Mg is upper diagonal
3654         // when common axis is assumed
3655 
3656         final var gradientEstimator = new GradientEstimator(
3657                 params -> evaluateCommonAxisWithGDependentCrossBiases(i, params));
3658 
3659         final var initialM = Matrix.identity(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
3660         initialM.add(getInitialMg());
3661 
3662         // Force initial M to be upper diagonal
3663         initialM.setElementAt(1, 0, 0.0);
3664         initialM.setElementAt(2, 0, 0.0);
3665         initialM.setElementAt(2, 1, 0.0);
3666 
3667         final var invInitialM = Utils.inverse(initialM);
3668         final var initialBg = getInitialBiasAsMatrix();
3669         final var initialB = invInitialM.multiplyAndReturnNew(initialBg);
3670         final var initialGg = getInitialGg();
3671         final var initialG = invInitialM.multiplyAndReturnNew(initialGg);
3672 
3673         fitter.setFunctionEvaluator(new LevenbergMarquardtMultiDimensionFunctionEvaluator() {
3674             @Override
3675             public int getNumberOfDimensions() {
3676                 // Before and after normalized gravity versors
3677                 return 2 * BodyKinematics.COMPONENTS;
3678             }
3679 
3680             @Override
3681             public double[] createInitialParametersArray() {
3682                 final var initial = new double[COMMON_Z_AXIS_UNKNOWNS_AND_CROSS_BIASES];
3683 
3684                 // biases b
3685                 for (var i = 0; i < BodyKinematics.COMPONENTS; i++) {
3686                     initial[i] = initialB.getElementAtIndex(i);
3687                 }
3688 
3689                 // upper diagonal cross coupling errors M
3690                 var k = BodyKinematics.COMPONENTS;
3691                 for (var j = 0; j < BodyKinematics.COMPONENTS; j++) {
3692                     for (var i = 0; i < BodyKinematics.COMPONENTS; i++) {
3693                         if (i <= j) {
3694                             initial[k] = initialM.getElementAt(i, j);
3695                             k++;
3696                         }
3697                     }
3698                 }
3699 
3700                 // g-dependent cross biases G
3701                 final var num = BodyKinematics.COMPONENTS * BodyKinematics.COMPONENTS;
3702                 for (int i = 0, j = k; i < num; i++, j++) {
3703                     initial[j] = initialG.getElementAtIndex(i);
3704                 }
3705 
3706                 return initial;
3707             }
3708 
3709             @Override
3710             public double evaluate(
3711                     final int i, final double[] point, final double[] params, final double[] derivatives)
3712                     throws EvaluationException {
3713                 EasyGyroscopeCalibrator.this.i = i;
3714 
3715                 // point contains fixed gravity versor values for current
3716                 // sequence
3717                 EasyGyroscopeCalibrator.this.point = point;
3718 
3719                 gradientEstimator.gradient(params, derivatives);
3720 
3721                 return evaluateCommonAxisWithGDependentCrossBiases(i, params);
3722             }
3723         });
3724 
3725         setInputData();
3726 
3727         fitter.fit();
3728 
3729         final var result = fitter.getA();
3730 
3731         final var bx = result[0];
3732         final var by = result[1];
3733         final var bz = result[2];
3734 
3735         final var m11 = result[3];
3736 
3737         final var m12 = result[4];
3738         final var m22 = result[5];
3739 
3740         final var m13 = result[6];
3741         final var m23 = result[7];
3742         final var m33 = result[8];
3743 
3744         final var g11 = result[9];
3745         final var g21 = result[10];
3746         final var g31 = result[11];
3747 
3748         final var g12 = result[12];
3749         final var g22 = result[13];
3750         final var g32 = result[14];
3751 
3752         final var g13 = result[15];
3753         final var g23 = result[16];
3754         final var g33 = result[17];
3755 
3756         final var b = new Matrix(BodyKinematics.COMPONENTS, 1);
3757         b.setElementAtIndex(0, bx);
3758         b.setElementAtIndex(1, by);
3759         b.setElementAtIndex(2, bz);
3760 
3761         final var m = new Matrix(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
3762         m.setElementAtIndex(0, m11);
3763         m.setElementAtIndex(1, 0.0);
3764         m.setElementAtIndex(2, 0.0);
3765 
3766         m.setElementAtIndex(3, m12);
3767         m.setElementAtIndex(4, m22);
3768         m.setElementAtIndex(5, 0.0);
3769 
3770         m.setElementAtIndex(6, m13);
3771         m.setElementAtIndex(7, m23);
3772         m.setElementAtIndex(8, m33);
3773 
3774         final var g = new Matrix(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
3775         g.setElementAtIndex(0, g11);
3776         g.setElementAtIndex(1, g21);
3777         g.setElementAtIndex(2, g31);
3778 
3779         g.setElementAtIndex(3, g12);
3780         g.setElementAtIndex(4, g22);
3781         g.setElementAtIndex(5, g32);
3782 
3783         g.setElementAtIndex(6, g13);
3784         g.setElementAtIndex(7, g23);
3785         g.setElementAtIndex(8, g33);
3786 
3787         setResult(m, b, g);
3788 
3789         // at this point covariance is expressed in terms of b, M and G, and must
3790         // be expressed in terms of bg, Mg and Gg.
3791         // We know that:
3792         // bg = M * b
3793         // Mg = M - I
3794         // Gg = M * G
3795 
3796         // b = [bx]
3797         //     [by]
3798         //     [bz]
3799 
3800         // M = [m11  m12  m13]
3801         //     [0    m22  m23]
3802         //     [0    0    m33]
3803 
3804         // G = [g11  g12  g13]
3805         //     [g21  g22  g23]
3806         //     [g31  g32  g33]
3807 
3808         // bg = [m11  m12  m13][bx] = [m11 * bx + m12 * by + m13 * bz] = [bgx]
3809         //      [0    m22  m23][by]   [           m22 * by + m23 * bz]   [bgy]
3810         //      [0    0    m33][bz]   [                      m33 * bz]   [bgz]
3811 
3812         // Mg = [sx  mxy  mxz] = [m11 - 1    m12         m13     ]
3813         //      [myx sy	  myz]   [0          m22 - 1     m23     ]
3814         //      [mzx mzy  sz ]   [0          0           m33 - 1 ]
3815 
3816         // Gg = [gg11  gg12  gg13] = [m11  m12  m13][g11  g12  g13]
3817         //      [gg21  gg22  gg23]   [0    m22  m23][g21  g22  g23]
3818         //      [gg31  gg32  gg33]   [0    0    m33][g31  g32  g33]
3819 
3820         // Defining the linear application:
3821         // F(b, M, G) = F(bx, by, bz, m11, m12, m22, m13, m23, m33, g11, g21, g31, g12, g22, g32, g13, g23, g33)
3822         // as:
3823         // [bgx] =  [m11 * bx + m12 * by + m13 * bz]
3824         // [bgy]    [           m22 * by + m23 * bz]
3825         // [bgz]    [                      m33 * bz]
3826         // [sx]     [m11 - 1]
3827         // [sy]     [m22 - 1]
3828         // [sz]     [m33 - 1]
3829         // [mxy]    [m12]
3830         // [mxz]    [m13]
3831         // [myx]    [0]
3832         // [myz]    [m23]
3833         // [mzx]    [0]
3834         // [mzy]    [0]
3835         // [gg11]   [m11 * g11 + m12 * g21 + m13 * g31]
3836         // [gg21]   [            m22 * g21 + m23 * g31]
3837         // [gg31]   [                        m33 * g31]
3838         // [gg12]   [m11 * g12 + m12 * g22 + m13 * g32]
3839         // [gg22]   [            m22 * g22 + m23 * g32]
3840         // [gg32]   [                        m33 * g32]
3841         // [gg13]   [m11 * g13 + m12 * g23 + m13 * g33]
3842         // [gg23]   [            m22 * g23 + m23 * g33]
3843         // [gg33]   [                        m33 * g33]
3844 
3845         // Then the Jacobian of F(b, M, G) is:
3846         // J = [m11  m12  m13  bx   by   0    bz   0    0    0    0    0    0    0    0    0    0    0  ]
3847         //     [0    m22  m23  0    0    by   0    bz   0    0    0    0    0    0    0    0    0    0  ]
3848         //     [0    0    m33  0    0    0    0    0    bz   0    0    0    0    0    0    0    0    0  ]
3849         //     [0    0    0    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0  ]
3850         //     [0    0    0    0    0    1    0    0    0    0    0    0    0    0    0    0    0    0  ]
3851         //     [0    0    0    0    0    0    0    0    1    0    0    0    0    0    0    0    0    0  ]
3852         //     [0    0    0    0    1    0    0    0    0    0    0    0    0    0    0    0    0    0  ]
3853         //     [0    0    0    0    0    0    1    0    0    0    0    0    0    0    0    0    0    0  ]
3854         //     [0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0  ]
3855         //     [0    0    0    0    0    0    0    1    0    0    0    0    0    0    0    0    0    0  ]
3856         //     [0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0  ]
3857         //     [0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0  ]
3858         //     [0    0    0    g11  g21  0    g31  0    0    m11  m12  m13  0    0    0    0    0    0  ]
3859         //     [0    0    0    0    0    g21  0    g31  0    0    m22  m23  0    0    0    0    0    0  ]
3860         //     [0    0    0    0    0    0    0    0    g31  0    0    m33  0    0    0    0    0    0  ]
3861         //     [0    0    0    g12  g22  0    g32  0    0    0    0    0    m11  m12  m13  0    0    0  ]
3862         //     [0    0    0    0    0    g22  0    g32  0    0    0    0    0    m22  m23  0    0    0  ]
3863         //     [0    0    0    0    0    0    0    0    g32  0    0    0    0    0    m33  0    0    0  ]
3864         //     [0    0    0    g13  g23  0    g33  0    0    0    0    0    0    0    0    m11  m12  m13]
3865         //     [0    0    0    0    0    g23  0    g33  0    0    0    0    0    0    0    0    m22  m23]
3866         //     [0    0    0    0    0    0    0    0    g33  0    0    0    0    0    0    0    0    m33]
3867 
3868         // We know that the propagated covariance is J * Cov * J', hence:
3869         final var jacobian = new Matrix(GENERAL_UNKNOWNS_AND_CROSS_BIASES, COMMON_Z_AXIS_UNKNOWNS_AND_CROSS_BIASES);
3870 
3871         jacobian.setElementAt(0, 0, m11);
3872         jacobian.setElementAt(0, 1, m12);
3873         jacobian.setElementAt(0, 2, m13);
3874         jacobian.setElementAt(0, 3, bx);
3875         jacobian.setElementAt(0, 4, by);
3876         jacobian.setElementAt(0, 6, bz);
3877 
3878         jacobian.setElementAt(1, 1, m22);
3879         jacobian.setElementAt(1, 2, m23);
3880         jacobian.setElementAt(1, 5, by);
3881         jacobian.setElementAt(1, 7, bz);
3882 
3883         jacobian.setElementAt(2, 2, m33);
3884         jacobian.setElementAt(2, 8, bz);
3885 
3886         jacobian.setElementAt(3, 3, 1.0);
3887         jacobian.setElementAt(4, 5, 1.0);
3888         jacobian.setElementAt(5, 8, 1.0);
3889 
3890         jacobian.setElementAt(6, 4, 1.0);
3891         jacobian.setElementAt(7, 6, 1.0);
3892 
3893         jacobian.setElementAt(9, 7, 1.0);
3894 
3895         jacobian.setElementAt(12, 3, g11);
3896         jacobian.setElementAt(12, 4, g21);
3897         jacobian.setElementAt(12, 6, g31);
3898         jacobian.setElementAt(12, 9, m11);
3899         jacobian.setElementAt(12, 10, m12);
3900         jacobian.setElementAt(12, 11, m13);
3901 
3902         jacobian.setElementAt(13, 5, g21);
3903         jacobian.setElementAt(13, 7, g31);
3904         jacobian.setElementAt(13, 10, m22);
3905         jacobian.setElementAt(13, 11, m23);
3906 
3907         jacobian.setElementAt(14, 8, g31);
3908         jacobian.setElementAt(14, 11, m33);
3909 
3910         jacobian.setElementAt(15, 3, g12);
3911         jacobian.setElementAt(15, 4, g22);
3912         jacobian.setElementAt(15, 6, g32);
3913         jacobian.setElementAt(15, 12, m11);
3914         jacobian.setElementAt(15, 13, m12);
3915         jacobian.setElementAt(15, 14, m13);
3916 
3917         jacobian.setElementAt(16, 5, g22);
3918         jacobian.setElementAt(16, 7, g32);
3919         jacobian.setElementAt(16, 13, m22);
3920         jacobian.setElementAt(16, 14, m23);
3921 
3922         jacobian.setElementAt(17, 8, g32);
3923         jacobian.setElementAt(17, 14, m33);
3924 
3925         jacobian.setElementAt(18, 3, g13);
3926         jacobian.setElementAt(18, 4, g23);
3927         jacobian.setElementAt(18, 6, g33);
3928         jacobian.setElementAt(18, 15, m11);
3929         jacobian.setElementAt(18, 16, m12);
3930         jacobian.setElementAt(18, 17, m13);
3931 
3932         jacobian.setElementAt(19, 5, g23);
3933         jacobian.setElementAt(19, 7, g33);
3934         jacobian.setElementAt(19, 16, m22);
3935         jacobian.setElementAt(19, 17, m23);
3936 
3937         jacobian.setElementAt(20, 8, g33);
3938         jacobian.setElementAt(20, 17, m33);
3939 
3940         final var jacobianTrans = jacobian.transposeAndReturnNew();
3941         jacobian.multiply(estimatedCovariance);
3942         jacobian.multiply(jacobianTrans);
3943         estimatedCovariance = jacobian;
3944     }
3945 
3946     /**
3947      * Internal method to perform general calibration when G-dependent cross
3948      * biases are being estimated.
3949      *
3950      * @throws AlgebraException                         if accelerometer parameters prevent fixing
3951      *                                                  measured accelerometer values due to numerical instabilities.
3952      * @throws FittingException                         if no convergence to solution is found.
3953      * @throws com.irurueta.numerical.NotReadyException if fitter is not ready.
3954      */
3955     private void calibrateGeneralAndGDependentCrossBiases() throws AlgebraException, FittingException,
3956             com.irurueta.numerical.NotReadyException {
3957         // The gyroscope model is
3958         // Ωmeas = bg + (I + Mg) * Ωtrue + Gg * ftrue + w
3959 
3960         // Ideally a least squares solution tries to minimize noise component, so:
3961         // Ωmeas = bg + (I + Mg) * Ωtrue + Gg * ftrue
3962 
3963         // For convergence purposes of the Levenberg-Marquardt algorithm, we
3964         // take common factor M = I + Mg
3965 
3966         // and the gyroscope model can be better expressed as:
3967 
3968         // Ωmeas = M*(Ωtrue + b + G * ftrue)
3969 
3970         // where:
3971         // bg = M*b --> b = M^-1*bg
3972         // Gg = M*G --> G = M^-1*Gg
3973 
3974         // Hence
3975         // Ωmeas - M*b - M*G*ftrue = M*Ωtrue
3976         // M^-1 * (Ωmeas - M*b - M*G*ftrue) = Ωtrue
3977 
3978         final var gradientEstimator = new GradientEstimator(
3979                 params -> evaluateGeneralWithGDependentCrossBiases(i, params));
3980 
3981         final var initialM = Matrix.identity(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
3982         initialM.add(getInitialMg());
3983 
3984         final var invInitialM = Utils.inverse(initialM);
3985         final var initialBg = getInitialBiasAsMatrix();
3986         final var initialB = invInitialM.multiplyAndReturnNew(initialBg);
3987         final var initialGg = getInitialGg();
3988         final var initialG = invInitialM.multiplyAndReturnNew(initialGg);
3989 
3990         fitter.setFunctionEvaluator(new LevenbergMarquardtMultiDimensionFunctionEvaluator() {
3991             @Override
3992             public int getNumberOfDimensions() {
3993                 // Before and after normalized gravity versors
3994                 return 2 * BodyKinematics.COMPONENTS;
3995             }
3996 
3997             @Override
3998             public double[] createInitialParametersArray() {
3999                 final var initial = new double[GENERAL_UNKNOWNS_AND_CROSS_BIASES];
4000 
4001                 // biases b
4002                 for (var i = 0; i < BodyKinematics.COMPONENTS; i++) {
4003                     initial[i] = initialB.getElementAtIndex(i);
4004                 }
4005 
4006                 // cross coupling errors M
4007                 final var num = BodyKinematics.COMPONENTS * BodyKinematics.COMPONENTS;
4008                 for (int i = 0, j = BodyKinematics.COMPONENTS; i < num; i++, j++) {
4009                     initial[j] = initialM.getElementAtIndex(i);
4010                 }
4011 
4012                 // g-dependent cross biases G
4013                 for (int i = 0, j = BodyKinematics.COMPONENTS + num; i < num; i++, j++) {
4014                     initial[j] = initialG.getElementAtIndex(i);
4015                 }
4016 
4017                 return initial;
4018             }
4019 
4020             @Override
4021             public double evaluate(
4022                     final int i, final double[] point, final double[] params, final double[] derivatives)
4023                     throws EvaluationException {
4024                 EasyGyroscopeCalibrator.this.i = i;
4025 
4026                 // point contains fixed gravity versor values for current
4027                 // sequence
4028                 EasyGyroscopeCalibrator.this.point = point;
4029 
4030                 gradientEstimator.gradient(params, derivatives);
4031 
4032                 return evaluateGeneralWithGDependentCrossBiases(i, params);
4033             }
4034         });
4035 
4036         setInputData();
4037 
4038         fitter.fit();
4039 
4040         final var result = fitter.getA();
4041 
4042         final var bx = result[0];
4043         final var by = result[1];
4044         final var bz = result[2];
4045 
4046         final var m11 = result[3];
4047         final var m21 = result[4];
4048         final var m31 = result[5];
4049 
4050         final var m12 = result[6];
4051         final var m22 = result[7];
4052         final var m32 = result[8];
4053 
4054         final var m13 = result[9];
4055         final var m23 = result[10];
4056         final var m33 = result[11];
4057 
4058         final var g11 = result[12];
4059         final var g21 = result[13];
4060         final var g31 = result[14];
4061 
4062         final var g12 = result[15];
4063         final var g22 = result[16];
4064         final var g32 = result[17];
4065 
4066         final var g13 = result[18];
4067         final var g23 = result[19];
4068         final var g33 = result[20];
4069 
4070         final var b = new Matrix(BodyKinematics.COMPONENTS, 1);
4071         b.setElementAtIndex(0, bx);
4072         b.setElementAtIndex(1, by);
4073         b.setElementAtIndex(2, bz);
4074 
4075         final var m = new Matrix(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
4076         m.setElementAtIndex(0, m11);
4077         m.setElementAtIndex(1, m21);
4078         m.setElementAtIndex(2, m31);
4079 
4080         m.setElementAtIndex(3, m12);
4081         m.setElementAtIndex(4, m22);
4082         m.setElementAtIndex(5, m32);
4083 
4084         m.setElementAtIndex(6, m13);
4085         m.setElementAtIndex(7, m23);
4086         m.setElementAtIndex(8, m33);
4087 
4088         final var g = new Matrix(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
4089         g.setElementAtIndex(0, g11);
4090         g.setElementAtIndex(1, g21);
4091         g.setElementAtIndex(2, g31);
4092 
4093         g.setElementAtIndex(3, g12);
4094         g.setElementAtIndex(4, g22);
4095         g.setElementAtIndex(5, g32);
4096 
4097         g.setElementAtIndex(6, g13);
4098         g.setElementAtIndex(7, g23);
4099         g.setElementAtIndex(8, g33);
4100 
4101         setResult(m, b, g);
4102 
4103         // at this point covariance is expressed in terms of b, M and G, and must
4104         // be expressed in terms of bg, Mg and Gg.
4105         // We know that:
4106         // bg = M * b
4107         // Mg = M - I
4108         // Gg = M * G
4109 
4110         // b = [bx]
4111         //     [by]
4112         //     [bz]
4113 
4114         // M = [m11  m12  m13]
4115         //     [m21  m22  m23]
4116         //     [m31  m32  m33]
4117 
4118         // G = [g11  g12  g13]
4119         //     [g21  g22  g23]
4120         //     [g31  g32  g33]
4121 
4122         // bg = [m11  m12  m13][bx] = [m11 * bx + m12 * by + m13 * bz] = [bgx]
4123         //      [m21  m22  m23][by]   [m21 * bx + m22 * by + m23 * bz]   [bgy]
4124         //      [m31  m32  m33][bz]   [m31 * bx + m32 * by + m33 * bz]   [bgz]
4125 
4126         // Mg = [sx  mxy  mxz] = [m11 - 1    m12         m13     ]
4127         //      [myx	sy	 myz]   [m21        m22 - 1     m23     ]
4128         //      [mzx	mzy  sz ]   [m31        m32         m33 - 1 ]
4129 
4130         // Gg = [gg11  gg12  gg13] = [m11  m12  m13][g11  g12  g13]
4131         //      [gg21  gg22  gg23]   [m21  m22  m23][g21  g22  g23]
4132         //      [gg31  gg32  gg33]   [m31  m32  m33][g31  g32  g33]
4133 
4134         // Defining the linear application:
4135         // F(b, M, G) = F(bx, by, bz, m11, m21, m31, m12, m22, m32, m13, m23, m33, g11, g21, g31, g12, g22, g32, g13, g23, g33)
4136         // as:
4137         // [bgx] =  [m11 * bx + m12 * by + m13 * bz]
4138         // [bgy]    [m21 * bx + m22 * by + m23 * bz]
4139         // [bgz]    [m31 * bx + m32 * by + m33 * bz]
4140         // [sx]     [m11 - 1]
4141         // [sy]     [m22 - 1]
4142         // [sz]     [m33 - 1]
4143         // [mxy]    [m12]
4144         // [mxz]    [m13]
4145         // [myx]    [m21]
4146         // [myz]    [m23]
4147         // [mzx]    [m31]
4148         // [mzy]    [m32]
4149         // [gg11]   [m11 * g11 + m12 * g21 + m13 * g31]
4150         // [gg21]   [m21 * g11 + m22 * g21 + m23 * g31]
4151         // [gg31]   [m31 * g11 + m32 * g21 + m33 * g31]
4152         // [gg12]   [m11 * g12 + m12 * g22 + m13 * g32]
4153         // [gg22]   [m21 * g12 + m22 * g22 + m23 * g32]
4154         // [gg32]   [m31 * g12 + m32 * g22 + m33 * g32]
4155         // [gg13]   [m11 * g13 + m12 * g23 + m13 * g33]
4156         // [gg23]   [m21 * g13 + m22 * g23 + m23 * g33]
4157         // [gg33]   [m31 * g13 + m32 * g23 + m33 * g33]
4158 
4159         // Then the Jacobian of F(b, M, G) is:
4160         // J = [m11  m12  m13  bx   0    0    by   0    0    bz   0    0    0    0    0    0    0    0    0    0    0  ]
4161         //     [m21  m22  m23  0    bx   0    0    by   0    0    bz   0    0    0    0    0    0    0    0    0    0  ]
4162         //     [m31  m32  m33  0    0    bx   0    0    by   0    0    bz   0    0    0    0    0    0    0    0    0  ]
4163         //     [0    0    0    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0  ]
4164         //     [0    0    0    0    0    0    0    1    0    0    0    0    0    0    0    0    0    0    0    0    0  ]
4165         //     [0    0    0    0    0    0    0    0    0    0    0    1    0    0    0    0    0    0    0    0    0  ]
4166         //     [0    0    0    0    0    0    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0  ]
4167         //     [0    0    0    0    0    0    0    0    0    1    0    0    0    0    0    0    0    0    0    0    0  ]
4168         //     [0    0    0    0    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0  ]
4169         //     [0    0    0    0    0    0    0    0    0    0    1    0    0    0    0    0    0    0    0    0    0  ]
4170         //     [0    0    0    0    0    1    0    0    0    0    0    0    0    0    0    0    0    0    0    0    0  ]
4171         //     [0    0    0    0    0    0    0    0    1    0    0    0    0    0    0    0    0    0    0    0    0  ]
4172         //     [0    0    0    g11  0    0    g21  0    0    g31  0    0    m11  m12  m13  0    0    0    0    0    0  ]
4173         //     [0    0    0    0    g11  0    0    g21  0    0    g31  0    m21  m22  m23  0    0    0    0    0    0  ]
4174         //     [0    0    0    0    0    g11  0    0    g21  0    0    g31  m31  m32  m33  0    0    0    0    0    0  ]
4175         //     [0    0    0    g12  0    0    g22  0    0    g32  0    0    0    0    0    m11  m12  m13  0    0    0  ]
4176         //     [0    0    0    0    g12  0    0    g22  0    0    g32  0    0    0    0    m21  m22  m23  0    0    0  ]
4177         //     [0    0    0    0    0    g12  0    0    g22  0    0    g32  0    0    0    m31  m32  m33  0    0    0  ]
4178         //     [0    0    0    g13  0    0    g23  0    0    g33  0    0    0    0    0    0    0    0    m11  m12  m13]
4179         //     [0    0    0    0    g13  0    0    g23  0    0    g33  0    0    0    0    0    0    0    m21  m22  m23]
4180         //     [0    0    0    0    0    g13  0    0    g23  0    0    g33  0    0    0    0    0    0    m31  m32  m33]
4181 
4182         // We know that the propagated covariance is J * Cov * J', hence:
4183         final var jacobian = new Matrix(GENERAL_UNKNOWNS_AND_CROSS_BIASES, GENERAL_UNKNOWNS_AND_CROSS_BIASES);
4184 
4185         jacobian.setElementAt(0, 0, m11);
4186         jacobian.setElementAt(0, 1, m12);
4187         jacobian.setElementAt(0, 2, m13);
4188         jacobian.setElementAt(0, 3, bx);
4189         jacobian.setElementAt(0, 6, by);
4190         jacobian.setElementAt(0, 9, bz);
4191 
4192         jacobian.setElementAt(1, 0, m21);
4193         jacobian.setElementAt(1, 1, m22);
4194         jacobian.setElementAt(1, 2, m23);
4195         jacobian.setElementAt(1, 4, bx);
4196         jacobian.setElementAt(1, 7, by);
4197         jacobian.setElementAt(1, 10, bz);
4198 
4199         jacobian.setElementAt(2, 0, m31);
4200         jacobian.setElementAt(2, 1, m32);
4201         jacobian.setElementAt(2, 2, m33);
4202         jacobian.setElementAt(2, 5, bx);
4203         jacobian.setElementAt(2, 8, by);
4204         jacobian.setElementAt(2, 11, bz);
4205 
4206         jacobian.setElementAt(3, 3, 1.0);
4207         jacobian.setElementAt(4, 7, 1.0);
4208         jacobian.setElementAt(5, 11, 1.0);
4209 
4210         jacobian.setElementAt(6, 6, 1.0);
4211         jacobian.setElementAt(7, 9, 1.0);
4212         jacobian.setElementAt(8, 4, 1.0);
4213 
4214         jacobian.setElementAt(9, 10, 1.0);
4215         jacobian.setElementAt(10, 5, 1.0);
4216         jacobian.setElementAt(11, 8, 1.9);
4217 
4218         jacobian.setElementAt(12, 3, g11);
4219         jacobian.setElementAt(12, 6, g21);
4220         jacobian.setElementAt(12, 9, g31);
4221         jacobian.setElementAt(12, 12, m11);
4222         jacobian.setElementAt(12, 13, m12);
4223         jacobian.setElementAt(12, 14, m13);
4224 
4225         jacobian.setElementAt(13, 4, g11);
4226         jacobian.setElementAt(13, 7, g21);
4227         jacobian.setElementAt(13, 10, g31);
4228         jacobian.setElementAt(13, 12, m21);
4229         jacobian.setElementAt(13, 13, m22);
4230         jacobian.setElementAt(13, 14, m23);
4231 
4232         jacobian.setElementAt(14, 5, g11);
4233         jacobian.setElementAt(14, 8, g21);
4234         jacobian.setElementAt(14, 11, g31);
4235         jacobian.setElementAt(14, 12, m31);
4236         jacobian.setElementAt(14, 13, m32);
4237         jacobian.setElementAt(14, 14, m33);
4238 
4239         jacobian.setElementAt(15, 3, g12);
4240         jacobian.setElementAt(15, 6, g22);
4241         jacobian.setElementAt(15, 9, g32);
4242         jacobian.setElementAt(15, 15, m11);
4243         jacobian.setElementAt(15, 16, m12);
4244         jacobian.setElementAt(15, 17, m13);
4245 
4246         jacobian.setElementAt(16, 4, g12);
4247         jacobian.setElementAt(16, 7, g22);
4248         jacobian.setElementAt(16, 10, g32);
4249         jacobian.setElementAt(16, 15, m21);
4250         jacobian.setElementAt(16, 16, m22);
4251         jacobian.setElementAt(16, 17, m23);
4252 
4253         jacobian.setElementAt(17, 5, g12);
4254         jacobian.setElementAt(17, 8, g22);
4255         jacobian.setElementAt(17, 11, g32);
4256         jacobian.setElementAt(17, 15, m31);
4257         jacobian.setElementAt(17, 16, m32);
4258         jacobian.setElementAt(17, 17, m33);
4259 
4260         jacobian.setElementAt(18, 3, g13);
4261         jacobian.setElementAt(18, 6, g23);
4262         jacobian.setElementAt(18, 9, g33);
4263         jacobian.setElementAt(18, 18, m11);
4264         jacobian.setElementAt(18, 19, m12);
4265         jacobian.setElementAt(18, 20, m13);
4266 
4267         jacobian.setElementAt(19, 4, g13);
4268         jacobian.setElementAt(19, 7, g23);
4269         jacobian.setElementAt(19, 10, g33);
4270         jacobian.setElementAt(19, 18, m21);
4271         jacobian.setElementAt(19, 19, m22);
4272         jacobian.setElementAt(19, 20, m23);
4273 
4274         jacobian.setElementAt(20, 5, g13);
4275         jacobian.setElementAt(20, 8, g23);
4276         jacobian.setElementAt(20, 11, g33);
4277         jacobian.setElementAt(20, 18, m31);
4278         jacobian.setElementAt(20, 19, m32);
4279         jacobian.setElementAt(20, 20, m33);
4280 
4281         final var jacobianTrans = jacobian.transposeAndReturnNew();
4282         jacobian.multiply(estimatedCovariance);
4283         jacobian.multiply(jacobianTrans);
4284         estimatedCovariance = jacobian;
4285     }
4286 
4287     /**
4288      * Internal method to perform calibration when common z-axis is assumed
4289      * for both the accelerometer and gyroscope and G-dependent cross biases
4290      * are ignored.
4291      *
4292      * @throws AlgebraException                         if accelerometer parameters prevent fixing
4293      *                                                  measured accelerometer values due to numerical instabilities.
4294      * @throws FittingException                         if no convergence to solution is found.
4295      * @throws com.irurueta.numerical.NotReadyException if fitter is not ready.
4296      */
4297     private void calibrateCommonAxis() throws AlgebraException, FittingException,
4298             com.irurueta.numerical.NotReadyException {
4299         // The gyroscope model is
4300         // Ωmeas = bg + (I + Mg) * Ωtrue + Gg * ftrue + w
4301 
4302         // Ideally a least squares solution tries to minimize noise component, so:
4303         // Ωmeas = bg + (I + Mg) * Ωtrue + Gg * ftrue
4304 
4305         // Where Mg is upper triangular
4306 
4307         // Since G-dependent cross biases are ignored, we can assume that Gg = 0
4308 
4309         // Hence:
4310         // Ωmeas = bg + (I + Mg) * Ωtrue
4311 
4312         // For convergence purposes of the Levenberg-Marquardt algorithm,
4313         // the gyroscope model can be better expressed as:
4314         // Ωmeas = T*K*(Ωtrue + b)
4315         // Ωmeas = M*(Ωtrue + b)
4316         // Ωmeas = M*Ωtrue + M*b
4317 
4318         // Hence
4319         // Ωmeas - M*b = M*Ωtrue
4320 
4321         // M^-1 * (Ωmeas - M*b) = Ωtrue
4322 
4323         // where:
4324         // M = I + Mg
4325         // bg = M*b = (I + Mg)*b --> b = M^-1*bg
4326 
4327         // Notice that M is upper diagonal because Mg is upper diagonal
4328         // when common axis is assumed
4329 
4330         final var gradientEstimator = new GradientEstimator(params -> evaluateCommonAxis(i, params));
4331 
4332         final var initialM = Matrix.identity(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
4333         initialM.add(getInitialMg());
4334 
4335         // Force initial M to be upper diagonal
4336         initialM.setElementAt(1, 0, 0.0);
4337         initialM.setElementAt(2, 0, 0.0);
4338         initialM.setElementAt(2, 1, 0.0);
4339 
4340         final var invInitialM = Utils.inverse(initialM);
4341         final var initialBg = getInitialBiasAsMatrix();
4342         final var initialB = invInitialM.multiplyAndReturnNew(initialBg);
4343 
4344         fitter.setFunctionEvaluator(new LevenbergMarquardtMultiDimensionFunctionEvaluator() {
4345             @Override
4346             public int getNumberOfDimensions() {
4347                 // Before and after normalized gravity versors
4348                 return 2 * BodyKinematics.COMPONENTS;
4349             }
4350 
4351             @Override
4352             public double[] createInitialParametersArray() {
4353                 final var initial = new double[COMMON_Z_AXIS_UNKNOWNS];
4354 
4355                 // biases b
4356                 for (var i = 0; i < BodyKinematics.COMPONENTS; i++) {
4357                     initial[i] = initialB.getElementAtIndex(i);
4358                 }
4359 
4360                 // upper diagonal cross coupling errors M
4361                 var k = BodyKinematics.COMPONENTS;
4362                 for (var j = 0; j < BodyKinematics.COMPONENTS; j++) {
4363                     for (var i = 0; i < BodyKinematics.COMPONENTS; i++) {
4364                         if (i <= j) {
4365                             initial[k] = initialM.getElementAt(i, j);
4366                             k++;
4367                         }
4368                     }
4369                 }
4370 
4371                 return initial;
4372             }
4373 
4374             @Override
4375             public double evaluate(
4376                     final int i, final double[] point, final double[] params, final double[] derivatives)
4377                     throws EvaluationException {
4378                 EasyGyroscopeCalibrator.this.i = i;
4379 
4380                 // point contains fixed gravity versor values for current
4381                 // sequence
4382                 EasyGyroscopeCalibrator.this.point = point;
4383 
4384                 gradientEstimator.gradient(params, derivatives);
4385 
4386                 return evaluateCommonAxis(i, params);
4387             }
4388         });
4389 
4390         setInputData();
4391 
4392         fitter.fit();
4393 
4394         final var result = fitter.getA();
4395 
4396         final var bx = result[0];
4397         final var by = result[1];
4398         final var bz = result[2];
4399 
4400         final var m11 = result[3];
4401 
4402         final var m12 = result[4];
4403         final var m22 = result[5];
4404 
4405         final var m13 = result[6];
4406         final var m23 = result[7];
4407         final var m33 = result[8];
4408 
4409         final var b = new Matrix(BodyKinematics.COMPONENTS, 1);
4410         b.setElementAtIndex(0, bx);
4411         b.setElementAtIndex(1, by);
4412         b.setElementAtIndex(2, bz);
4413 
4414         final var m = new Matrix(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
4415         m.setElementAtIndex(0, m11);
4416         m.setElementAtIndex(1, 0.0);
4417         m.setElementAtIndex(2, 0.0);
4418 
4419         m.setElementAtIndex(3, m12);
4420         m.setElementAtIndex(4, m22);
4421         m.setElementAtIndex(5, 0.0);
4422 
4423         m.setElementAtIndex(6, m13);
4424         m.setElementAtIndex(7, m23);
4425         m.setElementAtIndex(8, m33);
4426 
4427         setResult(m, b);
4428 
4429         // at this point covariance is expressed in terms of b, M and G, and must
4430         // be expressed in terms of bg, Mg and Gg.
4431         // We know that:
4432         // bg = M * b
4433         // Mg = M - I
4434         // Gg = M * G = 0
4435 
4436         // b = [bx]
4437         //     [by]
4438         //     [bz]
4439 
4440         // M = [m11  m12  m13]
4441         //     [0    m22  m23]
4442         //     [0    0    m33]
4443 
4444         // G = [g11  g12  g13] = 0
4445         //     [g21  g22  g23]
4446         //     [g31  g32  g33]
4447 
4448         // bg = [m11  m12  m13][bx] = [m11 * bx + m12 * by + m13 * bz] = [bgx]
4449         //      [0    m22  m23][by]   [           m22 * by + m23 * bz]   [bgy]
4450         //      [0    0    m33][bz]   [                      m33 * bz]   [bgz]
4451 
4452         // Mg = [sx  mxy  mxz] = [m11 - 1    m12         m13     ]
4453         //      [myx sy	  myz]   [0          m22 - 1     m23     ]
4454         //      [mzx mzy  sz ]   [0          0           m33 - 1 ]
4455 
4456         // Gg = [gg11  gg12  gg13] = 0
4457         //      [gg21  gg22  gg23]
4458         //      [gg31  gg32  gg33]
4459 
4460         // Defining the linear application:
4461         // F(b, M) = F(bx, by, bz, m11, m12, m22, m13, m23, m33)
4462         // as:
4463         // [bgx] =  [m11 * bx + m12 * by + m13 * bz]
4464         // [bgy]    [           m22 * by + m23 * bz]
4465         // [bgz]    [                      m33 * bz]
4466         // [sx]     [m11 - 1]
4467         // [sy]     [m22 - 1]
4468         // [sz]     [m33 - 1]
4469         // [mxy]    [m12]
4470         // [mxz]    [m13]
4471         // [myx]    [0]
4472         // [myz]    [m23]
4473         // [mzx]    [0]
4474         // [mzy]    [0]
4475         // [gg11]   [0]
4476         // [gg21]   [0]
4477         // [gg31]   [0]
4478         // [gg12]   [0]
4479         // [gg22]   [0]
4480         // [gg32]   [0]
4481         // [gg13]   [0]
4482         // [gg23]   [0]
4483         // [gg33]   [0]
4484 
4485         // Then the Jacobian of F(b, M) is:
4486         // J = [m11  m12  m13  bx   by   0    bz   0    0  ]
4487         //     [0    m22  m23  0    0    by   0    bz   0  ]
4488         //     [0    0    m33  0    0    0    0    0    bz ]
4489         //     [0    0    0    1    0    0    0    0    0  ]
4490         //     [0    0    0    0    0    1    0    0    0  ]
4491         //     [0    0    0    0    0    0    0    0    1  ]
4492         //     [0    0    0    0    1    0    0    0    0  ]
4493         //     [0    0    0    0    0    0    1    0    0  ]
4494         //     [0    0    0    0    0    0    0    0    0  ]
4495         //     [0    0    0    0    0    0    0    1    0  ]
4496         //     [0    0    0    0    0    0    0    0    0  ]
4497         //     [0    0    0    0    0    0    0    0    0  ]
4498         //     [0    0    0    0    0    0    0    0    0  ]
4499         //     [0    0    0    0    0    0    0    0    0  ]
4500         //     [0    0    0    0    0    0    0    0    0  ]
4501         //     [0    0    0    0    0    0    0    0    0  ]
4502         //     [0    0    0    0    0    0    0    0    0  ]
4503         //     [0    0    0    0    0    0    0    0    0  ]
4504         //     [0    0    0    0    0    0    0    0    0  ]
4505         //     [0    0    0    0    0    0    0    0    0  ]
4506         //     [0    0    0    0    0    0    0    0    0  ]
4507 
4508         // We know that the propagated covariance is J * Cov * J', hence:
4509         final var jacobian = new Matrix(GENERAL_UNKNOWNS_AND_CROSS_BIASES, COMMON_Z_AXIS_UNKNOWNS);
4510 
4511         jacobian.setElementAt(0, 0, m11);
4512         jacobian.setElementAt(0, 1, m12);
4513         jacobian.setElementAt(0, 2, m13);
4514         jacobian.setElementAt(0, 3, bx);
4515         jacobian.setElementAt(0, 4, by);
4516         jacobian.setElementAt(0, 6, bz);
4517 
4518         jacobian.setElementAt(1, 1, m22);
4519         jacobian.setElementAt(1, 2, m23);
4520         jacobian.setElementAt(1, 5, by);
4521         jacobian.setElementAt(1, 7, bz);
4522 
4523         jacobian.setElementAt(2, 2, m33);
4524         jacobian.setElementAt(2, 8, bz);
4525 
4526         jacobian.setElementAt(3, 3, 1.0);
4527         jacobian.setElementAt(4, 5, 1.0);
4528         jacobian.setElementAt(5, 8, 1.0);
4529 
4530         jacobian.setElementAt(6, 4, 1.0);
4531         jacobian.setElementAt(7, 6, 1.0);
4532 
4533         jacobian.setElementAt(9, 7, 1.0);
4534 
4535         final var jacobianTrans = jacobian.transposeAndReturnNew();
4536         jacobian.multiply(estimatedCovariance);
4537         jacobian.multiply(jacobianTrans);
4538         estimatedCovariance = jacobian;
4539     }
4540 
4541     /**
4542      * Internal method to perform general calibration when G-dependant cross
4543      * biases are ignored.
4544      *
4545      * @throws AlgebraException                         if accelerometer parameters prevent fixing
4546      *                                                  measured accelerometer values due to numerical instabilities.
4547      * @throws FittingException                         if no convergence to solution is found.
4548      * @throws com.irurueta.numerical.NotReadyException if fitter is not ready.
4549      */
4550     private void calibrateGeneral() throws AlgebraException, FittingException,
4551             com.irurueta.numerical.NotReadyException {
4552         // The gyroscope model is
4553         // Ωmeas = bg + (I + Mg) * Ωtrue + Gg * ftrue + w
4554 
4555         // Ideally a least squares solution tries to minimize noise component, so:
4556         // Ωmeas = bg + (I + Mg) * Ωtrue + Gg * ftrue
4557 
4558         // Since G-dependent cross biases are ignored, we can assume that Gg = 0
4559 
4560         // Hence:
4561         // Ωmeas = bg + (I + Mg) * Ωtrue
4562 
4563         // For convergence purposes of the Levenberg-Marquardt algorithm, the
4564         // gyroscope model can be better expressed as:
4565         // Ωmeas = T*K*(Ωtrue + b)
4566         // Ωmeas = M*(Ωtrue + b)
4567         // Ωmeas = M*Ωtrue + M*b
4568 
4569         // where:
4570         // M = I + Mg
4571         // bg = M*b = (I + Mg)*b --> b = M^-1*bg
4572 
4573         // Hence
4574         // Ωmeas - M*b = M*Ωtrue
4575 
4576         // M^-1 * (Ωmeas - M*b) = Ωtrue
4577 
4578         final var gradientEstimator = new GradientEstimator(params -> evaluateGeneral(i, params));
4579 
4580         final var initialM = Matrix.identity(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
4581         initialM.add(getInitialMg());
4582 
4583         final var invInitialM = Utils.inverse(initialM);
4584         final var initialBg = getInitialBiasAsMatrix();
4585         final var initialB = invInitialM.multiplyAndReturnNew(initialBg);
4586 
4587         fitter.setFunctionEvaluator(new LevenbergMarquardtMultiDimensionFunctionEvaluator() {
4588             @Override
4589             public int getNumberOfDimensions() {
4590                 // Before and after normalized gravity versors
4591                 return 2 * BodyKinematics.COMPONENTS;
4592             }
4593 
4594             @Override
4595             public double[] createInitialParametersArray() {
4596                 final var initial = new double[GENERAL_UNKNOWNS];
4597 
4598                 // biases b
4599                 for (var i = 0; i < BodyKinematics.COMPONENTS; i++) {
4600                     initial[i] = initialB.getElementAtIndex(i);
4601                 }
4602 
4603                 // cross coupling errors M
4604                 final var num = BodyKinematics.COMPONENTS * BodyKinematics.COMPONENTS;
4605                 for (int i = 0, j = BodyKinematics.COMPONENTS; i < num; i++, j++) {
4606                     initial[j] = initialM.getElementAtIndex(i);
4607                 }
4608 
4609                 return initial;
4610             }
4611 
4612             @Override
4613             public double evaluate(
4614                     final int i, final double[] point, final double[] params, final double[] derivatives)
4615                     throws EvaluationException {
4616                 EasyGyroscopeCalibrator.this.i = i;
4617 
4618                 // point contains fixed gravity versor values for current
4619                 // sequence
4620                 EasyGyroscopeCalibrator.this.point = point;
4621 
4622                 gradientEstimator.gradient(params, derivatives);
4623 
4624                 return evaluateGeneral(i, params);
4625             }
4626         });
4627 
4628         setInputData();
4629 
4630         fitter.fit();
4631 
4632         final var result = fitter.getA();
4633 
4634         final var bx = result[0];
4635         final var by = result[1];
4636         final var bz = result[2];
4637 
4638         final var m11 = result[3];
4639         final var m21 = result[4];
4640         final var m31 = result[5];
4641 
4642         final var m12 = result[6];
4643         final var m22 = result[7];
4644         final var m32 = result[8];
4645 
4646         final var m13 = result[9];
4647         final var m23 = result[10];
4648         final var m33 = result[11];
4649 
4650         final var b = new Matrix(BodyKinematics.COMPONENTS, 1);
4651         b.setElementAtIndex(0, bx);
4652         b.setElementAtIndex(1, by);
4653         b.setElementAtIndex(2, bz);
4654 
4655         final var m = new Matrix(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
4656         m.setElementAtIndex(0, m11);
4657         m.setElementAtIndex(1, m21);
4658         m.setElementAtIndex(2, m31);
4659 
4660         m.setElementAtIndex(3, m12);
4661         m.setElementAtIndex(4, m22);
4662         m.setElementAtIndex(5, m32);
4663 
4664         m.setElementAtIndex(6, m13);
4665         m.setElementAtIndex(7, m23);
4666         m.setElementAtIndex(8, m33);
4667 
4668         setResult(m, b);
4669 
4670         // at this point covariance is expressed in terms of b, M and G, and must
4671         // be expressed in terms of bg, Mg and Gg.
4672         // We know that:
4673         // bg = M * b
4674         // Mg = M - I
4675         // Gg = M * G = 0
4676 
4677         // b = [bx]
4678         //     [by]
4679         //     [bz]
4680 
4681         // M = [m11  m12  m13]
4682         //     [m21  m22  m23]
4683         //     [m31  m32  m33]
4684 
4685         // G = [g11  g12  g13] = 0
4686         //     [g21  g22  g23]
4687         //     [g31  g32  g33]
4688 
4689         // bg = [m11  m12  m13][bx] = [m11 * bx + m12 * by + m13 * bz] = [bgx]
4690         //      [m21  m22  m23][by]   [m21 * bx + m22 * by + m23 * bz]   [bgy]
4691         //      [m31  m32  m33][bz]   [m31 * bx + m32 * by + m33 * bz]   [bgz]
4692 
4693         // Mg = [sx  mxy  mxz] = [m11 - 1    m12         m13     ]
4694         //      [myx sy	  myz]   [m21        m22 - 1     m23     ]
4695         //      [mzx mzy  sz ]   [m31        m32         m33 - 1 ]
4696 
4697         // Gg = [gg11  gg12  gg13] = 0
4698         //      [gg21  gg22  gg23]
4699         //      [gg31  gg32  gg33]
4700 
4701         // Defining the linear application:
4702         // F(b, M) = F(bx, by, bz, m11, m21, m31, m12, m22, m32, m13, m23, m33)
4703         // as:
4704         // [bgx] =  [m11 * bx + m12 * by + m13 * bz]
4705         // [bgy]    [m21 * bx + m22 * by + m23 * bz]
4706         // [bgz]    [m31 * bx + m32 * by + m33 * bz]
4707         // [sx]     [m11 - 1]
4708         // [sy]     [m22 - 1]
4709         // [sz]     [m33 - 1]
4710         // [mxy]    [m12]
4711         // [mxz]    [m13]
4712         // [myx]    [m21]
4713         // [myz]    [m23]
4714         // [mzx]    [m31]
4715         // [mzy]    [m32]
4716         // [gg11]   [0]
4717         // [gg21]   [0]
4718         // [gg31]   [0]
4719         // [gg12]   [0]
4720         // [gg22]   [0]
4721         // [gg32]   [0]
4722         // [gg13]   [0]
4723         // [gg23]   [0]
4724         // [gg33]   [0]
4725 
4726         // Then the Jacobian of F(b, M) is:
4727         // J = [m11  m12  m13  bx   0    0    by   0    0    bz   0    0 ]
4728         //     [m21  m22  m23  0    bx   0    0    by   0    0    bz   0 ]
4729         //     [m31  m32  m33  0    0    bx   0    0    by   0    0    bz]
4730         //     [0    0    0    1    0    0    0    0    0    0    0    0 ]
4731         //     [0    0    0    0    0    0    0    1    0    0    0    0 ]
4732         //     [0    0    0    0    0    0    0    0    0    0    0    1 ]
4733         //     [0    0    0    0    0    0    1    0    0    0    0    0 ]
4734         //     [0    0    0    0    0    0    0    0    0    1    0    0 ]
4735         //     [0    0    0    0    1    0    0    0    0    0    0    0 ]
4736         //     [0    0    0    0    0    0    0    0    0    0    1    0 ]
4737         //     [0    0    0    0    0    1    0    0    0    0    0    0 ]
4738         //     [0    0    0    0    0    0    0    0    1    0    0    0 ]
4739         //     [0    0    0    0    0    0    0    0    0    0    0    0 ]
4740         //     [0    0    0    0    0    0    0    0    0    0    0    0 ]
4741         //     [0    0    0    0    0    0    0    0    0    0    0    0 ]
4742         //     [0    0    0    0    0    0    0    0    0    0    0    0 ]
4743         //     [0    0    0    0    0    0    0    0    0    0    0    0 ]
4744         //     [0    0    0    0    0    0    0    0    0    0    0    0 ]
4745         //     [0    0    0    0    0    0    0    0    0    0    0    0 ]
4746         //     [0    0    0    0    0    0    0    0    0    0    0    0 ]
4747         //     [0    0    0    0    0    0    0    0    0    0    0    0 ]
4748 
4749         // We know that the propagated covariance is J * Cov * J', hence:
4750         final var jacobian = new Matrix(GENERAL_UNKNOWNS_AND_CROSS_BIASES, GENERAL_UNKNOWNS);
4751 
4752         jacobian.setElementAt(0, 0, m11);
4753         jacobian.setElementAt(0, 1, m12);
4754         jacobian.setElementAt(0, 2, m13);
4755         jacobian.setElementAt(0, 3, bx);
4756         jacobian.setElementAt(0, 6, by);
4757         jacobian.setElementAt(0, 9, bz);
4758 
4759         jacobian.setElementAt(1, 0, m21);
4760         jacobian.setElementAt(1, 1, m22);
4761         jacobian.setElementAt(1, 2, m23);
4762         jacobian.setElementAt(1, 4, bx);
4763         jacobian.setElementAt(1, 7, by);
4764         jacobian.setElementAt(1, 10, bz);
4765 
4766         jacobian.setElementAt(2, 0, m31);
4767         jacobian.setElementAt(2, 1, m32);
4768         jacobian.setElementAt(2, 2, m33);
4769         jacobian.setElementAt(2, 5, bx);
4770         jacobian.setElementAt(2, 8, by);
4771         jacobian.setElementAt(2, 11, bz);
4772 
4773         jacobian.setElementAt(3, 3, 1.0);
4774         jacobian.setElementAt(4, 7, 1.0);
4775         jacobian.setElementAt(5, 11, 1.0);
4776 
4777         jacobian.setElementAt(6, 6, 1.0);
4778         jacobian.setElementAt(7, 9, 1.0);
4779         jacobian.setElementAt(8, 4, 1.0);
4780 
4781         jacobian.setElementAt(9, 10, 1.0);
4782         jacobian.setElementAt(10, 5, 1.0);
4783         jacobian.setElementAt(11, 8, 1.9);
4784 
4785         final var jacobianTrans = jacobian.transposeAndReturnNew();
4786         jacobian.multiply(estimatedCovariance);
4787         jacobian.multiply(jacobianTrans);
4788         estimatedCovariance = jacobian;
4789     }
4790 
4791     /**
4792      * Sets input data into Levenberg-Marquardt fitter.
4793      *
4794      * @throws AlgebraException if there are numerical instabilities.
4795      */
4796     private void setInputData() throws AlgebraException {
4797 
4798         final var ba = getAccelerometerBiasAsMatrix();
4799         final var ma = getAccelerometerMa();
4800 
4801         final var measuredBeforeF = new double[BodyKinematics.COMPONENTS];
4802         final var fixedBeforeF = new double[BodyKinematics.COMPONENTS];
4803 
4804         final var measuredAfterF = new double[BodyKinematics.COMPONENTS];
4805         final var fixedAfterF = new double[BodyKinematics.COMPONENTS];
4806 
4807         final var numSequences = sequences.size();
4808         final var x = new Matrix(numSequences, 2 * BodyKinematics.COMPONENTS);
4809         final var y = new double[numSequences];
4810         final var standardDeviations = new double[numSequences];
4811 
4812         // make a copy of input sequences that will be used to update
4813         // kinematics measurements with fixed values for memory efficiency
4814 
4815         fixedSequences = new ArrayList<>();
4816         for (final var sequence : sequences) {
4817             fixedSequences.add(new BodyKinematicsSequence<>(sequence));
4818         }
4819 
4820         accelerationFixer.setBias(ba);
4821         accelerationFixer.setCrossCouplingErrors(ma);
4822 
4823         var i = 0;
4824         for (final var sequence : sequences) {
4825             // sequence mean accelerometer samples of previous static
4826             // period will need to be fixed using accelerometer calibration
4827             // parameters
4828             measuredBeforeF[0] = sequence.getBeforeMeanFx();
4829             measuredBeforeF[1] = sequence.getBeforeMeanFy();
4830             measuredBeforeF[2] = sequence.getBeforeMeanFz();
4831             accelerationFixer.fix(measuredBeforeF, fixedBeforeF);
4832 
4833             measuredAfterF[0] = sequence.getAfterMeanFx();
4834             measuredAfterF[1] = sequence.getAfterMeanFy();
4835             measuredAfterF[2] = sequence.getAfterMeanFz();
4836             accelerationFixer.fix(measuredAfterF, fixedAfterF);
4837 
4838             // because we are only interested in gravity direction, we
4839             // normalize these vectors, so that gravity becomes independent
4840             // of current Earth position.
4841             ArrayUtils.normalize(fixedBeforeF);
4842             ArrayUtils.normalize(fixedAfterF);
4843 
4844             x.setSubmatrix(i, 0, i, 2, fixedBeforeF);
4845             x.setSubmatrix(i, 3, i, 5, fixedAfterF);
4846 
4847             y[i] = 0.0;
4848 
4849             standardDeviations[i] = computeAverageAngularRateStandardDeviation(sequence);
4850             i++;
4851         }
4852 
4853         fitter.setInputData(x, y, standardDeviations);
4854     }
4855 
4856     /**
4857      * Computes average angular rate standard deviation for measurements
4858      * in provided sequence.
4859      *
4860      * @param sequence a sequence.
4861      * @return average angular rate standard deviation expressed in radians
4862      * per second (rad/s).
4863      */
4864     private static double computeAverageAngularRateStandardDeviation(
4865             final BodyKinematicsSequence<StandardDeviationTimedBodyKinematics> sequence) {
4866         final var items = sequence.getSortedItems();
4867         final var size = items.size();
4868 
4869         var result = 0.0;
4870         for (final var item : items) {
4871             result += item.getAngularRateStandardDeviation() / size;
4872         }
4873 
4874         return result;
4875     }
4876 
4877     /**
4878      * Converts acceleration instance to meters per squared second.
4879      *
4880      * @param acceleration acceleration instance to be converted.
4881      * @return converted value.
4882      */
4883     private static double convertAcceleration(final Acceleration acceleration) {
4884         return AccelerationConverter.convert(acceleration.getValue().doubleValue(), acceleration.getUnit(),
4885                 AccelerationUnit.METERS_PER_SQUARED_SECOND);
4886     }
4887 
4888     /**
4889      * Converts angular speed value and unit to radians per second.
4890      *
4891      * @param value angular speed value.
4892      * @param unit  unit of angular speed value.
4893      * @return converted value.
4894      */
4895     private static double convertAngularSpeed(final double value, final AngularSpeedUnit unit) {
4896         return AngularSpeedConverter.convert(value, unit, AngularSpeedUnit.RADIANS_PER_SECOND);
4897     }
4898 
4899     /**
4900      * Converts angular speed instance to radians per second.
4901      *
4902      * @param angularSpeed angular speed instance to be converted.
4903      * @return converted value.
4904      */
4905     private static double convertAngularSpeed(final AngularSpeed angularSpeed) {
4906         return convertAngularSpeed(angularSpeed.getValue().doubleValue(), angularSpeed.getUnit());
4907     }
4908 
4909     /**
4910      * Makes proper conversion of internal cross-coupling, bias and g-dependent
4911      * cross bias matrices.
4912      *
4913      * @param m internal scaling and cross-coupling matrix.
4914      * @param b internal bias matrix.
4915      * @param g internal g-dependent cross bias matrix.
4916      * @throws AlgebraException if a numerical instability occurs.
4917      */
4918     private void setResult(final Matrix m, final Matrix b, final Matrix g) throws AlgebraException {
4919         setResult(m, b);
4920 
4921         // Gg = M*G
4922         m.multiply(g, estimatedGg);
4923     }
4924 
4925     /**
4926      * Makes proper conversion of internal cross-coupling and bias matrices.
4927      *
4928      * @param m internal scaling and cross-coupling matrix.
4929      * @param b internal bias matrix.
4930      * @throws AlgebraException if a numerical instability occurs.
4931      */
4932     private void setResult(final Matrix m, final Matrix b) throws AlgebraException {
4933         // Because:
4934         // M = I + Mg
4935         // b = M^-1*bg
4936 
4937         // Then:
4938         // Mg = M - I
4939         // bg = M*b
4940 
4941         if (estimatedBiases == null) {
4942             estimatedBiases = new double[BodyKinematics.COMPONENTS];
4943         }
4944 
4945         final var bg = m.multiplyAndReturnNew(b);
4946         bg.toArray(estimatedBiases);
4947 
4948         if (estimatedMg == null) {
4949             estimatedMg = m;
4950         } else {
4951             estimatedMg.copyFrom(m);
4952         }
4953 
4954         for (var i = 0; i < BodyKinematics.COMPONENTS; i++) {
4955             estimatedMg.setElementAt(i, i, estimatedMg.getElementAt(i, i) - 1.0);
4956         }
4957 
4958         if (estimatedGg == null) {
4959             estimatedGg = new Matrix(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
4960         } else {
4961             estimatedGg.initialize(0.0);
4962         }
4963 
4964         estimatedCovariance = fitter.getCovar();
4965         estimatedChiSq = fitter.getChisq();
4966         estimatedChiSqDegreesOfFreedom = fitter.getChisqDegreesOfFreedom();
4967         estimatedReducedChiSq = fitter.getReducedChisq();
4968         estimatedMse = fitter.getMse();
4969         try {
4970             estimatedP = fitter.getP();
4971             estimatedQ = fitter.getQ();
4972         } catch (final MaxIterationsExceededException ignore) {
4973             // if numerical instabilities arise, we assume worst case (no fit at all)
4974             // probability of finding a smaller chi square value is 1.0
4975             // quality of fit is 0.0
4976             estimatedP = 1.0;
4977             estimatedQ = 0.0;
4978         }
4979     }
4980 
4981     /**
4982      * Computes gravity versor error at the end of a sequence using provided
4983      * parameters.
4984      * This method is internally executed during gradient estimation and
4985      * Levenberg-Marquardt fitting needed for calibration computation.
4986      *
4987      * @param i      row position.
4988      * @param params array containing parameters for the general purpose case
4989      *               when G-dependent cross biases are taken into account. Must
4990      *               have length 18.
4991      * @return error between estimated and measured gravity versor.
4992      * @throws EvaluationException if there are numerical instabilities.
4993      */
4994     private double evaluateGeneralWithGDependentCrossBiases(final int i, final double[] params)
4995             throws EvaluationException {
4996 
4997         final var bx = params[0];
4998         final var by = params[1];
4999         final var bz = params[2];
5000 
5001         final var m11 = params[3];
5002         final var m21 = params[4];
5003         final var m31 = params[5];
5004 
5005         final var m12 = params[6];
5006         final var m22 = params[7];
5007         final var m32 = params[8];
5008 
5009         final var m13 = params[9];
5010         final var m23 = params[10];
5011         final var m33 = params[11];
5012 
5013         final var g11 = params[12];
5014         final var g21 = params[13];
5015         final var g31 = params[14];
5016 
5017         final var g12 = params[15];
5018         final var g22 = params[16];
5019         final var g32 = params[17];
5020 
5021         final var g13 = params[18];
5022         final var g23 = params[19];
5023         final var g33 = params[20];
5024 
5025         return evaluate(i, bx, by, bz, m11, m21, m31, m12, m22, m32, m13, m23, m33,
5026                 g11, g21, g31, g12, g22, g32, g13, g23, g33);
5027     }
5028 
5029     /**
5030      * Computes gravity versor error at the end of a sequence using provided
5031      * parameters.
5032      * This method is internally executed during gradient estimation and
5033      * Levenberg-Marquardt fitting needed for calibration computation.
5034      *
5035      * @param i      row position.
5036      * @param params array containing parameters for the general purpose case
5037      *               when G-dependent cross biases are taken into account. Must
5038      *               have length 15.
5039      * @return error between estimated and measured gravity versor.
5040      * @throws EvaluationException if there are numerical instabilities.
5041      */
5042     private double evaluateCommonAxisWithGDependentCrossBiases(final int i, final double[] params)
5043             throws EvaluationException {
5044 
5045         final var bx = params[0];
5046         final var by = params[1];
5047         final var bz = params[2];
5048 
5049         final var m11 = params[3];
5050 
5051         final var m12 = params[4];
5052         final var m22 = params[5];
5053 
5054         final var m13 = params[6];
5055         final var m23 = params[7];
5056         final var m33 = params[8];
5057 
5058         final var g11 = params[9];
5059         final var g21 = params[10];
5060         final var g31 = params[11];
5061 
5062         final var g12 = params[12];
5063         final var g22 = params[13];
5064         final var g32 = params[14];
5065 
5066         final var g13 = params[15];
5067         final var g23 = params[16];
5068         final var g33 = params[17];
5069 
5070         return evaluate(i, bx, by, bz, m11, 0.0, 0.0, m12, m22, 0.0, m13, m23, m33,
5071                 g11, g21, g31, g12, g22, g32, g13, g23, g33);
5072     }
5073 
5074     /**
5075      * Computes gravity versor error at the end of a sequence using provided
5076      * parameters.
5077      * This method is internally executed during gradient estimation and
5078      * Levenberg-Marquardt fitting needed for calibration computation.
5079      *
5080      * @param i      row position.
5081      * @param params array containing current parameters for the general purpose case
5082      *               when G-dependent cross biases are ignored. Must have length 9.
5083      * @return error between estimated and measured gravity versor.
5084      * @throws EvaluationException if there are numerical instabilities.
5085      */
5086     private double evaluateGeneral(final int i, final double[] params) throws EvaluationException {
5087 
5088         final var bx = params[0];
5089         final var by = params[1];
5090         final var bz = params[2];
5091 
5092         final var m11 = params[3];
5093         final var m21 = params[4];
5094         final var m31 = params[5];
5095 
5096         final var m12 = params[6];
5097         final var m22 = params[7];
5098         final var m32 = params[8];
5099 
5100         final var m13 = params[9];
5101         final var m23 = params[10];
5102         final var m33 = params[11];
5103 
5104         return evaluate(i, bx, by, bz, m11, m21, m31, m12, m22, m32, m13, m23, m33,
5105                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
5106     }
5107 
5108     /**
5109      * Computes gravity versor error at the end of a sequence using provided
5110      * parameters.
5111      * This method is internally executed during gradient estimation and
5112      * Levenberg-Marquardt fitting needed for calibration computation.
5113      *
5114      * @param i      row position.
5115      * @param params array containing current parameters for the common z-axis case
5116      *               when G-dependent cross biases are ignored. Must have length 9.
5117      * @return error between estimated and measured gravity versor.
5118      * @throws EvaluationException if there are numerical instabilities.
5119      */
5120     private double evaluateCommonAxis(final int i, final double[] params) throws EvaluationException {
5121         final var bx = params[0];
5122         final var by = params[1];
5123         final var bz = params[2];
5124 
5125         final var m11 = params[3];
5126 
5127         final var m12 = params[4];
5128         final var m22 = params[5];
5129 
5130         final var m13 = params[6];
5131         final var m23 = params[7];
5132         final var m33 = params[8];
5133 
5134         return evaluate(i, bx, by, bz, m11, 0.0, 0.0, m12, m22, 0.0, m13, m23, m33,
5135                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
5136     }
5137 
5138     /**
5139      * Computes gravity versor error at the end of a sequence using provided
5140      * parameters.
5141      * This method is internally executed during gradient estimation and
5142      * Levenberg-Marquardt fitting needed for calibration computation.
5143      *
5144      * @param i   row position.
5145      * @param bx  x coordinate of bias.
5146      * @param by  y coordinate of bias.
5147      * @param bz  z coordinate of bias.
5148      * @param m11 element 1,1 of cross-coupling error matrix.
5149      * @param m21 element 2,1 of cross-coupling error matrix.
5150      * @param m31 element 3,1 of cross-coupling error matrix.
5151      * @param m12 element 1,2 of cross-coupling error matrix.
5152      * @param m22 element 2,2 of cross-coupling error matrix.
5153      * @param m32 element 3,2 of cross-coupling error matrix.
5154      * @param m13 element 1,3 of cross-coupling error matrix.
5155      * @param m23 element 2,3 of cross-coupling error matrix.
5156      * @param m33 element 3,3 of cross-coupling error matrix.
5157      * @param g11 element 1,1 of g-dependent cross bias matrix.
5158      * @param g21 element 2,1 of g-dependent cross bias matrix.
5159      * @param g31 element 3,1 of g-dependent cross bias matrix.
5160      * @param g12 element 1,2 of g-dependent cross bias matrix.
5161      * @param g22 element 2,2 of g-dependent cross bias matrix.
5162      * @param g32 element 3,2 of g-dependent cross bias matrix.
5163      * @param g13 element 1,3 of g-dependent cross bias matrix.
5164      * @param g23 element 2,3 of g-dependent cross bias matrix.
5165      * @param g33 element 3,3 of g-dependent cross bias matrix.
5166      * @return error between estimated and measured gravity versor.
5167      * @throws EvaluationException if there are numerical instabilities.
5168      */
5169     private double evaluate(
5170             final int i, final double bx, final double by, final double bz,
5171             final double m11, final double m21, final double m31,
5172             final double m12, final double m22, final double m32,
5173             final double m13, final double m23, final double m33,
5174             final double g11, final double g21, final double g31,
5175             final double g12, final double g22, final double g32,
5176             final double g13, final double g23, final double g33) throws EvaluationException {
5177 
5178         try {
5179             final var measuredSequence = sequences.get(i);
5180             final var fixedSequence = fixedSequences.get(i);
5181 
5182             // generate new sequence using current parameters to fix angular rate measurements
5183             if (measuredSpecificForce == null) {
5184                 measuredSpecificForce = new Matrix(BodyKinematics.COMPONENTS, 1);
5185             }
5186             if (trueSpecificForce == null) {
5187                 trueSpecificForce = new Matrix(BodyKinematics.COMPONENTS, 1);
5188             }
5189 
5190             if (measuredAngularRate == null) {
5191                 measuredAngularRate = new Matrix(BodyKinematics.COMPONENTS, 1);
5192             }
5193             if (trueAngularRate == null) {
5194                 trueAngularRate = new Matrix(BodyKinematics.COMPONENTS, 1);
5195             }
5196 
5197             if (m == null) {
5198                 m = new Matrix(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
5199             }
5200             if (invM == null) {
5201                 invM = new Matrix(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
5202             }
5203             if (b == null) {
5204                 b = new Matrix(BodyKinematics.COMPONENTS, 1);
5205             }
5206             if (g == null) {
5207                 g = new Matrix(BodyKinematics.COMPONENTS, BodyKinematics.COMPONENTS);
5208             }
5209             if (tmp == null) {
5210                 tmp = new Matrix(BodyKinematics.COMPONENTS, 1);
5211             }
5212 
5213             m.setElementAt(0, 0, m11);
5214             m.setElementAt(1, 0, m21);
5215             m.setElementAt(2, 0, m31);
5216 
5217             m.setElementAt(0, 1, m12);
5218             m.setElementAt(1, 1, m22);
5219             m.setElementAt(2, 1, m32);
5220 
5221             m.setElementAt(0, 2, m13);
5222             m.setElementAt(1, 2, m23);
5223             m.setElementAt(2, 2, m33);
5224 
5225             Utils.inverse(m, invM);
5226 
5227             b.setElementAtIndex(0, bx);
5228             b.setElementAtIndex(1, by);
5229             b.setElementAtIndex(2, bz);
5230 
5231             g.setElementAt(0, 0, g11);
5232             g.setElementAt(1, 0, g21);
5233             g.setElementAt(2, 0, g31);
5234 
5235             g.setElementAt(0, 1, g12);
5236             g.setElementAt(1, 1, g22);
5237             g.setElementAt(2, 1, g32);
5238 
5239             g.setElementAt(0, 2, g13);
5240             g.setElementAt(1, 2, g23);
5241             g.setElementAt(2, 2, g33);
5242 
5243             // fix kinematics
5244             final var numItems = measuredSequence.getItemsCount();
5245             final var measuredItems = measuredSequence.getSortedItems();
5246             final var fixedItems = fixedSequence.getSortedItems();
5247             for (var j = 0; j < numItems; j++) {
5248                 final var measuredItem = measuredItems.get(j);
5249                 final var fixedItem = fixedItems.get(j);
5250 
5251                 fixKinematics(measuredItem.getKinematics(), fixedItem.getKinematics());
5252             }
5253 
5254 
5255             // integrate fixed sequence to obtain attitude change
5256             QuaternionIntegrator.integrateGyroSequence(fixedSequence, QuaternionStepIntegratorType.RUNGE_KUTTA, q);
5257 
5258             startPoint.setInhomogeneousCoordinates(point[0], point[1], point[2]);
5259             q.inverse();
5260             q.rotate(startPoint, endPoint);
5261 
5262             expectedEndPoint.setInhomogeneousCoordinates(point[3], point[4], point[5]);
5263 
5264             return expectedEndPoint.distanceTo(endPoint);
5265 
5266         } catch (final AlgebraException | RotationException e) {
5267             throw new EvaluationException(e);
5268         }
5269     }
5270 
5271     /**
5272      * Fixes provided kinematics with provided accelerometer parameters and
5273      * current gyroscope parameters.
5274      *
5275      * @param kinematics kinematics to be fixed with current values.
5276      * @param result     kinematics where result will be stored.
5277      * @throws AlgebraException if for some reason kinematics
5278      */
5279     private void fixKinematics(final BodyKinematics kinematics, final BodyKinematics result) throws AlgebraException {
5280 
5281         // Ωmeas = bg + (I + Mg) * Ωtrue + Gg * ftrue
5282         // Ωmeas = M*(Ωtrue + b + G * ftrue)
5283 
5284         // M = I + Mg
5285         // bg = M*b --> b = M^-1*bg
5286         // Gg = M*G --> G = M^-1*Gg
5287 
5288         // Ωtrue = M^-1 * Ωmeas - b - G*ftrue
5289 
5290         // fix specific force
5291         measuredSpecificForce.setElementAtIndex(0, kinematics.getFx());
5292         measuredSpecificForce.setElementAtIndex(1, kinematics.getFy());
5293         measuredSpecificForce.setElementAtIndex(2, kinematics.getFz());
5294 
5295         accelerationFixer.fix(measuredSpecificForce, trueSpecificForce);
5296 
5297         // fix angular rate
5298         measuredAngularRate.setElementAtIndex(0, kinematics.getAngularRateX());
5299         measuredAngularRate.setElementAtIndex(1, kinematics.getAngularRateY());
5300         measuredAngularRate.setElementAtIndex(2, kinematics.getAngularRateZ());
5301 
5302         g.multiply(trueSpecificForce, tmp);
5303         invM.multiply(measuredAngularRate, trueAngularRate);
5304         trueAngularRate.subtract(b);
5305         trueAngularRate.subtract(tmp);
5306 
5307         result.setSpecificForceCoordinates(
5308                 trueSpecificForce.getElementAtIndex(0),
5309                 trueSpecificForce.getElementAtIndex(1),
5310                 trueSpecificForce.getElementAtIndex(2));
5311 
5312         result.setAngularRateCoordinates(
5313                 trueAngularRate.getElementAtIndex(0),
5314                 trueAngularRate.getElementAtIndex(1),
5315                 trueAngularRate.getElementAtIndex(2));
5316     }
5317 
5318     /**
5319      * Resets any previous estimations.
5320      */
5321     private void reset() {
5322         estimatedBiases = null;
5323         estimatedMg = null;
5324         estimatedGg = null;
5325         estimatedCovariance = null;
5326         estimatedChiSq = 0.0;
5327         estimatedChiSqDegreesOfFreedom = 0;
5328         estimatedReducedChiSq = 0.0;
5329         estimatedMse = 0.0;
5330         estimatedP = 1.0;
5331         estimatedQ = 0.0;
5332         i = 0;
5333         fixedSequences = null;
5334         measuredSpecificForce = null;
5335         trueSpecificForce = null;
5336         measuredAngularRate = null;
5337         trueAngularRate = null;
5338         m = null;
5339         invM = null;
5340         b = null;
5341         g = null;
5342         tmp = null;
5343         point = null;
5344     }
5345 }