# PicSat's Enduring Legacy. Probing the Flight of a Small Astronomical Satellite

Cauê G. Menegaldo, Fábio de O. Fialho, Eduardo Janot-Pacheco, Felipe M. Pait, Vincent Lapeyrère

► To cite this version:

Cauê G. Menegaldo, Fábio de O. Fialho, Eduardo Janot-Pacheco, Felipe M. Pait, Vincent Lapeyrère. PicSat's Enduring Legacy. Probing the Flight of a Small Astronomical Satellite. Publications of the Astronomical Society of the Pacific, 2022, 134, 10.1088/1538-3873/ac4e1d. insu-03713332

HAL Id: insu-03713332

https://insu.hal.science/insu-03713332

Submitted on 4 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# PicSat's Enduring Legacy. Probing the Flight of a Small Astronomical Satellite

Cauê G. Menegaldo<fnref n="1" /> ID, Fábio de O. Fialho<fnref n="1" />,<fnref n="2" /> ID, Eduardo Janot-Pacheco³ ID, Felipe M. Pait<fnref n="1" /> ID, and Vincent Lapeyrère<fnref n="2" /> ID

Received 2021 November 24; accepted 2022 January 17; published 2022 March 14

## Abstract

Observation of planetary transits and other cutting-edge scientific missions can take advantage of affordable nanosatellites to probe interesting stellar targets. PicSat, a CubeSat dedicated to observing the Beta Pictoris star system, was designed to provide high-precision star pointing, a critical requirement for planetary transit detection. PicSat's Attitude Determination and Control System, responsible for delivering high-accuracy spacecraft pointing, requires dedicated development based on dynamic simulators. This paper presents a dynamic attitude and orbit propagation simulator for CubeSats in low Earth orbit, as well as for its de-tumbling mode. Validation has been performed through PicSat's in-flight data. High-precision dynamic models have been obtained for both attitude and orbit. Such models are well suited to the different mission phases, from spacecraft design to data exploitation. It is, therefore, a crucial tool to minimize the chance of failure of both the platform and the payload, especially in satellites such as PicSat, whose pointing depends on both. PicSat left an enduring legacy: its platform data allow us to obtain flight models that will be valuable for future missions.

*Unified Astronomy Thesaurus concepts: Artificial satellites (68); Astronomical instrumentation (799)*

## 1. Introduction

PicSat (Nowak et al. 2018) was a 3U (3 unit) CubeSat launched in 2018 and designed to observe the Beta Pictoris star system. The system comprises $\beta$ Pictoris, a young ($\sim$12 Myr) and nearby (20 pc) star, a circumstellar dust disk, and at least two young super-Jupiter planets: $\beta$ Pictoris b with $10–11 M_{Jup}$ and $\beta$ Pictoris c with $7.8 \pm 0.4 M_{Jup}$ (Lagrange et al. 2020). An artistic representation of the system is shown in Figure 1.

A CubeSat is a small satellite composed of unit cubes of $10 \text{ cm} \times 10 \text{ cm} \times 10 \text{ cm}$, and no more than 1.33 kg ($\sim$2.9 lb) per unit. The CubeSat concept was first proposed by California Polytechnic State University and Stanford University in 1999 Puig-Suari et al. (2001). A Cubesat often employs commercial off-the-shelf (COTS) components and its development does not need to follow space standards such as those from European Cooperation for Space Standardization (ECSS). This permits fast and inexpensive development, which is particularly suitable for university space missions in Low Earth Orbit (LEO).

The early history of CubeSats is marked by technological missions. After the advent of three-axis stabilization by means of star trackers in the early 2010s, a wide variety of

applications for these spacecraft has appeared, including high pointing accuracy scientific missions. According to the Nanosats Database (Kulu 2020), until 2021 January 1, 1357 CubeSats with a wealth of purposes were launched. Due to the short development life cycle and low cost, CubeSats are becoming more popular in space science and technology. Low cost can also justify higher risks and scientific missions with CubeSats are being progressively used to test new theories and techniques.

Villela et al. (2019) analyzed 855 CubeSat missions between 2017 and 2018 and found that 23.7% of launches carried out failed during commissioning or during the beginning of life. This condition, called infant mortality, occurs mainly with university satellites due to the lack of proper ground tests and adequate simulations (Langer et al. 2017). In those cases, a successful mission may require many launches, in a “fly-learn-refly” procedure. However, state-of-the-art scientific missions involve unique instruments and a trial and error procedure is unfeasible.

In astrophysics, successful CubeSat missions are still uncommon (Section 2 presents some statistics), although their use can clearly be of great advantage not only for preparing the next generation of space instrument scientists and engineers but also for doing high-quality science. To increase satellite reliability, a compromise between the practicality of CubeSats and the complexity of space standards such as ECSS is necessary.

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Department of Telecommunications and Control Engineering, University of São Paulo, Av. Prof. Luciano Gualberto, 158, SP, 05508-080, Brazil
caue.menegaldo@gmail.com</footnote>

<footnote n="2">² LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne University, UPMC Univ. Paris 06, Univ. Paris Diderot, Sorbonne Paris City, France
³Department of Astronomy, University of São Paulo, Rua do Matão, 1226, SP, 05508-090, Brazil</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Artistic representation of the Beta Pictoris system by Nicholas W. Beeson based on an image whose credit is: Axel Quetz/MPIA Graphics Department. Licensed under CC BY-SA 4.0.</figcaption></figure>

Although mostly built with COTS components and conceived to require shorter development cycles then bigger satellites, CubeSats are composed of several subsystems that need to be specified, evaluated and tested to assure mission success. The most complex and critical platform subsystem is the Attitude Determination and Control System (ADCS), which is in charge of stabilizing and pointing the satellite during the mission. ADCS is crucial for mission success. It comprises sensors that provide relative and absolute measurements, actuators and digital filters used for state estimation and control.

Nakasuka et al. (2018) present a discussion of micro/nano/ pico-satellites survivability based on their flight experience with university microsatellites. They discuss several ADCS anomalies that may arise and some strategies to mitigate them. It is advantageous to employ some techniques used in the engineering of larger satellites to the development of this subsystem. Although less complex than a state-of-the-art Attitude and Orbit Control System for large satellites, CubeSat ADCS requires dedicated development, specially in the case of scientific missions requiring high-accuracy pointing. It must include an ADCS high-level simulator followed by an On- board In-Loop Simulator (OILS) and finally a Hardware In the Loop Simulator (HILS). Those simulators are necessary not only for designing the mission, but also for preparing and supporting post-launch activities like commissioning, calibra- tion, and mission exploitation, including ADCS flight software parameterization, filtering strategies changes, and code update.

In spite of not having accomplished its astronomical goal, PicSat left an enduring legacy: data from the platform can be used to build a flight model that will be valuable for future missions. This paper presents a dynamic attitude and orbit propagation simulator for CubeSats in Low-Earth orbit as well as its de-tumbling mode, which are part of the ADCS high- level simulator. Validating such a simulator requires in-orbit data from previous missions. In our case, the simulator is validated using in-flight data from the PicSat mission.

This paper is organized as follows. Section 2 presents CubeSat background. Section 3 provides a description of PicSat as well as some detail on its unique fine pointing system.

Attitude dynamics and orbit propagation modeling are devel-
oped in Section 4. Model identification and validation based on
PicSat in-orbit data are detailed in Section 5. ADCS de-
tumbling mode and validation based on PicSat data are
explored in Section 6. Finally, Section 7 states the conclusions.

## 2. CubeSat Background

We present hereafter a short bibliographical review about Cubesat missions, ADCS and simulators.

## 2.1. CubeSat Missions

The time is ripe for the world science community to consider that small satellites may play important roles, specially in astrophysics, planetary science, heliophysics, and Earth science (Millan et al. 2019). Indeed, the last Decadal Surveys lead by the U.S. National Academies and published by NASA<fnref n="4" /> recommended increasing the level of investments in small satellite programs, mainly in the above mentioned areas. ESA has also had many initiatives concerning small satellites in the last years.<fnref n="5" />

CubeSats do not replace larger missions. They can be used advantageously to supply complementary measurements and explore targeted science goals. An important characteristic concerning the small satellites is the “fly-learn-refly” procedure, which is often adopted for when the budget allows it. A first flight model is launched and, if any serious problem appears, a second model is launched with the proper modifications. This approach increases significantly the number of successful missions. CubeSats thus bring versatility to space exploration. In astrophysics, in spite of the small dimensions which limit optical instrument aperture, long-term monitoring can be quite useful for stellar variability and exoplanetary studies. These were the main goals of PicSat.

Constellations of CubeSats with a variety of purposes, including long-term Earth monitoring and solar and space physics studies (e.g., Liddle et al. 2020, and references therein) have been launched during the last years (Kulu 2020).

## 2.2. CubeSat ADCS and Simulators

Attitude and Orbit Determination and Control System is responsible for estimating and controlling nanosatellite attitude and orbit. It is divided into ADCS itself, and Orbit Determination and Control System. These subsystems are subdivided into more specific ones: Attitude Determination System (ADS), Attitude Control System (ACS), Orbit Determination System (ODS), and Orbit Control System (OCS). Such a division separates measurements and state estimation algorithm from control law.

<footnote n="4">⁴ See https://science.nasa.gov/about-us/science-strategy/decadal-surveys.</footnote>

<footnote n="5">⁵ See https://www.esa.int/Enabling_Support/Preparing_for_the_Future/Discovery_and_Prepparation/CubeSats.</footnote>
Since each CubeSat unit weighs no more than 1.33 kg, technology miniaturization is required, orbit control is a challenge, and three-axis stabilization took time to become common place, while high accuracy remains unusual. A survey of 426 CubeSats launched until 2016 January (Polat et al. 2016) found that 56% of these CubeSats used magnetic actuators, 44% used reaction wheels, 13% used passive magnetic control, and 4% who did not use any type of control. Only 44%, the ones using reaction wheels, could implement three-axis control. Most of the CubeSats did without star trackers, limiting pointing accuracy to the order of degrees.

Except for MarCO-A and MarCO-B, all CubeSats were launched in LEO, of which 76% were in the 350–700 km range (Polat et al. 2016; Villela et al. 2019). As pointed out in Section 1, the development of high-accuracy ADCS requires simulators. Dynamic simulation of orbit and attitude allows assessing the orientation, trajectory and environmental conditions to which the spacecraft will be subject over time. For Sebestyen et al. (2018), in LEO, high-fidelity mathematical models must contain at least: (i) rigid body dynamics; (ii) orbital mechanics model; (iii) high-fidelity magnetic field model; (iv) high-fidelity atmospheric density model; (v) torque due to atmospheric drag; (vi) torque due to the pressure of solar radiation; (vii) torque due to the gravity gradient; (viii) torque due to the residual magnetic dipole; and (ix) position of the Sun.

Several academic simulators, when compared to the expected characteristics of a high-fidelity mathematical simulator, are deficient in at least one of the above items. This occurs both in simulators for specific applications, such as ACS and/or ADS for a particular nanosatellite Krogh & Schreder 2002; Gießelmann 2006; Francois-Lavet 2010; Jensen & Vinther 2010; Holst 2014; van Vuuren 2015; Rondão 2016; Thomsen & Nielsen 2016; Avanzini et al. 2019; Jonsson 2019), and in simulators with more comprehensive proposals to serve as a development platform (Menges et al. 1998; Lovera 2003; Nasirian et al. 2006; Naqvi & Raza 2007; Triharjanto et al. 2015; Habibkhah et al. 2017; Alshamy et al. 2019; Liu et al. 2019).

There exist open source simulators and toolboxes such as ODTBX, developed by NASA, which deals only with orbital dynamics;<fnref n="6" /> KPS available in C ++, MATLAB and Python, which has disturbance torques limitations (Omar 2017); PRO-PAT available as a MATLAB toolbox, also with disturbances limitations (Carrara 2015); SNAP available in Simulink and MATLAB, which has only attitude dynamics and limited disturbance modeling (Rawashdeh 2019); and OPEN-SES-SAME developed in C ++ (Turner 2003).

There exist proprietary software Systems Tool Kit and Orbit Determination Tool Kit simulators from AGI, for either standalone use or integrated with MATLAB;<fnref n="7" /> and Spacecraft Control Toolbox from Princeton Satellite System, available as a

MATLAB package.<fnref n="8" /> Given that the mathematical models of these simulators are not available, it is not possible to compare them to high-fidelity simulators criteria. However, all of them report being high-fidelity mathematical simulators.

Most of the cited simulators have no in-orbit data validation (Menges et al. 1998; Lovera 2003; Turner 2003; Nasirian et al. 2006; Carrara 2015; Triharjanto et al. 2015; Habibkhah et al. 2017; Omar 2017; Rawashdeh 2019). This is mainly due to the difficulty of obtaining real data. Alternatively, some studies compare their own simulation results with other academic (Liu et al. 2019) or commercial (Naqvi & Raza 2007; Alshamy et al. 2019) simulators.

In Sternberg et al. (2018), a team from NASA Jet Propulsion Laboratory uses flight data from the CubeSat MicroMAS-1 to validate their simulator. The same can be seen with CubeSats UWE-3 (Bangert et al. 2014), and ESTCube-1 (Sato et al. 2017). Some works are dedicated to validating the performance of one of their subsystems using in-flight data, which is the case for nanosatellites RAX-1 and RAX-2 (Slavinskis et al. 2016) and ESTCube-1 (Springmann & Cutler 2014). Finally, some other works use in-orbit data for an overall performance analysis (Taraba et al. 2009; Scholz et al. 2010; Xiang et al. 2012; Dechao et al. 2014; Sarda et al. 2014; Mason et al. 2017; Pong 2018).

ADCS simulators validated with in-flight data are altogether rare in the literature, showing the importance of using PicSat in-orbit data for this purpose, as we present in this paper.

### 3. PicSat

PicSat was a 3U CubeSat, developed by the High Angular Resolution Astronomy group of the Laboratory for Space Studies and Instrumentation in Astrophysics at Paris-Meudon Observatory. PicSat was first proposed in 2014, having its Preliminary Design Review in 2015 and its Critical Design Review in 2016 (Nowak et al. 2018), a noticeably short development timescale as expected for CubeSats.

The PicSat mission was intended to show that the diameter of β Pictoris b planet could be estimated with a CubeSat through the planetary transit detection method. This technique would also furnish precious information on the planet's close-by environment (Hill sphere). The PicSat mission would attain those scientific goals at an affordable cost. The project schedule was constrained by the predicted β Pictoris b transit window (between summer 2017 and summer 2018). The typical time frame and development costs of large (>1000 kg) or medium (500–1000 kg) class space-based observatories would most probably make unrealistic a dedicated space mission centered on the Beta Pictoris system (Nowak et al. 2018).

<footnote n="6">⁶ https://opensource.gsfc.nasa.gov/projects/ODTBX/</footnote>

<footnote n="7">⁷ https://www.agi.com/products</footnote>

<footnote n="8">⁸ http://www.psatellite.com/products/sct/</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** PicSat opto-mechanical payload. Credit: Lapeyrère et al. (2017).</figcaption></figure>

## 3.1. Payload and High-accuracy Pointing System

The retained in-house payload final design was based on an innovative technique of space-based interferometry: photometry was acquired by means of an optical system composed of two mirrors and a single-mode optical fiber, and a detector system: a solid-state photodetector Single-Photon Avalanche Diode (SPAD). Figure 2 shows the payload schematically.

The planetary transit detection method required high-precision photometry, which in turn required high-accuracy pointing and SPAD thermal stability. The mission required pointing the satellite toward Beta Pictoris with about 1" rms accuracy in the instrument boresight frame. Such accuracy could not be achieved using only COTS ADCS components. Coping with these requirements necessitated another innovation: the satellite fine pointing control was divided into two stages. The first one was a three-axis stabilization COTS ADCS using star tracker, reaction wheels and gyroscopes at platform level; and the second one an in-house developed active correction of the fiber position using piezoelectric actuators at payload level. The piezoelectric stage was mounted around the fiber to control its position regarding the target star and controlled in closed loop via a microprocessed system, enabling arc second accuracy on star pointing.

An embedded thermal control allowed SPAD noise to operate under 40 ppm hr⁻¹, adequate for the overall photometric error budget (110 ppm hr⁻¹) (Nowak et al. 2017).

## 3.2. ADCS Plus Piezo Stage

ADCS and piezo-controlled fiber formed a two-stage position control system. ADCS controlled satellite angular position up to 30" accuracy, and piezos controlled the linear fiber position through the observed target starlight, allowing 1" final accuracy.

Figure 3 presents a simplified control diagram of the system, indicating the relevant ACS and ADS components. Desired and actual states can be either angular position or angular velocity depending on the platform control mode (de-tumbling or target

pointing mode). The ADCS (model iADCS 100) was developed by Hyperion Technologies. Table 1 lists its actuators and Table 2 its sensors, according to its datasheet.

As iADCS-100 is a proprietary device, we do not have access to detailed information about its algorithms. Nevertheless, by the time of the mission (iADCS-100 software version 2.00), the target pointing mode was probably running a state estimator based on a linear Kalman filter and three Proportional-Integral-Derivative controllers (one for each axis) with anti-Windup (for the Integral term).

The ADCS residual pointing errors (30" accuracy) are the input of the Piezo-Stage Control System, whose objective is to correct them up to 1" accuracy. There is no feed feedback implemented between them. To do so, the fiber is mounted on a two-axis piezo actuator and positioned in the focal plane of the telescope. The piezo actuator is based on space-qualified components and is made by CEDRAT Technologies. Each axis is controlled independently. They are also equipped with strain gauges to measure the piezo excursion. The total range of each piezo is about 430 µm. The fiber position is modulated in the focal plane. After a first search to find the star in the field of view, multiple photometric measurements acquired at different positions along the modulation pattern are combined in an extended Kalman filter to derive a real-time estimate of the star position. The loop is closed using the best estimation of the star position as a reference position to the fiber (Nowak 2019).

## 3.3. The Satellite

Figure 4 presents platform and payload components and their locations inside the 3U envelope. The upper cube unit comprises the power subsystem (EPS), the communication subsystems (TRXVU and HDRM), the VHF antennas and the on-board computer (OBC). The intermediate unit comprises the UHF antennas, the ADCS and the electronic payload, including the photodiode. Optomechanical payloads (Piezo-stage and Telescope) and the star tracker are accommodated in the lower unit. Further details can be found in Nowak et al. (2018).

Payload subsystems were in-house developments and platform components were all COTS.

## 3.4. Mission Profile

PicSat was launched on 2018 January 12, by the PSLV-C40 (Polar Satellite Vehicle Launcher), operated by the Indian Space Research Organization. It was placed at a 505 km Sun-synchronous polar orbit (orbital period of about 95 minutes).

Two major failures prevented the mission from being fully successful: one in the ADCS pointing mode, which could never be activated, probably because of a failure in the star tracker's Lost in Space Mode; and another due to the electronic failure of a critical system by unknown reasons, which caused loss of communication with the CubeSat.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** PicSat's two-stage control: simplified ADCS and piezo stage diagrams.</figcaption></figure>

<table><caption>**Table 1**
iADCS 100's Actuators Components</caption><thead><tr><td>Component</td><td>Model</td><td>Quantity</td></tr></thead><tbody><tr><td>Magnetic actuator</td><td>MTQ 200.20</td><td>2 (for X and Y axes)</td></tr><tr><td>Magnetic actuator</td><td>MTQ 200.10S</td><td>1 (for Z-axis)</td></tr><tr><td>Reaction Wheel</td><td>RW 210.60</td><td>3</td></tr></tbody></table>

<table><caption>**Table 2**
iADCS 100's Sensors Components</caption><thead><tr><td>Component</td><td>Model</td><td>Quantity</td></tr></thead><tbody><tr><td>Star sensor</td><td>ST 200</td><td>1</td></tr><tr><td>3-axis magnetometer</td><td>...</td><td>1</td></tr><tr><td>3-axis gyrometer</td><td>...</td><td>1</td></tr></tbody></table>

After ten weeks of operations, on 2018 April 5 the mission was officially ended. Although unfortunately almost no payload data besides payload temperatures were obtained, sparse but still valuable platform data were acquired Nowak et al. (2018).

## 4. Attitude Dynamics and Orbit Propagation Modeling

PicSat had neither ODS nor OCS. Therefore, this paper deals only with ADCS. Nevertheless, for ADCS development and validation purposes an orbit propagator is necessary. Here we use a ground simulator that computes satellite orbital position based on given initial position and velocity at initial time. Based on the PicSat mission profile and iADCS 100 components, we model spacecraft dynamics, including actuators and sensors, disturbances, and orbit propagation. Modeled components, disturbances, and level of detail of models will depend on each

<figure><img src="image_6.png" /><figcaption>**Figure 4.** PicSat components. Credit: Lapeyrère et al. (2017).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5. Reference frames illustration.</figcaption></figure>

mission. In this paper such models are sufficient for PicSat needs, as will be seen in the validation section (Section 5).

## 4.1. Reference Frames

A crucial step before modeling is to define the set of reference frames necessary for describing attitude and orbit kinematics (positions and velocities) as well as dynamics. We adopt the frame notation presented in Craig (1989), to avoid ambiguity about the reference frame in which a quantity is computed and the one in which it is expressed.

Consider three reference frames {$A$}, {$B$}, and {$C$}, with {$B$} rotating with velocity $\Omega$. The velocity of frame {$B$} with respect to frame {$A$} is denoted $^A\Omega_B$, and $^C(^A\Omega_B)$ denotes the velocity of frame {$B$} with respect to frame {$A$} expressed in frame {$C$}. If {$A$} is an inertial frame, we simplify notation by using lowercase letters. Thus, the velocity of frame {$B$} with respect to the inertial frame {$A$} is $\omega_B$, and $^C\omega_B$ is the velocity of frame {$B$} with respect to the inertial frame {$A$} expressed in frame {$C$}. Figure 5 illustrates the four reference frames necessary for the modeling presented in this section.

1. Earth-Centered Inertial Reference Frame (ECI), denoted by {$I$}. The most important reference frame is the inertial one. All the physical quantities are computed with respect to this frame, and only afterwards expressed with respect to another useful frame. See Appendix A.1.1 for axes definition.

2. Earth-Centered Earth-Fixed Reference Frame (ECEF), denoted by {$E$}, is a non-inertial frame. See Appendix A.1.2 for axes definition.

3. Spacecraft Reference Frame (SCRF), denoted by {$S$}, is useful to represent some elements of the spacecraft, such as actuators, and is fixed with respect to the satellite body. See Appendix A.1.3 for axes definition.

4. Attitude Control Reference Frame (ACRF), denoted by {$A$}, is useful for attitude control. The spacecraft dynamics is simplified if represented in relation to its principal axes of inertia. See Appendix A.1.4 for axes definition.

Appendix A.2 presents the transformations between these reference frames. The next subsections present the key expressions for attitude and orbit dynamics in the reference frames. Derivations and further details are presented in Appendix.

## 4.2. Attitude Dynamics

The overall satellite attitude dynamic model is shown in Equation (1). Its derivation can be found in Appendices A.3–A.6

$$ \begin{aligned} & ^A\mathbf{T}_{\text{mtq}} - ^A\dot{\mathbf{H}}_{\text{rw}} + ^A\mathbf{T}_{\text{dst}} \\ & = ^A\mathbf{I}_{\text{sc}} \, ^A\dot{\omega}_{\text{A}} + ^A\omega_{\text{A}} \times (^A\mathbf{I}_{\text{sc}} \, ^A\omega_{\text{A}} + ^A\mathbf{H}_{\text{rw}}), \end{aligned} \quad (1) $$

where $^A\mathbf{T}_{\text{mtq}}$ is the total magnetic actuators torque, $^A\dot{\mathbf{H}}_{\text{rw}}$ is the controllable reaction wheels torque, $^A\mathbf{T}_{\text{dst}}$ is the torque due to all disturbances, $^A\mathbf{I}_{\text{sc}}$ is the spacecraft inertia matrix, $^A\dot{\omega}_{\text{A}}$ and $^A\omega_{\text{A}}$ are the satellite angular acceleration and velocity, respectively. Here $^A\omega_{\text{A}}$ is the angular velocity of attitude satellite frame {$A$} with respect to inertial frame {$I$} and expressed in terms of frame {$A$}. Finally, $^A\mathbf{H}_{\text{rw}}$ is the angular momentum of the reaction wheels.

### 4.2.1. Attitude Disturbance

Satellites are subject to disturbances that depend on a series of factors including mechanical characteristics of the spacecraft, mission profile and solar activity. In LEO, the main torques are due to the atmospheric (aerodynamic) drag, $T_d$, the solar radiation pressure, $T_{srp}$, the gradient of gravity, $T_{gg}$, and
the residual magnetic dipole, $T_{mag}$. The torque due to all disturbances, $T_{dst}$, is the sum of all these torques expressed in terms of frame {$S$}:

$$ ^{S}\mathbf{T}_{\text{dst}} = ^{S}\mathbf{T}_{\text{d}} + ^{S}\mathbf{T}_{\text{srp}} + ^{S}\mathbf{T}_{\text{gg}} + ^{S}\mathbf{T}_{\text{mag}}. \quad (2) $$

Appendix A.7.1 details these terms.

### 4.3. Orbit Propagation

One way to obtain orbit propagation is through Cowell's Formulation (Vallado 2013), a Keplerian orbit modified model considering a disturbance acceleration $a_{dst}$:

$$ ^{I}\mathbf{a}_{\text{sc}} = -\frac{\mu_{\oplus}}{r_{\text{sc}}^{3}} ^{I}\mathbf{r}_{\text{sc}} + ^{I}\mathbf{a}_{\text{dst}}, \quad (3) $$

where $\mu_{\oplus}$ is the orbital parameter of the Earth equal to 398,600 km$^3$ s$^{-2}$, $r_{sc}$ is the position vector of the spacecraft relative to the Earth and $a_{dst}$ is the acceleration due to all orbital disturbances written in {$I$}.

This modeling choice provides a high-accuracy orbit propagator as will be verified is Section 5.2.

#### 4.3.1. Orbital Disturbances

The most important orbital disturbances in LEO are the atmospheric drag, $a_d$, solar radiation pressure, $a_{srp}$, Earth's oblateness, $a_{ob}$, and third body influence (Sun and Moon), $a_{th}$. The acceleration due to all disturbances, $a_{dst}$, is given by adding the disturbance terms:

$$ ^{I}\mathbf{a}_{\text{dst}} = ^{I}\mathbf{a}_{\text{d}} + ^{I}\mathbf{a}_{\text{srp}} + ^{I}\mathbf{a}_{\text{ob}} + ^{I}\mathbf{a}_{\text{th}}. \quad (4) $$

Appendix A.7.2 details these terms.

### 4.4. Sensors

Sensor models consider measurement of noise and bias. Gyrometer measurements are expressed by

$$ \omega = \omega_{\text{actual}} + \mathbf{b}_{\omega} + \sigma_{\omega}, \quad (5) $$

where $\omega_{actual}$ is the actual angular velocity, $\mathbf{b}_{\omega}$ is the bias and $\sigma_{\omega}$ is a white noise.

Similarly, magnetometer measurements are expressed by

$$ \mathbf{B} = \mathbf{B}_{\text{actual}} + \mathbf{b}_{\text{B}} + \sigma_{\text{B}}, \quad (6) $$

where $\mathbf{B}_{actual}$ is the actual Earth's magnetic field orientation, $\mathbf{b}_{B}$ is the bias and $\sigma_{B}$ is a white noise.

## 5. Attitude Dynamics and Orbit Propagation: Model Identification and Validation Based on In-orbit Data

This section presents the identification and validation of models described in the previous section. As initial values for model coefficients identification, the following values were considered:

1. PicSat overall weight ($m_{sc}$), inertia matrix ($^{S}\mathbf{I}_{sc}$) and center of mass ($^{S}\mathbf{CM}$) with respect to the geometric center

<table><caption>**Table 3**
MTQ 200.20 Magnetic Actuator Data-sheet Values</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Average coil area</td><td>90.0 mm<sup>2</sup></td></tr><tr><td>Maximum magnetic moment</td><td>0.2 A m<sup>2</sup></td></tr><tr><td>Height</td><td>80.0 mm</td></tr><tr><td>Weight</td><td>39.6 mg</td></tr></tbody></table>

Note. Source:Hyperion Technologies.

<table><caption>**Table 4**
MTQ 200.10S Magnetic Actuator Data-sheet Values</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Average coil area</td><td>283.5 mm<sup>2</sup></td></tr><tr><td>Maximum magnetic moment</td><td>0.1 A m<sup>2</sup></td></tr><tr><td>Height</td><td>25.0 mm</td></tr><tr><td>Weight</td><td>44.8 mg</td></tr></tbody></table>

Note. Source:Hyperion Technologies.

(satellite reference frame {$S$}) were extracted from its CAD model.<fnref n="9" />

2. Parameter values for the magnetic actuators are shown in Tables 3 and 4, and for the reaction wheels in Table 5.

In-orbit data for model identification and validation are presented in Figure 6, showing PicSat in-orbit angular velocity measured by its gyrometer. All mission data available, ranging from 2018 January 22 to 2018 March 20, are plotted. Note that the z-axis angular velocity presents offset at several moments, including during the Beginning of Life (BoL) and during on-ground pre-launch data acquirement for calibration purposes. Therefore, pre-launch data was used by PicSat mission control to calibrate the gyrometer, and Z-axis offset was estimated by averaging the points prior to the launch. The value obtained was 0.1694 rad s$^{-1}$. Gyrometer bias is not completely corrected by this process. The actual value of angular velocity remains unknown and Equation (5) applies.

### 5.1. Attitude Dynamics Identification and Validation

Since PicSat had no sensors to measure torques provoked by external disturbances, attitude model identification must be performed in four steps.

**Step 1:** identification of the the open-loop satellite dynamics considering external disturbances as negligible.

<footnote n="9">⁹ The CAD model version made available was prior to final assembly, without cables and connections. In addition, those values were omitted here for confidentiality reasons.</footnote>
<table><caption>**Table 5**
RW 210.60 Reaction Wheel Data-sheet Values</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Maximum torque</td><td>0.1 × 10<sup>-3</sup> N m</td></tr><tr><td>Maximum angular momentum</td><td>6.0 × 10<sup>-3</sup> N m s</td></tr><tr><td>Weight</td><td>48.0 mg</td></tr></tbody></table>

Note. Source:Hyperion Technologies.

Step 2: validation of model obtained in Step 1.

Step 3: identification of the open-loop external disturbances, given the validated satellite dynamics.

Step 4: validation of model obtained in Step 3.

Figure 7 presents these steps in a flow diagram. One can refer to this figure during steps explanation.

We use the satellite dynamic model presented in Section 4.2 to determine the time window regarding PicSat open-loop data for which we can consider the external disturbances negligible. Isolating the derivative of the spacecraft angular velocity in Equation (1) gives

$$ ^{A}\dot{\omega}_{A} = (^{A}\mathbf{I}_{\text{sc}})^{-1} [^{A}\mathbf{T}_{\text{mtq}} - ^{A}\dot{\mathbf{H}}_{\text{rw}} + ^{A}\mathbf{T}_{\text{dst}} - ^{A}\omega_{A} \times (^{A}\mathbf{I}_{\text{sc}}^{A}\omega_{A} + ^{A}\mathbf{H}_{\text{rw}})]. \quad (7) $$

Equation (7) is the dynamic simulation form of Equation (1).

Considering the case in which magnetorquers are off ($T_{\text{mtq}} = 0$) and reaction wheels stopped ($H_{\text{rw}} = 0$), Equation (7) becomes

$$ ^{A}\dot{\omega}_{A} = (^{A}\mathbf{I}_{\text{sc}})^{-1} [^{A}\mathbf{T}_{\text{dst}} - ^{A}\omega_{A} \times (^{A}\mathbf{I}_{\text{sc}}^{A}\omega_{A})]. \quad (8) $$

Omitting the sub and superscribed reference frames and taking the integration of Equation (8) between times $t_i$ and $t_f$, we obtain

$$ \omega(t_f) = \omega(t_i) - \mathbf{I}_{\text{sc}}^{-1} \int_{t_i}^{t_f} \omega \times (\mathbf{I}_{\text{sc}}\omega) dt + \mathbf{I}_{\text{sc}}^{-1} \int_{t_i}^{t_f} \mathbf{T}_{\text{dst}} dt. \quad (9) $$

The first term in Equation (9) is the angular velocity of free motion, and depends only on the satellite initial velocity and its own gyroscopic torque:

$$ \omega_{\text{free}}(t_f) = \omega(t_i) - \mathbf{I}_{\text{sc}}^{-1} \int_{t_i}^{t_f} \omega \times (\mathbf{I}_{\text{sc}}\omega) dt. \quad (10) $$

In Equation (10), if $||\omega(t_i)|| \gg ||\mathbf{I}_{\text{sc}}^{-1} \int_{t_i}^{t_f} \omega \times (\mathbf{I}_{\text{sc}}\omega) dt||$ then $||\omega_{\text{free}}(t_f)|| \approx ||\omega(t_i)||$. The second term in Equation (9) is the angular velocity due to external disturbances:

$$ \omega_{\text{dst}}(t_f) = \mathbf{I}_{\text{sc}}^{-1} \int_{t_i}^{t_f} \mathbf{T}_{\text{dst}} dt. \quad (11) $$

The worst case for Equation (11) occurs when $T_{\text{dst}}$ is maximum and constant, that is,

$$ \omega_{\text{dstmax}}(t_f) = \mathbf{I}_{\text{sc}}^{-1} \mathbf{T}_{\text{dstmax}}(t_f - t_i). $$

<figure><img src="image_8.png" /><figcaption>**Figure 6.** PicSat pre-launch (before 2018 January 22) and in-orbit angular velocity (from 2018 January 22 to 2018 March 20).</figcaption></figure>

Thus, the ratio between the free motion and the maximum disturbance in terms of norm of angular velocity is

$$ \frac{\|\omega_{\text{free}}(t_f)\|}{\|\omega_{\text{dstmax}}(t_f)\|} \approx \frac{\|\omega(t_i)\|}{\|\mathbf{I}_{\text{sc}}^{-1}\mathbf{T}_{\text{dstmax}}\|(t_f - t_i)} $$

In attitude modeling for closed-loop control, we may neglect disturbances of less than one tenth of the open-loop response amplitude, that is, $||\omega_{\text{free}}||/||\omega_{\text{dst}}|| \ge 10$. Therefore

$$ (t_f - t_i) \le \frac{\|\omega(t_i)\|}{10 \|\mathbf{I}_{\text{sc}}^{-1} \mathbf{T}_{\text{dstmax}}\|} $$

Considering the suitable in-orbit data subsets for such evaluation, the lowest norm of the angular velocity was $\omega(t_i) = 0.028 \text{ rad s}^{-1}$.

Taking $\mathbf{I}_{\text{sc}}^{-1}$ from the PicSat CAD model as the best initial approximation of the true $\mathbf{I}_{\text{sc}}^{-1}$, and $T_{\text{dstmax}}$ estimated via the attitude dynamic model using initial values for parameters taken from satellite components data sheets, we find $||\mathbf{I}_{\text{sc}}^{-1} \mathbf{T}_{\text{dstmax}}|| = 4.15 \times 10^{-5} \text{ s}^{-2}$. Thus

$$ \Delta t \le \frac{0.028}{10.415 \times 10^{-5}} = 67.5 \text{ s.} $$

Therefore Step 1 must be applied to time windows up to 67.5 s. Since $T_{\text{dstmax}}$ is a very conservative estimation and a better value can be determined by iterating modeling results against in-orbit data, we will relax this constraint to 100 s allowing the use of more in-orbit data and improving model accuracy.

In order to measure $\omega(t_i)$ and $T_{\text{dstmax}}$ in open-loop attitude response with magnetorquers off and reaction wheels stopped, analysis of in-orbit angular velocity data allowed us to identify 10 data stretches with at least 10 consecutive samples having the following lengths: 290, 680, 570, 300, 770, 600, 690, 750, 720 and 620 s. To carry out Steps 1 and 2, these stretches were divided into subsets of 100 s. When necessary to complete the last subset, data from the previous one was used. There were 63 subsets altogether.
<figure><img src="image_9.png" /><figcaption>Figure 7. Attitude model identification and validation steps diagram.</figcaption></figure>

### 5.1.1. Step 1—Satellite Dynamics Identification

For Step 1, 32 from the original 63 subsets were selected. An arbitrary rule was applied: the odd index subsets were chosen, that is, the first, the third, the fifth, and so on. PicSat inertia matrix, ADCS gyrometers' bias and initial velocity were optimized/estimated based on a nonlinear least squares approach using those subsets, that is, they were all free parameters of the optimization algorithm. Principal inertia moments varied from about 1% in X and Y axes to 10% in Z-axis comparing CAD guess and identified values. Furthermore, the gyrometer biases were in the order of $10^{-3}$ rad s$^{-1}$ and their respective values were subtracted from each axis.

In order to quantify the difference between in-orbit data and the identified attitude model, we employ the metric

$$ \int_{t_i}^{t_f} \mathbf{diff} \approx \sum_{k=1}^{N} \| \omega_{\text{in-orbit}}(k) - \omega_{\text{model}}(k) \|_2, \quad (12) $$

where diff is a vector containing the $ℓ^2$-norm of the angular velocity difference between in-orbit and modeled satellite angular velocity, N is the number of measurements acquired between $t_i$ and $t_f$. As Equation (12) is a summation of modeling residues, its value increases with time for each velocity axis. More classical metrics would lead us to deal with in-orbit data

noise filtering, which is an extra work without any big advantage in terms of modeling accuracy regarding closed-loop control.

Figure 8 shows the identification results for the angular velocity in terms of Equation (12). In this plot each dot is the result of summation of 100 s of $ℓ^2$-norm data.

The average differences are distributed around the abscissa axis, indicating residual randomness. Summation of norm is typically around $7 \times 10^{-3}$ rad s$^{-1}$, which means that identification residuals are very small. We show them in more details in Step 2. Only for Z-axis, norm summation is as high as $2 \times 10^{-2}$ rad s$^{-1}$, which is due to the gyrometer noise.

### 5.1.2. Step 2—Satellite Dynamics Validation

The validation was performed with the other 31 subsets. In this case, even indexes were used, that is, the second, the fourth, the sixth and so on. An example of result for one subset is shown in Figure 9.

One can see an excellent match between in-orbit and modeled angular velocity, except for Z-axis, which is not due to the dynamic model. Instead, it is due to high gyrometer noise. In this time frame (100 s), the satellite behaves as a rigid body in free motion. Thus, we have kinetic energy being
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Step 1 results—summation of the $ℓ^2$-norm of the difference between in-orbit and modeled satellite angular velocity for all subsets.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(a) Model validation - in-orbit (real) and modeled (sim) satellite angular velocity</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Model validation - $ℓ^2$-norm of the difference between in-orbit and modeled satellite angular velocity</figcaption></figure>

<caption>**Figure 9.** Step 2 results for one subset.</caption>

exchanged among its three rotation axes, which means that angular velocities are coupled and are constrained by the rigid body dynamics. Thus, it is physically impossible to have Z-axis high-frequency velocity while X and Y axes move slowly. One can observe this behavior all around the subsets in Figure 8 as well. Therefore, for modeling and validation purposes, X and Y axes raw data are more reliable them Z-axis data. If we apply a low-pass filter matching the frequency domain of the satellite dynamics, they will be more consistently reliable. We have not designed and implemented such a filter, because a proper development requires detailed

knowledge about the flight sensors while the raw data accuracy is already enough for dynamic modeling in terms of ADCS.

Figure 10 shows the validation results for the angular velocity in terms of Equation (12). In this plot each dot is the result of summation of 100 s of $ℓ^2$-norm data. In other words, Figure 9 becomes three dots in Figure 10, one for each axis.

As for the identification results, the average differences appear to be random. This time, Figure 9 helps to understand what the amplitudes mean. The summation of the norm for the X-axis is $0.5 \times 10^{-2} \text{ rad s}^{-1}$, for the Y-axis is $1 \times 10^{-2} \text{ rad s}^{-1}$ and for the
<figure><img src="image_13.png" /><figcaption>Figure 10. Step 2 results—summation of the ℓ²-norm of the difference between in-orbit and modeled satellite angular velocity for all subsets.</figcaption></figure>

Z-axis is $2 \times 10^{-2}$ rad s⁻¹. Therefore, $1 \times 10^{-2}$ rad s⁻¹ or less indicates a really good match between real data and model.

### 5.1.3. Step 3—External Disturbances Identification

To accommodate uncertainties on the disturbance model parameters, and even data time stamping bias, a disturbance gain, $K_{dst}$, was introduced in Equation (2), leading to:

$$ ^{S}\mathbf{T}_{dst} = ^{S}K_{dst}(^{S}\mathbf{T}_{d} + ^{S}\mathbf{T}_{srp} + ^{S}\mathbf{T}_{gg} + ^{S}\mathbf{m}_{mag} \times ^{S}\mathbf{B}), \quad (13) $$

where the residual magnetic dipole ($\mathbf{m}_{mag}$) will be also a degree of freedom, as long as it is the most important disturbance in the specific case of PicSat.

As in Step 1, to perform Step 3 we selected 5 subsets alternately from the original 10 ones. The subsets are now longer than 100 s. The validated inertia matrix and gyro bias are obtained from Step 2. As previously, the procedure to identify $K_{dst}$ and $\mathbf{m}_{mag}$ utilized a nonlinear least squares approach. The values obtained were

$$ \begin{gathered} ^{S}\mathbf{m}_{mag} = [1.86 \pm 0.04 \quad 8.47 \pm 0.08 \quad 33.02 \pm 0.65]^{T} \times 10^{-3} \text{A m}^{2} \\ \text{and} \\ ^{S}\mathbf{K}_{dst} = [1.08 \pm 0.01 \quad 0.95 \pm 0.01 \quad 1.06 \pm 0.03]^{T}. \end{gathered} $$

Figure 11 shows the identification results for the angular velocity.

Step 3 utilizes a small number of subsets when compared to Steps 1 and 2. Another difference is that this time the duration is as long as 700 s. Even so, the summation of norm gives results in the order of $2 \times 4 \times 10^{-2}$ rad s⁻¹, which is remarkably low. Therefore, disturbances have been adequately modeled.

Feedback control can deal with disturbances which are unknown or known with some degree of accuracy. Step 3 results are sufficient for designing a compensator for PicSat.

### 5.1.4. Step 4—External Disturbances Validation

External disturbances validation was performed using the other 5 subsets. An example of result for one subset is shown in Figure 12.

<figure><img src="image_14.png" /><figcaption>Figure 11. Step 3 results—summation of the ℓ²-norm of the difference between in-orbit and modeled satellite angular velocity for all subsets.</figcaption></figure>

There is an excellent match between in-orbit and modeled angular velocity even in the presence of external disturbances. The spectral density of differences is not concentrated in any specific frequency, which indicates that no systematic modeling error is present.

Figure 13 presents a zoom-in in the modeled external disturbances. The torque due to the residual magnetic dipole, $\mathbf{T}_{mag}$, is the most important for PicSat. It is partially a consequence of letting $\mathbf{m}_{mag}$ as a free parameter and partially because, indeed, the residual magnetic dipole is typically the most relevant external disturbance torque for CubeSats. Nevertheless, any other choice would have produced a distortion in terms of disturbance weight and would not be appropriate. Figure 13 also shows that the torque due to the solar radiation pressure, $\mathbf{T}_{srp}$, is the less important.

In Section 5.1, using theoretical and nominal values, we have determined that for time frames smaller than 100 s, the influence of the external disturbances would be ten times less than the free motion dynamics. To verify this premise, Figure 14 presents simulations comparing the presence and absence of disturbances in the satellite's dynamic model against in-orbit data. One can observe a clear increase in the modeling error after 150 s if disturbance is not taking into account in the model.

More precisely, when simulation reaches 200 s, the ratio between external torque and free motion reaches about 1/10.
<figure><img src="image_15.png" /><figcaption>(a) Attitude disturbance validation - in-orbit (real) and modeled (sim) satellite angular velocity</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Attitude disturbance validation - $ℓ^2$-norm of the difference between in-orbit and modeled satellite angular velocity</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(c) Attitude disturbance validation - power spectrum for the difference between in-orbit and modeled satellite angular velocity</figcaption></figure>

**Figure 12.** Attitude disturbance model validation results for one subset.

<figure><img src="image_18.png" /><figcaption>**Figure 13.** Validated attitude disturbances (absolute torques) for Subset 1.</figcaption></figure>

Figure 15 shows the ratio evolution regarding the time window. Such result is computed based on the validated model.

The time window of 100 s adopted is consistent with negligible disturbance. The whole set of validation results for angular velocity is presented in Figure 16. As in the estimation

<figure><img src="image_19.png" /><figcaption>**Figure 14.** Comparison between presence and absence of disturbances over time —Subset 1. Real data and simulation with (w/) and without (s/) disturbances.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 15.** Disturbance to free motion ratio regarding time window.</figcaption></figure>

results, the summation stays at very low values indicating the modeling accuracy achieved suffices for feedback control.

## 5.2. Orbit Propagation Identification and Validation

The orbital position of a satellite can be obtained by means of a GPS device on board the spacecraft, or by tracking radars on the ground. PicSat lacked either system. However, the U.S. Space Surveillance Network tracks all objects in Earth orbit (artificial satellites and debris) and maintains a cataloged database of them (Vallado 2013). Their data are published by the Joint Space Operations Center on the portal Space-Track.org<fnref n="10" /> in a format called Two-Line Element (TLE).

To estimate its current position, the ODS has an embedded orbit propagator,<fnref n="11" /> which receives the last available TLE from the ground station.

For a given object, the time between TLE updates is not regular. It can vary from a fraction of the orbit to several days. Figure 17 shows the histogram of the update time between the PicSat TLEs from 2018 January 12 through 2019 September 29.

Altogether there were 1166 updates during this time span. The longest interval between TLEs was 30 orbits, equivalent to 47.4 hr and 75% of updates occurred in less than a 11 hr interval.

<footnote n="10">¹⁰ See https://www.space-track.org/.</footnote>

<footnote n="11">¹¹ In such cases, the propagator is often the Simplified General Perturbations 4 (SGP4) (Vallado et al. 2006), as it is compatible with TLE data (Vallado 2013). However, as iADCS-100 is a proprietary device, we cannot state which propagator was indeed embedded.</footnote>
<figure><img src="image_21.png" /><figcaption>**Figure 16.** Step 4 results—summation of the $ℓ^2$-norm of the difference between in-orbit and modeled satellite angular velocity for all subsets.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 17.** Histogram of PicSat's TLEs updates available from 2018 January 12 through 2019 September 29.</figcaption></figure>

The published TLEs lack accuracy information. However, in Riesing (2015) the Planet CubeSats constellation Flock 1B was used to determine the accuracy of the TLEs of these spacecraft, once they were equipped with a GPS module. As with PicSat, they were 3Us CubeSats, launched from the ISS. Ten satellites were deployed during 2014 September and 634 TLE updates occurred during that time. The statistical error in position could be estimated with good statistical accuracy to be less than 2.01 km in 25% of cases (first quartile), less than 4.52 km in 50% of cases (median), and less than 10.60 km in 75% of cases (third quartile).

Although contact with PicSat was lost on 2018 March 20, TLEs will continue to be available until the satellite re-enters Earth's atmosphere. We performed the identification and validation of the simulator's orbital propagation using the available TLEs for PicSat and the statistical error obtained in Riesing (2015). For this purpose, the following procedure was carried out: using TLEs as initial conditions for the simulator, we propagate the states until the next available TLE, comparing them at the end. We repeat this process ranging from TLEs available between one orbit up to 16 orbits (i.e., about 1 day between two TLEs). The results of the propagation error are shown in Figure 18. The horizontal dashed lines represent the statistical errors obtained in Riesing (2015).

<figure><img src="image_23.png" /><figcaption>**Figure 18.** Evolution of the orbital propagation difference between TLEs and model according to the number of propagated orbits.</figcaption></figure>

Model error increases with the propagation time. After 15 orbits, the median error of the model is 7.1 km, which is a very small error considering that 648,047.7 km were covered, i.e., $1.1 \times 10^{-7}\%$ propagation error in 15 orbits. For an observer on the ground, it is equivalent to 3.5 uncertainty in pointing toward the spacecraft. For payloads that do not demand very high accuracy orbit determination, such propagator accuracy is sufficient.

Based on this model, the predicted re-enter date for PicSat is 2022 December 24 03:44:34 (UTC).

## 6. ADCS De-tumbling Mode

PicSat was placed in orbit by a QuadPack CubeSat Deployer. Such deployers eject CubeSats by means of compressed springs that generate an initial angular velocity on the satellites. PicSat in-orbit data shows an initial velocity of $1.78 \text{ deg s}^{-1}$.

After launch, the ADCS should be able to reduce the spacecraft angular velocity down to $0.2 \text{ deg s}^{-1}$ (standard limit for proper working of the star tracker) and then pointing the target star with a maximum error of $\pm 30'' \text{ rms}$.

The iADCS-100 system has five operation modes: Idle, Safe, Measurement, De-tumbling and Target Pointing. Only the last two modes are closed-loop controllers. The De-tumbling Mode stabilizes the satellite angular velocity, reducing it to the requirement and the Target Pointing Mode points the payload toward the target star within the required accuracy.

The De-tumbling Mode uses only magnetorquers as actuators, and the Target Pointing Mode uses reaction wheels for pointing the satellite and magnetorquers for desaturating reaction wheels. By reducing the angular velocity, the De-tumbling Mode allows image acquisition from the star tracker's Lost in Space Mode, which provides the line of sight, preparing the ADCS to switch to Target Pointing Mode.

The De-tumbling Mode uses the simplest closed-loop controller. The exact structure of the proprietary de-tumbling controller used by PicSat's ADCS is unknown. Nevertheless, the Interface Control Document furnished with iADCS-100 clearly states that a B-dot controller is used for de-tumbling
<figure><img src="image_24.png" /><figcaption>Figure 19. PicSat in-orbit angular velocity norm.</figcaption></figure>

with two possible working configurations: proportional gain and maximum gain. First presented by Stickler & Alfriend (1976), the B-dot technique follows the principle of gradually reducing the spacecraft's rotational kinetic energy with the aid of the Earth's magnetic field. The control law polarizes the magnetorquers in the opposite direction of the local magnetic field derivative, $\dot{\mathbf{B}}$, from where the technique's name is derived.

The control law is expressed by:

$$^{A}\mathbf{m}_{\mathrm{m}}=-\mathbf{k}_{\mathrm{det}}^{A} \dot{\mathbf{B}}, \quad (14)$$

where $^{A}\mathbf{m}_{\mathrm{m}}$ is the dipole moment and $\mathbf{k}_{\mathrm{det}}$ is the controller gain (diagonal matrix), with $\mathbf{k}_{\mathrm{det}} > 0$.

## 6.1. De-tumbling Mode: Model Identification and In-orbit Data Validation

Figure 19 shows the PicSat angular velocity norm during the entire mission. Start and stop commands for de-tumbling mode are signaled with dashed and dotted vertical lines, respectively. Three of the occasions, where an effective reduction in angular velocity can be seen, are marked by numbers 1–3.

Considering the B-dot control structure, $\mathbf{k}_{\mathrm{det}}$ gain and initial velocity were estimated in a way similar to the dynamics validation performed in Section 5.1, using Equation (12). Free parameters of the algorithm were the initial velocity and the controller gain. Once again, the procedure was performed using a nonlinear least squares method. Data stretch number 2 was used to do the estimation and data stretches number 1 and 3 were used for the validation. B-dot gain identification and validation results are shown in Figure 20.

The identified controller (Equation (14)) presents a satisfactory result for the norm of angular velocities. This indicates that the kinetic energy reduction is achieved in a similar way for both flight and on-ground modeled controllers. However, if the in-orbit and modeled angular velocities are compared independently, axis by axis, the results are unsatisfactory, suggesting some inconsistency between them. Such an issue will be investigated and an advanced identification strategy as discussed in Romano & Pait (2017) will be applied to precisely determine what is missing in Equation (14).

<figure><img src="image_25.png" /><figcaption>(a) De-tumbling model identification result - in-orbit (real) and modeled (sim) satellite angular velocity norm</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(b) De-tumbling model validation result - in-orbit (real) and modeled (sim) satellite angular velocity norm</figcaption></figure>

<caption>Figure 20. De-tumbling control loop identification and validation results.</caption>

## 7. Conclusions

In this paper, we have presented and validated a dynamic attitude and orbit propagation simulator for PicSat as well as its de-tumbling mode. High-precision dynamic models have been obtained for both attitude and orbit. Norm of the difference between the attitude model and in-orbit angular velocity has been determined to be typically lower than 0.04 rad s⁻¹ over 700 s, while the orbit propagation is as accurate as ±500 m rms in one orbit. We have also properly simulated the de-tumbling angular velocity norm time series compared to in-orbit data.

Despite not having accomplished its astronomical goal, PicSat left us an enduring legacy: its platform data allowed us to perform a flight modeling that will be quite valuable for future missions. Such models are necessary for missions aiming to achieve high-performance pointing requirements (of the order of the arcsecond rms), as is the case of PicSat. They are also broadly suited for different mission phases, from spacecraft design to data exploitation.

Based on these models, many studies can be made for PicSat and any other CubeSat mission. Concerning PicSat, one can use the orbital propagator to apply orbital position correction to the available data as well as time stamping correction. The attitude dynamic model allows improving the two-stage control strategy, studying performance, observability, controllability, and system stability regarding nominal and degraded operation. It can be used to refine the control mode switch strategy and simulate issues as faced by PicSat to modify and improve on-board algorithms. For future missions, it can be also used for
implementation of the ADCS embedded software in a re-programmable Real-Time Operation System, allowing for re-programmable in-flight control algorithms, a highly desirable improvement compared to PicSat's ADCS. It is essential to minimize the possibility of failure for both the platform and the payload, especially in satellites such as PicSat, whose pointing depends on both.

There is still a lot to learn from PicSat flight data. This work made it possible to prepare the bases for a more detailed understanding of its months of operation, which will allow a remarkable increase of reliability for upcoming CubeSat missions, both in the preparatory and exploitation mission phases. In this sense, we are already working in a C++ version of the models and control loop presented here to produce On-board In-Loop (OIL) and Hardware In the Loop (HIL) simulators. We are also going to study in more detail the tumbling mode.

The models obtained are a novel result and can be applied to any CubeSat in low Earth orbit, requiring only minor adaptation to each spacecraft. The decision whether to develop models to support ADCS simulators is linked to the mission criticality. They are indispensable in scientific missions, therefore scientific space research laboratories worldwide need to employ CubeSat dynamic modeling and control to succeed in the most challenging endeavors.

The PicSat satellite was designed, integrated and operated by Observatoire de Paris—PSL/CNRS. We would like to thank the students, engineers and researchers who made this project happen. The core team was Sylvestre Lacour, Vincent Lapeyrere, Mathias Nowak, Lester David, Antoine Crouzier, Guillaume Schworer and Maarten Roos.

The PicSat team were supported by the European Union under ERC grant 639248 LITHIUM.

Part of the work presented in this paper was supported by Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) under grants 2016/13750-6 and 2018/19585-2.

# Appendix
## Detailed Spacecraft and Disturbance Modeling

This appendix details the spacecraft and disturbance modeling.

### A.1. Reference Frames

We start by clearly defining all necessary reference frames.

#### A.1.1. Earth-centered Inertial Reference Frame (ECI)

The ECI is denoted by {$I$}. It has its origin in the Earth's center of mass. $\hat{X}_I$ the X-axis unit vector is defined as the direction of the equinox of 2000 March, and represented by the

Greek letter $\gamma$, $\hat{Z}_I$ is taken at the north direction of the Earth's rotation axis and $\hat{Y}_I$ completes the right-handed Cartesian coordinate system.

#### A.1.2. Earth-centered Earth-fixed Reference Frame (ECEF)

The ECEF is a non-inertial frame denoted by {$E$}. It has its origin in the Earth's center of mass. $\hat{X}_E$ is defined as the direction of prime meridian, $\hat{Z}_E$ is taken at the north direction of the Earth's rotation axis and $\hat{Y}_E$ axis completes the right-handed Cartesian coordinate system and the ($\hat{X}_E, \hat{Y}_E$)-plane forms the equatorial plane.

#### A.1.3. Spacecraft Reference Frame (SCRF)

The SCRF, denoted by {$S$}, has its origin $O$ in the geometric center of the spacecraft, $\hat{Z}_S$ axis is normal to the side that contains the antennas, $\hat{X}_S$ axis is perpendicular to the satellite side and $\hat{Y}_S$ axis completes the right-handed coordinate system.

#### A.1.4. Attitude Control Reference Frame (ACRF)

The ACRF, denoted by {$A$}, has its origin is the spacecraft's center of mass $CM$, $\hat{Z}_A$ axis being in the direction of the biggest moment of inertia, $\hat{X}_A$ axis in the direction of the smaller moment of inertia and $\hat{Y}_A$ axis completes the right-handed coordinate system.

### A.2. Transformation between Reference Frames

Reference frames can be related to each other by frame transformations.

#### A.2.1. Transformation between ECI and ECEF

Once the ECI and ECEF ($\hat{X}\hat{Y}$)-planes are coplanar, the rotation between them comes only from the Earth's rotation movement around the $\hat{Z}$ axis. Thus, the rotation matrix and its representation in quaternions are expressed by

$$ {}^I\mathcal{Q}_E = \begin{bmatrix} \cos \theta_G & \sin \theta_G & 0 \\ -\sin \theta_G & \cos \theta_G & 0 \\ 0 & 0 & 1 \end{bmatrix} \quad (A1) $$

and

$$ {}^I q_E = \begin{bmatrix} \cos\left(\frac{\theta_G}{2}\right) & 0 & 0 & \sin\left(\frac{\theta_G}{2}\right) \end{bmatrix}^T, \quad (A2) $$

where $\theta_G$ is the Greenwich Sidereal Time, that is, the angle between the March equinox and the main meridian in degrees. As shown in Vallado (2013),

$$ \theta_G = \theta_{G0} + 360.98564724 \frac{t_U}{24}, \quad (A3) $$
where $t_U$ is the time in hours in UT1 scale, and $\theta_{G0}$ is the Greenwich Sidereal Time at 0 hr UT1. As presented in Seidelmann (2006), $\theta_{G0}$ is given by

$$ \theta_{G0} = 100.4606184 + 36000.77004 T_0 + 0,000387933 T_0^2 - 2.583 \times 10^{-8} T_0^3, \quad (A4) $$

where $T_0 = d_U/36525$, and $d_U$ is the time in Julian days since J2000.0.

### A.2.2. Transformation between SCRF and ACRF

The transformation of the spacecraft inertia matrix $I_{sc}$, from SCRF or just $\{S\}$ to ACRF or just $\{A\}$, is given by

$$ I = (AQ_S)(SI)(AQ_S)^T, \quad (A5) $$

where $^AI$ is the diagonal form of the spacecraft's inertia matrix, $^SI$ is the non-diagonal form of the spacecraft's inertia matrix and $AQ_S$ is the rotation matrix from $\{S\}$ to $\{A\}$.

As shown in Curtis (2014), the matrix $^AI$ is formed, on the main diagonal, by the eigenvalues of $^SI$ and the matrix $AQ_S$ is formed by the respective eigenvectors of $^SI$.

### A.3. Attitude Representation

We can represent the satellite orientation by means of rotation matrices (also called direction cosine matrices), Euler angles, or quaternions. Rotation matrices are hard to read and computationally expensive. Euler angles are easy to read but present singularities and require the use of rotation matrices when transforming reference frames. Quaternions are hard to read and computationally inexpensive. Therefore, from now on we present results using Euler angles, and compute transformations using quaternions, which are preferable for embedded ACS (Sidi 1997; Curtis 2014).

The Euler angles are three angles describing the orientation of a rigid body with respect to a fixed coordinate system. Here they are expressed by $(\phi, \theta, \psi)$ corresponding to yaw, pitch and roll angles in frame $\{A\}$.

A quaternion can be expressed as:

$$ \begin{aligned} q &= q_0 + q_1 \hat{i} + q_2 \hat{j} + q_3 \hat{k} \\ &= q_0 + \mathbf{q}; \end{aligned} \quad (A6) $$

where $q_0$ is the scalar part and $\mathbf{q}$ is the vector or imaginary part (Wertz 1978).

By restricting the quaternion to be unitary, $||q||=$ $\sqrt{q_0^2 + q_1^2 + q_2^2 + q_3^2} = 1$, it is possible to relate it to the symmetrical Euler parameters.

The transformation between them using the 3-2-1 sequence is given by:

$$ I_{q_A} = \begin{bmatrix} \cos(\frac{\psi}{2})\cos(\frac{\theta}{2})\cos(\frac{\phi}{2}) + \sin(\frac{\psi}{2})\sin(\frac{\theta}{2})\sin(\frac{\phi}{2}) \\ \sin(\frac{\psi}{2})\cos(\frac{\theta}{2})\cos(\frac{\phi}{2}) - \cos(\frac{\psi}{2})\sin(\frac{\theta}{2})\sin(\frac{\phi}{2}) \\ \cos(\frac{\psi}{2})\sin(\frac{\theta}{2})\cos(\frac{\phi}{2}) + \sin(\frac{\psi}{2})\cos(\frac{\theta}{2})\sin(\frac{\phi}{2}) \\ \cos(\frac{\psi}{2})\cos(\frac{\theta}{2})\sin(\frac{\phi}{2}) - \sin(\frac{\psi}{2})\sin(\frac{\theta}{2})\cos(\frac{\phi}{2}) \end{bmatrix}. \quad (A7) $$

### A.4. Satellite Kinematics

Concerning the satellite, kinematics describes its attitude and angular velocity time series without any concern about the torques that have created such movement.

Consider the spacecraft's Cartesian angular velocity described by

$$ \omega_A = [\omega_x \quad \omega_y \quad \omega_z]^T_A, \quad (A8) $$

where $\omega_x$, $\omega_y$ and $\omega_z$ are, respectively, the angular velocities around x, y, and z axes of frame $\{A\}$ w.r.t. frame $\{I\}$. As shown in Sidi (1997), the relationship between the time derivative of quaternions and the Cartesian angular velocity is given by:

$$ \frac{d}{dt} I_{q_A} = \frac{1}{2} W(\omega_A)^T I_{q_A}, \quad (A9) $$

where $W(\omega_A)$ is a skew-symmetric matrix given by

$$ W(\omega_A) = \begin{bmatrix} 0 & \omega_z & -\omega_y & \omega_x \\ -\omega_z & 0 & \omega_x & \omega_y \\ \omega_y & -\omega_x & 0 & \omega_z \\ -\omega_x & -\omega_y & -\omega_z & 0 \end{bmatrix}_A. \quad (A10) $$

As $\omega_x$, $\omega_y$ and $\omega_z$ are measured with respect to the inertial frame $\{I\}$, the satellite attitude is simply their integral. Knowing the satellite's angular velocity, it is straightforward to compute the equivalent quaternion derivative and the quaternion itself, which is useful for simulation purposes.

### A.5. Satellite Dynamics

Attitude dynamics describes the angular velocity as a function of the applied torque $T$.

Most CubeSats can be considered rigid bodies, because they do not have moving parts and orbit control systems, which usually employ of propellants. Therefore, mass distribution remains constant over time. Although PicSat has a piezoelectric positioning system as part of the payload (Piezo-stage in Figure 4), it can be considered a rigid body for attitude dynamics modeling. As a result, rigid body dynamics properties and equations apply. Considering $\omega_A$ the angular velocity
of {ACRF} or just {A} with respect to {ECI} or just {I}, the Euler equation is given by Curtis (2014); Sidi (1997):

$$ ^{A}\mathbf{T}_{\text{act}} + ^{A}\mathbf{T}_{\text{dst}} = ^{A}\mathbf{I}_{\text{sc}}{}^{A}\dot{\omega}_{\text{A}} + ^{A}\omega_{\text{A}} \times (^{A}\mathbf{I}_{\text{sc}}{}^{A}\omega_{\text{A}}) \quad (\text{A11}) $$

where $^{A}\mathbf{I}_{\text{sc}}$ and $^{A}\omega_{\text{A}}$ are the spacecraft inertia matrix and the angular velocity written in {A}, respectively, $^{A}\mathbf{T}_{\text{act}}$ is the torque provoked by actuators and $^{A}\mathbf{T}_{\text{dst}}$ is the torque provoked by disturbances.

### A.6. Actuators

According to Table 1, PicSat is equipped with magnetic actuators and reaction wheels. Therefore, the torque due to actuators is given by Equation (A12)

$$ ^{A}\mathbf{T}_{\text{act}} = ^{A}\mathbf{T}_{\text{mtq}} + ^{A}\mathbf{T}_{\text{rw}} \quad (\text{A12}) $$

where $^{A}\mathbf{T}_{\text{mtq}}$ is the total torque due to magnetorquers and $^{A}\mathbf{T}_{\text{rw}}$ is the total torque due to reaction wheels written in {A}.

We compute each of these torques in the next subsections.

#### A.6.1. Magnetorquers

In Low Earth Orbit (LEO), Earth's magnetic field intensity is not negligible (Wertz 1978). Its interaction with magnetorquers is capable of producing torque in the spacecraft. Physically, the actuator consists of a coil with a ferromagnetic core or air core. Therefore, it is possible, by Kirchhoff's Laws, to model it as a resistor and inductor series circuit. The magnetic moment is in the direction of the normal to the coil plane and the sense is given by the Right-Hand Rule. Its intensity is given by the Ampere's law. We neglect the demagnetization effect. The total intensity of the magnetic moment, $m_{\text{mtq}}$, is the sum of the intensities of the magnetic moment generated by the coil and by the ferromagnetic core. In vector form, the magnetorquer's torque is given by Equation (A13)

$$ ^{A}\mathbf{T}_{\text{mtq}} = \sum_{i=1}^{3} ^{A}\mathbf{m}_{\text{mtq}_i} \times ^{A}\mathbf{B} = ^{A}\mathbf{m}_{\text{mtq}} \times ^{A}\mathbf{B} \quad (\text{A13}) $$

where $m_{\text{mtq}_i}$ is the magnetic moment of each dipole, and $B$ the Earth's magnetic field.

#### A.6.2. Reaction Wheels

Reaction wheels (RW) allow more precise pointing control (Yost 2020). Considering $I_{\text{rw}}$ and $\omega_{\text{rw}}$, the moment of inertia and the angular velocity of reaction wheels, and applying the Euler characteristic equation, we obtain, for each reaction wheel i:

$$ ^{A}\mathbf{T}_{\text{rw}_i} = ^{A}\dot{\mathbf{H}}_{\text{rw}_i} + ^{A}\omega_{\text{A}} \times ^{A}\mathbf{H}_{\text{rw}_i} \quad (\text{A14}) $$

Only the first term of Equation (A14) ($^{A}\dot{\mathbf{H}}_{\text{rw}_i}$) can be used to control the spacecraft orientation. The second one ($^{A}\omega_{\text{A}} \times ^{A}\mathbf{H}_{\text{rw}_i}$) is a gyroscopic perturbation effect and it will be added to the satellite dynamics. Moreover, by action and reaction law, both terms must be negative. The final equation

is:

$$ ^{A}\mathbf{T}_{\text{mtq}} - ^{A}\dot{\mathbf{H}}_{\text{rw}} + ^{A}\mathbf{T}_{\text{dst}} = ^{A}\mathbf{I}_{\text{sc}}{}^{A}\dot{\omega}_{\text{A}} + ^{A}\omega_{\text{A}} \times (^{A}\mathbf{I}_{\text{sc}}{}^{A}\omega_{\text{A}} + ^{A}\mathbf{H}_{\text{rw}}). \quad (\text{A15}) $$

The simplicity of Equation (A15) is suitable for PicSat modeling purposes. Nonlinear effects such as RW's imbalance (microvibrations), Coulomb friction and magnetic modeling, magnetorquer's demagnetization could be included in additional terms, depending on the requirements for each mission.

### A.7. External Disturbances

This appendix section presents all relevant attitude and orbit external disturbances that PicSat faces during flight. Figure 21 presents a block diagram that summarizes them.

#### A.7.1. Attitude Disturbances

References Wertz (1978), Vallado (2013) and Curtis (2014) were used to implement attitude disturbances dynamic equations.

**Torque Due to Atmospheric Drag**—In orbit, the surface of a spacecraft collides with atmospheric particles, creating a force that will produce a torque in relation to the satellite's center of mass. The aerodynamic force on a surface element $dA$, written with respect to {S}, is given by

$$ ^{S}\mathbf{dF}_{\text{d}} = -\frac{1}{2}C_{\text{D}}\rho_{\text{rel}}^{2}(^{S}\hat{\mathbf{n}}_{\text{dA}} \cdot ^{S}\hat{\mathbf{v}}_{\text{rel}})^{2}^{S}\hat{\mathbf{v}}_{\text{rel}}dA, \quad (\text{A16}) $$

where $C_{\text{D}}$ is the drag coefficient, $\rho$ is the atmospheric density, $v_{\text{rel}}$ is the spacecraft's speed relative to the atmosphere, whose direction is represented by $\hat{\mathbf{v}}_{\text{rel}}$ and $\hat{\mathbf{n}}_{dA}$ is the unit vector of the surface element.

The torque due to atmospheric drag is expressed in {S} as

$$ ^{S}\mathbf{T}_{\text{d}} = \int (^{S}\mathbf{r}_{dA} \times ^{S}\mathbf{dF}_{\text{d}}) dA, \quad (\text{A17}) $$

where $\mathbf{r}_{dA}$ is the position vector of the spacecraft's center of mass up to the surface element $dA$ where the force $dF_d$ is being applied.

**Atmospheric Model**—Density ($\rho$) and other atmospheric properties vary depending on altitude, temperature, molecular composition, solar activity, etc. One of the most modern models used by the scientific community is the United States Naval Research Laboratory Mass Spectrometer and Incoherent Scatter Radar Exosphere (NRLMSISE-00) (Picone et al. 2002). It is an empirical model of the class Mass Spectrometer Incoherent Scatter Radar released in 2001 and valid up to an altitude of 1000 km. In this paper, the NRLMSISE-00 model integrated with the Simulink software is used.

**Torque Due to Solar Radiation Pressure**—The incidence of electromagnetic radiation on the spacecraft's surface generates disturbance torque. The main source of radiation is the direct radiation from the Sun. The intensity and direction of
<figure><img src="image_27.png" /><figcaption>Figure 21. Attitude and orbital external disturbances diagram.</figcaption></figure>

incidence, the shape and the optical properties of the satellite's surface are the main factors that influence the torque due to the solar radiation pressure.

Analogously to the torque due to atmospheric drag, we divide the spacecraft's surface into infinitesimal element and determine the force applied to each of them. Depending on the optical properties of the surface, radiation can be absorbed, reflected with the same incidence slope or reflected in a diffuse way. The differential radiation force due to each of these portions is given, respectively, by

$$^{S}d\mathbf{F}_{\text{ra}} = \frac{S}{c} [c_{\text{ra}}(-\cos \theta \, ^{S}\hat{\mathbf{n}}_{\text{dA}} + \sin \theta \, ^{S}\hat{\mathbf{s}}_{\text{dA}})] \cos \theta \, dA, \quad (\text{A18})$$

$$^{S}d\mathbf{F}_{\text{rr}} = \frac{S}{c} [(-1 + c_{\text{rr}}) \cos \theta \, ^{S}\hat{\mathbf{n}}_{\text{dA}} + (1 - c_{\text{rr}}) \sin \theta \, ^{S}\hat{\mathbf{s}}_{\text{dA}}] \cos \theta \, dA, \quad (\text{A19})$$

and

$$^{S}d\mathbf{F}_{\text{rd}} = \frac{S}{c} \left[ -\left(\cos \theta + \frac{2}{3}c_{\text{rd}}\right)^{S}\mathbf{n}_{\text{dA}} + \sin \theta \, ^{S}\mathbf{s}_{\text{dA}} \right] \cos \theta \, dA, \quad (\text{A20})$$

where S is the radiation intensity, c is the speed of light, $c_{\text{ra}}$, $c_{\text{rr}}$ and $c_{\text{rd}}$ are, respectively, the absorbed, reflected and diffuse radiation coefficients, constrained to $0 \le c_{\text{ra}} + c_{\text{rr}} + c_{\text{rd}} \le 1$. θ is the angle of incidence in relation to the normal, $\hat{\mathbf{n}}_{\text{dA}}$ and $\hat{\mathbf{s}}_{\text{dA}}$ are, respectively, the unit vector normal to the surface element and the unit vector parallel to the surface element and in the direction of radiation projection.

Thus the torque due to solar radiation pressure is given by

$$^{S}\mathbf{T}_{\text{rad}} = \int [({}^{S}\mathbf{r}_{\text{dA}} \times {}^{S}d\mathbf{F}_{\text{ra}}) + ({}^{S}\mathbf{r}_{\text{dA}} \times {}^{S}d\mathbf{F}_{\text{rr}}) + ({}^{S}\mathbf{r}_{\text{dA}} \times {}^{S}d\mathbf{F}_{\text{rd}})]dA, \quad (\text{A21})$$

where $\mathbf{r}_{\text{dA}}$ is the position vector of the spacecraft's center of mass up to the surface element $dA$ where forces $dF$ are being applied.

**Torque Due to Gravity Gradient**—The gravitational force depends on the distance to the object and its orientation. Therefore, each mass element of the spacecraft is subject to a different force. This force distribution produces torque when the mass of the satellite is not perfectly distributed. This phenomenon is referred to as the gravity gradient.

The torque due to the gravity gradient, $^{A}\mathbf{T}_{\text{gg}}$, when is written with respect to {A}, is given by

$$^{A}\mathbf{T}_{\text{gg}} = \frac{3\mu_{\oplus}}{r_{\text{sc}}^{3}} [^{A}\hat{\mathbf{r}}_{\text{sc}} \times (^{A}\mathbf{I}_{\text{sc}} \, ^{A}\hat{\mathbf{r}}_{\text{sc}})] \quad (\text{A22})$$

where $\mu_{\oplus}$ is the orbital parameter of the Earth, $\mathbf{r}_{\text{sc}}$ is the position vector of the spacecraft relative to the Earth and $\mathbf{I}_{\text{sc}}$ is the spacecraft inertia matrix.

**Torque Due to Residual Magnetic Dipole**—The interaction between the Earth's magnetic field and the spacecraft's magnetic dipole produces torque. This magnetic dipole may be desired, like the magnetic actuator presented in Appendix A.6.1, or undesired, like eddy current, hysteresis and other magnetic moments that may arise due to currents circulating in the electronic circuits of the spacecraft. Calling it the residual magnetic dipole, $\mathbf{m}_{\text{mag}}$, and using Equation (A13), the torque due to the residual magnetic dipole, $\mathbf{T}_{\text{mag}}$, is given by

$$^{A}\mathbf{T}_{\text{mag}} = ^{A}\mathbf{m}_{\text{mag}} \times ^{A}\mathbf{B}, \quad (\text{A23})$$
in which **B** is the Earth's magnetic field and its modeling is treated below.

**Model of the Earth's Magnetic Field**—With the help of satellites and observatories, mathematical models of the Earth's magnetic field have been developed since the 1960s. The models mostly used by the scientific community are the International Geomagnetic Reference Field (IGRF) (Thébault et al. 2015), developed by the International Association of Geomagnetism and Aeronomy, and the World Magnetic Model (WMM) (Chulliat et al. 2015) developed by the United State National Geophysical Data Center and British Geological Survey. Both models are updated every five years and are publicly available. Software such as MATLAB and Simulink have these models integrated. This paper uses the IGRF-12 version from 2015.

### A.7.2. Orbital Disturbances

References Vallado (2013), Curtis (2014) and Markley & Crassidis (2014) were used to implement orbital disturbances dynamic equations.

**Acceleration Due to Atmospheric Drag**—Atmospheric drag was initially explored in Appendix A.7.1. It interferes with the spacecraft's orbit by changing its acceleration. This disturbance in acceleration, **a**<sub>d</sub>, is given by

$$l\mathbf{a}_d = \frac{l\mathbf{F}_d}{m_{sc}}, \qquad (A24)$$

where $m_{sc}$ is the mass of the spacecraft and $F_d$ is the drag force as shown in Equation (A16).

**Acceleration Due to Solar Radiation Pressure**—Analogously to the drag, the solar radiation pressure interferes in the attitude and also in the orbit of the satellite.

The disturbance $a_{srp}$ is given by

$$l\mathbf{a}_{srp} = \frac{l\mathbf{F}_{ra} + l\mathbf{F}_{rr} + l\mathbf{F}_{rd}}{m_{sc}}, \qquad (A25)$$

where $m_{sc}$ is the mass of the spacecraft, $F_{ra}$ is the force due to the absorbed radiation (Equation (A18)), $F_{rr}$ is the force due to reflected radiation (Equation (A19)), and $F_{rd}$ is the force due to diffuse radiation (Equation (A20)).

**Acceleration Due to Earth's Oblateness**—The Earth's gravitational potential can be obtained through spherical harmonic expansion, generating an infinite series with associated Legendre polynomials.

To calculate the Earth's oblateness acceleration $a_{ob}$, it is necessary to truncate the series to a certain degree $m$ and order $n$, both parameters from Legendre polynomials. The most common implementation is performed with $m=0$ and $n=2$. However, as presented by Vallado (2013), for a high-fidelity mathematical model, it is necessary to consider truncation at least at $m=n=24$.

The most complete gravitational model is the Earth Gravitational Model 2008 (EGM2008) published by the U.S. National Geospatial-Intelligence Agency which allows truncation up to degree 2190 and order 2159 (Pavlis et al. 2012). In this paper, EGM2008 was used and the simulations were performed with truncation at $m=n=70$.

**Acceleration Due to Third-Body**—Kepler's Orbit takes into account the gravitational interaction between two bodies. Such a situation is often referred to as the Two-Body Problem. However, other celestial bodies, such as the Moon and the Sun, will also interact gravitationally with the spacecraft. This interaction can be considered as a disturbance and modeled as a three-body system. Applying Kepler's law and considering the Moon and the Sun as third-bodies, the acceleration is given by

$$l\mathbf{a}_{th} = \mu_m \left( \frac{l\mathbf{r}_{m/sc}}{r_{m/sc}^3} - \frac{l\mathbf{r}_m}{r_m^3} \right) + \mu_{\odot} \left( \frac{l\mathbf{r}_{\odot/sc}}{r_{\odot/sc}^3} - \frac{l\mathbf{r}_{\odot}}{r_{\odot}^3} \right), \qquad (A26)$$

where $\mu_m$ and $\mu_{\odot}$ are, respectively, the orbital parameters of the Moon and the Sun, $\mathbf{r}_{m/sc}$ and $\mathbf{r}_{\odot/sc}$ are, respectively, the position vectors of the Moon and the Sun with respect to the spacecraft, $r_m$ and $r_{\odot}$ are, respectively, the position vectors of the Moon and the Sun with respect to the Earth.

### ORCID iDs

Cauê G. Menegaldo https://orcid.org/0000-0003-1874-5396

Fábio de O. Fialho https://orcid.org/0000-0002-5463-5170

Eduardo Janot-Pacheco https://orcid.org/0000-0003-0079-3912

Felipe M. Pait https://orcid.org/0000-0001-8544-9306

Vincent Lapeyrère https://orcid.org/0000-0002-9651-8108

### References

Alshamy, H. M. I., Makled, A. E.-S., Elhalwagy, Y., & Hendy, H. 2019, IOP Conference Series: Materials Science and Engineering, 610, 012100

Avanzini, G., de Angelis, E. L., Giulietti, F., & Serrano, N. 2019, AcAau, 160, 625

Bangert, P., Busch, S., & Schilling, K. 2014, in Proc. 2nd IAA Conf. on Dynamics and Control of Space Systems (DYCOSS), 513, ed. F. Graziani, A. D. Guerman, & J. Contant (Rome: Univelt)

Carrara, V. 2015, in Proc. 17th Int. Symp. on Dynamic Problems of Mechanics, 1, ed. V. Steffen, D. A. Rade Jr, & W. M. Bessa (Natal RN: ABCM), http://abcm.org.br/anais/diname/2015/proceedings.html

Chulliat, A., Macmillan, S., Alken, P., et al. 2015, The US/UK World Magnetic Model for 2015–2020, Tech. Rep., National Geophysical Data Center, NOAA.

Craig, J. J. 1989, Introduction to Robotics: Mechanics and Control (Reading, Mass: Addison-Wesley)

Curtis, H. 2014, Orbital Mechanics for Engineering Students, Aerospace Engineering (3rd ed.; Oxford: Butterworth-Heinemann), 768

Dechao, R., Tao, S., Lu, C., Xiaoqian, C., & Yong, Z. 2014, ChJA, 27, 593

Francois-Lavet, V. 2010, Master's thesis, University of Liège, Liège, http://researchbank.rmit.edu.au/view/rmit:6277

Giebelmann, J. 2006, Master's thesis, RMIT University, Melbourne, http://www.raumfahrt.fh-aachen.de/compass-1/download/Development%20of%20an%20Active%20Magnetic%20Attitude%20Determination%20and%20Control%20System%20for%20Picosatellites%20on%20highly%20inclined%20circular%20Low%20Earth%20Orbits.pdf
Habibkhah, S., Arasi, J., & Bolandi, H. 2017, CSE, 19, 38

Holst, R. 2014, Master's thesis, Aalborg University, Aalborg, https://projekter.aau.dk/projekter/files/202321144/report.pdf

Jensen, K. F., & Vinther, K. 2010, Master's thesis, Aalborg University, Aalborg

Jonsson, L. 2019, Master's thesis, Luleå University of Technology, Space Technology, http://ltu.diva-portal.org/smash/record.jsf?pid=diva2:1298491&dswid=3150

Krogh, K., & Schreder, E. 2002, Attitude Determination for AAU CubeSat, Tech. Rep., No. IAS-1030, Aalborg University, Aalborg, http://studentspace.aau.dk/cubesat/dokumenter/ADC-report.pdf/

Kulu, E. 2020, Nanosats Database, https://www.nanosats.eu

Lagrange, A. M., Rubini, P., Nowak, M., et al. 2020, A&A, 642, A18

Langer, M., Weisgerber, M., Bouwmeester, J., & Hoehn, A. 2017, in Proc. 43rd IEEE Aerospace Conf., ed. M. S. Greco (Big Sky, MT: IEEE), 1

Lapeyrère, V., Lacour, S., David, L., et al. 2017, in Proc. 31st Annual AIAA/USU Conf. on Small Satellites, ed. P. Patterson (Logan: UT: AIAA/USU), 1, https://digitalcommons.usu.edu/smallsat/2017/all2017/86/

Liddle, J. D., Holt, A. P., Jason, S. J., O'Donnell, K. A., & Stevens, E. J. 2020, NatAs, 4, 1026

Liu, Y. H., Shangguan, D. S., Yao, J., Chen, L. P., & Ding, J. W. 2019, in Proc. 9th Asia Conf. on Mechanical and Aerospace Engineering 1215, ed. H. R. Li (Bristol: IOP Publishing), 105

Lovera, M. 2003, in Proc. 4th IMACS Symp. on Mathematical Modelling, ed. I. Troch & F. Breitenecker (Vienna: ARGESIM), 498

Markley, F. L., & Crassidis, J. L. 2014, Fundamentals of Spacecraft Attitude Determination and Control (New York: Springer)

Mason, J. P., Baumgart, M., Rogler, B., et al. 2017, JSSat, 6, 651

Menges, B., Guadiamos, C., & Pernicka, H. 1998, in Proc. 36th Aerospace Sciences Meeting and Exhibit (Reno, NV: AIAA)

Millan, R. M., von Steiger, R., Ariel, M., et al. 2019, AdSpR, 64, 1466

Nakasuka, S., Miyata, K., Tsuruda, Y., Aoyanagi, Y., & Matsumoto, T. 2018, AcAau, 145, 515

Naqvi, N. A., & Raza, M. 2007, in Proc. 3rd Int. Conf. on Recent Advances in Space Technologies, ed. C. W. Chen (Piscataway, NJ: IEEE), 448

Nasirian, M., Bolandi, H., Sedigh, A. K., & Khoogar, A. R. 2006, in Proc. 1st Int. Symp. on Systems and Control in Aerospace and Astronautics, ed. B. Filman (Harbin: IEEE), 160

Nowak, M. 2019, Theses, Université Paris Sciences et Lettres, Paris, FR, https://tel.archives-ouvertes.fr/tel-02890983

Nowak, M., Lacour, S., Crouzier, A., et al. 2018, Proc. SPIE, 10698, 1069821-1

Nowak, M., Lacour, S., Lapeyrère, V., et al. 2017, in Proc. 31st Annual AIAA/USU Conf. on Small Satellites, 1, ed. P Patterson (Logan, UT: AIAA/USU), https://digitalcommons.usu.edu/smallsat/2017/all2017/103/

Omar, K. H. 2017, in Proc. 55th AIAA Aerospace Sciences Meeting (Grapevine, TX: American Institute of Aeronautics and Astronautics), 1

Pavlis, N. K., Holmes, S. A., Kenyon, S. C., & Factor, J. K. 2012, JGRB, 117, B4

Picone, J. M., Hedin, A. E., Drob, D. P., & Aikin, A. C. 2002, JGRA, 107, A12

Polat, H. C., Virgili-Llop, J., & Romano, M. 2016, JSSat, 5, 513

Puig-Suari, J., Turner, C., & Twiggs, R. 2001, in Proc. 15th Annual AIAA/USU Conf. on Small Satellites, ed. P. Patterson (Logan, UT: AIAA/USU), 1, https://digitalcommons.usu.edu/smallsat/2015/all2015/58/

Rawashdeh, S. A. 2019, International Journal of Aerospace Engineering, 2019, 1

Riesing, K. 2015, in Proc. 29th AIAA/USU Conf. on Small Satellites, 1, ed. P. Patterson (Logan, UT: AIAA/USU), https://digitalcommons.usu.edu/smallsat/2015/all2015/58/

Romano, R. A., & Pait, F. 2017, ITAC, 62, 2180

Rondão, D. O. d. M. A. 2016, Master's thesis, Técnico Lisboa, https://fenix.tecnico.ulisboa.pt/cursos/meaer/dissertacao/846778572211765/

Sarda, K., Grant, C., Chaumont, M., et al. 2014, in Proc. 28th Annual AIAA/USU Conf. on Small Satellites, ed. P. Patterson (Logan, UT: (AIAA/USU), 1, https://digitalcommons.usu.edu/smallsat/2014/YearReview/4/

Sato, Y., Fujita, S., Kuwahara, T., et al. 2017, in Proc. 2017 IEEE/SICE Int. Symp. on System Integration (SII), ed. M. S. Greco (Taipei: IEEE), 686

Scholz, A., Ley, W., Dachwald, B., Miau, J. J., & Juang, J. C. 2010, AcAau, 67, 1289

Sebestyen, G., Fujikawa, S., Galassi, N., & Chuchra, A. 2018, Low Earth Orbit Satellite Design (Berlin: Springer)

Seidelmann, P. K. 2006, Explanatory Supplement to the Astronomical Almanac (Sausalito, CA: Univ. Science Books)

Sidi, M. J. 1997, Spacecraft Dynamics and Control (Cambridge: Cambridge Univ. Press), 432

Slavinskis, A., Ehrpais, H., Kuuste, H., et al. 2016, JAerE, 29, 1

Springmann, J. C., & Cutler, J. W. 2014, AcAau, 99, 201

Sternberg, D. C., Pong, C., Filipe, N., et al. 2018, JSpRo, 55, 322

Stickler, A. C., & Alfriend, K. T. 1976, JSpRo, 13, 282

Taraba, M., Rayburn, C., Tsuda, A., & MacGillivray, C. 2009, in Proc. 23rd Annual AIAA/USU Conf. on Small Satellites, 1, ed. P Patterson (Logan, UT: AIAA/USU), https://digitalcommons.usu.edu/smallsat/2009/all2009/64/

Thébault, E., Finlay, C. C., Beggan, C. D., et al. 2015, EP&S, 67, 79

Thomsen, B. G., & Nielsen, J. 2016, Master's thesis, Aalborg University, Aalborg, https://projekter.aau.dk/projekter/files/239482913/report_final.pdf

Triharjanto, R. H., Ubay, M. S. N., Utama, S., & Poetro, R. E. 2015, in Proc. 2015 IEEE Int. Conf. on Aerospace Electronics and Remote Sensing Technology (ICARES), ed. G. Setti (Bali: IEEE), 1

Turner, A. J. 2003, Master's thesis, Virginia Tech, Blacksburg, https://vtechworks.lib.vt.edu/handle/10919/34727

Vallado, D. 2013, Fundamentals of Astrodynamics and Applications (Hawthorne, CA: Microcosm Press)

Vallado, D., Crawford, P., Hujsak, R., & Kelso, T. S. 2006, in AIAA/AAS Astrodynamics Specialist Conf. and Exhibit, 1 (Reston, VA: American Institute of Aeronautics and Astronautics)

van Vuuren, G. H. J. 2015, Master's thesis, Stellenbosch University, Stellenbosch, https://scholar.sun.ac.za/handle/10019.1/96979

Villela, T., Costa, C. A., Brandão, A. M., Bueno, F. T., & Leonardi, R. 2019, International Journal of Aerospace Engineering, 2019, 1

Wertz, J. R. 1978, Spacecraft Attitude Determination and Control, Vol. 73 (Dordrecht: Reidel), 858

Xiang, T., Meng, T., Wang, H., Han, K., & Jin, Z.-H. 2012, AcAau, 77, 182

Yost, B. 2020, State of the Art of Small Spacecraft Technology, https://sst-soa.arc.nasa.gov