Modeling of Resident Space Object Light Curves with Blender Software

A key aspect of orbit determination and propagation method improvements lies in gaining better insights into the characteristics of each object, such as its shape, mass, and operational status. While some of this information may not directly avert collisions, it plays a crucial role in identifying potential threats and advancing our overall understanding of objects in near-Earth space. Estimation of some of these parameters can be performed by photometry research of resident space object shape.

For spacecraft, regardless of their orbital location, ground-based optical telescopes assigned to space surveillance often face limitations in resolving these objects. Consequently, they appear as spatially unresolved point sources, making it difficult to differentiate between them. The process of characterization involves delving deeper into object’s properties to distinguish it from others. For resident space objects, these characteristics include orientation, rate of change of orientation, physical shape, and surface material composition. Previous research on determining spacecraft characteristics using unresolved observations has predominantly focused on analyzing photometric light curves (Luu et al., 2003; Scott et al., 2008; Somers, 2011; Be’dard, 2013; Jolley, 2014). These light curves represent the brightness of the spacecraft as a function of time and have been successfully used to determine the spin rate and axis of uncontrolled spacecraft, as well as differentiate co-located satellites. Notably, the shape of RSOs holds significant importance as it allows for inferring their size, mass, and potentially their operational status.

Photometric measurements have gained popularity in both the space situational awareness (SSA) community and the asteroid community in recent decades for estimating the size and shape (Jah et al., 2007; Wetterer, 2009; Linares et al., 2010). Among the techniques employed, LC inversion (LCI) has proven to be a valuable tool, particularly in the asteroid community for generating accurate shapes of numerous asteroids (Kaasalainen et al., 2001a; Kaasalainen et al., 2001b; Kaasalainen et al., 2002; Torppa et al., 2003). As a good example, there is a “Voxteroid: Generating asteroid light curves using Threejs and Blender” project (Cullen and Paterson, 2021). This project is designed as a web interface for generating the light curves of various three-dimensional shapes including existing asteroids, comets, and other objects using the Threejs engine (Cabello, 2024) and Blender (Community B.O., 2024). After choosing the rotation axis, the sun’s angle, and the asteroid’s rotation speed, the user can generate its corresponding light curve by sending a request to a Blender instance. The project provides users with an accessible interface to find light curves for various objects and allows them to upload their own light curves.

However, the application of LCI to satellites has been limited, often assuming known attitude states and focusing on simplistic shapes like polyhedrons or cube-shaped objects (Fulcoly et al., 2012). Some studies in the SSA community have explored the use of various Kalman filters, including multiple-model adaptive estimation (MMAE), to estimate polyhedron and cube-shaped objects as well as attitude (Linares et al., 2010; Holzinger et al., 2012). While promising, these methods tend to overlook the more complex shapes typically associated with satellites.

When it comes to measuring the reflectance of a spacecraft, there are a few things to consider. The materials used in its construction all play a role in determining its reflectance, with the amount of each material affecting the overall result. Rather than relying on a lab test, which can be difficult to carry out, a computer-aided design (CAD) model of the spacecraft can be used along with samples of its composite materials to create a simulated model. This method is much simpler, as the samples can be easily manipulated and analyzed.

Since it is not always possible to physically examine the composition of a spacecraft, accurate modeling is crucial. By creating synthetic photometric observations, we can better interpret photometric measurements and learn more about unresolved objects in space.

This paper aims to assess the feasibility of leveraging Blender, a software primarily used for animation and modeling, to enhance our understanding of the shape and orientation of objects under study. By using Blender, the plausibility of our perceptions regarding these properties can be better evaluated and potentially offer valuable insights for future space situational awareness efforts.

To reproduce RSO’s pass parameters in the Blender environment, Python (Van Rossum & Drake, 2009) scripts can be used. Python scripts are a versatile way to extend Blender functionality. Most areas of Blender can be scripted, including animation, rendering, import and export, object creation, and automating repetitive tasks. To interact with Blender, scripts can use the tightly integrated API, which has good written documentation. The two most common ways to execute Python scripts in Blender are using the built-in text editor or entering commands in the Python console. Rather than manually configuring your spaces for Python development, you can use the Scripting workspace accessible from the Topbar tabs. From the text editor, you can open “.py” files or paste them from the clipboard, and then test the script using the “Run Script” command. The Python Console is typically used for typing in snippets and for testing to get immediate feedback, but it can also have entire scripts pasted into it. Scripts can also be executed from the command line with Blender.

All our scenes will be scaled (1:1000) because of the enormous distance from the Earth to the Sun (150 · 10⁶ km). We start by creating an Earth prototype as a sphere with a 6371 km radius and a Sun that can be created as a Blender light source – the Sun-like object with an energy of 1.38 · 10³ W/m². The angular dimension of the Sun is 0.5 degrees, and Blender allows us to set this property directly without additional computations. Earth is situated in the center of the coordinate system, and Sun’s location is calculated at the middle moment of the RSO pass. The observer position is calculated from latitude and longitude as an (x,y,z) vector for a geographic position on, above, or below the Earth’s surface in the ITRS reference frame. The RSO and Sun positions are also calculated in the ITRS frame.

The model of the RSO can be made in a separate Blander file or downloaded from NASA 3D Models website (NASA, 2024) if it is available. We load the RSO model to the invisible container that will move through the object’s orbit and provide rotation and precession of the object in the coordinate system that is parallel to the Earth’s surface. To reproduce the orbit of the RSO, skyfield Python package is used (Rhodes, 2019). We can calculate the RSO position according to available TLE elements at the needed moment with any time resolution. The camera is situated in the observer point on the surface of the Earth and is looking permanently at the RSO. The camera’s focal distance is variable, to prevent change of the RSO size in the frame. An example of the loaded Blender scene is presented in Figure 1.

Figure 1.

Rendering occurs with the use of the Cycles render engine. Cycles is Blender’s physically-based path tracer for production rendering. It is designed to provide physically based results out-of-the-box, with artistic control and flexible shading nodes. It can work with GPU and CPU as well. We also tried to use an alternate render engine Evee, which is much faster but also less accurate when it comes to objects’ details shading. As a result, we obtain a series of frames (or animated video file) that covers RSO pass with the needed resolution (FPS), each frame of which can be processed to measure the light that we obtain from the entire frame. After we measure the flux of all frames, we can calculate additional parameters such as air mass, elevation, and range to satellite according to available TLE. Operating with such data, we can calculate standard magnitude (m_st) according to the equations: (1) minst=−2.5⋅logF {m_{inst}} = - 2.5 \cdot \log \left( F \right) (2) mst=A+minst+mz+mr {m_{st}} = A + {m_{inst}} + {m_z} + {m_r} where m_inst is an instrumental magnitude, F is a registered flux, A is a system zero point, and m_z, m_r are corrections for air mass and correction for standard range, which is defined as 1000 km for a case of artificial satellites, m_z and m_r are defined by equations: (3) mz=k⋅1coszmr=−5⋅logR1000 \matrix{{{m_z} = k \cdot {1 \over {cos \left( z \right)}}} \cr {{m_r} = - 5 \cdot \log {R \over {1000}}} \cr } where z is the zenith angle of the satellite, k is the coefficient of extinction, and R is the distance from the observer to the satellite. This process is described in detail in the paper by Kudak et al. (2017). For comparison with real observed LCs, all magnitudes were normalized.

The light curves that we analyzed for this project were obtained in the Laboratory of Space Research at Uzhhorod National University in Ukraine. The observations were taken at the Derenivaka observation point (coordinates Lat: 48.563417 N and Lon: 22.453758 E). We used a refractor telescope with a diameter of 120 mm and an equivalent focus of 114 mm, mounted on the alt-azimuth mount of the TPL-1M telescope (1-meter class) that can be operated from a PC. Our setup included a QHY174M-GPS camera and a Johnson R photometric filter, providing a field of view of 2.82°×1.76° and a scale of 10.6 arcsec/pixel. Light curves were processed using ccd_phot scripts described in the paper by Kudak et al. (2022). We also conducted observations in Uzhhorod city (coordinates Lat: 48.633689 N and Lon: 22.298902 E) using an electro-photometer with photo-multiplier and B, V filters on a base of an AFU-75 satellite camera.

Material types of the RSO’s 3D model are very important for the simulation of synthetic light curves. This means that we must not only correctly choose mirror surfaces like solar panels and diffuse scattering surfaces like RSO’s body but also choose the right color and other parameters for these details. For these purposes, in Blender software, we mainly use two materials for our 3D models: Principled BSDF for diffuse scattering and Glossy BSDF for specular reflecting surfaces.

The Principled BSDF combines multiple layers into a single easy-to-use node. It is based on the Disney principled model described in Burley & Studios (2012). This shader includes multiple layers to create a wide variety of materials. The base layer is a user-controlled mix between diffuse, metal, subsurface scattering, and transmission. On top of that, there is a specular layer, a sheen layer, and a clear coat layer. The Principled BSDF material has more than 20 parameters, the main of which are as fallows: Color – diffuse or metal surface color; Subsurface IOR – index of refraction for Subsurface Scattering; Metallic – blends between a non-metallic and metallic material model. A value of 1.0 gives a fully specular reflection tinted with the base color, without diffuse reflection or transmission. At 0.0, the material consists of a diffuse or transmissive base layer, with a specular reflection layer on top; Specular – amount of dielectric specular reflection. Specifies facing (along normal) reflectivity in the most common 0 – 8% range; Roughness – specifies microfacet roughness of the surface for diffuse and specular reflection; Alpha – controls the transparency of the surface, with 1.0 fully opaque. Usually linked to the Alpha output of an Image Texture node. and others.

The Glossy BSDF node is used for reflection with microfacet distribution and used for materials such as metal or mirrors. This material’s main properties are as fallows: Color of the surface, or physically speaking, the probability that light is reflected for each wavelength; Roughness – input for the surface roughness resulting in sharp to blurry reflections.

For the simulation of the TOPEX/Poseidon (NORAD:22076, COSPAR:92052A) satellite’s LCs NASA 3D model of this RSO

https://nasa3d.arc.nasa.gov/detail/topex-poseidon

We can conclude that TOPEX/Poseidon is spinning around three axes. This motion consists of the satellite’s own rotation around its main axis with the period P, the precession of this axis around the direction “satellite – the center of the Earth” (Z axis) with the period P_pr, and its rotation around the Earth with the orbital period P_orb = 112.4 minutes. The axis of the satellite’s own rotation also has a precessional motion around the direction “satellite – Earth center” that is trying to describe a circle in space. According to the different publications, RSO’s spin rate is P = 9 ∼ 10 seconds, the precession angle is in a wide range of values 12 – 45°, and the period of precession P_pr is 140 seconds or is discussed as orbital period precession with a period of 173 days. TOPEX/Poseidon has an inclined solar panel with a pitch angle of 12° (Sagnieres et al., 2020; Koshkin et al., 2018; Kucharski et al., 2017; Kudak et al., 2017).

COSMOS 2502 (NORAD:40358, COSPAR:14086A) is a satellite to collect intelligence for the Russian government (N2YO, 2024). Russian authorities identified the payload as a communications satellite, but independent observers believe that it is the second platform in a new series of spacecraft designed to eavesdrop on radio transmissions. The satellite’s orbit is the same type of orbit achieved in a November 2009 launch of a new Lotos S military intelligence spacecraft.

For simulation of the COSMOS 2502 satellite we build our 3D model from open source images and information available online (Krebs, 2024) This object was not as well studied as the TOPEX/Poseidon satellite. The only information about this RSO is available in ATLAS of light curves of space objects published in Odessa Astronomical Observatory (Koshkin et al., 2021). In the photometry database of Uzhhorod Laboratory of Space Research, we have 167 LCs of this object starting from 2020. According to observed LCs, we can clearly define the spin period (P ∼ 40 sec) of this object. 3D models of selected RSOs are presented in Figure 2.

Figure 2.

We have manually selected the appropriate rotation parameters for these objects and simulated synthetic LCs with Blender software according to their real passes upon the observation point. The parameters of synthetic LCs are presented in Table 1. Comparison of synthetic LCs with observed LCs presented in Figure 3.

Figure 3.

For automatization of finding self-rotation RSO parameters, we try to use the Markov Chain Monte Carlo method (Goodman & Weare, 2010) and its implementation written in Python – ecmcee described in paper Foreman-Mackey et al. (2013). In this case, the goodness of the fit is evaluated by the logarithm of the likelihood function ln p. The likelihood is the probability of a dataset given the model parameters and is equivalent to describing the generative procedure for the data. Alongside the likelihood function, prior knowledge about the model parameters is also used in the process of accepting or rejecting the generated states of the sampler. Parameters with uniform priors are required to specify the upper and lower boundaries of the sampling. On the other hand, the parameters with normal prior distribution are drawn around the most probable value with the expected standard deviation. Additionally, the normal prior distribution can be bounded with the lower and upper boundaries similar to the uniform distribution to prevent the sampler from reaching the undesired areas of the parameter space.

In other words, we write a function that can take initial parameters like the period of RSO self-rotation (P), phase of rotation (P_phase), period of spin axis precession (P_pr), its phase (Pr_phase), and angle of precession (Pr_angle) and produce an LC according to the 3D model of RSO and its orbital motion predicted from TLE elements. We can also use other parameters to fit – the refraction coefficient of some RSO’s materials, the size of solar panels, even changing the spin axis, or other parameters. Of course, the fewer parameters we have to fit, the less time it takes to find them.

As the next step, an MCMC method runs a Sampler that goes over all possible variants of our parameters within the boundaries that we set for each of them. For example, the spin period can be estimated from LC by other methods, like the Lomb – Scargle periodogram analysis (Lomb, 1976; Scargle, 1998), and that is why we can set this parameter to be variable only in a small range, for example, 1%. Other parameters can have wider ranges. All limits of variable parameters can be seen on the corner plot (Figure 5).

The disadvantage of this method is the long time that must be spent to go over all possible variants of variable parameters and find a set that has the smallest residuals with observed LC (our likelihood function). If the model is simple and the dataset is small, we can use a large number of iterations which is good for the MCMC method. In our case, one model estimation takes ∼ 30 – 90 seconds (depending on LC duration). To speed up the process, we activate the Blender option to use a GPU card to render frames of our LC. As the testing PC, we use an Intel Core i7-8700 (6 cores at 3200 Mhz) CPU with 16GB of RAM and a Nvidia RTX 4060 Ti graphic card (2860 CUDA cores at 2860 Mhz). To fit 30 seconds long LC of the Topex/Poseidon satellite with 500 iterations and 20 walkers, we spend ∼ 20 hours of processing time on our testing PC. We also used parallelization in two threads. In this case, usage of CPU reached 100%, and usage of GPU was up to 20%. But, the results were worth it. We obtain a good fit of LC (see Figure 4), and corner plots of fitted parameters are presented in Figure 5. With the increasing number of iterations, we can achieve better precision for variable parameters.

Figure 4.

Figure 5.

In RSO orientation research, one of the main problems is proving the reliability of the parameters of the satellite’s own rotation proposed by the researcher. For such purposes, we test the Blender software. We use well-studied RSO TOPEX/Poseidon to model observed LCs. Rotation parameters of this RSO were studied by many authors, and at the same time, there is a certain discrepancy in the values of the rotation parameters in different publications. First, we use manual generation of synthetic LC for this RSO using different parameters of its rotation. Not all parameters from previous research available in publications correspond well to observed LC. The best results for manual fit for this RSO are presented in Table 1 and Figure 3. The main parameters that affect the shape of observed LC are spin period and precession angle. The angle of the solar panel (SP) also affects the shape of the synthetic LCs, but this value has been well-studied in previous studies, and the SP angle of 12 degrees compares well with observed LCs. To confirm our results of the manual fit of synthetic LC we use the MCMC method that allows us to choose the studied parameters of RSO rotation in a wide range. According to Table 1 with the MCMC method, we obtain similar main parameters (spin period and precession angle). The phases of the spin period and precession period are different because different parts of the same observed LC were used. As for the precession period, we can conclude that large errors are related to the fact that to determine this parameter, it is necessary to analyze long light curves samples. Our synthetic LCs do not perfectly fit the observed LC, this fact is because of using an not ideal 3D model of RSO and its surface parameters. But even with such basic models, we were able to obtain good enough results in the simulation of selected LCs.

If we compare the manually generated synthetic LC of TOPEX/Poseidon (Figure 3) with the MCMC generated LC (Figure 4), we can say that the MCMC light curve has smaller residuals (∼0.4 compared to ∼0.7 in manually created synthetic LC) and much better correspondence to real observed LC mainly when spike occurs. We can conclude that the manual fit of LC sometimes can be very painful for the user. It is a way of right and bad choices that are made in an almost blind way. The user must have an intuition based on the huge experience to tell what part of observed LC can correspond to some parts on the surface of RSO and also guess the orientation of this RSO. Finally, we must also emphasize that it is not always possible to perform manual modeling of an arbitrarily selected part of observed LC. Modeling synthetic LCs with the help of the MCMC method is a much easier way, but it demands more time to process. In both cases (manual and MCMC modeling of observed LCs), it makes sense only if you have good and fast photometry of the selected RSO. The higher the spin of the RSO the faster photometry is needed to be sure in the identification of some parts of the RSO body.

As another example of synthetic LC modeling, we use the COSMOS 2502 satellite. This object was not as well studied as the TOPEX/Poseidon satellite. We were able to build a 3D model of this RSO from available information. Starting from the year 2020, we have 167 LCs for this object in our photometric database. From these LCs, we determine the spin period of RSO, and other main parameters were determined during the manual fit process (see Table 1).

Solving the task of inverse LC modeling with the help of open-source software available for a wide range of scientists interested in RSO LC analysis can lead to the same effect as we can see for the eclipsing binary systems. There is open-source software like PHOEBE, Binary Maker, and ELISa (Prša et al., 2016; Bradstreet, 2005; Čokina et al., 2021) capable of determining eclipsing binaries system’s physical parameters according to precisely observed LCs and radial velocities, and a lot of scientists use these options to expand general knowledge about variable stars and their parameters and thus understand their evolution. Despite the greater complexity of the inverse photometry task for RSO, even partially solving it, we can better understand the physics of RSO rotation and its behavior in near-Earth space and make in-depth conclusions about the behavior of space debris and working spacecraft in orbit.

To summarize the above, we can say that Blender software is a good choice for modeling RSO light curves. It has useful features like Python scripting and GPU support benefits. Modeling the observed LC can be done manually if the operator knows at least approximately the RSO rotation parameters or using an MCMC method that fits all the necessary parameters if you have an RSO model that is close to the true RSO shape. In theory, it is also possible to specify the parameters of the 3D model if all the rotation parameters are well-known. However, this situation is extremely rare for this type of research.