Thermal Analysis of Satellite Libertad 2: a Guide to CubeSat Temperature Prediction

ABSTRACT:

We present a model for predicting the temperature of three-unit CubeSat on a low Earth orbit, which supposes a single temperature common to all satellite components. Our exposition includes a detailed, to a large extent analytical, computation of the external heat fluxes for a particular orbit and spacecraft assumptions based on the features foreseen for satellite Libertad 2 under development at Universidad Sergio Arboleda. Moreover, supported by specialized thermal analysis software, we compute the heat fluxes and their associated temperature for all possible orbital orientations, and combine these results with a description of the satellite orbital plane rotation (nodal regression) and the solar motion on the ecliptic, to determine the minima and maxima of the orbital temperature oscillation for a mission lifetime of a year. We find that, for feasible model parameters, the temperature extremes are mostly within the operating temperature range of the most sensitive satellite component, 0 ºC ≤ T ≤ 60 ºC, suggesting mission viability. Finally, we discuss possible model improvements which would allow testing of satellite design upgrades. In this regard, it is worth noting that the calculation of the external heat fluxes here described can be carried over, almost unchanged, to a more accurate model describing heat transfer between satellite parts with different temperatures.

KEYWORDS:
CubeSat; Low Earth orbit; Thermal analysis; Nodal regression; Beta angle; Numerical simulation; Linearization

INTRODUCTION

The Universidad Sergio Arboleda of Bogotá, Colombia, is building a nano-satellite (for a classification of satellites according to their mass see Fortescue et al. 2003) with the main goal of capturing photographs of Earth's surface, with potential use in precision agriculture. The satellite, named Libertad 2, will have the shape of a rectangular parallelepided of sides 30 cm × 10 cm × 10 cm (Fig. 1), thus conforming to the CubeSat standard as a three-unit (3U) spacecraft (a description of this standard can be found on The Cubesat Program 2015). Libertad 2 constitutes the continuation of Universidad Sergio Arboleda's satellite development program started by pico-satellite Libertad 1 (1U, NORAD Catalog Number 31128), launched on April 17, 2007 (Joya 2007; NASA 2014).

Figure 1
Views of Libertad 2 with the body fixed coordinate system {x,y,z}. The faces are labeled by a numerical index and a word. The front face has dimensions 10 cm × 10 cm and the bottom face is 10 cm × 30 cm.The dark gray hexagons represent solar cells while the light gray circle represents the camera lens opening.

The CubeSat standard was developed in 1999 with the objective, among others, of facilitating the involvement of universities in the aerospace industry. This goal could be evaluated as "fulfilled" on the basis of a review by Swartwout (2013), that counted 77 university-led CubeSat-class missions out of a total of 112, between 2000 and 2012. However, the same review reports that university missions have been plagued by a high rate of unsuccessful missions, contributing 27 out of the total 34 failures. Many university missions were also found at the lowest level of mission impact, that of "beepsats", able only to send basic telemetry. A cause of this situation is lack of training of the teams undertaking CubeSat missions at universities, since the "easy path" to space may have attracted many newcomers. Hence, to channel fruitfully this interest, it is urgent to develop resources for a rapid and low-cost acquisition of the pertinent aerospace education. Failing to do so, we may be defeating one of the purposes of the CubeSat standard by accepting that relevant CubeSat missions are the sole privilege of well-established institutions (universities, government organizations, or businesses) of the aerospace sector. An example of a CubeSat mission by experienced developers is the Dynamic Ionosphere CubeSat Experiment (DICE) (Fish et al. 2014); a growing interest of commercial developers in nano-satellites was reported by Buchen (2014).

Acknowledging this need, this paper intends to be a resource for quick learning of the concepts and methods involved in predicting the temperature of a CubeSat. Necessary background material is included to make the presentation self-contained, model assumptions are given a clear mathematical formulation, and the computations feature intermediate steps to facilitate their reproduction, thereby making the paper accessible to anyone with a solid quantitative training.

Prediction of a satellite's temperature is justified by the need for all satellite components to function within their operating temperature ranges - defined as "the maximum and minimum temperature limits between which components successfully and reliably meet their specified operating requirements" (Garzon 2012) -, otherwise risking malfunctioning or damage. To forecast the effect of design choices, like the materials on the external surfaces, on the temperature, it is necessary to create a model of the heat transferred between the satellite and its surroundings and between the satellite parts.

The modeling techniques can be found in the available literature on thermal analysis of nano-satellites. Dinh (2012) studied the temperature of internal electronics in a 1U CubeSat using the software packages Thermal Desktop and ANSYS. Jacques (2009) conducted the thermal analysis for the OUFTI-1U CubeSat relying on the software ESATAN-TMS, whereas Garzon (2012) simulated the temperatures of the OSIRIS-3U CubeSat employing the multiphysics software COMSOL. Bulut and Sozbir (2015) investigated the temperature for different solar panel combinations in a 1U CubeSat. Optimization of thermal design parameters through genetic algorithms was pursued by Escobar et al. (2016). Kang and Oh (2016) performed an experimental validation of their model predictions using a thermal vacuum chamber. Active thermal control with phase change materials was investigated by Shinde et al. (2017) and Kang and Oh (2016). A simulation code implemented in MATLAB was developed by Corpino et al. (2015). Mason et al. (2018) compared a thermal model of the Miniature X-Ray Solar Spectrometer 3U CubeSat with actual on-orbit temperature measurements, finding agreement within a few degrees Celsius. Remarkably, their thermal design allowed for the payload to stay at an almost constant on-orbit temperature of -40.91 ºC (standard deviation of 0.19 ºC), isolated from the much wider temperature oscillations of the solar arrays, roughly from -40 ºC to 50 ºC.

In the present study, a model for predicting the temperature of satellite Libertad 2 is examined. For simplicity, the internal heat transfer is supposed to happen instantaneously, resulting in a single temperature shared by all satellite components. The ensuing single temperature model (also known as a single node model), serves as a stepping stone toward the development of a more accurate, multiple node model (with different parts allowing different temperatures), which will be the subject of a future work. In particular, the irradiances associated to the external heat sources, here presented, can be passed unchanged to the multiple node model.

This study is structured as follows. In the section Design of Libertad 2, the satellite components and their operating temperature ranges are presented. The following section introduces the single node model for a special orbit and shows the model-predicted temperature (model solutions) together with some verification criteria. Then, the variation of the external heat fluxes over a mission life of a year is considered; the resulting effect on the satellite temperature is computed with the help of the thermal analysis software ESATAN-TMS release 5. Finally, the last section summarizes our results and discusses possible future work.

DESIGN OF LIBERTAD 2

The satellite will carry several components in addition to the photographic camera. A system of wheels and electromagnets will control the satellite orientation for the purpose of pointing the camera toward a particular target. These actuators together with orientation sensors constitute what is known as the Attitude Determination and Control System (ADCS). Instructions will be sent to the satellite via Very High Frequency (VHF) and Ultra High Frequency (UHF) radio waves. S-Band microwaves will be used to transmit the photographs as digital data. Hence, the satellite will have transmitter and receiver equipment on the above mentioned telecommunication bands (named, respectively, VHF-UHF transceiver and S-band transceiver). Solar cells will provide the satellite energy, part of which will be stored in a battery (to provide energy in the periods of low or nil solar radiation, for instance, when the satellite traverses Earth's shadow), that will be distributed between the previously referred components by a special circuitry, the Electric Power System (EPS). The processing of instructions and the coordination of all the components tasks will be performed by an On Board Computer (OBC).

Table 1 shows the operating temperature range (T^O_min, T^O_max), for each satellite subsystem. Operation outside these ranges should be avoided.

Table 1
Operating minimum and maximum temperature values (T^O_min and T^O_max, respectively) for the satellite subsystems.

Subsystem	TOmin (ºC)	TOmax (ºC)
Photographic camera	0	60
Batteries	0	85
ADCS	-20	85
OBC	-40	85
EPS	-40	85
VHF-UHF transceiver	-40	85
S-Band transceiver	-40	85
Solar panels	-40	150

SINGLE NODE MODEL

In this section, the thermal model is developed from first principles. The variation of the satellite temperature T is governed by the first law of thermodynamics for the case in which no work is performed by (or over) the system. Accordingly, the instantaneous rate of change of the satellite's internal energy CdT/dt, with C the satellite's heat capacity, is equal to the difference between the heat fluxes entering the satellite and those leaving it. This results in the differential equation for T (Eq. 1), as explained next.

Satellite Libertad 2 will be put on an orbit with altitude in the range 600 - 900 km, which fits into the definition of a Low Earth Orbit (LEO), according to Vallado (1997). At this altitude, the most significant heat fluxes entering the i-th satellite face (cf. Fig. 1 for a labeling of the faces) are Q_S,i(t), Q_alb,i(t), and Q_E,i(t), whose sources are, respectively, the direct solar radiation, the solar radiation reflected on the surface of Earth (albedo radiation), and the infrared radiation emitted by Earth (Fortescue 2003; Gilmore 2002). For the outgoing heat flux, we assume that the i-th satellite face emits infrared radiation as a gray body, with radiosity ε_iσT⁴, where ε_i is the face's emissivity, 0 < ε_i < 1, a property of the face's surface finish, and σ is the Stefan-Boltzmann's constant. Combining the previous assumptions leads to the following evolution equation for the temperature T (Eq. 1),

(1) $\mathit{C}\frac{\mathit{dT}}{\mathit{dt}}=-\sigma {\mathit{A}}_{\mathit{eff}}{\mathit{T}}^4+{\mathit{Q}}_{\mathit{tot}}\left(\mathit{t}\right),$

where (Eq. 2),

(2) ${\mathit{A}}_{\mathit{eff}}=\sum _{\mathit{i}=1}^6{\mathit{A}}_{\mathit{i}}{\epsilon }_{\mathit{i}}$

with A_i the area of the i-th face, and Q_tot (t) is the total external heat flux (Eq. 3),

(3) $\begin{array}{c}{\mathit{Q}}_{\mathrm{tot}}\left(\mathit{t}\right)={\mathit{Q}}_s\left(\mathit{t}\right)+{\mathit{Q}}_{\mathrm{alb}}\left(\mathit{t}\right)+{\mathit{Q}}_E\left(\mathit{t}\right)\\ \mathrm{with}{\mathit{Q}}_{\mathit{S}}\left(\mathit{t}\right)={\sum }_{\mathit{i}=1}^6{\mathit{Q}}_{\mathit{S},\mathit{i}}\left(\mathit{t}\right),{\mathit{Q}}_{\mathit{alb}}\left(\mathit{t}\right)={\sum }_{\mathit{i}=1}^6{\mathit{Q}}_{\mathit{alb},\mathit{i}}\left(\mathit{t}\right),\mathrm{and}{\mathit{Q}}_{\mathit{E}}\left(\mathit{t}\right)={\sum }_{\mathit{i}=1}^6{\mathit{Q}}_{\mathit{E},\mathit{i}}\left(\mathit{t}\right).\end{array}$

The heat fluxes Q_S,i (t), Q_alb,i (t), and Q_E,i (t) depend on the satellite orbit and orientation, with Q_S,i (t) and Q_alb,i (t) further affected by the direction of the sunrays, denoted by the (unit) solar vector ŝ, which points toward the sun. In the following sections we show, in a special case, how these factors enter the computation of the heat fluxes.

SATELLITE ORBIT, SATELLITE ATTITUDE, AND SOLAR VECTOR

The satellite orbit is described in relation to the Geocentric Equatorial Coordinate System (GECS) (Vallado1997) with origin at Earth's center and axes {X, Y, Z}, with unit vectors {û₁, û₂, û₃}, oriented as follows. û₃ points toward the North Pole (Fig. 2), while û₁ and û₂ are contained in the plane of Earth's equator and have fixed directions (described in the next section) in relation to the distant stars.

Figure 2
Satellite orbit and orientation in relation to the Geocentric Equatorial Coordinate System {û₁, û₂, û₃}.

Throughout the paper, we assume for Libertad 2 a perfectly circular orbit of radius R = 7110 km (similar to the orbit of Libertad 1, discussed in the next section). In this section, the orbit is supposedly in the û₁ - û₃ plane (Fig. 2) and the solar vector is taken as ŝ = û₁. General orbital planes and ŝ are considered in the following sections.

The satellite's orientation (or attitude) is described in terms of a body fixed coordinate system {x, y, z}, with unit vectors {ê₁, ê₂, ê₃}, oriented with respect to the satellite as shown in Fig. 1. The camera's optical axis will be parallel to ê₁ and the lens will look through an orifice on face 2. A satellite on a LEO is subject to atmospheric drag that reduces the satellite's kinetic energy causing it to eventually fall on Earth. To minimize atmospheric friction, the satellite's longer edges, parallel to ê₃, will be aligned with the satellite's velocity. The vector -ê₁ will point toward the center of Earth (-ê₁ parallel to nadir) to enable the camera to photograph the globe's surface. This orientation scenario (Fig. 2), which we label as velocity-nadir attitude, could be maintained throughout the entire orbit, despite external torques, thanks to the ADCS.

For quick face identification, we find it convenient to use a word labeling (Fig. 1), alternative to the index i, based on the mnemonic of imagining the satellite as a miniature land vehicle with the vector from the rear to the windshield parallel to the velocity and its bottom looking toward Earth.

Having specified the orbit and satellite orientation, we proceed to the calculation of the incoming heat fluxes.

CALCULATION OF DIRECT SOLAR RADIATION

Figure 3 shows the satellite on its orbit as viewed from the -û₂ direction. An arc of the orbit, with angle measure 2ξ, given by (Eq. 4),

(4) $\mathit{\xi }={\sin }^{-1}\left({\mathit{R}}_{\mathit{E}}/\mathit{R}\right)\approx 63.8°$

Figure 3
Satellite orbit as viewed from the -û₂ direction (see Fig. 2 for the coordinate system definition). The dark gray region represents the orbit eclipse zone. ξ is half the angle measure of the eclipse zone and v is the true anomaly.

with R_E = 6378 km (Bate et al. 1971) the radius of Earth, passes through the shadow projected by Earth (eclipse zone) where the direct solar radiation heat fluxes (solar heat fluxes, for short) Q_S,i are zero. Outside the eclipse zone, Q_S,i will depend on the angle b_i between the outward unit vector normal to the face, nˆ_i, and the solar vector ŝ (Eq. 5)

(5) $\cos {\mathit{b}}_{\mathit{i}}={\hat{\mathit{n}}}_{\mathit{i}}·\hat{\mathit{s}}.$

cos b_i ≤ 0 means the sunrays do not reach the face, so, also in this case, Q_S,i = 0. When cos b_i > 0, the face will project, onto a plane perpendicular to the solar rays, an area A cos b_i, absorbing a heat flux (Eq. 6)

(6) ${\mathit{Q}}_{\mathit{S},\mathit{i}}={\mathit{\alpha }}_{\mathit{i}}{\mathit{A}}_{\mathit{i}}\cos {\mathit{b}}_{\mathit{i}}{\mathit{I}}_{\mathit{S}},$

where α_i is the face's absorptivity, 0 < α_i < 1 which, as the emissivity, is a property of the face's external surface finish, and I_S is the solar irradiance (the energy transported by the solar radiation per unit of area and unit of time). The solar irradiance I_S depends on the satellite-sun distance, implying that, for a LEO, I_S is highest (lowest) at perihelion (aphelion). In this work we use the recommended mean value, known as the solar constant, I_S = 1367 W/m² (Anderson et al. 2001), thus neglecting the variation of I_S around this mean (±3.4%).

The angle b_i will be determined by the satellite position on the orbit, denoted by the true anomaly v, measured from the û₁ direction (Fig. 3). To calculate the product nˆ_i · ŝ we express both nˆ_i and ŝ in the {û₁, û₂, û₃} basis. Therefore, for the top face, we have nˆ₄ = eˆ₁ = cos v uˆ₁ + sin v uˆ₃. So, cos b₄ = nˆ₄ · ŝ = cos v. Similar calculations apply to the remaining faces.

Faces left, right and space will be partially covered by six solar cells (Fig. 1), each cell with an area of 30.18 cm² (AZUR SPACE 2016). A fraction of the solar power absorbed by a cell, equal to the cell's efficiency η = 28% (AZUR SPACE 2016), becomes electric energy rather than heat. Hence, the cell's effective absorptivity is given by α_eff = α_c - η, where α_c = 0.91 is the cell's nominal absorptivity. As a result, each one of the faces possessing solar cells will exhibit an average absorptivity (Eq. 7),

(7) ${\mathit{\alpha }}_{\mathit{avg}}={\mathit{f}}_{\mathit{c}}{\mathit{\alpha }}_{\mathit{eff}}+\left(1-{\mathit{f}}_{\mathit{c}}\right){\mathit{\alpha }}_{\mathit{u}},$

where f_c = 60.36% is the fraction of the face's area covered by solar cells, and α_u = 0.5 is the absorptivity of the non-solar-cell surfaces, yielding α_avg = 0.578.

Figure 4, shows the total solar heat flux as function of the true anomaly v.

Figure 4
Total solar heat flux Q_S (solid line) and total albedo heat flux Q_alb (dashed line) as functions of the true anomaly v. Both Q_S and Q_alb fall to zero inside the eclipse zone, 180º - ξ ≤ v ≤ 180º + ξ (see Fig. 3).

By considering the variation in time of v, we can express the heat fluxes as functions of time, Q_S,i (t) = Q_S,i [v(t)]. For a circular orbit, the true anomaly changes with constant angular velocity ω (Eq. 8),

(8) $\mathit{v}\left(\mathit{t}\right)={\mathit{v}}_0+\mathit{\omega }\mathit{t},$

where the initial value is chosen for convenience as v₀ = 0 and ω is given by Kepler's third law (Vallado 1997) (Eq. 9),

(9) $\mathit{\omega }=\sqrt{\frac{\mathit{\mu }}{{\mathit{R}}^3}},$

where µ = 3.986004415 × 10⁵ km³/s² is the Earth's gravitational parameter (Vallado 1997). This gives an orbital period (Eq. 10),

(10) $\mathit{P}=\frac{2\mathit{\pi }}{\mathit{\omega }}\approx 99.4\min .$

CALCULATION OF INFRARED RADIATION

The Earth emits thermal radiation with highest values of the spectral distribution inside the infrared wavelength band (Anderson et al. 2001; Anderson and Smith 1994). A patch on the surface of Earth is approximated as a perfectly diffuse radiator (Palmer and Grant 2010; Wolfe 1998) with radiosity I_E, defined as the radiant flux leaving the patch per unit area. I_E is higher (lower) in warmer (colder) areas of the globe and clouds absorb infrared radiation, decreasing the value of I_E perceived by a spacecraft on orbit (Anderson et al. 2001; Anderson and Smith 1994). As a result, I_E changes both in space and time, I_E = I_E (θ, ϕ, t), where θ and ϕ are angular coordinates describing a point on the surface of Earth. However, for the sake of simplicity, I_E may be modeled as uniform over the globe, I_E = I_E(t), with the time-variation taking into account the change of the part of the globe viewed by the satellite as it traverses its orbit, as well as the intrinsic time-fluctuation of I_E.

During the Earth Radiation Budget Experiment (ERBE) the infrared radiation falling over three LEO satellites was recorded every 16 s during 28 months (Anderson et al. 2001; Anderson and Smith 1994). The infrared intensity was reported as the radiosity I_E(t) of a spherical surface at an altitude of 30 km above the surface of the Earth, termed Top of Atmosphere (TOA). The measured I_E(t) changed erratically and was treated as a random time series. For orbits of inclination above 60º (see the next section for the definition of inclination) a time-averaged I_E -value of 211 W/m² was found (Anderson et al. 2001). In this work, we assume a spatially uniform time-constant radiosity I_E equal to this mean, scaled following conservation of energy from TOA to the surface of Earth by the factor (R_E + 30 km)²/, yielding the final value of I_E = 213 W/m².

Due to the uniformity and time-invariance of I_E, the Earth heat fluxes Q_E,i are independent of the satellite position on its orbit, thus constant in time. Conforming to Lambert's law, a surface element on the globe of area dS will emit, in the direction of the satellite, a radiant intensity (power per solid angle) I_E cos a_E dS/π, where a_E is the angle between the unit vector pointing from the surface element toward the satellite ρˆ and the unit vector normal to the surface element nˆ_E (see Fig. 5). The i-th face will intercept the radiant energy traveling in a solid angle A_i cos a_i/ρ², with a_i the angle between - ρˆ and nˆ_i. Of this energy a fraction ε_i will be absorbed. As a result, the contribution of the surface element to the heat flux Q_E,i is (Eq. 11),

(11) ${\mathit{dQ}}_{\mathit{E},\mathit{i}}={\epsilon }_{\mathit{i}}{\mathit{A}}_{\mathit{i}}{\mathit{I}}_{\mathit{E}}\frac{\cos {\mathit{a}}_{\mathit{E}}\cos {\mathit{a}}_{\mathit{i}}}{{\mathit{\pi }\mathit{\rho }}^2}\mathit{dS}.$

Figure 5
Quantities involved in the calculation of the infrared heat flux over the i-th satellite face, Q_iE (Eqs. 12 and 13).

Integrating over all surface elements yields (Eq. 12),

(12) ${\mathit{Q}}_{\mathit{E},\mathit{i}}={\epsilon }_{\mathit{i}}{\mathit{A}}_{\mathit{i}}{\mathit{I}}_{\mathit{E}}{\mathit{F}}_{\mathit{iE}},$

where F_iE is the radiative view factor (Eq. 13),

(13) ${\mathit{F}}_{\mathit{iE}}=\int \frac{\cos {\mathit{a}}_{\mathit{E}}\cos {\mathit{a}}_{\mathit{i}}}{{\mathit{\pi }\mathit{\rho }}^2}\mathit{dS}.$

To evaluate the integral in (Eq. 13), we express it in spherical coordinates (r, θ, ϕ), with û₃ acting as the z-axis. We begin recalling that cos a_E = ρˆ · nˆ_E and cos a_i = - ρˆ · nˆ_i. Next, we define ρ = r_sat - r_E, where r_sat is the satellite position and r_E is the surface element position (Eq. 14),

(14) ${\mathit{r}}_{\mathit{E}}={\mathit{R}}_{\mathit{E}}\left(\sin \mathit{\theta }\cos \varphi {\hat{\mathit{u}}}_1+\sin \mathit{\theta }\sin \varphi {\hat{\mathit{u}}}_2+\cos \mathit{\theta }{\hat{\mathit{u}}}_3\right)$

Therefore, ρ = ||ρ||, ρˆ = ρ/ρ, and nˆ_E = r_E /R_E. The view factor F_iE is independent of the true anomaly v. Hence, for ease, we chose r_sat = Rû₃. As an example, we compute F_iE for the bottom face, F_2E, in which case, nˆ_i = -û₃. Substituting the previous vector relationships into (Eq. 13) leads to (Eq. 15),

(15) ${\mathit{F}}_{2\mathit{E}}=\frac{1}{\mathit{\pi }}{\int }_0^{{\mathit{\theta }}_{\mathit{M}}}{\int }_0^{2\mathit{\pi }}\frac{\left(1-\mathit{h}\cos \mathit{\theta }\right)\left(\cos \mathit{\theta }-\mathit{h}\right)}{{\left(1+{\mathit{h}}^2-2\mathit{h}\cos \mathit{\theta }\right)}^2}\sin \mathit{\theta }\mathit{d}\varphi \mathit{\theta },$

where h = R/R_E and the limit of integration θ < θ_M, with θ_M = cos^-1 (R_E/R), is determined by the condition cos a_E > 0, which states that infrared rays do not traverse the Earth. Evaluation of Eq. 15 yields (Eq. 16),

(16) ${\mathit{F}}_{2\mathit{E}}=\frac{1}{{\mathit{h}}^2}\approx 0.80.$

For faces front, rear, left and right (i = 1,3,5,6), evaluation of Eq. 13 yields (Eq. 17)

(17) ${\mathit{F}}_{\mathit{iE}}=-\frac{\sqrt{{\mathit{h}}^2-1}}{\mathit{\pi }{\mathit{h}}^2}+\frac{1}{\mathit{\pi }}{\tan }^{-1}\left(\frac{1}{\sqrt{{\mathit{h}}^2-1}}\right),\approx 0.23.$

We assume an emissivity ε_u = 0.05, corresponding to a shiny metallic surface finish (Gilmore 2002), for all surfaces not covered by solar cells (this includes the totality of the faces devoid of solar cells). Hence, the faces that do possess solar cells will display average emissivity (Eq. 18),

(18) ${\epsilon }_{\mathit{avg}}={\mathit{f}}_{\mathit{c}}{\epsilon }_{\mathit{c}}+\left(1-{\mathit{f}}_{\mathit{c}}\right){\epsilon }_{\mathit{u}}$

where ε_C = 0.89 is the solar cells' emissivity (AZUR SPACE 2014), so ε_avg = 0.557.

Table 2 presents the value of Q_E,i calculated for each face.

Table 2
Infrared heat fluxes Q_E,i computed using Eq. 12 and the program ESATAN. The table total corresponds to Q_E = ∑⁶_{i = 1} Q_E,i

Face	QE,i(W)
	Eq. 12	ESATAN
Front	0.0243	0.0242
Bottom	0.2571	0.2585
Left	0.8119	0.8006
Top	0	0
Right	0.8119	0.8404
Rear	0.0243	0.0247
Total	1.9294	1.9484

CALCULATION OF REFLECTED SOLAR RADIATION

We suppose that a patch on the surface of Earth reflects light diffusely and equally in all directions (Lambertian reflectance) with diffuse reflectance (also known as albedo) γ, defined as the reflected fraction of the incident radiant flux. More accurate models use the bidirectional reflectance distribution function (Schaepman-Strub et al. 2006), employed, for instance, in processing the data provided by NASA's MODIS Instruments on board the Aqua and Terra satellites (NASA 2015; Strahler et al. 1999). The albedo varies considerably across the globe: it is usually higher over areas covered with ice and snow, deserts, and cloudy regions, and generally lower over the oceans, in the absence of clouds (Anderson et al. 2001; Anderson and Smith 1994). Therefore, the albedo depends on place and time, γ = γ(θ, ϕ, t), but, for simplicity, it is modeled as if it were uniform, γ = γ(t). The albedo γ(t) of the TOA measured during the ERBE behaved like a random time series (as did I_E(t)) with a mean value of 0.23 for orbits of inclination superior to 60º (Anderson et al. 2001). In this work, we assume for γ a time-constant value equal to this mean plus the correction for zero β angle (see the next section for the definition of the β angle) of 0.04 found in Anderson et al. (2001), scaled from the TOA to the surface of Earth, for a final value of γ = 0.273.

For the calculation of the albedo radiation heat fluxes, we follow Fitz et al. (1963). An element on the surface of Earth, of area dS, diffusely reflecting solar radiation may be considered as a diffuse radiator with radiosity cos (b_E) γI_S, where b_E is the angle between the normal to the surface element nˆ_E and the solar vector ŝ (Eq. 19),

(19) $\cos \left({\mathit{b}}_{\mathit{E}}\right)={\hat{\mathit{n}}}_{\mathit{E}}·\hat{\mathit{s}}.$

Hence, the albedo heat flux contributed by the surface element can be obtained by replacing in Eq. 11I_E by cos (b_E) γI_S and ε_i by α_i, yielding (Eq. 20),

(20) ${\mathit{dQ}}_{\mathit{alb},\mathit{i}}={\mathit{\alpha }}_{\mathit{i}}{\mathit{A}}_{\mathit{i}}\mathit{\gamma }{\mathit{I}}_{\mathit{S}}\frac{\cos {\mathit{a}}_{\mathit{E}}\cos {\mathit{a}}_{\mathit{i}}\cos {\mathit{b}}_{\mathit{E}}}{{\mathit{\pi }\mathit{\rho }}^2}\mathit{dS}.$

Integrating over the surface of Earth, we have (Eq. 21),

(21) ${\mathit{Q}}_{\mathit{alb},\mathit{i}}={\mathit{\alpha }}_{\mathit{i}}{\mathit{A}}_{\mathit{i}}\mathit{\gamma }{\mathit{I}}_{\mathit{S}}{\tilde{\mathit{F}}}_{\mathit{iE}},$

where the albedo view factor F̃_iE is given by (Eq. 22),

(22) ${\tilde{\mathit{F}}}_{\mathit{iE}}=\int \frac{\cos {\mathit{a}}_{\mathit{E}}\cos {\mathit{a}}_{\mathit{i}}\cos {\mathit{b}}_{\mathit{E}}}{{\mathit{\pi }\mathit{\rho }}^2}\mathit{dS}.$

For the evaluation of the integral in Eq. 22, we use a rotating system of coordinates, {ê'₁, ê'₂, ê'₃}, with origin at the center of the Earth and ê'₃-axis aligned with the position of the satellite r_sat. This system agrees with the GECS, {û₁, û₂, û₃}, when v = 0. In the rotating system the solar vector is given by (Eq. 23),

(23) $\hat{\mathit{s}}=\sin \mathit{v}{{\hat{\mathit{e}}}^{\prime }}_1+\cos \mathit{v}{{\hat{\mathit{e}}}^{\prime }}_3.$

In a way analogous to the computation of the infrared view factors presented before, the cosines in Eq. 22 are written as dot products of vectors expressed as linear combinations of the basis {ê'₁, ê'₂, ê'₃} with coefficients in terms of spherical coordinates (θ, ϕ, r), with ê'₃ acting as the polar axis. Applying this procedure, for the bottom face, Eq. 22 is transformed into Eq. 24,

(24) ${\tilde{\mathit{F}}}_{2\mathit{E}}\left(\mathit{V}\right)=\frac{1}{\mathit{\pi }}{\int }_{{\mathit{\theta }}_{\mathit{m}}}^{{\mathit{\theta }}_{\mathit{M}}}{\int }_{\mathit{I}\varphi \left(\mathit{\theta }\right)}\frac{\left(1-\mathit{h}\cos \mathit{\theta }\right)\left(\cos \mathit{\theta }-\mathit{h}\right)}{{\left(1+{\mathit{h}}^2-2\mathit{h}\cos \mathit{\theta }\right)}^2}\times \left(\sin \mathit{\theta }\cos \varphi \sin \mathit{v}+\cos \mathit{\theta }\cos \mathit{v}\right)\sin \mathit{\theta }\mathit{d}\varphi \mathit{d}\mathit{\theta },$

where the integration limit θ_m and the integration range over ϕ, I_ϕ(θ), are discussed next. The region of integration is the intersection of the spherical region R = {(θ, ϕ) such that θ ≤ θ_M} (θ_M has been defined previously) and the lit hemisphere, bounded by the great circle normal to the solar vector ŝ. For v ∈ [0, π/2 - θ_M] and v ∈ [3π/2 + θ_M, 2π], the region R is contained into the lit hemisphere, so the region of integration is R itself, θ_m = 0, I_ϕ(θ) = [0, 2π]. For v ∈ [π/2 + θ_M, 3π/2 - θ_M], the intersection is empty, so F̃_2E(v) = 0 (in general F̃(v)_iE = 0). For π/2 - θ_M < v < π/2 + θ_M, the part of the region R invaded by the dark hemisphere must be excluded from the integration region, as depicted by Fig. 6.

Figure 6
Region of integration of the albedo view factor (Eq. 24), for π/2 - θ_M < v < π/2 (light gray). The dark gray zone represents the part of the region R (definition in text) invaded by the dark hemisphere, bounded by the curves ϕ_a(θ) = cos^-1 (- cot v cot θ) and ϕ_b(θ) = 2π - ϕ_a(θ). For π/2 - v < θ < θ_M, the range of integration in ϕ, I_ϕ(θ), is the union of the intervals [0, ϕ_a(θ)] and [ϕ_b(θ), 2π].

Similar considerations apply to the cases v ∈ [π/2, π/2 + θ_M] and v [3π/2 -θ_M, 3π/2 + θ_M].

The integral (Eq. 24) was computed numerically with the MATLAB function quad2d for different values of v to obtain the albedo view factor F̃_2E shown in Fig. 7. Figure 4 displays the total albedo heat flux with Q_alb,i = ∑⁶_{i= 1} Q_alb,i computed using Eq. 21.

Figure 7
Albedo view factor F̃_2E(v).

SOLUTION AND TESTING CRITERIA (LINEARIZATION OF THE TEMPERATURE EQUATION)

We now present the temperature predicted by the model. To the satellite was assigned the heat capacity of one kilogram of aluminum, C = 921.6 J/ºC (Gilmore 2002). For the Stefan-Boltzmann constant, we used the 2014 CODATA recommended value σ = 5.670373 × 10^-8 (NIST 2015). The total external heat power Q_tot = Q_S + Q_alb + Q_E is portrayed in Fig. 9a.

We solved numerically Eq. 1 using the function ode45 of MATLAB (relative and absolute tolerances set to 10^-6 and 10^-4, respectively). To avoid the possible adverse effect of the discontinuities of the total external heat flux Q_tot (t) and of its derivative on the solution's accuracy, we divided the interval [0, P] in the following five subintervals, [0, 90º/ω], [90º/ω, (180º - ξ)/ω], [(180º - ξ)/ω, (180º + ξ)/ω)], and [270º/ω, 360º/ω]. The intervals [nP, (n + 1)P], n = 1,2,..., were divided analogously. Using the solution computed at the end of a given subinterval as initial condition for the next subinterval, we prevent the differential equation solver from stepping across discontinuities. The solutions calculated for different initial temperatures are shown in Fig. 8. They approach the same periodic oscillation regardless of the initial value.

Figure 8
Predicted satellite temperature obtained from the numerical solution of Eq. 1 for different initial conditions.

Figure 9
(a) Total external heat flux Q_tot (solid line) and three-term partial sum of the Fourier series expansion of Q_tot, Q₂(t) = Q₀ + a₁ cos(v) + a₂ cos(2v) (dashed line). (b) Numerically computed asymptotic periodic temperature T(t) (solid line) and three-Fourier-mode approximation T₂(t) = T₀ + A₁cos (ωt - ψ₁) + A₂cos (2ωt - ψ₁)(dashed line).

A checking criterion is available for the numerical solution: the long term mean value of T⁴(t) can be known without solving Eq. 1, as follows. Integrating Eq. 1 over one period gives (Eq. 25),

(25) $\mathit{C}{\int }_{\mathit{nP}}^{\left(\mathit{n}+1\right)\mathit{P}}\frac{\mathit{dT}}{\mathit{dt}}\mathit{dt}=-{\mathit{A}}_{\mathit{eff}}\mathit{\sigma }{\int }_{\mathit{nP}}^{\left(\mathit{n}+1\right)\mathit{P}}{\mathit{T}}^4\mathit{dt}+{\int }_{\mathit{nP}}^{\left(\mathit{n}+1\right)}{\mathit{Q}}_{\mathit{tot}}\left(\mathit{t}\right)\mathit{dt}.$

For n sufficiently large, T(t) has reached the asymptotic periodic behavior, T(nP) = T[(n + 1)P], and so the integral on the left-hand side of Eq. 25 vanishes, giving the relation (Eq. 26),

(26) ${\mathit{A}}_{\mathit{eff}}\mathit{\sigma }{\int }_{\mathit{nP}}^{\left(\mathit{n}+1\right)\mathit{P}}{\mathit{T}}^4\mathit{dt}={\int }_0^{\left(\mathit{n}+1\right)\mathit{P}}{\mathit{Q}}_{\mathit{tot}}\left(\mathit{t}\right)\mathit{dt},$

where we have used the periodicity of Q_tot (t) to change the limits of integration on the right-hand side. Defining the mean temperature T₀ as (Eq. 27),

(27) ${\mathit{T}}_0^4=\frac{1}{\mathit{P}}{\int }_{\mathit{nP}}^{\left(\mathit{n}+1\right)\mathit{P}}{\mathit{T}}^4\left(\mathit{t}\right)\mathit{dt},$

from Eq. 26, we obtain (Eq. 28),

(28) ${\mathit{T}}_0^4=\frac{{\mathit{Q}}_0}{{\mathit{A}}_{\mathit{eff}}\mathit{\sigma }}$

where Q₀ is the mean heat flux (Eq. 29),

(29) ${\mathit{Q}}_0=\frac{1}{\mathit{P}}{\int }_0^{\mathit{P}}{\mathit{Q}}_{\mathit{tot}}\left(\mathit{t}\right)\mathit{dt}.$

We computed T₀ by two alternative ways: (i) applying the definition Eq. 26 to the numerically computed T(t) with n = 50; and (ii) using Eq. 27. The two T₀-values found agreed within the first six significant figures, T₀ = 270.210 K ≈ -2.9 ºC, supporting the validity of the numerical solution.

Additional features of the solution T(t) can be understood through linearization of Eq. 1 around T₀. We define the deviation from T₀ as T^* (t) = T(t) - T₀. Substituting T^* into Eq. 1, we have (Eq. 30)

(30) $\mathit{C}\frac{{\mathit{dT}}^{*}}{\mathit{dt}}=-{\mathit{A}}_{\mathit{eff}}\mathit{\sigma }{\left({\mathit{T}}_0+{\mathit{T}}^{*}\right)}^4+{\mathit{Q}}_0+{\mathit{Q}}^{*}\left(\mathit{t}\right),$

where Q^* (t) = Q_tot - Q₀. Assuming |T^*| << T₀ allows the use of the approximation (1 + T^*/T₀)⁴ ≈ 1 + 4T^*/T₀ which, together with Eq. 28, turns Eq. 30 into its linearization (Eq. 31),

(31) $\frac{{\mathit{dT}}^{*}}{\mathit{dt}}=-\frac{{\mathit{T}}^{*}}{\mathit{\tau }}+\frac{{\mathit{Q}}^{*}\left(\mathit{t}\right)}{\mathit{C}},$

where the time constant τ is given by (Eq. 32),

(32) $\mathit{\tau }=\frac{\mathit{C}}{4{\mathit{A}}_{\mathit{eff}}\mathit{\sigma }{\mathit{T}}_0^3}\approx 65.2\min .$

To solve Eq. 30, we rely on the Fourier series expansion of Q^* (t) (Eq. 33),

(33) ${\mathit{Q}}^{*}\left(\mathit{t}\right)=\sum _{\mathit{n}=1}^{\infty }{\mathit{a}}_{\mathit{n}}\cos \left(\mathit{n}\mathit{\omega }\mathit{t}\right),$

where (Eq. 34),

(34) ${\mathit{a}}_{\mathit{n}}=\frac{2}{\mathit{P}}{\int }_0^{\mathit{P}}{\mathit{Q}}^{*}\left(\mathit{t}\right)\cos \left(\mathit{n}\mathit{\omega }\mathit{t}\right)\mathit{dt}.$

Fig. 9a shows the approximation to Q_tot given by the three-term partial sum Q₂(t) = Q₀ + a₁ cos(ωt) + a₂ cos(2ωt). The steady state solution of Eq. 31 can be written as (Eq. 35),

(35) ${\mathit{T}}^{*}\left(\mathit{t}\right)={\Sigma }_{\mathit{n}=1}^{\infty }{\mathit{p}}_{\mathit{n}}\left(\mathit{t}\right)$

where p_n (t) is the steady state solution of (Eq. 36),

(36) $\frac{{\mathit{dp}}_{\mathit{n}}}{\mathit{dt}}=\frac{{\mathit{p}}_{\mathit{n}}}{\mathit{\tau }}+\frac{{\mathit{a}}_{\mathit{n}}}{\mathit{C}}\cos \left(\mathit{n}\mathit{\omega }\mathit{t}\right),$

given by (Eq. 37),

(37) ${\mathit{p}}_{\mathit{n}}\left(\mathit{t}\right)={𝒜}_{\mathit{n}}\cos \left(\mathit{n}\mathit{\omega }\mathit{t}-{\mathit{\psi }}_{\mathit{n}}\right)$

with (Eq. 38),

(38) ${𝒜}_{\mathit{n}}=\frac{{\mathit{a}}_{\mathit{n}}\mathit{\tau }}{\mathit{C}\sqrt{1+{\left(\mathit{n}\mathit{\tau }\mathit{w}\right)}^2}}$

and (Eq. 39),

(39) ${\mathit{\psi }}_{\mathit{n}}={\tan }^{-1}\left(\mathit{n}\mathit{\omega }\mathit{\tau }\right)$

Figure 9b shows the approximation T₂(t) to the solution T(t) given by the first three Fourier modes, . Note how T₂(t) is able to capture the amplitude of the oscillations. The analysis based on linearization, here shown, although approximate, has the advantage of providing simple relations (Eqs. 28, 32 and 38) between the model parameters and attributes of the solution. To appreciate the usefulness of Eq. 38, the reader is invited to reckon a₁ from visual inspection of Fig. 9a, substitute the approximate a₁-value on Eq. 38, and compare the temperature oscillation amplitude so calculated with the actual one shown in Fig. 9b.

We end this section explaining a detail of Fig. 9b: T₂(t) is slightly above T(t). This is due to T₂(t) - T₀ and T^*(t) having different time-averages. 1/P ∫^P₀ [T₂(t) - T₀]dt = 0, while δ = 1/P ∫^P₀ T^*(t)dt < 0. To see this, integrate T^* (t) = T(t) - T₀ over one period to obtain (Eq. 40),

(40) $\mathit{\delta }=\overline{\mathit{T}}-{\mathit{T}}_0$

where T = 1/P ∫^P₀ T^*(t)dt. Jensen's inequality (Rudin 1970) implies that, for non-constant T(t), T < T₀. Hence, δ < 0.

YEARLY VARIATION OF THE RADIATION SCENARIO

In this section, we discuss important changes in the solar and albedo heat fluxes that take place over a time scale of several months for general satellite orbits. They are the result of variations in the angle between the satellite orbital plane and the solar vector ŝ (the minimum of the angles between all orbital positions r_sat and ŝ), known as the β angle. To see how β affects the solar heat flux Q_S falling on the satellite, contrast the time-periodic Q_S found for the orbit in the previous section (β = 0º) to the time-constant Q_S expected for velocity-nadir attitude and a circular orbit on a plane perpendicular to ŝ (β = 90º). The orientation of the orbital plane is given by the normal vector lˆ computed as the direction of the orbital angular momentum l = r_sat × v_sat, with v_sat the satellite velocity, lˆ = l/||l||. Moreover, from the definition of β, we have (Eq. 41)

(41) $\cos \left(90°-\mathit{\beta }\right)=\hat{\mathit{s}}·\hat{\mathit{l}}$

During the year, both lˆ and ŝ rotate with respect to the GECS. We will express both vectors on the {û₁, û₂, û₃} basis, and, using Eq. 41, obtain β. But first, we introduce the standard way of presenting a general spacecraft orbit.

KEPLERIAN ORBITAL ELEMENTS

On a first approximation, the satellite orbit is described by a solution of Kepler problem (Vallado 1997), namely, an ellipse with one of its foci at the origin of the GECS. The orbit shape is described by the semimajor axis a and the eccentricity e = c/a, where c is half the distance between the foci. As e → 0, the ellipse approaches a circle (e = 0).

lˆ is determined by the inclination i and the right ascension of the ascending node (RAAN) Ω through (Eq. 42) (see Fig. 10),

(42) $\hat{\mathit{l}}=\sin \mathit{i}\sin \Omega {\hat{\mathit{u}}}_1-\sin \mathit{i}\cos \Omega {\hat{\mathit{u}}}_2+\cos \mathit{i}{\hat{\mathit{u}}}_3.$

Figure 10
Four of the six Keplerian orbital elements. i: inclination; Ω: right ascension of the ascending node; g: argument of perigee; and v: true anomaly.

The ascending node is the orbital point on the equatorial plane at which the satellite moves from south to north. The orientation of the semimajor axis within the orbital plane is given by the argument of perigee g, defined as the angle between the position vectors of the ascending node and the periapsis, the location at which the satellite is closest to the origin (Fig. 10). Finally, the satellite position on the ellipse is set by the true anomaly v, measured from the periapsis. The Keplerian orbital elements a, , Ω, i, g, and v completely specify a solution of Kepler problem (Vallado 1997; Bate et al. 1971).

The orbital elements of Libertad 1 are available online from Space-Track (a service of the United States Department of Defense (JFCCS 2015)) in the two line element set (TLE) format. Elements , Ω, i, and g can be retrieved straightforwardly from a TLE. On the other hand, a needs to be computed from the mean motion (mean angular velocity) using Eq. 9 with R substituted by a. Similarly, ν is obtained from the mean anomaly (Vallado 1997).

Table 3 shows the orbital elements of Libertad 1 on April 18, 2007, 21 h 54 min 52.4 s Coordinated Universal Time (according to Portilla (2012), the orbital elements of Libertad 1 disclosed for prior times are possibly inaccurate). Note the low value of the eccentricity, which justifies our assumption of a circular orbit for Libertad 2. Observe also the closeness between a and the value of R assumed in the previous section.

Table 3
Orbital elements of satellite Libertad 1 on April 18, 2007, 21 h 54 min 52.4 s (JFCCS 2015). a: semi-major axis; e: eccentricity; v: true anomaly; i: inclination; Ω: Right Ascension of the Ascending Node (RAAN); g: argument of perigee (see Fig. 10).

a (km)	e	i (º)	Ω (º)	g (º)	v (º)
7097.87	0.0102	98.08	184.09	207.97	154.07

The orbital elements of a spacecraft on a LEO do not stay constant. Of relevance to the thermal analysis is the change of the RAAN Ω, known as nodal regression and determined next.

NODAL REGRESSION

The assumptions of Kepler problem are satisfied only when the planet has a perfectly spherical mass distribution. In actual fact, Earth presents a protrusion around the Equator, which can be thought of as a ring added to a sphere. This ring applies a gravitational torque that rotates the satellite orbital angular momentum, causing a drift of Ω (Fortescue et al. 2003) (Eq. 43),

(43) $\Omega ={\Omega }_0-\dot{\Omega }\mathit{t}$

with a constant rate of nodal regression Ω̇. The protrusion is measured by the coefficient J₂ = 1.0826 × 10^-3 of the expansion in spherical harmonics of the gravitational potential generated by Earth (Fortescue et al. 2003; Vallado 1997). J₂ determines the rate of nodal regression through the relation (Fortescue et al. 2003; Vallado 1997) (Eq. 44),

(44) $\dot{\Omega }=\frac{3}{2}\frac{{\mathit{J}}_2{\mathit{R}}_{\mathit{E}}^2}{{\mathit{R}}^2}\mathit{\omega }\cos \mathit{i}$

where terms O (J²₂) were neglected and a circular orbit of radius R and angular velocity ω, given by Eq. 9, was assumed.

VARIATION OF THE SOLAR VECTOR Sˆ and of β

For a satellite on a LEO, the solar vector ŝ is essentially the unit vector in the direction of the Sun's position with respect to the GECS. Due to the tilt of the Earth's rotational axis relative to the plane of its orbit, from the geocentric point of view, the Sun traverses an orbit on a plane (the ecliptic) inclined roughly 23º with respect to the Equatorial plane (Fig. 11). In fact, by definition, the û₁ vector points toward the ascending node of the Sun's orbit, or point of Aries () (Fig. 11). ŝ is determined by the Sun's longitude Ω_S and declination δ_S with respect to the GECS (Fig. 11) (Eq. 45),

(45) $\hat{\mathit{s}}=\cos {\mathit{\delta }}_{\mathit{s}}\cos {\Omega }_{\mathit{s}}{\hat{\mathit{u}}}_1+\cos {\mathit{\delta }}_{\mathit{s}}\sin {\Omega }_{\mathit{s}}{\hat{\mathit{u}}}_2+\sin {\mathit{\delta }}_{\mathit{s}}{\hat{\mathit{u}}}_3.$

Figure 11
Solar orbit on the Ecliptic in relation to the GECS. The figure also depicts the obliquity of the Ecliptic E, the ecliptic longitude of the Sun λ_ecl, the Sun's longitude Ω_S, and the Sun's declination δ_S.

The evolution of Ω_S and δ_S can be computed by algorithms available in Vallado (1997), Nautical Almanac Office (2008), and Fitzpatrick (2010) or it can be obtained from online resources (Jubier 2004), such as NASA's HORIZONS system (JPL 2015a). To produce an initial β angle similar to that of Libertad 1, we assume Libertad 2 will be deployed into orbit on April 19, 2019 at 0 h 0 min 0 s UTC. Figure 12 shows Ω_S and δ_S from the time of insertion into orbit, computed using the algorithm in Fitzpatrick (2010) (presented in Appendix A). Our values agree with HORIZONS apparent ephemeris within 0.01º.

Figure 12
Right ascension, Ω_S, and declination, δ_S, of the Sun in the days following the hypothetical time of deployment into orbit of Libertad 2, April 19, 2019 0 h 0 min 0 s UTC (t = 0).

For the computation of β, we assume that Libertad 2 will be inserted into a circular orbit with the same radius R used in the previous section (7110 km), and values of inclination and RAAN similar to those in Table 3, i = 98º, and Ω₀ = 184º. For these, Eq. 44 yields Ω̇ ≈ -0.948º/day.

Figure 13 displays β(t) obtained from Eq. 41 with Ω(t) given by Eq. 43, i assumed constant (Portilla 2012), and Ω_S(t) and δ_S(t) calculated as described previously.

Figure 13
β angle in the days following the time of Libertad 2 insertion into orbit.

YEARLY VARIATION OF THE TEMPERATURE EXTREMES

Here, we determine the changes of the minima and maxima (extremes) of the temperature oscillation in response to the variation of the β angle throughout the year. A conceivable way of accomplishing this would be to compute the year-long solution T(t) of Eq. 1 with Q_tot(t) now including the dependence on β (in addition to the true anomaly v), Q_tot(t) = Q_tot [v(t), β(t)]. However, a more efficient and insightful approach is possible. Since β changes much more slowly (over intervals of several months) than ν (over a period of about 100 min), for an interval of a few orbital periods, the solution for Q_tot [v(t), β(t)] can be approximated by the solution for Q_tot [v(t), β] with β fixed in time at some representative value. Let T_β (t) denote the long term periodic solution for constant β. We can then approximate the year-long solution T(t) using the solutions T_β (t) by means of T(t) ≈ T_β(t)(t); the accuracy of this approximation, which we did not appraise, might be estimated by the perturbative technique of multiple scales (Hinch 1991).

Next, we look for a convenient way of computing T_β(t). To this end, it is pertinent noting that circular orbits with identical radius R but different orientation (i, Ω) are equivalent for the purposes of thermal analysis if they share the same β-value. For velocity-nadir attitude, orbits with identical β are subject to the same variation of the solar heat flux, Q_S(v - ψ_S), with possible differences only in the phase angle ψ_S. Furthermore, under the model of uniform time-constant albedo γ, they also present equal profiles of albedo heat flux change, Q_alb(v - ψ_alb), again, except for eventual differences in ψ_alb.

For polar orbits, those with i = 90º, the simple relationship β = Ω holds when ŝ = û₁ (Fig. 14). Considering the straightforward geometrical configuration involved, we use polar orbits to compute T_β(t). Note that, because any circular orbit is equivalent to a polar one, all possible radiation scenarios are encompassed when considering polar orbits within the range 0 ≤ Ω ≤ π/2.

Figure 14
Polar orbits (thick lines) seen from the û₃ direction. The dark-gray region represents the eclipse zone of Earth (light-gray disk). Values of β = Ω displayed: (i) 0º, (ii) 30º, and (iii) 90º.

For velocity-nadir attitude and polar orbits, the computation of the albedo view factor F_iE through the integration in Eq. 22 becomes impractically complicated for Ω-values other than 0º or 90º. An alternative for F_iE calculation is the Monte Carlo ray tracing method (MCRTM), which follows the paths of rays with random directions, send from random points on the surfaces exchanging radiation (Walker et al. 2010). To establish T_β (t), we used the thermal analysis software ESATAN-TMS release 5 (ITP Engines UK Ltd. 2013a; 2013b), which implements the MCRTM and solves the associated differential equations for the temperatures. In operating ESATAN, we approximated the single-node model by a multiple-node model as explained below.

USE OF ESATAN FOR TEMPERATURE COMPUTATION

The calculation of a satellite's temperature using ESATAN proceeds through three stages (ITP Engines UK Ltd. 2013a; 2013b): (i) the geometrical model, (ii) the radiative case, and (iii) the analysis case. During the first stage, the satellite's geometrical model is assembled from shells with predetermined shapes. For our case, we built a parallelepiped-shaped empty box using six rectangular shells of thickness 2 mm. The inwardly looking shell faces were declared thermally inactive, while the outer faces were given thermo-optical properties (ε_avg, α_avg) or (ε_u, α_u). To each shell we assigned a node (thus, an individual temperature), producing a six-node satellite model. To approximate a single-node model the shells' bulk material was given an unrealistically high thermal conductivity of 10⁶ W/(m⋅K). The density and specific heat were set to 3.291429 × 10⁶ kg/m³ and 2 J/(kg K), respectively, so that the box total heat capacity was equal to the value of C used previously.

In the radiative case, the thermal environment, orbit, and satellite attitude are defined. We gave the solar irradiance I_S, the albedo γ, and the Earth radiosity I_E the values employed in the previous section. I_S was set using the solar constant override option, while I_E was assigned indirectly through the planet emissivity ε_E and the planet temperature T_E following the relationship (Eq. 46),

(46) ${\mathit{I}}_{\mathit{E}}={\epsilon }_{\mathit{E}}\mathit{\sigma }{\mathit{T}}_{\mathit{E}}^4.$

In particular, we used ε_E = 1 and T_E = 247.567 K. The Sun's right ascension and declination were set to Ω_S = δ_S = 0, yielding ŝ = û₁.

The satellite's altitudes of apogee and perigee were both assigned values of 732 km, thus producing a circular orbit with the same radius R used in the previous section. We set the inclination to i = 90º and used 52 different values of 0 ≤ Ω ≤ 90º. For all orbits, the satellite was put in velocity-nadir attitude. The planet's gravitational parameter µ was set indirectly using the acceleration of gravity g_a = 9.79871 m/s², which, through the relationship (Eq. 47),

(47) ${\mathit{g}}_{\mathit{a}}=\frac{\mathit{\mu }}{{\mathit{R}}_{\mathit{E}}^2}$

produces the value of µ used in the previous section. We directed ESATAN to calculate the heat fluxes for a hundred different values of the true anomaly v (orbital positions) with default parameters used for the MCRTM.

In the analysis case, the Stefan-Boltzmann constant was set, through the variable STEFAN, to the value referred previously. We chose the routine SLCRNC for the numerical solution of the temperature differential equations with time step of 0.01 s. The six nodal temperatures were computed for an interval of about 20 orbital periods, long enough for the temperature oscillations to reach the periodic regime.

Each Ω-value required an independent ESATAN temperature simulation. For Ω = 0, the simulation parameters were set up using ESATAN's graphical user interface ("workbench"), which saves logfiles with the configuration instructions. For the remaining Ω-values, we found it more expeditious to use a Cygwin (Cygwin 2017) shell script that updated Ω in the logfiles and fed them to ESATAN processes running in batch mode.

TEMPERATURE EXTREMES

For Ω = 0, the six nodal temperatures computed with ESATAN at no time differed by more than 0.03 ºC from the nodal average, indicating close approximation to a single node model. On the periodic regime, ESATAN's nodal average temperature was at all times roughly 0.2 ºC above the single node model temperature computed previously (Fig. 15). This small difference is attributed to a less exact estimation of the heat fluxes by the MCRTM (see Table 2) compared to the integrations in the previous section.

Figure 15
Comparison of temperatures computed with the procedure of the third section (solid gray line) and with ESATAN (black dashed line). Both computations start at T (0) = 0 ºC.

Let the maximum and minimum temperatures be denoted by T_max(β) = max_t T_β(t) and T_min(β) = min_t T_β(t). Figure 16 shows T_max(β) and T_min(β) computed with ESATAN for three different values of the emissivity ε_u, keeping all other parameters values.

Figure 16
Maximum, T_max, (upper branch) and minimum, T_min, (lower branch) of the temperature oscillation as functions of the β angle, for ε_u = 0.05(solid line), 0.1 (dashed line), and 0.5 (dot-dashed line).

In the three cases, T_max(β) increases from β = 0 to β = ξ (see Eq. 4 for the definition of ξ), where it reaches a maximum with a kink, to then decrease. Similar behavior is observed for T_min(β). In this context, ξ has a special meaning: when β < ξ, orbits have an eclipse zone, while when β ≥ ξ, they do not. For β = 90º, the external heat fluxes are constant, and so T_max(β) = T_min(β). Additionally, for fixed β, the higher ε_u is, the lower both T_max(β) and T_min(β) are. Figure 17 shows the temperature extremes T_max and T_min as functions of the time t since deployment into orbit. T_max(t) is obtained by composition of the function T_max(β) (Fig. 16) with the function β(t) (see Fig. 13), T_max(t) = T_max[β(t)].

Figure 17
Maximum, T_max, (upper branch) and minimum, T_min, (lower branch) of the temperature oscillation as functions of the time t since deployment into orbit, for ε_u = 0.05 (solid line), 0.1 (dashed line), and 0.5 (dot-dashed line).

For the single node model, the acceptable satellite temperature range is the intersection of the operating temperature intervals of all the components. Hence, according to Table 1, for Libertad 2, this range is 0 ºC ≤ T ≤ 60 ºC. For ε_u = 0.05, T_max(t) is within the acceptable range but T_min(t) can be below the lower limit by as much as 9 ºC. Similar behavior is observed for ε_u = 0.1, with the lower limit being exceeded by a higher amount, 13 ºC, in the worst case. ε_u = 0.5, in contrast, corresponds to T_min(t) substantially below the lower limit, implying that an external cover with this emissivity must be ruled out.

The problem of T_min(t) exceeding the lower limit for ε_u = 0.05 and 0.1 could be solved in at least two ways. First, an orbit exposed to higher β angles (Fig. 16) could be chosen. Second, the satellite could be wrapped on a thermally insulating material, such as a multilayer insulation blanket (MLIB) (Gilmore 2002). An MLIB can withstand a high temperature difference between the elements enclosed within it, on one side, and its outer cover, on the other. So while the outer cover's temperature may fluctuate wildly, due to direct exposition to the external heat fluxes, the internal components' temperature would oscillate with much smaller amplitude. As a result, T_min would be drawn up, closer to the oscillation average, which is determined by the outer cover's absorptivity and emissivity. A quantitative prediction of this effect, however, requires a model with at least two nodes (outer cover and internal components, respectively), which is beyond the scope of the present work.

CONCLUDING REMARKS

We carried out a preliminary thermal analysis of a three-unit CubeSat with the features envisioned for satellite Libertad 2, under development at Universidad Sergio Arboleda (Bogota, Colombia).

In particular, we thoroughly showed the computation of the most significant heat fluxes affecting the CubeSat's temperature on a particular low-Earth, circular, polar orbit. These are the fluxes due to the direct solar radiation, the solar radiation reflected on Earth (albedo), and the infrared radiation released by Earth. The time-periodic fluxes were included into a differential equation model describing the evolution of a single temperature, assumed shared by all satellite components (single-node model). We found that the model solutions approach a unique periodic temperature oscillation, independently of the initial condition.

Additionally, we considered circular orbits with arbitrary orientation. For these, the heat fluxes are determined by the angle between the orbital plane and the direction of the sunrays, or β angle, making any orbit equivalent to a polar one, for thermal analysis purposes. With the help of the thermal analysis software ESATAN-TMS, we calculated the temperature T_β(t) for polar orbits with 0º ≤ β ≤ 90º. Based on data from the previous CubeSat mission Libertad 1, we established a hypothetical evolution of β(t) for a mission life of a year, tracking the solar position on the ecliptic and including the nodal regression of the satellite orbit. Combining T_β(t) and β(t), we determined the minima and maxima (extremes) of the temperature oscillation during the mission life. For plausible values of the relevant satellite parameters, the temperature extremes were mostly within the operating range of the most sensitive satellite component, 0º ≤ T ≤ 60º, suggesting viability of the mission.

The thermal model considered can be improved in ways that build upon the work presented here. The most important enhancement would consist of assigning distinct temperatures to different satellite parts, producing a multiple node model able to handle heat exchange between satellite components (as described in Karam (1998), for instance). For this model, the irradiances associated to the external heat fluxes would be identical to the ones presented here, as a result of they being determined solely by the external geometrical configuration. Since the solar cells' effective absorptivity is influenced by the electrical power drained from them, it would be pertinent to relate the thermal analysis to the satellite's electrical energy consumption (a preliminary analysis of electrical energy utilization for Libertad 2 is provided in Sanchez-Sanjuan et al. 2016). Coupling of the thermal and electrical analyses would be also necessary when studying the heat release from electrical circuits, be it unintended or deliberate, as in the use of heaters for active thermal control. As priorly mentioned, the multiple node model would allow the representation of passive control mechanisms (thermal insulation) or additional active control methods (such as louvers) (Gilmore 2002).

The temperature analysis could additionally include the fluctuations of the thermal environment, represented by the infrared radiosity I_E and the albedo γ. This would involve solving a statistical ensemble of thermal models, each with a particular realization of uniform and time invariant I_E and γ drawn from the joint probability distribution of these two parameters found by the ERBE (Anderson et al. 2001). Of special interest would be the minima and maxima of I_E and γ, as they would supply the temperature extremes. Alternatively, if the desired degree of accuracy calls for it, the parameters could be introduced in the model as stochastic signals I_E(t) and γ(t).

Other phenomena to consider in the model are the following: the yearly variation of the solar irradiance I_S depending upon the Earth-Sun distance (Anderson et al. 2001), the correction of the albedo coefficient of reflection γ following the change of the orbit's β angle (Anderson et al. 2001), and the possible effect of the deposition of space dust on the external surfaces' emissivity and absorptivity (Gilmore 2002).

ACKNOWLEDGMENTS

This work was partially supported by the Young Research Award (YAV) of Universidad Sergio Arboleda, Colombia. We thank ITP Engines UK for kindly sponsoring the Misión Libertad 2 project with the licenses for the analysis and simulation software ESATAN-TMS.

APPENDIX - A COMPUTATION OF THE SUN'S LONGITUDE AND DECLINATION

Here, we explain the procedure followed for calculation of the Sun's longitude and declination relative to the GECS, Ω_S and δ_S, respectively. Basic geometry shows that Ω_S and δ_S are related to the Sun's longitude on the ecliptic plane λ_ecl (depicted in Fig. 11 of the main text) by (Eqs. 1 and 2),

(1) ${\Omega }_{\mathit{s}}={\tan }^{-1}\left(\cos \mathit{E}\tan {\mathit{\lambda }}_{\mathit{ecl}}\right),$

(2) ${\mathit{\delta }}_{\mathit{S}}={\sin }^{-1}\left(\sin \mathit{E}\sin {\mathit{\lambda }}_{\mathit{ecl}}\right),$

where E = 23.439º (Vallado 1997) is the ecliptic obliquity (see Fig. 11 of the main text) and Ω_S is in the same quadrant as λ_ecl.

To determine λ_ecl, we used the algorithm presented in Fitzpatrick (2010) (similar algorithms can be found in Vallado (1997) and in Nautical Almanac Office (2008)). As shown by Fig. 1, λ_ecl is related to the Sun's true anomaly v_S and the Sun's orbit argument of periapsis g_S by (Eq. 3),

(3) ${\mathit{\lambda }}_{\mathit{ecl}}={\mathit{g}}_{\mathit{S}}+{\mathit{v}}_{\mathit{S}}.$

Figure 1
The Sun's orbit seen from a direction perpendicular to it. Definitions of v_S, g_S, and λ_ecl shown. The Earth sits at one of the ellipse's foci; the symbol × marks the other one. The orbit's eccentricity e_E (actual value shown in Table 4) has been exaggerated for clarity.

v_S can be found with the help of the Sun's mean anomaly M, which changes in time at a constant rate Ṁ (value presented in Table 1) (Eq. 4),

(4) $\mathit{M}={\mathit{M}}_{\mathit{ep}}+\dot{\mathit{M}}\left({\mathit{t}}^{\prime }-{{\mathit{t}}^{\prime }}_{\mathit{ep}}\right),$

where t' is the time expressed as the Julian Date (Vallado 1997) and M_ep (Table 1) is the mean anomaly at the reference time t'_ep = 2451545.0 (Julian Date), corresponding to the 1^st January 2000 at 12:00 UT, known as the J2000 epoch.

To second order in the Sun's orbit eccentricity e_E (Table 1), the relationship between v_S and M is (Eq. 5),

(5) ${\mathit{v}}_{\mathit{S}}=\mathit{M}+2{\mathit{e}}_{\mathit{E}}\sin \mathit{M}+\frac{5}{4}{\mathit{e}}_{\mathit{E}}^2\sin 2\mathit{M}.$

Substitution of Eq. 5 into Eq. 3 leads to (Eq. 6),

(6) ${\mathit{\lambda }}_{\mathit{ecl}}=\mathit{L}+2{\mathit{e}}_{\mathit{E}}\sin \mathit{M}+\frac{5}{4}{\mathit{e}}_{\mathit{E}}^2\sin 2\mathit{M},$

where L is the Sun's mean ecliptic longitude, L = g_S + M, that changes linearly with time (Eq. 7),

(7) $\mathit{L}={\mathit{L}}_{\mathit{ep}}+\dot{\mathit{L}}\left({\mathit{t}}^{\prime }-{{\mathit{t}}^{\prime }}_{\mathit{ep}}\right)$

with values of L_ep and L̇ displayed in Table 1.

Table 1
Parameter values used in the calculation of the Sun's ecliptic longitude λ_ecl (Fitzpatrick 2010).

Parameter	Value
Mep	357.588º
M.	0.98560025º/day
Lep	280.458º
L.	0.98564735º/day
eE	0.016711

t' is related to the time elapsed since orbital insertion t (in units of day) used in Eq. 43 (main text) by t' = t'₀ + t, where t'₀ is the Julian Date of the time of orbital insertion, determined from the Gregorian calendar date (April 19, 2019, 0 h 0 min 0 s) with the function juliandate of Matlab, an algorithm for the calculation of the Julian Date from the Gregorian calendar date is explained in Vallado (1997). Internet resources for the same purpose are available in JPL (2015b) and Jubier (2004).

REFERENCES

• Anderson BJ, Justus CG, Batts GW (2001) Guidelines for the selection of near-Earth thermal environment parameters for spacecraft design. (TM-2001-211221). NASA Technical Report.
• Anderson BJ, Smith, RB (1994) Natural orbital environment definition guidelines for use in aerospace vehicle development. (TM-4527). NASA Technical Report.
• AZUR SPACE Solar Power GmbH (2014) Personal communication.
• AZUR SPACE Solar Power GmbH (2016) 28% Triple Junction GaAs Solar Cell 3G28C datasheet; [accessed 2018 May 9]. http://www.azurspace.com/images/products/0002490-00-04_DB_GA_175.pdf
  » http://www.azurspace.com/images/products/0002490-00-04_DB_GA_175.pdf
• Bate R, Mueller D, White J (1971) Fundamentals of astrodynamics. New York: Dover Publications.
• Buchen E (2014) SpaceWorks' 2014 nano/microsatellite market assesssment. Presented at: AIAA/USU Conference on Small Satellites; Logan, USA.
• Bulut M, Sozbir N (2015) Analytical investigation of a nanosatellite panel surface temperatures for different altitudes and panel combinations. Applied Thermal Engineering 75:1076-1083. doi: 10.1016/j.applthermaleng.2014.10.059
  » https://doi.org/10.1016/j.applthermaleng.2014.10.059
• Corpino S, Caldera M, Nichele F, Masoero M, Viola N (2015) Thermal design and analysis of a nanosatellite in low Earth orbit. Acta Astronautica 115:247-261. doi: 10.1016/j.actaastro.2015.05.012
  » https://doi.org/10.1016/j.actaastro.2015.05.012
• Cygwin (2017) Cygwin; [accessed 2017 May 1]. https://www.cygwin.com
  » https://www.cygwin.com
• Dinh D (2012) Thermal modeling of nanosat (Master's Thesis). San Jose: San Jose State University.
• Escobar E, Diaz M, Zagal JC (2016) Evolutionary design of a satellite thermal control system: real experiments for a CubeSat mission. Applied Thermal Engineering 105:490-500. doi: 10.1016/j.applthermaleng.2016.03.024
  » https://doi.org/10.1016/j.applthermaleng.2016.03.024
• Fish CS, Swenson CM, Crowley G, Barjatya A, Neilsen T, Gunther J, Azeem I, Pilinski M, Wilder R, Allen D, et al. (2014). Design, development, implementation, and on-orbit performance of the dynamic ionosphere CubeSat experiment mission. Space Science Reviews 181(1-4):61-120. doi: 10.1007/s11214-014-0034-x
  » https://doi.org/10.1007/s11214-014-0034-x
• Fitz CD, Lukasik SJ, Mayer EA, Simon DR (1963). Albedo and planet radiation intercepted by an Earth satellite. (TDR-63-92). AEDC Technical Report.
• Fitzpatrick R (2010) Determination of ecliptic longitude; [accessed 2015 June 16]. http://farside.ph.utexas.edu/Books/Syntaxis/Almagest/node34.html
  » http://farside.ph.utexas.edu/Books/Syntaxis/Almagest/node34.html
• Fortescue P, Stark J, Swinerd G, editors (2003) Spacecraft systems engineering. 3rd ed. Chichester: John Wiley & Sons.
• Garzon MM (2012) Development and analysis of the thermal design for the OSIRIS-3U CubeSat (Master's Thesis). Pennsylvania: Pennsylvania State University.
• Gilmore DG (2002) Spacecraft thermal control handbook: fundamental technologies. El Segundo: The Aerospace Press.
• Hinch EJ (1991) Perturbation Methods. Cambridge: Cambridge University Press.
• ITP Engines UK Ltd. (2013a) ESATAN-TMS Workbench Getting Started Guide. Whetstone: ITP Engines UK Ltd.
• ITP Engines UK Ltd. (2013b) ESATAN-TMSWorkbench User Manual. Whetstone: ITP Engines UK Ltd.
• Jacques L (2009) Thermal design of the OUFTI-1 nanosatellite (Master's Thesis). Liège: Université de Liège.
• [JFCCS] Joint Functional Component Command for Space (2015) Space-Track; [accessed 2017 May 1]. https://www.space-track.org
  » https://www.space-track.org
• Joya R (2007) Libertad 1, primer satélite colombiano en el espacio. Revista Innovación y Ciencia 14:16-23.
• [JPL] Jet Propulsion Laboratory (2015a) JPL Solar System Dynamics - horizons system; [accessed 2015 June 16]. http://ssd.jpl.nasa.gov/?horizons
  » http://ssd.jpl.nasa.gov/?horizons
• [JPL] Jet Propulsion Laboratory (2015b) JPL Solar System Dynamics - time conversion tool; [accessed 2015 June 17]. http://ssd.jpl.nasa.gov/tc.cgi
  » http://ssd.jpl.nasa.gov/tc.cgi
• Jubier XM (2004) Solar, Lunar, and Planets Ephemerides; [accessed 2015 June 17]. http://xjubier.free.fr/en/site_pages/astronomy/ephemerides.html
  » http://xjubier.free.fr/en/site_pages/astronomy/ephemerides.html
• Kang S-J, Oh H-U (2016) On-orbit thermal design and validation of 1U standardized CubeSat of STEP CubeLab. International Journal of Aerospace Engineering 2016:4213189. doi: 10.1155/2016/4213189
  » https://doi.org/10.1155/2016/4213189
• Karam RD (1998) Satellite thermal control for systems engineers. Reston: American Institute of Aeronautics and Astronautics. doi: 10.2514/4.866524
  » https://doi.org/10.2514/4.866524
• Mason JP, Lamprecht B, Woods TN, Downs C (2018) CubeSat on-orbit temperature comparison to thermal-balance-tuned-model predictions. Journal of Thermophysics and Heat Transfer 32:237-255. doi: 10.2514/1.T5169
  » https://doi.org/10.2514/1.T5169
• NASA (2015) Moderate Resolution Imaging Spectroradiometer; [accessed 06-July-2015]. http://modis.gsfc.nasa.gov/
  » http://modis.gsfc.nasa.gov/
• NASA (2014) Space Science Data Coordinated Archive Master Catalog Spacecraft Query; [accessed 2015 July 27]. http://nssdc.gsfc.nasa. gov/nmc/SpacecraftQuery.jsp
  » http://nssdc.gsfc.nasa.gov/nmc/SpacecraftQuery.jsp
• [NIST] National Institute of Standards and Technology (2015) Fundamental Physical constants; [accessed 2015 June 16]. http://physics.nist.gov/cuu/Constants/index.html
  » http://physics.nist.gov/cuu/Constants/index.html
• Palmer JM, Grant BG (2010) The art of radiometry. Bellingham: SPIE Press.
• Portilla JG (2012) La órbita del satélite Libertad 1. Rev Acad Colomb Cienc 36:491-500.
• Rudin W (1970) Real and complex analysis. New York: McGraw-Hill.
• Sanchez-Sanjuan S, Gonzalez-Llorente J, Hurtado-Velasco R (2016) Comparison of the incident solar energy and battery storage in a 3U CubeSat satellite for different orientation scenarios. J Aerosp Technol Manag 8(1):91-102. doi: 10.5028/jatm.v8i1.531
  » https://doi.org/10.5028/jatm.v8i1.531
• Schaepman-Strub G, Schaepman ME, Painter TH, Dangel S, Martonchik JV (2006) Reflectance quantities in optical remote sensing-definitions and case studies. Remote Sensing of Environment, 103(1):27-42. doi: 10.1016/j.rse.2006.03.002
  » https://doi.org/10.1016/j.rse.2006.03.002
• Shinde P, Fernandez A, Tansel I, Tosunoglu S (2017) Active thermal control system for CubeSat. Presented at: 30th Florida Conference on Recent Advances in Robotics; Boca Raton, USA.
• Strahler AH, Muller J-P, MODIS Science Team Members (1999) MODIS BRDF/Albedo Product: Algorithm Theoretical Basis Document Version 5.0; [accessed 2015 July 6]. http://modis.gsfc.nasa.gov/data/atbd/atbd_mod09.pdf
  » http://modis.gsfc.nasa.gov/data/atbd/atbd_mod09.pdf
• Swartwout M (2013) The first one hundred cubesats: a statistical look. Journal of Small Sattelites 2(2):213-233.
• The CubeSat Program (2015) CubeSat Design Specification Rev. 13; [accessed 2015 July 27]. http://www.cubesat.org/resources/
  » http://www.cubesat.org/resources/
• Nautical Almanac Office (2008) The Astronomical Almanac for the Year 2015. Washington: TSO.
• Vallado DA (1997) Fundamentals of astrodynamics and applications. New York: McGraw-Hill.
• Walker T, Xue S-C, Barton GW (2010) Numerical determination of radiative view factors using ray tracing. Journal of Heat Transfer 132(7):072702. doi: 10.1115/1.4000974
  » https://doi.org/10.1115/1.4000974
• Wolfe WL (1998) Introduction to Radiometry. Bellingham: SPIE Press.

Edited by

Section Editor: Marcia Mantelli

Publication Dates

• Publication in this collection
  08 Nov 2018
• Date of issue
  2018

History

• Received
  17 Jan 2018
• Accepted
  24 May 2018