Requirements remove numpy version

https://github.com/nliaudat/magnetometer_calibration

.idea/Python-PS2000B.iml

@@ -1,7 +1,9 @@
 <?xml version="1.0" encoding="UTF-8"?>
 <module type="PYTHON_MODULE" version="4">
   <component name="NewModuleRootManager">
-    <content url="file://$MODULE_DIR$" />
+    <content url="file://$MODULE_DIR$">
+      <excludeFolder url="file://$MODULE_DIR$/venv" />
+    </content>
     <orderEntry type="jdk" jdkName="Python 3.9 (Helmholtz_Test_Bench)" jdkType="Python SDK" />
     <orderEntry type="sourceFolder" forTests="false" />
   </component>

Documentation/Abbreviation.tex

@@ -173,4 +173,5 @@
 \newacronym{pla}{PLA}{Polylactic Acid, a thermoplastic polyester}
 \newacronym{pc}{PC}{Personal Computer}
 \newacronym{tbd}{TBD}{To Be Determined}
-\newacronym{mpi}{MPI}{Message Passing Interface}
+\newacronym{mpi}{MPI}{Message Passing Interface}
+\newacronym{dut}{DUT}{Device Under Test}

Documentation/Chapters/03_Software.tex

@@ -165,7 +165,7 @@ More details on the execution thread is provided in Section \ref{sec:csv_exec}.
 
 \myparagraph{Magnetometer Calibration Mode Class \code{CalibrateMagnetometer}}
 This class constructs magnetometer calibration interface. 
-It is placed in the application's main area, its layout is shown in Figure TODO.
+It is placed in the application's main area, its layout is shown in Figure \ref{fig:magcalibrationpure}.
 This class creates and manages a separate execution thread defined in \code{calibration.py} to preform the calibration without blocking the \gls{ui}.
 Hardware control is acquired upon creation of the thread, or an exception is returned otherwise.
 More details on the calibration thread is provided in Section \ref{sec:calibration_processes}.
@@ -236,7 +236,7 @@ The line plot is modified to visually reflect the discrete and nearly instantane
 
 \subsection{Calibration Procedures \code{calibration.py}}
 \label{sec:calibration_processes}
-This file contains the worker thread objects \code{AmbientFieldCalibration}, \code{CoilConstantCalibration}, and \code{MagnetometerCalibration}.
+This file contains the worker thread objects \code{AmbientFieldCalibration}, \code{CoilConstantCalibration}, \code{MagnetometerCalibrationSimple}, and \code{MagnetometerCalibrationComplete}.
 All of these threads start by checking for and acquiring hardware during instantiation, while still in the main thread.
 Then upon running, the main code in \code{calibration\_procedure} is executed; This function may freely block and make use of sleep commands.
 Typically, calibration procedures should save two data detail levels: Processed data that is displayed as a final result, and detailed, raw data points that can be exported by the user to apply custom algorithms or verify the applications function.
@@ -337,7 +337,8 @@ The magnetometer object mirrors the state of the TCP connection, which is presen
 \section{Conversion to Executable}
 The program is compiled to a .exe executable file to enable simple use without a Python installation. For this, the "Auto Py to Exe" tool developed by B. Vollebregt \cite{auto_py_to_exe} is used. The resulting files are distributed in a separate repository on the \gls{irs} git server.\footnote{\url{https://egit.irs.uni-stuttgart.de/zietzm/Helmholtz_Test_Bench_Releases}} When a new software version is ready for publication, the following procedure should be used to release it:
 \begin{enumerate}
-	\item Run "auto-py-to-exe.exe" (provided in main development repository)
+	\item \code{pip install auto-py-to-exe} in the project's virtual environment.
+	\item Execute "auto-py-to-exe" in a console. To make sure the right version is executed, don't install auto-py-to-exe in your system environment.
 	\item In the Auto Py to Exe \gls{ui}, select \code{main.py} file as "Script Location"
 	\item Select options "One File" and "Window Based"
 	\item Select file "Helmholtz.ico" as the "Icon"

Documentation/Chapters/04_Users Guide.tex

@@ -118,8 +118,8 @@ These instructions assume a disassembly according to Section \ref{sec:disassembl
 \section{Software Users Guide}\label{sec:software_guide}
 \subsection{Installation}
 \begin{enumerate}
-	\item Download latest release: \url{https://egit.irs.uni-stuttgart.de/zietzm/Helmholtz_Test_Bench_Releases/releases}
-	\item Unpack ZIP-folder and run "Helmholtz Cage Control.exe"
+	\item Download latest release: \url{https://egit.irs.uni-stuttgart.de/eive/Helmholtz_Test_Bench/releases}
+	\item Unpack ZIP-folder and run "Release{\textbackslash}Helmholtz Control.exe"
 	\item Setup hardware and program according to Section \ref{sec:software_init}
 \end{enumerate}
 
@@ -232,7 +232,7 @@ The manual input mode is used to set static currents or magnetic fields on the t
 	\item For magnetic fields, choose whether ambient field should be compensated by (un)ticking the checkbox
 	\item Press "Execute" button, devices will now implement set values
 	\item Check console output to see if any errors occurred
-	\item Monitor behaviour in status display and on devices
+	\item Monitor behavior in status display and on devices
 	\item When finished, press "Power Down All" button to remove currents from the test bench
 \end{enumerate}
 
@@ -240,9 +240,9 @@ The manual input mode is used to set static currents or magnetic fields on the t
 This mode is used to run timed sequences of magnetic fields. These have to be defined in a \gls{csv} file of the following form:
 \begin{itemize}
 	\item \textit{Column separator:} Semicolon (;)
-	\item \textit{Decimal:} Comma (,)
+	\item \textit{Decimal:} Comma (,) or Period (.)
 	\item \textit{Line terminator:} Tested with Windows standard (\code{\textbackslash r\textbackslash n}), other options may work as well
-	\item \textit{Columns:} Time in seconds; X-axis, Y-axis and Z-axis  flux density in Tesla
+	\item \textit{Columns:} Time in seconds; X-axis, Y-axis and Z-axis flux density in Tesla
 \end{itemize}
 An example for the \gls{csv} file structure is given below:
 %[caption=Example \gls{csv} file]
@@ -284,11 +284,124 @@ The \gls{ui} layout is shown in Figure \ref{fig:csvmodepure}, its main elements
 	\item Monitor execution in console, status display and on devices
 \end{enumerate}
 
+\subsubsection*{Ambient Field and Coil Constant Calibration}
+\label{sec:ambient_field_calibration}
+The application offers calibration tools to determine the exact ambient field with the intention of cancelling it, as well for measuring the test bench's coil constants.
+These are both integrated into one application view, depicted in Figure \ref{fig:ambientcalibrationpure}.
+All the calibration tools require access to the complete Helmholtz cage hardware, as well as a magnetometer.
+These must be connected before starting the test.
+It is important that the magnetometer is centered as well as possible to achieve accurate results.
+For these tests, where the magnetometer itself is not the \gls{dut}, it is recommended to use the IRS's FGM3D reference magnetometer.
+Further, an adapter script (\code{fgm3d\_adapter.py}) already exists for this sensor.
+For more details on writing adapter scripts and connecting magnetometers, please refer to Section \ref{sec:tcp_api}.
+
+\begin{figure}[h]
+	\centering
+	\includegraphics[width=\linewidth]{media/ambient_calibration_pure}
+	\caption{Ambient field and coil constant calibration view.}
+	\label{fig:ambientcalibrationpure}
+\end{figure}
+
+After setting up the magnetometer, operation is simple: Click on the buttons to ``Calibrate Ambient Field'' or ``Calibrate Coil Constants''.
+This will cause the calibration procedure to start.
+The current status will be shown in the progress bar underneath.
+After completing, the results of the calibration procedure will be saved into the data fields on the right hand side of Figure \ref{fig:ambientcalibrationpure}.
+At this point, the buttons underneath the respective will become available.
+These two calibrations are primarily used to determine important application config parameters, namely the ambient field and coil constant as their names imply.
+To automatically apply the newly measured values, click ``Save and apply''.
+The ``Copy to clipboard'' will put the results table as shown in the \gls{ui} into the system clipboard, from which it can be pasted into software such as Microsoft Excel and LibreOffice Calc.
+If semi-raw experiment data is required to verify the program functioning or to apply custom algorithms, it can be exported with the ``Export raw CSV'' button.
+
+\vspace{5mm}
+\textbf{Ambient Field Calibration Method}
+This calibration uses a P-controller (implemented as PI-controller with $I=0$) to attempt to reach zero as a set point for the magnetometer.
+Updates are preformed approximately every half second by default and the calibration executes for 45 seconds.
+The procedure consistently achieves zero offsets below \SI{10}{\nano\tesla}.
+The ambient field exported by this calibration has been processed by multiplying the current that was required with the current coil constants, thus it is strongly recommended to first preform the coil constant calibration if both are to be executed.
+
+\vspace{5mm}
+\textbf{Coil Constant Calibration Method}
+The coil consant calibration uses a linear distribution of currents in each axis individually to collect coil consant data.
+Each setpoint is held for 3 seconds before being measured, and by default 8 points distributed across \SI{-3}{\ampere} to \SI{3}{\ampere}.
+To calculate the individual setpoints corresponding to each setpoint, first the field magnitude compared to the inital conditions is calculated, second the sign is estimated and reapplied, third the field is divided by the applied current.
+The final result is the average of all constants.
+In addition to the primary coil constant measurement, the maximum setpoints are also sampled once for each axis to calculate the angles between the coils.
+This is simply done by calculating the angle between the measured vectors.
+
+To understand the calibration methods in all detail, it is recommended to look at the \code{calibration.py} source file.
+
+\subsubsection*{Magnetometer Calibration}
+\label{sec:magnetometer_calibration}
+
+The helmholtz control software supports calibrating magnetometers.
+In preparation, a calibration of the ambient field using the reference magnetometer (FGM3D) should be conducted beforehand to achieve accurate results.
+Afterwards, the magnetometer must be replaced with the \gls{dut}.
+Since the software only supports one magnetometer at once, the reference magnetometer (meaning: its adapter script) should be disabled beforehand.
+Tip: To mount CubeSat magnetometers a purpose built PC104 holder can be found in accompanyment of the Helmholtz cage.
+
+As the coordinate system of the Helmholtz Cage and the magnetometer may not match, there are two options that available.
+First, the calibration can be conducted as-is, and the calibration results will reflect the unexpected coordinate system by often yielding negative sensitivities and untypical sensor axis angles.
+These values are not invalid, but describe a sensor correction that would transform it into the Helmholtz coordinate system.
+With some work, the desired calibration parameters could be extracted mathematically by changing signs and adding or subtracting \SI{90}{\degree} increments.
+The second option would be to change the sensors orientation in software in the adapter script that is used.
+This yields a more expected set of calibration parameters, but these now do not describe the initial axes anymore.
+In this case, the parameters must be correlated back to their initial axes and angles must be negated if the axis was flipped.
+This common issue should be solved robustly in the control/calibration software in the future.
+
+\begin{figure}[h]
+	\centering
+	\includegraphics[width=\linewidth]{media/magnetometer_calibration_pure}
+	\caption{Magnetometer calibration view.}
+	\label{fig:magcalibrationpure}
+\end{figure}
+The procedure to start the magnetometer can be done as follows:
+\begin{itemize}
+\item Activate the FGMM3D box, start the \textit{FGM3D TD Application} software and check if the device was recognissed automatically (Port: FGM3D TD (000125).
+ If the device was not recognized, click on \textit{File -> Rescan} devices, else restart the computer.
+\item Click on \textit{connct} and \textit{start}. Telemetry should now be displayed in the application. Use the \textit{Live Streaming} console and stream to \textbf{COM10}
+\item Start the \textit{Helmholtz Control.exe} software and navigate to e.g. \textit{Menu -> Ambient Filed Calibration}
+Start the \textit{fgm3dcadapter.py} script that forwards the data live stream from \textit{COM10} on the Sensys app to \textbf{COM11} of the  Helmholtz application.
+%\item Start the data live stream in the FGM3D software by first starting the data acquisition and then live streaming function. Check if the console of the applocation says "Magnetometer State: active"
+\item The test bench should be ready to run e.g. the Magnetometer Calibrtation script.
+\end{itemize}
+
+
+The calibration method used (Zikmund \cite{ref:calibration_procedure_magnetometer_helmholtz_cage}) relies on a simplified error model and the Helmholtz test bench.
+After either measuring the local geomagnetic field or cancelling it, the magnetometer-under-test runs through a sequence of magnetic fields generated by the Helmholtz coils, which supplies sufficient data to solve a system of equations containing the coefficients of interest.
+The non-linear system is constructed with the equation below on a per-axis basis, with one row for every sample.
+
+\begin{equation*}
+    B_{meas} = S \left( B_E \sin{\alpha_E} + B_x \cos{\alpha} \cos{\beta} + B_y \cos{\alpha} \sin{\beta} + B_z \sin{\alpha} \right)
+\end{equation*}
+
+The calibration procedure makes use of nearly equidistantly distributed vectors in all directions as test points.
+The number of these points, as well as the settle time before taking a measurement, can be set in the \gls{ui}.
+Generally, a high number of points, such as larger than 8, is recommended to achieve both an accurate result and also to get useful residual data.
+The meaning of the individual calibration coefficients can be derived from Zikmund \cite{ref:calibration_procedure_magnetometer_helmholtz_cage} or also the thesis accompanying the Helmholtz test bench \cite{ref:leons_test_bench}.
+
+After completing, the results of the calibration procedure will be saved into the data fields on the right hand side of Figure \ref{fig:magcalibrationpure}.
+This also enables the save buttons underneath.
+The ``Copy to clipboard'' will put the results table as shown in the \gls{ui} into the system clipboard, from which it can be pasted into software such as Microsoft Excel and LibreOffice Calc.
+If semi-raw experiment data is required to verify the program functioning or to apply custom algorithms, it can be exported with the ``Export raw CSV'' button.
+
 \subsubsection*{Data Logging Configuration Page}\label{sec:logging_guide}
-The application has the ability to log test bench data to a \gls{csv} file. The data is temporarily stored internally and must be saved to an external file by user request.\\
-An example of a log file is given in Appendix \ref{app:example_files}. The first three columns are time stamps: date, system time and time since the start of logging in seconds. The other columns contain the data, as selected by the user. All dynamic values from the status display can be logged, see Table \ref{tab:status_contents} for explanations. Each type of data is logged for all three axes. The units for numerical values are Volt, Ampere and Tesla.\\
-There are two options as for when and how often data is logged. The first option logs a row of data in a regular time interval specified by the user. The second option logs, whenever a significant command is sent to the test bench, for example when a new field vector is commanded. Both options can be used simultaneously.\\
-The logging configuration \gls{ui} is shown in Figure \ref{fig:loggingpure}. Its elements are listed below.
+The application has the ability to log test bench data to a \gls{csv} file.
+The data is temporarily stored internally and must be saved to an external file by user request.
+The logging output is highly customizable using the options shown in Figure \ref{fig:loggingpure}.
+
+An example of a log file is given in Appendix \ref{app:example_files}.
+Unless explicitly disabled, the first three columns are time stamps: date, system time and time since the start of logging in seconds.
+All dynamic values from the status display as well as the magnetometer status can be logged, see Table \ref{tab:status_contents} for explanations.
+Each type of data is logged for all three axes by default, but this can also be specified.
+The units for numerical values are Volt, Ampere and Tesla.
+
+There are two options as for when and how often data is logged.
+The first option logs a row of data in a regular time interval specified by the user.
+The second option logs, whenever a significant command is sent to the test bench, for example when a new field vector is commanded.
+Both options can be used simultaneously.
+
+The logging configuration \gls{ui} is shown in Figure \ref{fig:loggingpure}.
+Its elements are listed below.
 
 \begin{figure}[h]
 	\centering
@@ -309,7 +422,8 @@ The logging configuration \gls{ui} is shown in Figure \ref{fig:loggingpure}. Its
 	\item \textbf{"Datapoints logged" counter:} Displays the current number of logged data rows
 	\item \textbf{"Log in regular intervals" controls:} Enable checkbox to periodically log data, set interval (in seconds) in the entry field to the right
 	\item \textbf{"Log whenever test bench is commanded" checkbox:} Enable, to log data on significant changes to the test bench (e.g. a new field vector is commanded)
-	\item \textbf{Data selection checkboxes:} Select what data to log, explanations are given in Table \ref{tab:status_contents}
+	\item \textbf{Data selection checkboxes:} Select what data to log, explanations for most options are given in Table \ref{tab:status_contents}.
+	In addition, the Log X/Y/Z-Axis Data checkboxes specify whether the selected logging variables will be included for the respective axis.
 \end{itemize}
 \textbf{To collect and save log data:}
 \begin{enumerate}
@@ -351,11 +465,10 @@ Figure \ref{fig:settingspure} shows a screenshot of the \gls{ui} layout with the
 		\item Opens dialogue to let user choose new file path and name (must be *.ini)
 		\item Reinitializes test bench devices with current settings
 	\end{itemize}
-	\item \textbf{"\gls{psu} Serial Port" entries:} Input \gls{com} ports for both \gls{psu}s here
+	\item \textbf{"Serial Port" entries:} Input \gls{com} ports for both \gls{psu}s and the switch box here
 	\begin{itemize}
-		\item Use Windows device manager to find correct port names (connect \gls{psu}s separately to differentiate between devices)
+		\item Use Windows device manager to find correct port names (connect devices separately to differentiate between devices)
 		\item Test bench X- and Y-axes need to be connected to channel 1 and 2 of one \gls{psu}, Z-axis to channel 1 of the other
-		\item \textit{Note: Switch box Arduino should be found automatically}
 	\end{itemize}
 	\item \textbf{Program constant entry fields:} Set constants here, details are listed in Table \ref{tab:settings_entries}
 	\item \textbf{"Update and Reinitialize" button:} Implements any changed settings in the program and on test bench devices, needs to be pressed for changes to take effect
@@ -419,7 +532,7 @@ Figure \ref{fig:settingspure} shows a screenshot of the \gls{ui} layout with the
 		\item If program has not been configured before:
 		\begin{enumerate}
 			\setcounter{enumii}{2}
-			\item Use Windows device manager to find correct serial \gls{com} ports for \gls{psu}s (connect /disconnect in turn to differentiate between devices)
+			\item Use Windows device manager to find correct serial \gls{com} ports for \gls{psu}s and Arduino (connect /disconnect in turn to differentiate between devices)
 			\item Enter \gls{com} port names in application (switch box should be found automatically)
 			\item Press "Update and Reinitialize" button
 		\end{enumerate}
@@ -435,6 +548,7 @@ Figure \ref{fig:settingspure} shows a screenshot of the \gls{ui} layout with the
 		\setcounter{enumii}{5}
 		\item Check console print to see if all devices were found, otherwise check physical connections and \gls{com} port settings
 	\end{enumerate}
+	\item If using a magnetometer, the listening port is open and the adapter script can now be started
 	\item Test configuration
 	\begin{enumerate}
 		\item Go to manual mode (Menu$\rightarrow$Static Manual Input)
@@ -444,6 +558,7 @@ Figure \ref{fig:settingspure} shows a screenshot of the \gls{ui} layout with the
 			\item Current should be activated on correct \gls{psu} channel
 			\item For negative currents, corresponding status \gls{led} on switch box should light up and relay actuation be audible as clicking sound
 		\end{itemize}
+		\item If using a magnetometer: Check the device output in either of the calibration-views
 	\end{enumerate}
 	\item Go back to the settings page (Menu$\rightarrow$Settings...)
 	\item Change program constants as needed (e.g. enter measured ambient field), see Section \ref{sec:settings_guide}
@@ -457,4 +572,50 @@ Figure \ref{fig:settingspure} shows a screenshot of the \gls{ui} layout with the
 \subsection{TCP Remote Control}
 \label{sec:tcp_api}
 
-TODO.
+The Helmholtz Control Software offers a TCP automation interface that is opened upon application start-up.
+To find out which port is being listened on, look for the corresponding information in the application console, but typically, it will be port 6677.
+The commands that are accepted by the TCP interface are documented in Table \ref{tab:tcp_comamnds} and also in the \code{socket\_control.py} source file.
+All field commands will implicitly preform the usual safety checks to ensure safe operation.
+The commands that are shown must all be terminated with a single {\textbackslash}n (newline) char.
+Commands may be split across multiple packets if desired.
+Important Note: Before useful commands can be sent, \code{declare\_api\_version} must be called.
+
+The \code{tools} folder in the application git repository contains an example implementation (\code{fgm3d\_adapter.py}) of a magnetometer interface using the TCP socket.
+
+%\small
+\begin{longtable}{lp{8.5cm}}
+	\caption{TCP Remote Control Commands}\\
+	\hline
+	\textbf{Command} & \textbf{Description}\\ \hline
+	\code{set\_raw\_field} [X] [Y] [Z] & Returns: 0 or 1 for success.
+	Accepts decimal point formatted floats, with or without scientific notation. The float() cast must understand it.
+	The field units are Tesla.
+	This causes an additional field of the given strength to be generated, without regard for the pre-existing geomagnetic/external fields. \\ \hline
+
+	\code{set\_compensated\_field} [X] [Y] [Z] & Returns: 0 or 1 for success.
+	Accepts decimal point formatted floats, with or without scientific notation. The float() cast must understand it.
+	The field units are Tesla.
+	This causes a field of exactly the given magnitude to be generated by compensating external factors such as the
+	geomagnetic field. \\ \hline
+
+	\code{set\_coil\_currents} [X] [Y] [Z] & Returns: 0 or 1 for success.
+	Accepts decimal point formatted floats, with or without scientific notation. The float() cast must understand it.
+	The field units are Ampere.
+    This establishes the requested current in the individual coils. \\ \hline
+
+	\code{magnetometer\_field} [X] [Y] [Z] & Returns: 1.
+	Accepts decimal point formatted floats, with or without scientific notation. The float() cast must understand it.
+	The field units are Tesla.
+	Sets the state of a virtual magnetometer object which mirrors a physical sensor by means of
+	this command. \\ \hline
+
+	\code{get\_api\_version} & Returns: a string uniquely identifying each API version.
+    This function can be called before \code{declare\_api\_version}. \\ \hline
+
+	\code{declare\_api\_version} [version] & Returns: 0 or 1.
+	Declare the API version the client application was programmed for.
+	It must be compatible with the current API version. This prevents unexpected behavior by forcing programmers to specify which API they are expecting.
+	This function must be called before sending HW commands. \\ \hline
+	\label{tab:tcp_comamnds}
+\end{longtable}
+%\normalsize

Documentation/Preamble/TitlePage.tex

@@ -14,6 +14,7 @@
 \begin{center}
 Author:\\
 Martin Zietz\\
+Leon Teichröb\\
 
 Supervisors:\\
 M.Sc. Markus T. Koller\\

Documentation/References.bib

@@ -172,4 +172,23 @@
 	author={Stier, Annika and Schweigert, Robin and Galla, Daniel and Lengowski, Michael and Klinkner, Sabine},
 	year={2020},
 	publisher={University of Leicester}
+}
+
+@inproceedings{ref:calibration_procedure_magnetometer_helmholtz_cage,
+    author = {Zikmund, A. and Janosek, Michal},
+    year = {2014},
+    month = {05},
+    pages = {473-476},
+    title = {Calibration procedure for triaxial magnetometers without a compensating system or moving parts},
+    isbn = {9781467363860},
+    journal = {Conference Record - IEEE Instrumentation and Measurement Technology Conference},
+    doi = {10.1109/I2MTC.2014.6860790}
+}
+
+@thesis{ref:leons_test_bench,
+	author = {Leon Teichröb},
+	title = {Mapping and Calibration of a Helmholtz Magnetic Field Cage and Test of the EIVE Attitude Control System},
+	date = {2021},
+	institution = {Institute for Space Systems (IRS), University of Stuttgart},
+	location = {Stuttgart},
 }

requirements.txt

@@ -3,7 +3,7 @@ cycler==0.10.0
 future==0.18.2
 kiwisolver==1.3.2
 matplotlib==3.3.4
-numpy==1.19.3
+# numpy==1.19.3  ## do not include versioning to avoid versioning conflict
 pandas==1.1.5
 Pillow==8.4.0
 pyparsing==2.4.7
@@ -12,3 +12,4 @@ python-dateutil==2.8.2
 pytz==2021.3
 scipy==1.7.1
 six==1.16.0
+screeninfo~=0.8.1

src/calibration.py

@@ -3,6 +3,9 @@ import time
 from datetime import datetime
 from threading import Thread
 import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.backends.backend_tkagg import FigureCanvasTkAgg
+from tkinter import LabelFrame
 import scipy.optimize
 
 from src.utility import ui_print
@@ -21,6 +24,7 @@ class AmbientFieldCalibration(Thread):
     P_CONTROL = -7e3  # 0.2 A/s slew-rate at 40uT
     I_CONTROL = 0  # -1e4  # 0.01A/s slew-rate for 1uTs
     I_LIMIT = 1e-7  # uTs, Limit I to 0.025 A/s slew-rate to prevent wind-up
+
     # D_CONTROL = Not implemented for now
 
     def __init__(self, view_queue):
@@ -63,7 +67,7 @@ class AmbientFieldCalibration(Thread):
         target_time = 0
         current_time = datetime.now()
         while (current_time - start_time).seconds < self.SETTLE_TIME:
-            # Each axis runs its own PID controller. They are slightly coupled by unorthogonality, which should
+            # Each axis runs its own PID controller. They are slightly coupled by non-orthogonality, which should
             # hopefully not destabilize the feedback loop
             for i in range(3):
                 # Error in tesla
@@ -71,11 +75,11 @@ class AmbientFieldCalibration(Thread):
                 e = g.MAGNETOMETER.field[i]
                 # Change in control current
                 du = e * self.P_CONTROL + self.error_integral[i] * self.I_CONTROL
-                self.axis_currents[i] += du*dt
+                self.axis_currents[i] += du * dt
 
                 # Update integral
                 # Add increment
-                self.error_integral[i] += e*dt
+                self.error_integral[i] += e * dt
                 # Clamp range
                 self.error_integral = np.clip(self.error_integral, -self.I_LIMIT, self.I_LIMIT)
 
@@ -149,7 +153,7 @@ class CoilConstantCalibration(Thread):
         # All generated fields will be compared to this using a simple difference method
         ambient_field = g.MAGNETOMETER.field
 
-        # This generates linearly spaced current setpoints and excludes zero
+        # This generates linearly spaced current set points and excludes zero
         currents = np.linspace(-self.MEASUREMENT_RANGE, self.MEASUREMENT_RANGE, self.MEASUREMENT_POINTS * 2 + 1)
         currents = np.delete(currents, self.MEASUREMENT_POINTS)
 
@@ -234,14 +238,16 @@ class CoilConstantCalibration(Thread):
         self.view_queue.put({'cmd': command, 'arg': arg})
 
 
-class MagnetometerCalibration(Thread):
-    TEST_VECTOR_MAGNITUDE = 100e-6  # In Tesla. Chosen so it can be achieved with a 3A PSU.
+class MagnetometerCalibrationSimple(Thread):
 
-    def __init__(self, view_queue, calibration_points, calibration_interval):
+    def __init__(self, view_queue, calibration_points, calibration_interval, calibration_mag_field,
+                 mgm_to_helmholtz_cos_trans):
         Thread.__init__(self)
         self.view_queue = view_queue
         self.calibration_points = calibration_points
+        self.calibration_mag_field = calibration_mag_field
         self.calibration_interval = calibration_interval
+        self.matrix_trans_mgm_to_hh = [[x.get() for x in row] for row in mgm_to_helmholtz_cos_trans]
 
         # Hardware checks are done in the init method to allow for exception handling in main thread
         # This means the run method should/must be called directly after Thread object creation.
@@ -281,8 +287,9 @@ class MagnetometerCalibration(Thread):
         self.cage_dev.set_field_compensated([0, 0, 0])
         # Sleep for a certain duration to allow psu to stabilize output and magnetometer to supply readings
         time.sleep(self.calibration_interval)
+        matrix_trans_mgm_to_hh_np = np.array(self.matrix_trans_mgm_to_hh)
         # The offsets can easily be read from the magnetometer
-        offsets = g.MAGNETOMETER.field
+        offsets = matrix_trans_mgm_to_hh_np.dot(g.MAGNETOMETER.field)
         # Save data point to raw_data list
         raw_data.append({'applied_x': 0, 'applied_y': 0, 'applied_z': 0,
                          'measured_x': offsets[0], 'measured_y': offsets[1], 'measured_z': offsets[2]})
@@ -299,14 +306,14 @@ class MagnetometerCalibration(Thread):
         # Collect sensor data for each test vector
         for vec_idx, test_vec in enumerate(test_vectors):
             # Command output
-            applied_vec = test_vec * self.TEST_VECTOR_MAGNITUDE
+            applied_vec = test_vec * self.calibration_mag_field
             self.cage_dev.set_field_raw(applied_vec)
 
             # Sleep for a certain duration to allow psu to stabilize output and magnetometer to supply readings
             time.sleep(self.calibration_interval)
 
             # Read output and save to array for solver later
-            raw_reading = g.MAGNETOMETER.field
+            raw_reading = matrix_trans_mgm_to_hh_np.dot(g.MAGNETOMETER.field)
             reading = raw_reading - offsets
             for i in range(3):
                 row = {'m': reading[i], 'b_x': applied_vec[0], 'b_y': applied_vec[1], 'b_z': applied_vec[2]}
@@ -338,12 +345,13 @@ class MagnetometerCalibration(Thread):
         b_e_x = g.CAGE_DEVICE.axes[0].ambient_field
         b_e_y = g.CAGE_DEVICE.axes[1].ambient_field
         b_e_z = g.CAGE_DEVICE.axes[2].ambient_field
-        b_e = sqrt(b_e_x**2 + b_e_y**2 + b_e_z**2)
+        b_e = sqrt(b_e_x ** 2 + b_e_y ** 2 + b_e_z ** 2)
 
         # Perform least squares optimization on all magnetometer axes
         sensor_parameters = []
         for axis, axis_samples in enumerate(samples):
-            result = scipy.optimize.least_squares(self.residual_function, (1.0, pi/4, pi/4, pi/4), args=(b_e, axis_samples), gtol=1e-13)
+            result = scipy.optimize.least_squares(self.residual_function, (1.0, pi / 4, pi / 4, pi / 4),
+                                                  args=(b_e, axis_samples), gtol=1e-13)
             s, alpha_e, alpha, beta = result.x
             residual = np.max(np.abs(result.fun))
             sensor_parameters.append({'sensitivity': s,
@@ -368,7 +376,8 @@ class MagnetometerCalibration(Thread):
             b_x = sample['b_x']
             b_y = sample['b_y']
             b_z = sample['b_z']
-            res.append(m - s * (b_e*sin(alpha_e) + b_x*cos(alpha)*cos(beta) + b_y*cos(alpha)*sin(beta) + b_z*sin(alpha)))
+            res.append(m - s * (b_e * sin(alpha_e) + b_x * cos(alpha) * cos(beta)
+                                + b_y * cos(alpha) * sin(beta) + b_z * sin(alpha)))
         return res
 
     @staticmethod
@@ -393,4 +402,686 @@ class MagnetometerCalibration(Thread):
         return points
 
     def put_message(self, command, arg):
-        self.view_queue.put({'cmd': command, 'arg': arg})
+        self.view_queue.put({'cmd': command, 'arg': arg})
+
+
+class MagnetometerCalibrationComplete(Thread):
+
+    def __init__(self, view_queue, calibration_points, calibration_interval, calibration_mag_field,
+                 mgm_to_helmholtz_cos_trans, right_column):
+        Thread.__init__(self)
+        self.view_queue = view_queue
+        self.calibration_points = calibration_points
+        self.calibration_interval = calibration_interval
+        self.calibration_mag_field = calibration_mag_field
+        self.matrix_trans_mgm_to_hh = [[x.get() for x in row] for row in mgm_to_helmholtz_cos_trans]
+        self.right_column = right_column
+
+        # Hardware checks are done in the init method to allow for exception handling in main thread
+        # This means the run method should/must be called directly after Thread object creation.
+
+        # Make sure we really have magnetometer data
+        if not g.MAGNETOMETER.connected:
+            ui_print("\nError: The magnetometer is not connected. Required for ambient field calibration.")
+            raise DeviceAccessError("The magnetometer is not connected. Required for ambient field calibration.")
+
+        # Acquire cage device. This resource will only be released after the thread is ended.
+        try:
+            self.cage_dev = g.CAGE_DEVICE.request_proxy()
+        except DeviceBusy:
+            ui_print("\nError: Failed to acquire coil control. Required for ambient field calibration.")
+            raise DeviceAccessError("Failed to acquire coil control. Required for ambient field calibration.")
+
+    def run(self):
+        try:
+            raw_data = self.calibration_procedure()
+            self.put_message('finished', None)
+            return raw_data
+        except Exception as e:
+            self.put_message('failed', e)
+        finally:
+            self.cage_dev.close()
+
+    def calibration_procedure(self):
+        # According to method outlined in:
+        # Chi, C. & Janosek, Lv, J-W., Wang, D (2018). Calibration of triaxial magnetometer with ellipsoid
+        # fitting method, DOI:10.1088/1755-1315/237/3/032015 and https://github.com/nliaudat/magnetometer_calibration
+
+        # This contains the raw experiment data for exporting
+        # Each row is a dict containing the applied vector and measured vector
+        raw_data = []
+
+        # Find sensor offsets. They must be found prior to applying the chosen calibration algorithm
+        # This will be accurate if the cage was recently calibrated
+        self.cage_dev.set_field_compensated([0, 0, 0])
+        # Sleep for a certain duration to allow psu to stabilize output and magnetometer to supply readings
+        time.sleep(self.calibration_interval)
+        # The offsets can easily be read from the magnetometer
+        offsets = g.MAGNETOMETER.field
+        # Save data point to raw_data list
+        raw_data.append({'applied_x': 0, 'applied_y': 0, 'applied_z': 0,
+                         'measured_x': offsets[0], 'measured_y': offsets[1], 'measured_z': offsets[2]})
+        # Set new progress indicator for UI
+        self.set_progress(True, 0)
+
+        # Generate our set of test vectors
+        test_vectors = self.fibonacci_sphere(self.calibration_points)
+
+        # Holds the known variables for each row of our system of equations. These are M, B_x, B_y, B_z
+        # (B_E is constant for the test and not stored in the array)
+        # Each sensor axis has its own independent system of equations
+        samples = [[], [], []]
+        # Collect sensor data for each test vector
+        for vec_idx, test_vec in enumerate(test_vectors):
+            # Command output
+            applied_vec = test_vec * self.calibration_mag_field
+            self.cage_dev.set_field_raw(applied_vec)
+
+            # Sleep for a certain duration to allow psu to stabilize output and magnetometer to supply readings
+            time.sleep(self.calibration_interval)
+
+            # Read output and save to array for solver later
+            raw_reading = g.MAGNETOMETER.field
+            reading = raw_reading - offsets
+            for i in range(3):
+                row = {'m': reading[i], 'b_x': applied_vec[0], 'b_y': applied_vec[1], 'b_z': applied_vec[2]}
+                # self.put_message("[Axis {}] {}".format(i, row))
+                samples[i].append(row)
+
+            # Save data point to raw_data list
+            raw_data.append({'applied_x': applied_vec[0], 'applied_y': applied_vec[1], 'applied_z': applied_vec[2],
+                             'measured_x': raw_reading[0], 'measured_y': raw_reading[1], 'measured_z': raw_reading[2]})
+
+            # Set new progress indicator for UI
+            self.set_progress(True, vec_idx + 1)
+
+        # Put device into an off and ready state
+        self.cage_dev.idle()
+        return raw_data
+
+    def set_progress(self, offset_complete, test_vec_index):
+        progress = int(offset_complete) * 0.2 + (test_vec_index / self.calibration_points) * 0.8
+        self.put_message('progress', progress)
+
+    @staticmethod
+    def solve_system(raw_data, matrix_trans_mgm_to_hh_tk):
+        u_tesla = 10 ** 6
+        # Unpack data:
+        mag_x_set = np.zeros(len(raw_data))
+        mag_y_set = np.zeros(len(raw_data))
+        mag_z_set = np.zeros(len(raw_data))
+        mag_x_m = np.zeros(len(raw_data))
+        mag_y_m = np.zeros(len(raw_data))
+        mag_z_m = np.zeros(len(raw_data))
+        for i in range(len(raw_data)):
+            mag_x_set[i] = raw_data[i]['applied_x']
+            mag_y_set[i] = raw_data[i]['applied_y']
+            mag_z_set[i] = raw_data[i]['applied_z']
+            mag_x_m[i] = raw_data[i]['measured_x']
+            mag_y_m[i] = raw_data[i]['measured_y']
+            mag_z_m[i] = raw_data[i]['measured_z']
+
+        # Apply coordinate transformation from Helmholtz system to MGM system
+        matrix_trans_mgm_to_hh_np = np.zeros((3, 3))
+        for row in range(3):
+            for col in range(3):
+                matrix_trans_mgm_to_hh_np[row][col] = matrix_trans_mgm_to_hh_tk[row][col]
+        # matrix_trans_mgm_to_hh_np = [[-1, 0, 0], [0, 1, 0], [0, 0, -1]]  # hardcoded for MGM 1 / 3
+        # matrix_trans_mgm_to_hh_np = [[0, 1, 0], [-1, 0, 0], [0, 0, 1]]  # hardcoded for MGM 0 / 2
+        ui_print("Applying transformation matrix to data. h_{set,mgm} = mgm_T_hh * h_{set,hh}")
+        ui_print(matrix_trans_mgm_to_hh_np)
+        # Transform hh set magnetic field cos to sensor cos
+        mag_xyz_set = np.c_[mag_x_set, mag_y_set, mag_z_set]
+        mag_xyz_set_hh = np.matmul(matrix_trans_mgm_to_hh_np, mag_xyz_set.T).T
+        mag_x_set = mag_xyz_set_hh[:, 0]
+        mag_y_set = mag_xyz_set_hh[:, 1]
+        mag_z_set = mag_xyz_set_hh[:, 2]
+        # Calculate total error
+        err_x = 0
+        err_y = 0
+        err_z = 0
+        for i in range(1, len(mag_x_m)):
+            err_x += (mag_x_m[i] - mag_x_set[i]) ** 2
+            err_y += (mag_y_m[i] - mag_y_set[i]) ** 2
+            err_z += (mag_z_m[i] - mag_z_set[i]) ** 2
+        err_t = np.sqrt(err_x ** 2 + err_y ** 2 + err_z ** 2)
+        ui_print('Error in X/Y/Z/total: {:.2e} / {:.2e} / {:.2e} / {:.2e} [uT]'.format(err_x * u_tesla, err_y * u_tesla,
+                                                                                       err_z * u_tesla,
+                                                                                       err_t * u_tesla))
+
+        # Filter raw data
+        try:
+            mag_x_set, mag_y_set, mag_z_set, mag_x_m, mag_y_m, mag_z_m = \
+                MagnetometerCalibrationComplete.filter_magnetometer_data(mag_x_set, mag_y_set, mag_z_set,
+                                                                         mag_x_m, mag_y_m, mag_z_m)
+        except Warning as waring_filter_mgm_data:
+            ui_print(f"{waring_filter_mgm_data}")
+        # Get general magnitude of the set magnetic filed
+        mag_amp_set = np.sqrt(
+            mag_x_set[1:len(mag_x_set)] ** 2 + mag_y_set[1:len(mag_x_set)] ** 2 + mag_z_set[1:len(mag_x_set)] ** 2)
+        mag_amp_avg_set = np.mean(mag_amp_set)
+        mag_amp_m = np.sqrt(
+            mag_x_m[1:len(mag_x_m)] ** 2 + mag_y_m[1:len(mag_x_m)] ** 2 + mag_z_m[1:len(mag_x_m)] ** 2)
+        mag_amp_avg_m = np.mean(mag_amp_m)
+
+        # Calculate ellipsoid fit
+        try:
+            q_mat, n, d = MagnetometerCalibrationComplete.fit_ellipsoid(mag_x_m, mag_y_m, mag_z_m)
+            ui_print('Q matrix =')
+            ui_print(q_mat)
+            ui_print('n =  [T]')
+            ui_print(n)
+            # Retrieve calibration parameters
+            q_mat_inv = np.linalg.inv(q_mat)
+            b = -np.dot(q_mat_inv, n)
+            a_mat_inv = np.real(mag_amp_avg_set / np.sqrt(np.dot(n.T, np.dot(q_mat_inv, n)) - d) * scipy.linalg.sqrtm(q_mat))
+            a_mat = np.linalg.inv(a_mat_inv)
+            # Calculate error
+            cal_x = np.zeros(mag_x_m.shape)
+            cal_y = np.zeros(mag_y_m.shape)
+            cal_z = np.zeros(mag_z_m.shape)
+            total_error = 0
+            for i in range(len(mag_x_m)):
+                h = np.array([[mag_x_m[i], mag_y_m[i], mag_z_m[i]]]).T
+                h_hat = np.matmul(a_mat_inv, h - b)  # Scaling issue
+                cal_x[i] = h_hat[0]
+                cal_y[i] = h_hat[1]
+                cal_z[i] = h_hat[2]
+                mag = np.dot(h_hat.T, h_hat)
+                err = (mag[0][0] - 1) ** 2
+                total_error += err
+            # Scale unity solution back to original scaling
+            a_mat_inv = a_mat_inv
+            # Console output
+            ui_print("Average magnitude: mag_amp_avg_set = {:.4f} uT, mag_amp_avg_m = {:.4f} uT".format(
+                mag_amp_avg_set * 10 ** 6, mag_amp_avg_m * 10 ** 6))
+            ui_print('Soft-iron correction matrix a_mat_inv = ')
+            ui_print(a_mat_inv)
+            ui_print('Hard-iron offset b = [T]')
+            ui_print(b)
+            ui_print("Normalized total error E = {:4e} [0..1]]".format(total_error))
+            sensor_parameters = [{'a_mat': a_mat.tolist(),
+                                  'a_mat_inv': a_mat_inv.tolist(),
+                                  'b': b.tolist(),
+                                  'q_mat': q_mat.tolist(),
+                                  'n': n.tolist(),
+                                  'mag_amp_avg_set': mag_amp_avg_set.tolist(),
+                                  'total_error': total_error}]
+            return sensor_parameters, mag_x_set, mag_y_set, mag_z_set, mag_x_m, mag_y_m, mag_z_m, cal_x, cal_y, cal_z, mag_amp_avg_set
+        except Warning as warning_message:
+            ui_print('A_inv could not be calculated! A warning occurred.')
+            ui_print('Please check if transformation matrix is input correctly.')
+            ui_print(f"{warning_message}")
+        except Exception as exception_message:
+            ui_print('A_inv could not be calculated! An unknown error occurred.')
+            ui_print('Please check if transformation matrix is input correctly.')
+            ui_print(f"{exception_message}")
+
+    @staticmethod
+    def filter_magnetometer_data(mag_x_set, mag_y_set, mag_z_set, mag_x_m, mag_y_m, mag_z_m):
+        flag_validity = np.ones(len(mag_x_set))
+        # Issue 1: Sometimes, the old and new measurement are very similar.
+        #          The set magnetic field has not changed fast enough.
+
+        for i in range(1, len(mag_x_set)):
+            # Filter entry if all measurement results are very close to one another
+            if mag_x_m[i - 1] != 0:
+                div_x_m = mag_x_m[i] / mag_x_m[i - 1]
+            else:
+                div_x_m = 1
+            if mag_y_m[i - 1] != 0:
+                div_y_m = mag_y_m[i] / mag_y_m[i - 1]
+            else:
+                div_y_m = 1
+            if mag_z_m[i - 1] != 0:
+                div_z_m = mag_z_m[i] / mag_z_m[i - 1]
+            else:
+                div_z_m = 1
+            if mag_x_set[i - 1] != 0:
+                div_x_set = mag_x_set[i] / mag_x_set[i - 1]
+            else:
+                div_x_set = 1
+            if mag_y_set[i - 1] != 0:
+                div_y_set = mag_y_set[i] / mag_y_set[i - 1]
+            else:
+                div_y_set = 1
+            if mag_z_set[i - 1] != 0:
+                div_z_set = mag_z_set[i] / mag_z_set[i - 1]
+            else:
+                div_z_set = 1
+            bound_low = 0.9
+            bound_up = 1.1
+            if bound_low < div_x_m < bound_up and not bound_low < div_x_set < bound_up:
+                flag_validity[i] = 0
+            if bound_low < div_y_m < bound_up and not bound_low < div_y_set < bound_up:
+                flag_validity[i] = 0
+            if bound_low < div_z_m < bound_up and not bound_low < div_z_set < bound_up:
+                flag_validity[i] = 0
+            # Filter entry if signs of measurement result and set result are different
+            thresh = 70.0e-06  # [T] threshold value, difference in magnetic field component too large to be reasonable
+            if np.abs(mag_x_set[i] - mag_x_m[i]) > thresh or np.abs(mag_y_set[i] - mag_y_m[i]) > thresh or np.abs(
+                    mag_z_set[i] - mag_z_m[i]) > thresh:
+                flag_validity[i] = 0
+            # ui_print("Ix {} X_set: {:.2e}, X_m: {:.2e}, Y_set: {:.2e}, Y_m: {:.2e}, Z_set: {:.2e}, Z_m: {:.2e}"
+            #       .format(i, mag_x_set[i], mag_x_m[i], mag_y_set[i],  mag_y_m[i],  mag_z_set[i], mag_z_m[i]))
+            # ui_print("Ix {} X_set-X_m: {:.2e}, Y_set-Y_m: {:.2e}, Z_set-Z_m: {:.2e}"
+            #       .format(i, np.abs(mag_x_set[i]-mag_x_m[i]), np.abs(mag_y_set[i]-mag_y_m[i]), np.abs(mag_z_set[i]-mag_z_m[i])))
+        # Issue: first element is zero-field measurement and always zero -> might lead to issues during calculation
+        flag_validity[0] = 0
+        # Delete faulty measurements from arrays
+        n_zeros = np.count_nonzero(flag_validity == 0)
+        ix_del = np.zeros(n_zeros)
+        j = 0
+        for i in range(0, len(mag_x_set)):
+            if flag_validity[i] == 0:
+                ix_del[j] = int(i - 0)
+                j = j + 1
+        ui_print('{} flagged measurements with indices {} removed.'.format(len(ix_del), ix_del))
+
+        mag_x_set = np.delete(mag_x_set, ix_del.astype(int), 0)
+        mag_y_set = np.delete(mag_y_set, ix_del.astype(int), 0)
+        mag_z_set = np.delete(mag_z_set, ix_del.astype(int), 0)
+        mag_x_m = np.delete(mag_x_m, ix_del.astype(int), 0)
+        mag_y_m = np.delete(mag_y_m, ix_del.astype(int), 0)
+        mag_z_m = np.delete(mag_z_m, ix_del.astype(int), 0)
+
+        n_filtered = len(mag_x_set)
+        if n_filtered < 5:
+            raise Warning("Filtered vectors only contains less than 5 Elements, filter not applied!")
+        else:
+            return mag_x_set, mag_y_set, mag_z_set, mag_x_m, mag_y_m, mag_z_m
+
+    @staticmethod
+    def fit_ellipsoid_old(mag_x_m, mag_y_m, mag_z_m):
+        # Script is adapted version from ThePoorEngineer - Calibrating the magnetometer
+        # https://thepoorengineer.com/en/calibrating-the-magnetometer/#Calibration
+        a1 = mag_x_m ** 2
+        a2 = mag_y_m ** 2
+        a3 = mag_z_m ** 2
+        a4 = 2 * np.multiply(mag_y_m, mag_z_m)
+        a5 = 2 * np.multiply(mag_x_m, mag_z_m)
+        a6 = 2 * np.multiply(mag_x_m, mag_y_m)
+        a7 = 2 * mag_x_m
+        a8 = 2 * mag_y_m
+        a9 = 2 * mag_z_m
+        a10 = np.ones(len(mag_x_m)).T
+        d_mat = np.array([a1, a2, a3, a4, a5, a6, a7, a8, a9, a10])
+
+        # Equation 7, k = 4
+        c1 = np.array([[-1, 1, 1, 0, 0, 0],
+                       [1, -1, 1, 0, 0, 0],
+                       [1, 1, -1, 0, 0, 0],
+                       [0, 0, 0, -4, 0, 0],
+                       [0, 0, 0, 0, -4, 0],
+                       [0, 0, 0, 0, 0, -4]])
+
+        # Equation 11
+        s_mat = np.matmul(d_mat, d_mat.T)
+        s11 = s_mat[:6, :6]
+        s12 = s_mat[:6, 6:]
+        s21 = s_mat[6:, :6]
+        s22 = s_mat[6:, 6:]
+
+        # Equation 15, find eigenvalue and vector
+        # Since s_mat is symmetric, s12.T = s21
+        tmp = np.matmul(np.linalg.inv(c1), s11 - np.matmul(s12, np.matmul(np.linalg.inv(s22), s21)))
+        eigen_value, eigen_vector = np.linalg.eig(tmp)
+        u1 = eigen_vector[:, np.argmax(eigen_value)]
+
+        # Equation 13 solution
+        u2 = np.matmul(-np.matmul(np.linalg.inv(s22), s21), u1)
+
+        # Total solution
+        u = np.concatenate([u1, u2]).T
+
+        q_mat = np.array([[u[0], u[5], u[4]],
+                          [u[5], u[1], u[3]],
+                          [u[4], u[3], u[2]]])
+
+        n = np.array([[u[6]],
+                      [u[7]],
+                      [u[8]]])
+
+        d = u[9]
+
+        return q_mat, n, d
+
+    @staticmethod
+    def fit_ellipsoid(mag_x_m, mag_y_m, mag_z_m):
+        """ Estimate ellipsoid parameters from a set of points.
+            Parameters
+            ----------
+            mag_x_m, mag_y_m, mag_z_m : array_like, array_like, array_like
+              The samples (M,N) where M=3 (x,y,z) and N=number of samples.
+            Returns
+            -------
+            s : array_like
+              The samples (M,N) where M=3 (x,y,z) and N=number of samples.
+            Returns
+            -------
+            M, n, d : array_like, array_like, float
+              The ellipsoid parameters M, n, d.
+            References
+            ----------
+            .. [1] Qingde Li; Griffiths, J.G., "Least squares ellipsoid specific
+               fitting," in Geometric Modeling and Processing, 2004.
+               Proceedings, vol., no., pp.335-340, 2004
+            .. https://github.com/nliaudat/magnetometer_calibration/blob/main/calibrate.py
+        """
+
+        # Converts to samples (M,N) where M=3 (x,y,z) and N=number of samples.
+        s = np.array([mag_x_m, mag_y_m, mag_z_m])
+
+        # d (samples)
+        d = np.array([s[0] ** 2., s[1] ** 2., s[2] ** 2.,
+                      2. * s[1] * s[2], 2. * s[0] * s[2], 2. * s[0] * s[1],
+                      2. * s[0], 2. * s[1], 2. * s[2], np.ones_like(s[0])])
+
+        # s, s_11, s_12, s_21, s_22 (eq. 11)
+        s = np.dot(d, d.T)
+        s_11 = s[:6, :6]
+        s_12 = s[:6, 6:]
+        s_21 = s[6:, :6]
+        s_22 = s[6:, 6:]
+
+        # c (Eq. 8, k=4)
+        c = np.array([[-1, 1, 1, 0, 0, 0],
+                      [1, -1, 1, 0, 0, 0],
+                      [1, 1, -1, 0, 0, 0],
+                      [0, 0, 0, -4, 0, 0],
+                      [0, 0, 0, 0, -4, 0],
+                      [0, 0, 0, 0, 0, -4]])
+
+        # v_1 (eq. 15, solution)
+        e = np.dot(np.linalg.inv(c),
+                   s_11 - np.dot(s_12, np.dot(np.linalg.inv(s_22), s_21)))
+
+        e_w, e_v = np.linalg.eig(e)
+
+        v_1 = e_v[:, np.argmax(e_w)]
+        if v_1[0] < 0:
+            v_1 = -v_1
+
+        # v_2 (eq. 13, solution)
+        v_2 = np.dot(np.dot(-np.linalg.inv(s_22), s_21), v_1)
+
+        # Quadratic-form parameters
+        m = np.array([[v_1[0], v_1[3], v_1[4]],
+                      [v_1[3], v_1[1], v_1[5]],
+                      [v_1[4], v_1[5], v_1[2]]])
+        n = np.array([[v_2[0]],
+                      [v_2[1]],
+                      [v_2[2]]])
+        d = v_2[3]
+
+        return m, n, d
+
+    def fit_ellipsoid(mag_x_m, mag_y_m, mag_z_m):
+        """ Estimate ellipsoid parameters from a set of points.
+            Parameters
+            ----------
+            mag_x_m, mag_y_m, mag_z_m : array_like, array_like, array_like
+              The samples (M,N) where M=3 (x,y,z) and N=number of samples.
+            Returns
+            -------
+            s : array_like
+              The samples (M,N) where M=3 (x,y,z) and N=number of samples.
+            Returns
+            -------
+            M, n, d : array_like, array_like, float
+              The ellipsoid parameters M, n, d.
+            References
+            ----------
+            .. [1] Qingde Li; Griffiths, J.G., "Least squares ellipsoid specific
+               fitting," in Geometric Modeling and Processing, 2004.
+               Proceedings, vol., no., pp.335-340, 2004
+            .. https://github.com/nliaudat/magnetometer_calibration/blob/main/calibrate.py
+        """
+
+        # Converts to samples (M,N) where M=3 (x,y,z) and N=number of samples.
+        s = np.array([mag_x_m, mag_y_m, mag_z_m])
+
+        # d (samples)
+        d = np.array([s[0] ** 2., s[1] ** 2., s[2] ** 2.,
+                      2. * s[1] * s[2], 2. * s[0] * s[2], 2. * s[0] * s[1],
+                      2. * s[0], 2. * s[1], 2. * s[2], np.ones_like(s[0])])
+
+        # s, s_11, s_12, s_21, s_22 (eq. 11)
+        s = np.dot(d, d.T)
+        s_11 = s[:6, :6]
+        s_12 = s[:6, 6:]
+        s_21 = s[6:, :6]
+        s_22 = s[6:, 6:]
+
+        # c (Eq. 8, k=4)
+        c = np.array([[-1, 1, 1, 0, 0, 0],
+                      [1, -1, 1, 0, 0, 0],
+                      [1, 1, -1, 0, 0, 0],
+                      [0, 0, 0, -4, 0, 0],
+                      [0, 0, 0, 0, -4, 0],
+                      [0, 0, 0, 0, 0, -4]])
+
+        # v_1 (eq. 15, solution)
+        e = np.dot(np.linalg.inv(c),
+                   s_11 - np.dot(s_12, np.dot(np.linalg.inv(s_22), s_21)))
+
+        e_w, e_v = np.linalg.eig(e)
+
+        v_1 = e_v[:, np.argmax(e_w)]
+        if v_1[0] < 0:
+            v_1 = -v_1
+
+        # v_2 (eq. 13, solution)
+        v_2 = np.dot(np.dot(-np.linalg.inv(s_22), s_21), v_1)
+
+        # Quadratic-form parameters
+        m = np.array([[v_1[0], v_1[3], v_1[4]],
+                      [v_1[3], v_1[1], v_1[5]],
+                      [v_1[4], v_1[5], v_1[2]]])
+        n = np.array([[v_2[0]],
+                      [v_2[1]],
+                      [v_2[2]]])
+        d = v_2[3]
+
+        return m, n, d
+
+    @staticmethod
+    def plot_magnetometer_calibration(target_column, mag_x_set, mag_y_set, mag_z_set, mag_x_m, mag_y_m, mag_z_m,
+                                      cal_x, cal_y, cal_z, mag_amp_avg_set):
+        plot_fontsize = 5
+        ax_width = 0.2
+
+        # Plot frame (overwrite plotframe)
+        plot_frame = LabelFrame(target_column, text="Result plots:")
+        plot_frame.grid(row=1, column=0, padx=(100, 0), pady=20, sticky="nw")
+
+        # Plot calibrated results
+        fig1 = plt.figure('MGM_cal_complete_left', figsize=(2.5, 3), dpi=100)
+        fig1.clf()  # clear figure from previous use
+        canvas1 = FigureCanvasTkAgg(fig1, plot_frame)
+        u_tesla = 10 ** 6  # Tesla to microTesla conversion factor
+        meas_no = len(mag_x_set)  # Measurement number
+
+        # uncalibrated result plot
+        ax1 = fig1.add_subplot(221, projection='3d')
+        # ax1.clf()
+        ax1.set_xlabel(r'$B_x [{\mu}T]$', fontsize=plot_fontsize)
+        ax1.set_ylabel(r'$B_y [{\mu}T]$', fontsize=plot_fontsize)
+        ax1.set_zlabel(r'$B_z [{\mu}T]$', fontsize=plot_fontsize)
+        ax1.tick_params(axis='both', labelsize=plot_fontsize)
+        for axis in ['top', 'bottom', 'left', 'right']:
+            ax1.spines[axis].set_linewidth(ax_width)
+        ax1.xaxis.set_tick_params(width=ax_width)
+        ax1.yaxis.set_tick_params(width=ax_width)
+        ax1.zaxis.set_tick_params(width=ax_width)
+        ax1.grid(which='major', linewidth=ax_width)
+        ax1.grid(which='minor', linewidth=ax_width / 2)
+        # linking lines between measured and calibrated points
+        for i, j, k, l, m, n in zip(mag_x_set * u_tesla, mag_y_set * u_tesla, mag_z_set * u_tesla, mag_x_m * u_tesla,
+                                    mag_y_m * u_tesla,
+                                    mag_z_m * u_tesla):
+            ax1.plot([i, l], [j, m], [k, n], linewidth=0.2, color='b')
+        # set magnetic vectors
+        l1 = ax1.scatter(mag_x_set * u_tesla, mag_y_set * u_tesla, mag_z_set * u_tesla, s=1, color='k',
+                         label="$B_{set}$ Helmholtz field")
+        # measured values
+        l2 = ax1.scatter(mag_x_m * u_tesla, mag_y_m * u_tesla, mag_z_m * u_tesla, s=1, color='b',
+                         label="$B_{raw}$ uncalibrated")
+        # plot sphere with magnitude
+        u = np.linspace(0, 2 * np.pi, 100)
+        v = np.linspace(0, np.pi, 100)
+        x = np.outer(np.cos(u), np.sin(v)) * mag_amp_avg_set
+        y = np.outer(np.sin(u), np.sin(v)) * mag_amp_avg_set
+        z = np.outer(np.ones(np.size(u)), np.cos(v)) * mag_amp_avg_set
+        ax1.plot_wireframe(x * u_tesla, y * u_tesla, z * u_tesla, rstride=10, cstride=10, alpha=0.7, color='y',
+                           linewidth=0.1)
+        ax1.plot_surface(x * u_tesla, y * u_tesla, z * u_tesla, alpha=0.3, color='y', linewidth=0.1)
+        # ax1.legend(loc='upper right', fontsize=plot_fontsize)
+
+        # Calibrated result plot
+        ax2 = fig1.add_subplot(223, projection='3d')
+        ax2.set_xlabel(r'$B_x [\mu T]$', fontsize=plot_fontsize)
+        ax2.set_ylabel(r'$B_y [\mu T]$', fontsize=plot_fontsize)
+        ax2.set_zlabel(r'$B_z [\mu T]$', fontsize=plot_fontsize)
+        ax2.tick_params(axis='both', labelsize=plot_fontsize)
+        for axis in ['top', 'bottom', 'left', 'right']:
+            ax2.spines[axis].set_linewidth(ax_width)
+        ax2.xaxis.set_tick_params(width=ax_width)
+        ax2.yaxis.set_tick_params(width=ax_width)
+        ax2.zaxis.set_tick_params(width=ax_width)
+        ax2.grid(which='major', linewidth=ax_width)
+        ax2.grid(which='minor', linewidth=ax_width / 2)
+        # linking lines between measured and calibrated points
+        for i, j, k, l, m, n in zip(mag_x_set * u_tesla, mag_y_set * u_tesla, mag_z_set * u_tesla, cal_x * u_tesla,
+                                    cal_y * u_tesla,
+                                    cal_z * u_tesla):
+            ax2.plot([i, l], [j, m], [k, n], linewidth=0.2, color='r')
+        # set magnetic vectors
+        l3 = ax2.scatter(mag_x_set * u_tesla, mag_y_set * u_tesla, mag_z_set * u_tesla, s=1, color='k',
+                         label="$B_{set}$ Helmholtz field")
+        # calibrated values ellipsoid fit
+        l4 = ax2.scatter(cal_x * u_tesla, cal_y * u_tesla, cal_z * u_tesla, s=1, color='r',
+                         label="$B_{cal,el}$ ellipsoid fit")
+        # plot sphere with magnitude
+        u = np.linspace(0, 2 * np.pi, 100)
+        v = np.linspace(0, np.pi, 100)
+        x = np.outer(np.cos(u), np.sin(v)) * mag_amp_avg_set
+        y = np.outer(np.sin(u), np.sin(v)) * mag_amp_avg_set
+        z = np.outer(np.ones(np.size(u)), np.cos(v)) * mag_amp_avg_set
+        ax2.plot_wireframe(x * u_tesla, y * u_tesla, z * u_tesla, rstride=10, cstride=10, alpha=0.7, color='y',
+                           linewidth=0.1)
+        ax2.plot_surface(x * u_tesla, y * u_tesla, z * u_tesla, alpha=0.3, color='y', linewidth=0.1)
+        # ax2.legend(loc='upper right', fontsize=plot_fontsize)
+        ax3 = fig1.add_subplot(222)
+        ax3.axis('off')
+        ax3.legend([l1, l2], ["$B_{set}$", "$B_{raw}$"],
+                   loc='right', fontsize=plot_fontsize)
+        ax4 = fig1.add_subplot(224)
+        ax4.axis('off')
+        ax4.legend([l3, l4], ["$B_{set}$", "$B_{cal,el}$"],
+                   loc='right', fontsize=plot_fontsize)
+        fig1.subplots_adjust(bottom=0.15, left=0.05, right=0.95, top=1.0)
+        canvas1.draw()
+        canvas1.get_tk_widget().grid(row=0, column=0)
+
+        # 2d_math plots
+        fig2 = plt.figure('MGM_cal_complete_right', figsize=(4, 3), dpi=100)
+        fig2.clf()  # clear figure from previous use
+        canvas2 = FigureCanvasTkAgg(fig2, plot_frame)
+        # x panel
+        ax5 = fig2.add_subplot(311)
+        ax5.grid(which='major', linewidth=ax_width)
+        ax5.grid(which='minor', linewidth=ax_width / 2)
+        ax5.set_ylabel(r'$B_x [\mu T]$', fontsize=plot_fontsize)
+        ax5.plot(np.linspace(1, meas_no, meas_no), mag_x_set * u_tesla, linewidth=0.2, color='k', linestyle='solid',
+                 label=r'$B_{x,set}$')
+        ax5.plot(np.linspace(1, meas_no, meas_no), mag_x_m * u_tesla, linewidth=0.2, color='b', linestyle='solid',
+                 label=r'$B_{x,m}$')
+        ax5.plot(np.linspace(1, meas_no, meas_no), cal_x * u_tesla, linewidth=0.2, color='r', linestyle='solid',
+                 label=r'$B_{x,el}$')
+        ax5.legend(loc='upper right', fontsize=plot_fontsize)
+        ax5.tick_params(axis='both', labelsize=plot_fontsize)
+        ax5.axes.get_xaxis().set_visible(False)
+        for axis in ['top', 'bottom', 'left', 'right']:
+            ax5.spines[axis].set_linewidth(ax_width)
+        ax5.xaxis.set_tick_params(width=ax_width)
+        ax5.yaxis.set_tick_params(width=ax_width)
+        # y panel
+        ax6 = fig2.add_subplot(312)
+        ax6.grid(which='major', linewidth=ax_width)
+        ax6.grid(which='minor', linewidth=ax_width / 2)
+        ax6.set_ylabel(r'$B_y [\mu T]$', fontsize=plot_fontsize)
+        ax6.plot(np.linspace(1, meas_no, meas_no), mag_y_set * u_tesla, linewidth=0.2, color='k', linestyle='solid',
+                 label=r'$B_{y,set}$')
+        ax6.plot(np.linspace(1, meas_no, meas_no), mag_y_m * u_tesla, linewidth=0.2, color='b', linestyle='solid',
+                 label=r'$B_{y,m}$')
+        ax6.plot(np.linspace(1, meas_no, meas_no), cal_y * u_tesla, linewidth=0.2, color='r', linestyle='solid',
+                 label=r'$B_{y,el}$')
+        ax6.legend(loc='upper right', fontsize=plot_fontsize)
+        ax6.tick_params(axis='both', labelsize=plot_fontsize)
+        ax6.axes.get_xaxis().set_visible(False)
+        for axis in ['top', 'bottom', 'left', 'right']:
+            ax6.spines[axis].set_linewidth(ax_width)
+        ax6.xaxis.set_tick_params(width=ax_width)
+        ax6.yaxis.set_tick_params(width=ax_width)
+        # z panel
+        ax7 = fig2.add_subplot(313)
+        ax7.grid(which='major', linewidth=ax_width)
+        ax7.grid(which='minor', linewidth=ax_width / 2)
+        ax7.set_xlabel("Measurement number", fontsize=plot_fontsize)
+        ax7.set_ylabel(r'$B_z [\mu T]$', fontsize=plot_fontsize)
+        ax7.plot(np.linspace(1, meas_no, meas_no), mag_z_set * u_tesla, linewidth=0.2, color='k', linestyle='solid',
+                 label=r'$B_{z,set}$')
+        ax7.plot(np.linspace(1, meas_no, meas_no), mag_z_m * u_tesla, linewidth=0.2, color='b', linestyle='solid',
+                 label=r'$B_{z,m}$')
+        ax7.plot(np.linspace(1, meas_no, meas_no), cal_z * u_tesla, linewidth=0.2, color='r', linestyle='solid',
+                 label=r'$B_{z,el}$')
+        ax7.legend(loc='upper right', fontsize=plot_fontsize)
+        ax7.tick_params(axis='both', labelsize=plot_fontsize)
+        for axis in ['top', 'bottom', 'left', 'right']:
+            ax7.spines[axis].set_linewidth(ax_width)
+        ax7.xaxis.set_tick_params(width=ax_width)
+        ax7.yaxis.set_tick_params(width=ax_width)
+        canvas2.draw()
+        canvas2.get_tk_widget().grid(row=0, column=1)
+
+    # Function passed to scipy for the optimization
+    @staticmethod
+    def residual_function(x, b_e, samples):
+        # Unpack vector. These unknown parameters are described in the calibration paper
+        s, alpha_e, alpha, beta = x
+        # Residual vector
+        res = []
+        for sample in samples:
+            # Unpack row coefficients:
+            m = sample['m']
+            b_x = sample['b_x']
+            b_y = sample['b_y']
+            b_z = sample['b_z']
+            res.append(m - s * (b_e * sin(alpha_e) + b_x * cos(alpha) * cos(beta)
+                                + b_y * cos(alpha) * sin(beta) + b_z * sin(alpha)))
+        return res
+
+    @staticmethod
+    def fibonacci_sphere(samples):
+        """
+        Algorithm to generate roughly equally spaced points on a sphere
+        From https://stackoverflow.com/a/26127012"""
+        points = []
+        phi = pi * (3.0 - sqrt(5.0))  # golden angle in radians
+
+        for i in range(samples):
+            y = 1 - (i / float(samples - 1)) * 2  # y goes from 1 to -1
+            radius = sqrt(1 - y * y)  # radius at y
+
+            theta = phi * i  # golden angle increment
+
+            x = cos(theta) * radius
+            z = sin(theta) * radius
+
+            points.append(np.array([x, y, z]))
+
+        return points
+
+    def put_message(self, command, arg):
+        self.view_queue.put({'cmd': command, 'arg': arg})

src/exceptions.py

@@ -11,3 +11,8 @@ class DeviceAccessError(Exception):
 class DeviceBusy(DeviceAccessError):
     """Error thrown when the HW proxy (i.e. access) cannot be acquired"""
     pass
+
+
+class MagFieldOutOfBounds(Exception):
+    """Set magnetic field must be a positive number between 0 and 200µT!"""
+    pass

src/globals.py

@@ -50,3 +50,6 @@ default_psu_config = {
 # Configuration for socket interface
 SOCKET_PORT = 6677
 SOCKET_MAX_CONNECTIONS = 5
+
+# Hardware safety limits
+MAG_MAG_FIELD = 5*37*1e-6  # [µT=A*µT/A] min coil constant x PSU  current limit

src/utility.py

@@ -1,6 +1,6 @@
 import csv
 from tkinter import filedialog
-
+import numpy as np
 import src.globals as g
 
 
@@ -31,4 +31,21 @@ def save_dict_list_to_csv(filename, data, query_path=False):
 
         csv_writer.writeheader()
         for row in data:
-            csv_writer.writerow(row)
+            csv_writer.writerow(row)
+
+
+def load_dict_list_from_csv(filename, query_path=False):
+    """ Reads a csv file under the specified path containing one row for each dict in the list 'data'.
+    The file header containing the keys of the first dict entry is deleted upon reading.
+    Each dict should use the same keys."""
+    if query_path:
+        filename = filedialog.askopenfilename(initialfile=filename, title="Select csv file location...",
+                                              filetypes=(("CSV", "*.csv"),))
+    data = np.genfromtxt(filename, dtype=float, delimiter=',')
+    data = data[1:len(data), :]  # remove header
+    # Save data point to raw_data list
+    raw_data = []
+    for i in range(data.shape[0]):
+        raw_data.append({'applied_x': data[i][0], 'applied_y': data[i][1], 'applied_z': data[i][2],
+                         'measured_x': data[i][3], 'measured_y': data[i][4], 'measured_z': data[i][5]})
+    return data, raw_data, filename

tools/fgm3d_adapter.py

@@ -1,11 +1,15 @@
 import serial
+import time
 from datetime import datetime
 import socket
 
 # Lookup table for cage axis to magnetometer transformation
-# First entry: Magnetometer y axis to cage x (reversed)
-# axis_mapping = [[1, -1 ...]
-axis_mapping = [[2, -1], [0, 1], [1, -1]]
+# First entry: Magnetometer z axis to cage x (reversed) [1,1]
+# axis_mapping = [[corresponding axis (0=x,1=y,2=z), direction (+-1) ...]
+# axis_mapping = [[2, -1], [0, 1], [1, -1]] # Z_MGM = -X_HH, X_MGM = +Y_HH, Y_MGM = -Z_HH,
+axis_mapping = [[0, -1],  # X_MGM = -X_HH
+                [2, -1],  # Y_MGM = -Z_HH
+                [1, -1]]  # Z_MGM = -Y_HH,
 
 # Helmholtz control software tcp port
 s = socket.socket(socket.AF_INET, socket.SOCK_STREAM)
@@ -15,16 +19,25 @@ s.sendall("declare_api_version 2\n".encode())
 # FGM3D software virtual serial port
 ser = serial.Serial("COM11")
 
+# check whether day time saving is current or not
+if time.localtime().tm_isdst:
+    dst = 2
+else:
+    dst = 1
+
 line = b""
 ready = False
 while True:
     line += ser.read()
     if line[-2:] == b"\r\n":
         new_line = line[:-2].decode('ascii')
-        delta = datetime.now().timestamp()*1000 - int(new_line.split(';')[0])+(2*60*60*1000)
+
+        delta = datetime.now().timestamp()*1000 - int(new_line.split(';')[0])+(dst*60*60*1000)
         if delta < 500 and not ready:
             ready = True
-            print("Program ready!")
+            print("FGM3D adapter script ready!")
+        elif not ready:
+            print("FGM3D adapter script not ready!")
         if ready:
             # Data is not valid otherwise
             # Only use the x y and z values