*2.2. PiDOSE*

The piDOSE radiation detector is based on a PIN diode X100-7 coupled with a CsI:Tl scintillator with dimensions of 4 × 8 × 8 mm. The detector operates in particle counting mode. The energy deposition is not measured. The integration time is approximately 55 s. The data were converted to counts per minute and corrected for dead time. As shown in [36], the detector can register protons with energies higher than 30 MeV. Since the incident angle of particles and rotation of the satellite could not be determined, the counts were 4 π normalized. Hence, isotropical irradiation is assumed. Although the irradiation in the SAA is strongly directional, the satellite slowly rotates, which helps to mitigate the different directional sensitivities of piDOSE. A typical observed radiation field is shown in Figure 5.

**Figure 5.** Example of the measurement of the radiation field (number of registered particles per minute (CMP)).

#### *2.3. Position of the Measurement*

As the satellite is not equipped with an altitude determination and control system, the satellite slowly rotates. The rotation rate is determined from the current of individual

solar panels or from the periodic variation in the carrier to noise (C/N0) indicator of the GPS receiver. The rotation rate is not constant; the typical value is one turn per minute [34]. This rotation causes the GPS position to be unavailable some of the time [34]. The typical availability of PNT (position, navigation, and timing) solutions is 80% of the time. The problem was solved by extrapolation of the satellite position.

The following proposed interpolation algorithm is based on the satellite motion equation and its modification for the Earth-centered Earth-fixed (ECEF) coordinate system [37,38]:

$$\begin{aligned} \frac{d\mathbf{x}}{dt} &= V\_{\mathbf{x}}; \frac{d\mathbf{y}}{dt} = V\_{\mathbf{y}}; \frac{dz}{dt} = V\_{z} \\ \frac{d\mathbf{V\_{x}}}{dt} &= -\frac{\mu}{\mathbf{r}^{3}}\mathbf{x} - \frac{3}{2}\mathbf{J}\_{0}^{2}\frac{\mu\mathbf{x}\_{0}^{2}}{\mathbf{r}^{5}}\mathbf{x}\left(1 - \frac{5\mathbf{x}^{2}}{\mathbf{r}^{2}}\right) + \omega^{2}\mathbf{x} + 2\omega\mathbf{V\_{y}} \\ \frac{d\mathbf{V\_{y}}}{dt} &= -\frac{\mu}{\mathbf{r}^{3}}\mathbf{y} - \frac{3}{2}\mathbf{J}\_{0}^{2}\frac{\mu\mathbf{x}\_{0}^{2}}{\mathbf{r}^{5}}\mathbf{y}\left(1 - \frac{5\mathbf{x}^{2}}{\mathbf{r}^{2}}\right) + \omega^{2}\mathbf{y} + 2\omega\mathbf{V\_{x}} \\ \frac{d\mathbf{V\_{x}}}{dt} &= -\frac{\mu}{\mathbf{r}^{3}}\mathbf{z} - \frac{3}{2}\mathbf{J}\_{0}^{2}\frac{\mu\mathbf{x}\_{0}^{2}}{\mathbf{r}^{3}}\mathbf{z}\left(1 - \frac{5\mathbf{x}^{2}}{\mathbf{r}^{2}}\right) \end{aligned} \tag{1}$$

where (*x*, *y*, *z*) and *Vx*, *Vy*, *Vz* are the position and velocity vectors, respectively, of the satellite in ECEF coordinates provided by the onboard GPS receiver; *r* = x2 + y<sup>2</sup> + z<sup>2</sup> is a satellite radius; <sup>μ</sup> = 398600.4 × <sup>10</sup><sup>9</sup> <sup>m</sup>3/s2 is the standard gravitation parameter of the Earth; J 2 <sup>0</sup> = 1082625.7 × <sup>10</sup>−<sup>9</sup> is the second-harmonic coefficient of geopotential; <sup>ω</sup> = 7.292115 × <sup>10</sup>−<sup>5</sup> rad/s is the Earth rotation rate; and *ae* is the equatorial radius of the reference ellipsoid.

The satellite position is calculated by numerical integration of (1) using the fourthorder Runge–Kutta method [38].

The described algorithm enables the calculations of missing positions of the satellite and also enables the transformation of the position measurement from the end of the measurement cycle to its middle, as the satellite registers the position of the end of the measurement.

The ECEF position is, then, transformed to the geometrical coordinates LLH (longitude, latitude, and height) for further processing [38].

#### *2.4. Position of the South Atlantic Anomaly (SAA)*

The position of the SAA can be calculated as a centroid of the radiation field [13]. For this purpose, the data should be transformed to the latitude-longitude grid, and then, the centroid is calculated. The following formulas are presented based on [13]:

$$\begin{array}{l}\text{Latitude}\_{centroid} = \frac{\sum(\text{Interpolated Flux} \times \text{Flux's Littleule})}{\sum \text{Interpolated Flux}}\\\text{Longitude}\_{centroid} = \frac{\sum(\text{Interpolated Flux} \times \text{Flux's Longitud})}{\sum \text{Interpolated Flux}}\end{array}\tag{2}$$

The problem with Equation (2) is that it does not consider the grid area; therefore, the individual grid points have equal mass. The problem is that the grid area decreases with latitude. Moreover, the method does not take into account that a one-degree grid is a spherical surface, and data are processed in two-dimension (2D). The proposed method is based on centroids in the space of constant curvature defined in [39]. The mass of the individual grids is calculated as a product of the number of registered particles per minute and grid area. The coordinates of the grid are transformed from the geodetic LLA coordinates to the ECEF, then the ECEF centroid coordinates are calculated using (3) and the results are transformed back to the LLH. The adjusted formulas for calculation of the longitudinal ECEF (*xc*, *yc*, *zc*) coordinates of the centroid adjusted to the grid area size are as follows:

$$\begin{array}{c} \mathbf{x}\_{\mathbb{C}} = \frac{\frac{\sum(\text{CPM}(\mathbf{x}, \mathbf{y}, \mathbf{z}) \cdot \cos(\boldsymbol{\varrho}) \cdot \mathbf{x})}{\sum(\text{CPM}(\mathbf{x}, \mathbf{y}, \mathbf{z})}}{\sum(\text{CPM}(\mathbf{x}, \mathbf{y}, \mathbf{z}) \cdot \cos(\boldsymbol{\varrho}) \cdot \mathbf{y})}\\ y\_{\mathbb{C}} = \frac{\frac{\sum(\text{CPM}(\mathbf{x}, \mathbf{y}, \mathbf{z})}{\sum(\text{CPM}(\mathbf{x}, \mathbf{y}, \mathbf{z}) \cdot \cos(\boldsymbol{\varrho}) \cdot \mathbf{z})}{\sum(\text{CPM}(\mathbf{x}, \mathbf{y}, \mathbf{z})}} \end{array} \tag{3}$$

where *λ* and *ϕ* are coordinates of the grid in LLH, (*x*, *y*, *z*) are coordinates of the grid in ECEF, and *CPM*(*x*, *y*, *z*) is the number of counts per minute in the grid.

#### *2.5. Data Processing*

The registered scientific data are downloaded by a ground station for further processing. One dataset contains data from approximately 22 orbits. The data processing can be summarized as follows:


#### **3. Results**

The radiation measurement in the SAA region is presented in Figure 6. The figure displays the development of the radiation field expressed as the number of particles registered per minute (CPM counts per minute) over the measurement campaigns, including the positions of the radiation maximum and centroid. The SAA positions are summarized in Table 2.

**Table 2.** Development of the position of the SAA.


Figure 7 displays the north–south and east–west cross-sections of the SAA radiation field measured in the individual campaigns.

The positions of the SAA radiation field maximum and centroid positions are also displayed in Figure 8. It is evident that the position of the maximum features a much higher scatter than the position of the centroid in which the scatter is considerably lower. The details of the centroid position scatter are shown in Figure 9.

The development of the position of the centroid was interpolated by a first-order polynomial (straight line), as shown in Figure 10. The resulting polynomials for latitude *λ<sup>F</sup>* and longitude *ϕ<sup>F</sup>* are as follows:

$$
\lambda\_F = -0.00093d - 26.0229q\_F = -0.00062d - 48.65654\tag{4}
$$

where *d* is the time in days from the first measurement (30 August 2019).

The position of the centroid is, therefore, moving by 0.00093◦ to the west and 0.00062◦ to the south per day, which is 0.34◦ to the west and 0.23◦ to the south per year.

**Figure 6.** Development of the radiation field (CPM) in the SAA over measurements with the positions of the maximum of radiation (cross) and centroid (circle). (**a**) 30 August 2019; (**b**) 30 September 2019; (**c**) 27 March 2020; (**d**) 10 October 2020; (**e**) 1 November 2020; (**f**) 17 November 2020; (**g**) 28 December 2020; (**h**) 2 January 2021.

**Figure 7.** North–south and east–west cross-sections of the SAA radiation field.

**Figure 8.** Development of the SAA maximum and centroid.

**Figure 9.** Development of the SAA centroid.

**Figure 10.** Straight line interpolation of the centroid position.

#### **4. Discussion**

The position of the South Atlantic Anomaly can be defined by various methods. The most intuitive is the position of the maximum of the fluence function. Alternatively, the position can be understood as a centroid of this function. In this paper, we applied both methods.

The position of the fluence maximum observed by individual measurement campaigns fluctuates more than the fluence centroid position. The positions of the fluence maxima are more sensitive than the position of the centroid to measurement noise. Although the centroid method is relatively insensitive to the number of measurement points, it would benefit from denser sampling of the SAA region. Denser sampling can be obtained by selecting an orbit with higher inclination or by releasing a swarm of satellites.

Our results show the average westward movement of the SAA by 0.34◦/year, which is in good agreement with previously reported measurements summarized in Table 1. The data show an average southward movement rate of 0.23◦/year. Even though most of the publications show northward drift, southward drift was also observed in [2,3,21]. The inconsistency in the north–south direction can be explained by applying the improved centroid calculation method that correctly takes into account grid area, in contrast to the application of the less precise centroid position calculation method used in [13].

As shown in [2], the drift of the SAA in latitude is dependent on the energy range of the protons that are measured. Protons within the SAA with lower energy tend to drift towards the north, whereas the most energetic protons of the SAA move towards the south. In [34], piDOSE was estimated to be sensitive to protons with energies larger than 30 MeV due to the shielding of sensitive volume. Moreover, the latitudinal movement of the SAA depends on the solar cycle [21]. Often, there was a shift in direction during the solar cycle minima and maxima. Since our data captured the time of the solar cycle minimum, it can be expected to observe a similar shift in SAA movement. Another effect that influences the measurement of latitudinal movement is a systematic error caused by the fast latitudinal movement of a satellite. The particle counting time is about one minute. In this time, the satellite travels approximately four degrees in orbit which is projected mainly in the latitude. The maximal error in the determination of the latitude position of the fluence maximum within the counting time is two degrees. For the uniform distribution the mean systematic error is about one degree. The error in the east–west direction is considerably lower due to the slower movement of the satellite in longitude.
