Effects of rotation errors on goniometric measurements

Goniometric measurements are essential for the determination of many optical quantities, and quantifying the effects of errors in the rotation axes on these quantities is a complex task. In this paper, we show how a measurement model for a four-axis goniometric system can be developed to allow the effects of alignment and rotation errors to be included in the uncertainty of the measurement. We use three different computational methods to propagate the uncertainties due to several error sources through the model to the rotation angles and then to the measurement of bidirectional reflectance and integrated diffuse reflectance, a task that would otherwise be intractable. While all three methods give the same result, the GTC Python package is the simplest and intrinsically provides a full uncertainty budget, including all correlations between measurement parameters. We then demonstrate how the development of a measurement model and the use of GTC has improved our understanding of the system. As a consequence, taking advantage of negative correlations between measurements in different geometries allows us to minimise the total uncertainty in integrated diffuse reflectance, lowering the standard uncertainty from 0.0029 to 0.0015.


Introduction
Many spectrophotometric and photometric quantities depend critically on the geometry of measurement. For example, the bidirectional reflectance distribution function (BRDF) is a function of the polar and azimuthal angles of the incident light, θ i and φ i , and the polar and azimuthal angles of the detected light, θ d and φ d . Therefore, instruments designed to measure BRDF must contain at least four axes in order to sample the whole hemisphere in front of the sample.
Despite the best efforts of the instrument designers, each rotation around these axes will be subject to errors, including: • displacement of the centres of rotation of the axes; • misalignment (non-orthogonality) of the axes; • angular resolution of the axis motors; • angular accuracy of each axis; • positioning of the zero angle of each axis. * Author to whom any correspondence should be addressed.
To capture the effect of various errors on an estimate of a measurand, a measurement model must be constructed [1]. For quantities measured using instruments with multiple axes, the measurement model should account for the parameters that define the axes, and all the errors in these parameters. For such goniometric systems, this model is extremely complex, including many correlated terms. For example, if θ i and θ d are not controlled by independent rotation stages, their errors become correlated with each other. Once the measurement model has been developed, the errors can be propagated through the equations of the model to determine the effect of each error on the angles at which the measurement is made and then, subsequently, on the measured value. Not only does the development of a measurement model allow for errors in the system to be propagated through to the measured value, it also helps the user develop a more detailed understanding of the instrument.
Many different instruments have been designed for measuring BRDF over the past couple of decades [2][3][4][5][6][7][8][9][10], and a lot of work has been done to understand the effects of various errors on the measurements made using goniospectrophotometers [11][12][13]. However, because of the complexity of these systems, there has been relatively little work done on investigating the effects of rotation and misalignment errors. In fact, many of the previous analyses have assumed that the rotation stages are perfectly aligned and that the angles of incidence and detection are exactly equal to the angles set on the instrument.
There has been a small amount of work done considering alignment errors. For example, Jaanson [14] considered the effect of errors in φ i and sample offsets on BRDF measurements. Three different misalignment scenarios were considered: an error in the viewing polar angle θ d , and sample displacements in the x-and y-directions. In each of these cases, the effect on the measured flux was determined analytically. Similarly, in the measurement models for the Physikalisch-Technische Bundesanstalt gonioreflectometer [6] and the Spanish National Research Council gonioreflectometer [9], uncertainty in the positioning of the sample relative to the centre of rotation was considered. However, neither of these models propagate these errors through the rotations of the sample to the measurements of BRDF.
Li et al [15] considered the angular accuracy of the rotation stages in their gonioreflectometer, and estimated the size of the angular accuracy error. Earlier, Leloup et al [7] derived expressions for the actual positions of their sample, light source, and detector using rotation matrices, then converted these into spherical coordinates to calculate the BRDF. However, neither Li et al nor Leloup et al considered the effect of alignment errors, both assuming that the axes on their instruments were all perfectly aligned.
There has also been a lot of work done to characterise the errors in coordinate measuring machines, which are susceptible to similar misalignment errors. For example, Muralikrishnan et al [16] considered the effects of several rotation errors in x-ray computed tomography instruments. They used repeated measurements to judge the sensitivity of their instrument to these errors using a frequentist approach to estimate the effects of the various errors. Meanwhile, Gaska et al [17] suggested experiments to determine errors in five-axis coordinate systems, and developed a simulation model to correct for these [18]. He et al [19] described the various angular offset and displacement errors that five-axis machine tools are subject to.
In this paper, we develop a measurement model for a fouraxis goniometric system, taking into account all the rotational and misalignment errors. We then explain how such a measurement model can be useful in developing a better understanding of the instrument and improving uncertainties in the measurement. Section 2 introduces the measurement model, using the MSL goniospectrophotometer designed to measure diffuse reflectance as an example. Then, in section 3, we investigate how the uncertainties due to all the error sources propagate through our measurement model to the estimate of the measurand, using three methods of error propagation: Monte Carlo (MC) simulations, the analytical calculation of sensitivity coefficients, and the GTC Python package [20,21]. We show that GTC, which automatically accounts for all the correlations between the errors, is an invaluable tool for analysing such a complex system. In section 4, we discuss how the information from this measurement model can be used to get a better understanding of the instrument and to give robust uncertainty estimates.

Measurement model
A measurement model expresses a mathematical relation between the measurand and all the factors that can influence the measurement process, including all the errors. (Even if the best estimate of an error is zero, the actual value may not be. Therefore, uncertainty in the estimate must be propagated through the model to contribute to the total uncertainty in the measured value as an estimate of the measurand.) Some quantities, such as BRDF, need to be measured using instruments with multiple axes. Regardless of how much time and effort is invested in the alignment of such instruments, the alignment will never be perfect, and each movement will be subject to errors due to misalignment. These errors will influence the measured value and should be included in the measurement model. In this section, we will demonstrate the development of our measurement model for BRDF using the MSL goniospectrophotometer. Measurement models for systems with different geometries or a different number of axes can be developed in a similar manner.
The BRDF is a function of the polar and azimuthal angles of incidence and detection, θ i , φ i , θ d , φ d [22]: where Φ i and Φ d are the incident and detected fluxes, respectively, and Ω d is the solid angle of detection, which depends on both the distance between the sample and the detector, L, and the surface area of the detector, A. Note that the BRDF also depends on the wavelength of the incident flux, but this dependence has been omitted here for simplicity. The solid angle is given by (2) Figure 1 shows the geometry of the incidence and detection angles.
In the MSL goniospectrophotometer,the sample is mounted on three orthogonal rotation stages positioned in the centre of a slew ring on which the detector rotates (see figure 2). The system is controlled by varying d, the angle of the detector (the slew ring angle), and the pitch, u, yaw, v, and roll, w, of the sample via the angles of the three rotation stages. However, BRDF is generally defined in terms of θ i , φ i , θ d , and φ d . Since the rotations are controlled in a coordinate system different to that of the measurand, transformations between the two coordinate systems are required. The coordinate transformations These allow us to calculate the angles that d, u, v, and w need to be set to in order to achieve any given geometry in the BRDF space.
To keep track of the position of the sample as it rotates about the different axes, we define a Cartesian coordinate space with respect to the instrument. We define the z-axis to be along the incident beam, the y-axis to be vertical with respect to the slew ring, and the x-axis to be orthogonal to both the y-and z-axes.
The origin is the point where the incident beam meets the sample. Our rotation stages are positioned so that all three rotation axes intersect at the origin. Pitch is defined as rotation about the x-axis by an angle u, yaw is rotation about the y-axis by an angle v, and roll is rotation about the z-axis by an angle w (see figure 2). The detector rotates about the y-axis by an angle d. This rotation is independent of the yaw rotation.
Ideally, the rotation stages would rotate perfectly around the x-, y-, and z-axes. However, due to the inevitable errors in the alignment of the system, the rotations will not occur perfectly around these axes. Instead, we define two types of error to describe the alignment errors in the system: displacement errors and angular offset errors. The displacement errors describe the error in the position of the centre of rotation for each axis, while the angular offset errors describe the error in the tilt of the rotation axis.
Each axis of rotation is subject to two displacement errors and two angular offset errors-one in each of the directions orthogonal to the nominal axis of rotation (by convention, we choose there to be no error in the direction of the nominal axis of rotation). The centre of rotation, Δu, and the axis of rotation, u, are given by, for example, for the pitch rotation (here, bold Roman type is used to denote points or vectors): where E Δu y and E Δu z are the displacement errors of the rotation axis in the y-and z-directions, E δu y and E δu z are the angular offset errors, describing the tilt of the rotation axis in the y-and z-directions, andî is the unit vector pointing along the nominal pitch axis. The angles of rotation are also subject to errors due to the finite resolution and accuracy of the encoders on each axis and the errors in identifying the zero angle for each axis. The measurement equations for the angles of the stages are given by, for example, for the pitch axis: where u is the angle that the axis was set to (the nominal value), the E U quantities are the unknown errors in the angular resolution, angular accuracy, and zero position for the axis, respectively, and U is the unknown actual pitch angle. Similar equations to (7)-(9) hold for the other three axes. Rotation matrices can be used to describe the rotation of each axis of the instrument. The rotation matrix about a unit vector defining the axis of rotation, for examplê u = u x u y u z T , by an angle u is given by [23]: Note that equation (10) rotates points about the origin. Therefore, in order to achieve rotation about a displaced axis, we first need to translate the point by the displacement, apply the rotation matrix in equation (10) to the point we want to rotate, then translate the rotated point back by the displacement error. Thus, we can find the position of a point a after rotating about a displaced and misaligned axis using an operator R(û, Δu, u) defined by: where Δu is the displacement error of the axis of rotation, as given in equation (7), and R u,û is the rotation matrix given in equation (10).
To find the position of a point after rotation about three different axes, care needs to be taken to ensure that the measurement equations describe the rotations as the instrument performs them. For example, in the MSL goniospectrophotometer, the rotation stages have been designed in such a way that the roll stage is contained within the pitch stage, which, in turn, is contained within the yaw stage (see figure 2). So, for example, when we rotate the pitch stage, it carries the roll stage with it and, therefore, the orientation of the roll axis changes. Thus, the rotation operators need to be applied in the right order. First, the point needs to be rotated about the roll axis, then the pitch axis, and finally the yaw axis. (If the calculation is done in a different order, then some of the axes of rotation would need to be rotated as well.) The vector a obtained by rotating a vector a about all three axes in the MSL goniospectrophotometer can be found by rotating the points at each end of the vector and calculating the difference: where o = R(v, Δv, v)R(û, Δu, u)R(ŵ, Δw, w)o is the position of the centre of the sample (the rotated origin) after rotating about the three axes.
The orientation of the sample after its rotations can be determined by findingî ,ĵ , andk , the rotated unit vectors. For example: Figure 3 shows the positions of these vectors before and after rotation.
Because the centre of the sample moves as the sample rotates, the point where the incident beam intersects the sample will change. To find the point where the incident beam meets the sample after the rotations, we can treat the sample as a plane. Since the incident beam is defined to be along the zaxis, the point where the incident beam hits the sample is the point where the sample plane intersects the z-axis. Using vector calculus, the offset of the sample position in the z-direction, Δz, is given by: wherek z is the z-component of thek vector, as shown in figure 3. Similarly to the rotation stages, we can use a rotation matrix to describe the rotation of the detector. The final position of the detector can be found by applying the rotation matrix for the detector axis to the initial position of the detector: where R = [0, 0, R] is the initial position of the detector (corresponding to d = 0 • ), which depends on the distance between the sample and the detector, R. The distance between the sample and detector after rotations, L, can then be calculated by finding the vector between the point where the beam meets the sample and the rotated position of the detector: Once all of the sample vectors have been rotated, we can calculate the values of θ i , φ i , θ d , φ d using the following: whereL is L/L, a unit vector pointing from the sample in the direction of the detector. Equations (16)- (20) are what we need to calculate the uncertainty in the BRDF (see equations (1) and (2)) due to rotation and alignment errors. Standard uncertainties and assigned probability distributions for the alignment and angle errors. The standard uncertainties in both displacement errors for each axis and for both angular offset errors in the pitch and yaw axes are the same, so only one value has been included in the table for clarity. The roll zero position is arbitrary, so there is no zero position error for the roll axis. The length, L, between the sample and the detector has been assigned a Gaussian distribution with a standard uncertainty of 0.10 mm. Each of the errors described above needs to be characterised, either experimentally or through knowledge of the system. For the MSL system, each of the alignment and rotational errors were characterised as follows: • Displacement errors (e.g., E Δu y ): two alignment lasers were positioned outside the slew ring and aligned to cross in the centre of the slew ring, which we aligned to the origin. A target was mounted on the rotation stages so that in the zero position, the alignment lasers met in the centre of the target. Then, each rotation stage was rotated by a large angle (at least 40 • ) in turn, and the offset of the crossing point of the alignment lasers on the target was used to assess the possible range of the displacement errors to enable an uncertainty to be assigned. • Angular offset errors (e.g., E δu y ): one of the alignment lasers was positioned so that it reflected off a mirror mounted on the sample holder onto a target. Each rotation stage was rotated, one at a time, and the off-axis displacement of the reflected beam was measured as the angle changed. From the variability of these measurements, we could get a good estimate of the variability of the angular offset errors. • Resolution errors (e.g., E U resolution ): the number of steps per degree was found from the specifications for the stepper motors in each of the rotation axes. These values were used to calculate the smallest step size for each axis. We expect the resolution error to follow a uniform distribution with the width given by the smallest step size. • Accuracy errors (e.g., E U accuracy ): the angular accuracy of each axis was checked using a regular polygon and an autocollimator. For each axis, the autocollimator was aligned to retroreflect off one of the faces of the polygon. The axis was then rotated so the autocollimator was incident on the next face of the polygon, and by determining the variability in the angle that the stage had to be set to in order for the autocollimator to retroreflect, we could obtain an estimate for the uncertainty associated with the accuracy error in each axis.
• Zero position errors (e.g., E U zero ): a mirror was mounted on the sample holder and the incident beam was reflected back off the mirror, with each rotation stage set to 0 • . The zero position errors were assessed by how well we could determine that the beam was retroreflecting.
The alignment has been carefully adjusted until the best estimate for each of the displacement, angular offset, accuracy, resolution, and zero position errors is 0. Note, however, that there is still a standard uncertainty associated with each error that will contribute to the uncertainty in the measured value. We cannot know the size of the remaining bias, but the standard uncertainty associated with each error describes the likely magnitude of this bias. The standard uncertainties for the MSL goniospectrophotometer and their assigned probability distributions are given in table 1. For the errors assigned a uniform distribution, there are infinite degrees of freedom. The number of degrees of freedom for the errors assigned Gaussian distributions varies between 5 and 50.
While the analysis here is specific to the MSL goniospectrophotometer, equations (7)-(11) describe the position and rotation of any axis, taking into account the errors discussed. These can be used in any system and propagated through to uncertainty in BRDF or other measured quantity.

Propagating uncertainty
Equations (7)- (20) are the measurement model for θ i , φ i , θ d , and φ d . By propagating the uncertainties in all of the input parameters through these equations, we can calculate the uncertainty in θ i , φ i , θ d , and φ d . We have used three different methods of uncertainty propagation and compared the results. Firstly, the sensitivity coefficients of the outputs to each of the error contributions were derived by analytical differentiation using Mathematica (e.g., ∂θ i /∂E U accuracy ), then the methods specified in the GUM [1] were used to calculate the uncertainties. Secondly, probability distributions were assigned to each of the sources of error (see table 1) and MC simulations were carried out by randomly selecting values for each error from those distributions and calculating the standard deviation of the resulting calculated values of the angles. Thirdly, the GTC Python tool [20,21] for the automatic propagation of uncertainty was used, which takes inputs called uncertain numbers that consist of a value with the corresponding standard uncertainty and number of degrees of freedom. GTC then uses automatic differentiation to apply the GUM methods to calculate the propagated uncertainties, while carrying the information about all contributions (including correlations) through the calculations.
We found that all three methods produce the same propagated uncertainties. For example, figure 4 shows the resulting uncertainties in θ i , φ i , θ d , and φ d as a function of the yaw angle, using each of the three methods, when the pitch angle was set to u = 30 • , the roll angle was set to w = 0 • , and the detector angle was set to d = 30 • . The two analytic solutions, using derived sensitivity coefficients and using GTC, give the exact same solutions, so are plotted directly on top of one another. The uncertainties calculated using MC simulations vary randomly around the analytic solution, where the deviation from the analytic solution is dependent on the number of simulations run. The figure shows the results with 10 000 simulations. The standard uncertainties in each of the input parameters are given in table 1. Similar plots could be made for different combinations of incident geometries, or different values for the alignment and rotational uncertainties.
While all three methods of uncertainty propagation gave the same results, we found that GTC was the easiest to implement, could give us the correlations between quantities (e.g., u(θ i , φ i )) without any additional effort, and could produce an uncertainty budget for the final estimate. When integrating the BRDF to calculate diffuse reflectance, the covariance between the BRDF measured at one angle and that measured at another angle must be known for a robust estimate of the uncertainty. In this case, the GTC package becomes invaluable, as all the information about each uncertainty contribution is carried through the calculations intrinsically, so all correlations are automatically taken into account with no additional effort.
On the other hand, the MC simulations required a large number of runs to obtain a satisfactory approximation for the distribution of measured values, which needed a lot of computational time to complete. Furthermore, determining the independent contribution of each error component on the final value for diffuse reflectance would need to be done independently for each input error. Similarly, when analytically evaluating sensitivity coefficients and using the GUM methods to propagate uncertainty, the determination of correlations between components for any further processing would require a lot of care and additional work. Therefore, all uncertainties in the rest of this paper have been calculated using the GTC package. The advantages of being able to look at the correlations between quantities and the contributions of each individual error on the final uncertainty are expanded on in section 4.

Application
A measurement model and method for propagating uncertainties through the model, like that described above, gives a complete and robust uncertainty budget for our measured values. It can also be a useful tool in developing a deeper understanding of the instrument. There are two 'poses' that the MSL goniospectrophotometer can take for each combination of θ i , φ i , θ d , and φ d . For example, for in-plane geometries, if we move the detector to the other side of the incident beam and roll the sample by 180 • , the same geometry is realised. The original design of the MSL instrument minimises the travel around the detector axis when selecting the pose. We can examine the impact of this default decision by looking at the uncertainties.
As an example, figures 5(a) and (b) show the uncertainty in θ i and θ d , respectively, for the MSL goniospectrophotometer, as θ d and φ d vary over the whole hemisphere, with the incident angles fixed at θ i = 45 • and φ i = 0 • . Note that in order to achieve any given θ d , φ d geometry, the u, v, w, and d angles may all change. Since the uncertainty in θ i depends on u, v, w, and d, its uncertainty changes even though its value remains constant in figure 5(a).
There are a couple of main points to note when looking at the uncertainty in θ i in figure 5(a). First, the area of low uncertainty (the lighter area) occurs when θ d is larger than θ i . Only the detector needs to move to achieve these geometries, so there is no change in the uncertainty in θ i . However, if θ d is less than θ i , or if φ d > 0 • , our system rolls the sample by default, to minimise movement of the detector. Figure 5(a) clearly shows an increase in uncertainty in this case, which arises from the angular offset uncertainty in the x-direction of the roll axis. Similarly, figure 5(b) shows the uncertainty in θ d for the same geometries. Again, we can see the increase in uncertainty as we introduce a roll to the sample.
Through analysis like this, we can improve the measurement process so that the uncertainties are minimised. For example, from figures 5(a) and (b), we can identify the dominance of the roll angular offset error on the uncertainty in angles achieved by the system. Therefore, instead of picking the default pose, we can choose the pose of the system that gives the minimum uncertainty, thus allowing us to achieve lower uncertainties in our final measurements. Figure 6 shows the standard uncertainties in θ i and θ d when choosing the pose of the system to minimise the roll, and we can see that there is a significant drop in the uncertainty for large parts of the hemisphere for both θ i and θ d .
Ultimately, the main point of this analysis is to determine the contribution of the rotation and alignment errors to the final measured value of, in our case, the BRDF or the integral of some part of the BRDF (e.g., to determine hemispherical reflectance). Equation (1) gives the measurement equation for BRDF, and by propagating the uncertainties of the angles and the distance to the detector through to BRDF, we can look at how the alignment errors in the rotation stages contribute to uncertainty in BRDF. Note that the dependence of the measured flux, Φ d , on θ i , φ i , θ d , and φ d must be determined empirically for each sample.
As an example, figure 7 shows the BRDF shape and the relative uncertainty in the BRDF for a perfect Lambertian material and a material with a sharp peak in its BRDF, both measured with 0 • : x • (θ i = 0 • , θ d = x • ) geometry. The solid curve and the dotted curve show the BRDFs and the two dashed curves show the relative standard uncertainty in BRDF. For the Lambertian material, the flux, Φ d , is independent of the incident  angle, and has a cosine dependence on θ d , so we can see that the propagated uncertainty increases as θ d increases. For these measurements, only the detector angle changes; however, the uncertainty in both θ d and the length L, as well as the sensitivity coefficient ∂Φ d ∂θ d , increase with the detector angle, d. The largest contributions from rotational and alignment errors are due to the zero position errors in the detector and yaw axes and the accuracy error in the detector axis.
For the non-Lambertian material, the signal varies with θ d as cos(θ d ) ∂β ∂θ d and with θ i as ∂β ∂θ i (not shown). Note that the errors in θ d and θ i are correlated because both depend on the settings of the instrument in u, v, w, d space. Accounting for this (and all other correlations between angles) appears almost intractable, but the development of the measurement model and the use of the GTC package enable us to solve the problem relatively easily. We find for the non-Lambertian material that the final uncertainty in BRDF at small angles is large because ∂β ∂θ d is large. At larger θ d , the uncertainty increases again as the uncertainties in θ d and L grow with the detector angle.
Analysing the different contributions of the rotation and alignment errors to the uncertainty in the measured value also helps find methodical ways to lower the final uncertainty. The first thing is to be able to quantify the effect of rotation and alignment errors relative to other sources of error, such as sample uniformity, temperature sensitivity, wavelength, or angular resolution. If the rotation contributions are large, then it is worth investing in improvements to the stages, control gear, and alignment capability; otherwise other contributions can be addressed first. Specific aspects of the rotation stages can be worked on as it can easily be identified which error, on which axis, needs to be better characterised or corrected.
Secondly, measurement schemes that take advantage of knowledge about correlations between errors can be employed to reduce the final uncertainty. For example, when calculating the hemispherical reflectance of a sample, we need to integrate BRDF values measured over the full hemisphere. We showed earlier that certain poses give lower uncertainty than other poses, so it might seem that selecting those poses would be the best approach. However, some errors are negatively correlated between the two possible poses, so by alternating the measurement steps between one pose selection and the other (e.g., making measurements at θ d = 5 • , 15 • , . . . in one pose and measurements at θ d = 10 • , 20 • , . . . in the other pose), we can achieve a lower uncertainty than by making all the measurements in the 'lower uncertainty' pose. Using the complete measurement model, and having access to the uncertainty budget of the final result (easily accessed when using GTC), we can quantify and track these effects and select the best measurement scheme.
In the case of a Lambertian sample, the standard uncertainty in the diffuse reflectance due to alignment and rotational errors, with the same input errors described in table 1, is 0.0029 if all measurements of BRDF are made on one side of the instrument. This uncertainty is dominated again by the zero position and accuracy errors in the yaw and detector axes. If, instead, we were to measure the BRDF by taking measurements alternately on both sides, using the same number of total measurement points, and then integrate, the standard uncertainty decreases to 0.0015, as the zero position errors are negatively correlated on the two sides of our instrument. This example also demonstrates the value of being able to carry information about correlations between measurements that share common errors through the calculations and, therefore, produce robust estimates of uncertainty.

Conclusion
Alignment errors of instruments with multiple rotation axes contribute to the total uncertainty in the measured values of radiometric and photometric quantities. This work has shown that even quite small uncertainties in the rotation and alignment of rotation stages can accumulate, their effects can vary strongly with the configuration of the stages, and they can propagate through the measurements to have significant contributions in the final uncertainty budget. A complete understanding of this was achieved by building a measurement model that explicitly relates the measured value to the physical sources of error in the rotation stages. Here we have considered displacement, angular offset, accuracy, resolution, and zero offset errors. We have used the GTC Python package to input and track the uncertainties through the measurement equations of the model so that the contribution of each can be identified clearly in the uncertainty budget of the measured value. This obviates the need to derive complicated sensitivity coefficients.
We have also demonstrated how understanding the ways in which the various errors influence the measured value can give a much better understanding of an instrument. Using the MSL goniospectrophotometer as an example, we showed how the construction of the measurement model allowed us to improve our uncertainties in BRDF and hemispherical reflectance. By looking at the effects of each input error, we were able to select poses and measurement schemes to optimise the uncertainty achieved. Moreover, the uncertainty analysis of all future measurements has become very simple, including accounting for all correlations between the measured angles and distances, and their combined effect on the measured value.

Acknowledgments
This work was funded by the New Zealand government and has been done in the frame of the project 18SIB03 BxDiff, that has received funding from the EMPIR programme co-financed by the Participating States and from the European Union's Horizon 2020 research and innovation programme.