© The Authors 2026

1 Introduction

Companies operating low Earth orbit satellite (LEO) constellations plan to launch tens of thousands of communications satellites over the next decade. Some of these companies are ONEWEB, Starlink, Amazon Kuiper, Guowang, and StarNet (Boley et al. 2022). The number of satellite launches is continually increasing, with the number of satellites in LEO expected to exceed one hundred thousand by the end of the decade (Venkatesan et al. 2020) and over ∼560 000 based on requests submitted to the U.S. Federal Communications Commission (FCC) and the International Telecommunication Union (ITU) filings¹.

Low Earth orbit satellites are generally deployed at altitudes between 350 km and 2000 km. Depending on the orbital height, the satellites may reflect the sunlight and become visible to observers beyond twilight and into the night. For this reason, satellites are most visible for one to three hours after sunset or before sunrise. Additionally, during summer at high latitudes, LEOsats are visible almost all night (McDowell 2020).

Due to the number of current and planned LEOsats in orbit, and their ability to reflect sunlight, they negatively impact several areas of observational astronomy, especially the observation of rare transients that require scanning large portions of the sky. Examples of such cases include the search for and tracking of near-Earth objects (NEOs), inner solar system asteroids, as well as deep multi-object spectroscopic surveys such as those performed by the Dark Energy Spectroscopic Instrument (DESI) wide-field spectrograph installed on the Nicholas U. Mayall 4 m telescope at Kitt Peak National Observatory (Flaugher & Bebek 2014) and the Dark Energy Camera (DECam) installed on the 4 m Victor M. Blanco telescope at the Cerro Tololo Inter-American Observatory in Chile (Diehl et al. 2012). Satellite brightness might also affect wide-field exoplanet transit surveys such as the Hungarian-made Automated Telescope Network (HATNet), whit its 10.4^◦ × 10.4 ^◦; field of view (FoV) and the Next Generation Transit Survey (NGTS), consisting of 12 telescopes each with an 8 deg² FoV, providing a total coverage of 96 deg² (Wheatley et al. 2013). Low Earth orbit also reflect in both the near-infrared (NIR) and infrared (IR) regimes; thus, surveys using these wavelengths, including the wide-field NIR camera at the 4 m United Kingdom Infra-Red Telescope on Mauna Kea (Casali et al. 2007), are affected.

The most noticeable impact from LEOsats on astronomical observations, in both the optical and the NIR, is bright streaks across the detector. These streaks affect amateur and professional astronomers (Walker et al. 2020). Different groups in the observational astronomy community have reported the reflective brightness of Starlink (a subsidiary of SpaceX) satellites (e.g., Tregloan-Reed et al. 2020, 2021; Tyson et al. 2020; Boley et al. 2022; Mróz et al. 2022; Horiuchi et al. 2023), while attempting to quantify the risks that LEOsat constellations pose to ground-based astronomy. Recently, an international observational campaign led by the IAU Centre for the Protection of the Dark and Quiet Sky from Satellite Constellation Interference (CPS)² has been successful in measuring the reflective brightness of the AST SpaceMobile Bluewalker 3, a prototype 4G/5G satellite, which at the time of launch was the largest communications satellite in LEO with a 64 m² array. The reported results found that BlueWalker 3 shows a stationary V magnitude of 0.4 ± 0.1 (Nandakumar et al. 2023). These results support the importance of the development of national and international legal frameworks, along with recommendations and mitigation strategies, for both LEOsat operators and observatories (Walker et al. 2020, Dark & quiet skies).

The impact of individual LEOsats on astronomy depends on how bright they appear. The reflected sunlight from satellites comprises two specular and diffuse components. Specular reflection (also called regular reflection) corresponds to light scattered in a preferential direction (mirror-like reflection), whereas diffuse reflection occurs when light is scattered in several directions. The type of scattering depends on the detailed shape of the satellites.

Both the materials used for the construction of the satellites and the roughness of the surface materials affect the type of scattering. For instance, phased antennas produce a more diffuse scattering effect and less of a specular reflection. This is one reason why the first-generation of Starlink LEOsats, which lacked reflective brightness mitigations, had an on-station visible stationary magnitude of V ∼ 5, making them visible to the naked eye (Tregloan-Reed et al. 2020, 2021; Walker et al. 2020; Tyson et al. 2020; Krantz et al. 2023; Boley et al. 2022). During the satellites’ deployment phase, their orientation in space also matters in how they reflect sunlight. Krantz et al. (2023), and Boley et al. (2022) studied the effects of the Sun phase angle and range on the apparent magnitude of Starlink satellites using a brightness model: a Lambertian diffuse-specular sphere by McCue et al. (1971); Williams & McCue (1966); Pradhan et al. (2019). Another reason why LEOsats are much brighter than geostationary satellites is because of their low orbits, with altitudes between 350 and 1200 km. When considering the inverse square law, a satellite at 500 km appears 3600 times brighter (∼8.9 mag brighter) than a satellite at 30 000 km.

To date, several studies have been done to measure the effectiveness of different mitigation designs aimed at reducing the reflective brightness of LEOsats (e.g., Hainaut & Williams 2020; McDowell 2020; Tyson et al. 2020; Tregloan-Reed et al. 2020, 2021). Studies report that for the first Starlink mitigation design, Darksat, the optical magnitude was reduced from ∼5 to ∼6 mag while in the NIR it went from ∼2 to ∼3 mag (Tregloan-Reed et al. 2020, 2021; Tyson et al. 2020). The impact of reflected sunlight from LEOsats on astronomy will be even greater for telescopes with instruments designed for long exposures on ultra-wide fields, due to the combined effect of the telescope light-collecting area and its angular FoV (Hainaut & Williams 2020; McDowell 2020). A notable example of this is the Charles Simonyi telescope’s camera detectors at the National Science Foundation’s Vera C. Rubin Observatory. Tyson et al. (2020) analyzed the impact of LEOsats on this telescope and camera. They concluded that an optical magnitude of g ≥ 7 is necessary for satellites located at an orbital height of 550 km to enable corrections of the trails left by LEOsats to a level below the background noise. For orbital heights above 550 km, the LEOsats must be fainter because their lower angular velocity means they illuminate they illuminate the detector pixels for longer. The first workshop on the issue of satellite constellations and their impact on astronomy (Walker et al. 2020) found that the magnitude limit m_lim, for which the above-mentioned correction is possible, is related to the orbital height h by $$m_{\text{lim}}=\{\begin{array}{ll}7\hfill & \text{when}h\le 550\text{\hspace{0.17em}km}\hfill \\ 7+2.5{\log }_{10}(h/550\text{\hspace{0.17em}km})\hfill & \text{when}h>550\text{\hspace{0.17em}km}\hfill \end{array}.$$(1)

For example, satellite constellations with orbital heights of 1200 km, such as ONEWEB, must not exceed g = 7.9 mag.

The characterization of the reflected brightness of an object can be represented by different mathematical models. For example, the McCue et al. (1971) model offers an approximation of the reflected brightness of an object with different geometries (spherical, cylindrical, and flat). On the other hand, the Minnaert (1941) model studied how light is reflected on an opaque surface. Both models have been used to characterize observations of Starlink LEOsats by Pradhan et al. (2019); Mallama (2020); Tregloan-Reed et al. (2020, 2021); Lawler et al. (2022); Boley et al. (2022); Krantz et al. (2023).

The aim of this work is the characterization of the stationary magnitude of the ONEWEB satellite constellation using two different sphere models by McCue et al. (1971) which consider the phase angle, range, effective area, cross section, albedo and mixing coefficient. In addition, we study the reflective opacity of LEOsats using the approach of Tregloan-Reed et al. (2020), based on the model of Minnaert (1941), which considers the incidence and observer angles and edge darkening parameter. All models are applied to a group of over 300 satellites.

A satellite constellation consists of a group of artificial satellites with identical circular orbits, height and inclination, grouped in a series of orbital planes with nodes located along the equatorial plane, as defined by Walker (1984). The ONEWEB constellation consists of a set of 720 satellites designed and manufactured by Airbus Defence and Space and ONEWEB Ltd³. The satellites have a box-wing design with a central box-shaped bus and two solar panel arrays. The design also includes two extended antenna reflectors (Krantz et al. 2023). The constellation has 18 circular polar-orbital planes at an altitude of 1200 km and inclined is 87^◦, providing global coverage. The main aim of this constellation is to enable high-speed, low-latency inter-net access worldwide, including remote parts of the world (del Portillo et al. 2019).

The ONEWEB constellation has three different stage orbits (i.e., orbital parked, orbital rise, and on-station), which depend on their activity. Orbital parked is when the satellites have an altitude of less than ∼630 km, and they are parked at the start of their service life while waiting to perform an orbital rise or complete service. Orbital rise is when the satellite has an altitude between ∼630 to ∼1100 km. In this stage satellites are maneuvered to their final orbit while a increasing altitude. On-station is when the satellites have reached their operational height (t∼ 1200 km) and are ready to start operations. A plot of the range versus altitude of the data studied in this paper is shown in Figure C.1 in Appendix C. As of 2024, all the ONEWEB constellation satellites occupy their final on-station position.

2 Observation and data reduction

To observe LEOsats, it is necessary to know where to point the telescope. This information is encoded in the two-line elements (TLE), which is a standardized format used to describe the orbit elements of a satellite. Two-line elements consist of a set of data that includes orbital parameters, such as the satellite’s position, velocity, and other relevant information at a specific epoch or time (Croitoru & Oancea 2016). Two-line elements were developed by the North American Aerospace Defense Command (NORAD). NORAD is an organization established by the United States and Canada to monitor North American airspace for objects in orbit between ∼100 and 40 000 km altitude. All TLEs are displayed on a website where they are available to the public. They can be used for multiple technical studies when the position of a satellite is needed at some instant in time.

The task of transforming the data provided by the TLEs into celestial coordinates was carried out using the Satellite-Tracking⁴ code. To use the code, users must select the observatory and the observation time (evening or morning twilight). Then, the Satellite-Tracking code downloads the latest available TLE from the Celestrak⁵ website, outputs the results in an easy-to-read text file as a function of time, allowing observers to scroll down the list, selecting the targets to observe, as seen in Table 1.

In general, when observing LEOsats, the wait-and-capture method is used, which consists of tracking background stars at sidereal speed. This type of observation produces an image where the stars are dots and the LEOsats are represented by a line crossing the image. The angular velocity of LEOsats depends on their orbital height and sky position (their angular velocity decreases at higher elevations), and it can range from a few arcmin s⁻¹ to about 2^◦ s⁻¹. This effectively means that LEOsats can transverse the FoV in less than a second (depending on the angular velocity and FoV size), which requires precise timing of the opening and closing of the camera shutter. During exposures of a few seconds, the light collected from the LEOsat must be scaled to the exposure time, as sky and stellar light are effectively collected throughout, whereas the satellite’s full trail may partially fall outside the FoV. This is expressed in Equation (2).

Data reduction and astrometric calibration were performed using the LEOSatpy⁶ code. The code was written in Python and uses Astropy for flat-field, bias, and dark subtraction. The astro-metric calibration begins with the alignment of images with the same background, both with and without satellites, followed by the identification of references stars from the Gaia Early Data Release 3 catalog (Lanzafame et al. 2023).

The next step is to compute the stationary magnitude m_s and position of the satellite at the midpoint of the line. To this end, the magnitude of each LEOsat is extrapolated using the equation by Boley et al. (2022): $$m_s=m_0+m_{\text{ins}}-2.5{\log }_{10}\left(\frac{\varphi }{L}{t}_{\text{exp}}\right),$$(2)

where m₀ is the zero-point magnitude of the stars, m_ins is the instrumental aperture magnitude, ϕ is the angular velocity of the satellite in the sky, L is the length of the line left by the LEOsat, and t_exp is the exposure time. The stationary magnitudes distribution of our sample of LEOsats is shown in Figure 1.

After that, the range-corrected magnitude H¹²⁰⁰ is set as the brightness that the LEOsat would have at a range of 1200 km, the nominal altitude of the ONEWEB constellation, as follows: $$H^{1200}=m_s-5{\log }_{10}\left(\frac{r}{1200\text{\hspace{0.17em}km}}\right),$$(3)

where r is the range of the satellite in km. The observations used in this work were made using a 60 cm-aperture telescope operated by the Chungbuk National University Observatory in Jincheon (CbNUOJ) in South Korea (Kim et al. 2014; Wang et al. 2019). The telescope used an SBIG STX-16083 charge-coupled device (CCD) camera of 4096 × 4096 pixels, allowing a FoV of 72′ × 72′ and a pixel scale of 1.05 arcsec/pixel. In addition, Johnson/Cousins UBVRI filters can be used through an electronic filter wheel system. This study uses images taken with V and R filters with a binning of 2 × 2 pixels to compensate for the altitude of the ONEWEB satellites by collecting more light.

The telescope control system uses an internet program that connects to a Yb (Ytterbium) optical lattice clock at KRISS (Korea Research Institute of Standards and Science, Heo et al. 2022) to synchronize the time each minute. The precision of the time synchronization is approximately ± 0.01 s. In addition, the precision of CCD shutter operation is measured to ± 0.005 s.

The data presented in this work consist of 626 observations of 322 individual ONEWEB satellites, observed over 41 nights (from October 2021 to September 2022) during evening twilight. A total of 175 observations were conducted for off-station satellites (579–1100 km) whose average altitude is 684.94±108.79 km and average range is 914.37±182.01 km, and 451 observations of on-station satellites (1100–1229 km) whose average altitude is 1211.40±9.05 km and average range is 1542.20±201.87 km, as shown in Figure C.1.

The data were studied as a whole and also separated into four samples according to the filter used (V or R) and their on- or off-station position to study how the brightness models behave for satellites at different altitudes. More details about the samples per filter can be seen in Table A.1. Appendix D contains a glossary of all the variables and notations used in this work.

Fig. 1 Distribution of stationary magnitudes of all ONEWEB satellites studied. The sample is divided by filter and orbit type. There are 451 on-station satellites, of which 410 were observed with the V filter and 41 with the R filter. There are 175 off-station satellites observed 158 with the V filter and 17 with the R filter. The average µ and standard deviation σ of ms are shown in the legend.

3 Modeling of the satellite brightness

Satellites shine because they reflect sunlight, i.e., sunlight hits the satellites and their parts (the main body, antennas, and solar panels), which reflect a fraction of this sunlight toward our direction. This sunlight reflected by the satellite is what we can observe from Earth and measure through astronomical observations by the stationary magnitude m_s. The satellite does not reflect all the light it receives, and the amount of light we observe depends on several factors:

• (i)

  the physical characteristics of the satellite: its shape and those of its parts (main body, antennas, and solar panels), the type of scattering (diffuse or specular), the materials they are made of, and the roughness of the surfaces;

• (ii)

  the geometry between the Sun, the satellite, and the observer;

• (iii)

  the orientation of the satellite itself, its solar panels, and antennas with respect to the observer;

• (iv)

  the distance between the satellite and the observer.

In practice, no function can describe all these factors, but it is necessary to find models that can estimate most of them. The approach is to assume several hypotheses that allow to simulate the general trends and brightness distributions of satellites.

In this work, we consider two types of models: the Lambertian spherical models which allow the observed magnitude to be reproduced as a function of the Sun-satellite-observer geometry and type of scattering (diffuse or specular), and the relative reflectance (RR) model which compares the magnitudes observed from two different orientations. These models are generally referred to in the literature as bidirectional reflectance distribution functions or BRDFs (Shepard 2017, p. 26).

3.1 Lambertian spherical models

In this model, we consider the satellites as having a uniform spherical shape. In addition, we consider two cases of dispersion: one in which the uniform sphere is diffuse, and another in which the uniform sphere is both diffuse and specular. The specular reflection corresponds to light scattered in one direction as in a mirror, whereas diffuse reflection is when light is scattered in several directions.

The stationary magnitude m_s of LEOsats can be reproduced using the Lambertian sphere model described by Williams & McCue (1966) as follows: $$m_s=m_{\odot }-2.5\cdot {\log }_{10}[\zeta \cdot F_0(\Phi )]+5\cdot {\log }_{10}(r),$$(4)

where m_⊙ = −26.77 is the Sun’s apparent magnitude in the Johnson V filter (Willmer 2018), ζ is the so-called effective area of the sphere reflecting light from the Sun, F₀(Φ) is the so-called phase function which describes the type of scattering mentioned above as a function of the Sun phase angle Φ, and r is the range of the satellite (distance between it and the observer). ζ = ρ · A, where ρ is the albedo of the reflecting sphere and A = πR² is its cross-area. The Sun phase angle Φ is defined as the angle between the Sun, the satellite, and the observer (see Fig. 2). Φ is measured on the plane formed by two vectors: the one going from the satellite to the Sun and the one from the satellite to the observer. Of the several factors affecting the observed brightness and itemized at the beginning of Section 3, in Equation (4), Φ accounts for the effects of item (ii), ζ likewise for item (i), while the log₁₀(r) term measures the effect of item (iv). Appendix B briefly describes the derivation of this equation.

3.1.1 Lambertian diffuse spherical model

The diffuse spherical (DS) model assumes that a LEOsat reflects light like a uniform diffuse sphere, in which case the phase function is given by $$F_0(\text{Φ})=\frac{2}{3\cdot {\pi }^{p+1}}\cdot {[}^{(}.$$(5)

Boley et al. (2022) has used this model.

Fig. 2 Graphical representation of the position of the satellite and definitions of the Sun incidence angle θ0 (in orange), the observer angle θ1 (in dark blue), the Sun phase angle Φ (in green), and the solar elongation angle γ (in brown). Notice that θ0, θ1 and Φ are not contained in the same plane but form a triangular pyramid. The distance is not to scale.

3.1.2 Lambertian diffuse and specular spherical model

The diffuse and specular spherical (DSS) model assumes that a LEOsat reflects light like a uniform diffuse and specular sphere, in which case the phase function is given by $$F_0(\text{Φ})=\beta F_1(\text{Φ})+(1-\beta )F_2(\text{Φ}),$$(6)

with 0 ≤ β ≤ 1 being the mixing coefficient of diffuse-specular reflection. The function F₁(Φ) represents the diffuse part of the phase function, while F₂(Φ) represents the specular portion of it, given respectively by $$F_1(\text{Φ})=\frac{2}{3\cdot {\pi }^2}[(\pi -\text{Φ})\cdot \cos (\text{Φ})+\sin (\text{Φ})],$$(7) $$F_2(\text{Φ})=\frac{1}{4\pi }.$$(8)

When β = 0, F₀(Φ) is only specular, and when β = 1, F₀(Φ) is only diffuse and equivalent to Equation (5) for p = 1. Williams & McCue (1966); McCue et al. (1971); Lawler et al. (2022); Krantz et al. (2023) have used this model.

In the Lambertian spherical models, both A and ρ shift the model up or down (larger area or albedo, brighter magnitudes, and reciprocally). For the DS model, p changes the slope of the function, which is always concave down as a function of Φ. On the other hand, in the DSS model, β changes both the slope and concavity of the model. The larger the specular portion, the larger the change in the observed magnitude with Φ, in the sense that a mirror looks brighter than an opaque surface. Also, at low values of Φ the sphere becomes brighter (sphere is in “full phase”) and at large values of Φ, it becomes fainter (sphere is in “new phase”). In Figures C.2 and C.3 in Appendix C, we show the behavior of m_s for both models, considering different values of ρ, A, and p (DS); and ρ, A and β (DSS), for r = 1200 km. Since in both cases, ρ and A have the same functional form, they cannot be fitted separately; therefore, we fit ζ = ρ · A.

3.1.3 Effective albedo for the DS and DSS models

The DS and DSS models are a good approximation for the apparent brightness of satellites, but they do not take into account the actual geometry of the satellite. The apparent brightness of satellites is not only the result of the combination of the reflection of their diffuse and specular parts, but also of small changes in the geometry of the illuminated surface due to the orientation of the satellite with respect to the observer Krantz et al. (2023). The orientation of the satellite only allows it to reflect light with some of its parts (Williams & McCue 1966; McCue et al. 1971).

These changes in the reflected light due to the orientation changes can be represented by a changing albedo for a fixed cross section area A = 1 square meter. Such albedo for A = 1 is called effective albedo and does not represent a real albedo. In fact, it can be larger than 1, but it represents the relative change between the different geometries of the incident light on the satellite and the reflection produced by its surfaces, i.e., item (iii) listed in Section 3.

For the DS model, with p = 1, using Equations (4) and (5), the effective albedo is given by $${\rho }_{\text{eff,ds}}=\frac{3{\pi }^2\cdot r^2\cdot {10}^{\frac{m_{\odot }-m_s}{2.5}}}{2\cdot [(\pi -\text{Φ})\cdot \cos \text{Φ}+\sin \text{Φ}]},$$(9)

which has been studied by Krantz et al. (2023). In this work, we also examine the effective albedo of the DSS model, with β = 0.5, which from Equations (4, 6, 7) and (8) is given by $${\rho }_{\text{eff,dss}}=\frac{8\cdot 3\cdot {\pi }^2\cdot r^2\cdot {10}^{\frac{m_{\odot }-m_s}{2.5}}}{4\cdot 2\cdot [(\pi -\text{Φ})\cdot \cos \text{Φ}+\sin \text{Φ}]+3\pi }.$$(10)

3.2 Relative reflectance model

The RR model allows us to compare the magnitudes of a satellite observed at the same range from two different orientations, following Tregloan-Reed et al. (2021). Each orientation is parameterized by the angles (θ₀, θ₁) and $({\widehat{\theta }}_0,{\widehat{\theta }}_1)$, where θ₀, and ${\widehat{\theta }}_0$ are the Sun incidence angles and θ₁ and ${\widehat{\theta }}_1$ are the observer angles in each orientation (see Fig. 2). By definition, the difference in stationary magnitudes between two different orientations can be obtained from the ratio of the reflected fluxes by $$m_s({\theta }_0,{\theta }_1)-m_s({\widehat{\theta }}_0,{\widehat{\theta }}_1)=-2.5{\log }_{10}\left(\frac{F_{{\theta }_0,{\theta }_1}(r)}{F_{{\widehat{\theta }}_0,{\widehat{\theta }}_1}(r)}\right).$$

In the RR model, such ratio is $$\frac{F_{{\theta }_0,{\theta }_1}(r)}{F_{{\widehat{\theta }}_0,{\widehat{\theta }}_1}(r)}=ℛ={\left(\frac{\cos {\theta }_0\times \cos {\theta }_1}{\cos {\widehat{\theta }}_0\times \cos {\widehat{\theta }}_1}\right)}^{k-1},$$

where 0 ≤ k ≤ 1 is known as the edge darkening parameter, which describes the angular distribution of the brightness of the reflected light. When k = 1, reflectance in all directions is the same, that is the model describes the reflectance of the surface as uniform, i.e., the surface is diffuse (Tregloan-Reed et al. 2021).

The ratio ℛ is the same regardless of the satellite’s properties and its range. In other words, the specific values of the stationary magnitudes being compared will change with r, the RR model only describes their difference. It is reasonable to assume that a satellite has the same stationary magnitude at a fixed reference orientation and range; therefore, we can convert any observed m_s to that fixed value. Consequently, with the RR model, we reproduce the observed range-corrected magnitude H¹²⁰⁰, using the following equation: $$H_{({\theta }_0,{\theta }_1)}^{1200}=H_{(70,0)}^{1200}-2.5{\log }_{10}\left[{\left(\frac{\cos {\theta }_0\times \cos {\theta }_1}{\cos {70}^{\circ }\times \cos 0^{\circ }}\right)}^{k-1}\right],$$(11)

where we chose the reference orientation (70^◦, 0^◦) and $H_{(70,0)}^{1200}$ is an additional parameter of the model to estimate. The chosen ${\widehat{\theta }}_0={70}^{\circ }$ corresponds to the Sun at 30^◦ below the horizon and ${\widehat{\theta }}_1=0^{\circ }$ for the satellite at the observer’s zenith. This brightness model attempts to describe the relative effects of item (iii) in Section 3.

Figure 2 shows a graphical representation of the angles used in the brightness models: Sun phase angle Φ, Sun incidence angle θ₀, and observer angle θ₁. These angles are defined as follows: Φ is formed by vectors $\overrightarrow{(\text{SatSun})}$ and $\overrightarrow{(\text{SatObs})}$, θ₀ is formed by $\overrightarrow{(\text{SatSun})}$ and $\overrightarrow{({\text{SatNadir}}_{\text{Sat}})}$, and θ₁ is formed by $\overrightarrow{({\text{SatNadir}}_{\text{Sat}})}$ and $\overrightarrow{(\text{SatObs})}$. The Sun incidence angle is obtained from θ₀ = 180^◦ − z_⊙,sat, where z_⊙,sat is the Sun’s zenith distance at the Satellite, estimated from the expression: $$\begin{array}{ll}\cos (z_{\odot ,\text{sat}})=\hfill & \cos (\text{ΔGeog}(\text{obs,\hspace{0.17em}sat}))\cdot \cos (z_{\odot ,\text{sat}})\hfill \\ \hfill & +\sin (\text{ΔGeog}(\text{obs,\hspace{0.17em}sat}))\cdot \sin (z_{\odot ,\text{obs}})\hfill \\ \hfill & \cos (\text{ΔAz}(\odot ,\text{sat})),\hfill \end{array}$$

where ∆Geog(obs, sat) is the angular distance between observer’s longitude and latitude and the satellite’s longitude and latitude, ∆Az(⊙, sat) is the difference between the azimuth of the satellite and the azimuth of the Sun, and z_⊙,obs is the Sun’s zenith distance for the observer.

On the other hand, the observer angle θ₁ is obtained from the triangle △O_⊕O_satO_obs: $$\begin{array}{l}{R}_{\oplus }^2=r^2+{(R_{\oplus }+H_{\text{sat}})}^2-2\cdot r\cdot (R_{\oplus }+H_{\text{sat}})\cdot \cos {\theta }_1\Rightarrow \hfill \\ {\theta }_1=\arccos \left(\frac{r^2+2\cdot R_{\oplus }\cdot H_{\text{sat}}+H_{\text{sat}}^2}{2\cdot r\cdot (R_{\oplus }+H_{\text{sat}})}\right)\hfill \end{array}$$

where R_⊕ is the Earth’s radius, H_sat is the perpendicular height (above Earth) of the satellite, and r is the range of the LEOsat. The angles θ₀ and θ₁ are generally not in the same plane, due to a third direction or vector (from the satellite to its nadir) considered in their definitions.

The solar elongation angle γ is obtained from the law of cosines on the equatorial coordinates of the Sun (α_⊙, δ_⊙) and satellite (α_sat, δ_sat): $$\cos (\gamma )=\sin ({\delta }_{\odot })\sin ({\delta }_{\text{sat}})+\cos ({\delta }_{\odot })\cos ({\delta }_{\text{sat}})\cos ({\alpha }_{\odot }-{\alpha }_{\text{sat}})$$(12)

Finally, to compute the Sun phase angle Φ, we must notice that since the Sun can be taken at being at an infinite distance; therefore, the directions satellite-Sun and observer-Sun are parallel and Φ = 180^◦ − γ.

4 TLE-predicted versus actual position of satellites

Two-line elements have some limitations in terms of the accuracy of the data they provide. The Simplified General Perturbation 4 (SGP4) model used to define the orbital parameters includes the dynamic effects of atmospheric drag and solar radiation pressure, providing a fast propagation time but low accuracy (e.g., Acciarini et al. (2025), their Section 4.1). Furthermore, TLEs do not consider maneuvers made by operating companies to adjust the orbit. In addition, the TLE format does not include covariance matrices that allow error estimation from the orbital parameters (Vallado & Cefola 2012).

One way to estimate the accuracy of TLE positions is to use satellite observations to estimate these errors. Astronomical observations are independent of the data and the SGP4 model that generates the TLEs. The results of the Satellite-Tracking code tell us the time and equatorial coordinates at which we must point the telescope to see a satellite pass through the center of the image. However, in most cases, the satellite does not pass through the center of the image, and there is an angular difference σ_tle between the TLE-predicted and the actual observed coordinate.

In our sample, these differences are within 2 arc minutes (see Fig. 3), which is around the typical value observed by other authors with other LEOsats constellations, –for example Nandakumar et al. (2023, Blue Walker, 7.2 ± 3.1 arcmin), Krantz et al. (2023, ONEWEB constellation, between 1.8 and 2.4 arcmin), and Halferty et al. (2022, Starlink constellation, on average 0.576 deg).

Fig. 3 Histogram of σtle for all the ONEWEB satellites studied, showing the difference between the TLE-predicted and the observed equatorial coordinated of satellites.

5 Data analysis: Differential evolution Monte Carlo

We start by looking at the stationary magnitude distribution of the ONEWEB satellites, as obtained from Equation (2) and displayed in Figure 1. Our data fall within the range presented by Krantz et al. (2023), with 3.9 ≤ m_s ≤ 11.0 mag. The average magnitude for the whole sample is 7.8 mag, in which 74% are on-station satellites and 26% off-station satellites. The mean stationary magnitude error is 0.19 with a standard deviation of 0.06, and it ranges between 0.1 and 0.3 magnitudes for 98% of the data.

We used the differential evolution Markov chain Monte Carlo (DEMC) method to fit the data with the proposed brightness models. This method, which combines the Markov chain Monte Carlo and differential evolution algorithms, makes the search process more efficient by avoiding sampling local minimal approximations. The method adjusts multiple chains in parallel that learn from each other instead of running independently. New chains are accepted or rejected according to the Metropolis-Hastings criterion (Sherri et al. 2019).

A Python routine was created using the PyMC library to fit the DS and DSS models to the stationary magnitudes m_s, and the RR model to the range-corrected magnitudes H¹²⁰⁰, using DEMC. All brightness models were computed for all satellites, the on-station, and the off-station subsets.

In general, to use the PyMC library and DEMC, users must perform the following: assign the prior distributions of both the parameters and the noise of the observations; define the deterministic function to be fit, which depends on the parameters and the observations; define the likelihood of the observations, which in this case corresponds to a normal distribution of the observations around the function being fitted (see Table 2); and, lastly, define the step method and the parameters for iterations. In our case, everything was constant for the three models (DS, DSS, RR). The Metropolis-Hastings criterion (Metropolis, class of PyMC) was used, setting the tuning frequency (tune_interval) equal to 150 and the initial scale factor for proposal (scaling) equal to 0.01. The iteration parameters were set using the sample function of PyMC. Moreover, following parameters were set: the number of samples to draw (draws) equal to 25 000, number of iterations to tune (tune) equal to 25 000, number of chains to sample (chains) equal to 10, step function or collection of functions (step) equal to Metropolis, and random seed(s) used by the sampling steps (random_seed) equal to 42.

6 Results

6.1 Lambertian diffuse spherical model

The DS model in Equations (4) and (5) describes the stationary magnitude m_s as a function of the Sun phase angle Φ and fixed parameters ζ and p. The DEMC algorithm finds the values of ζ and p for which the DS model reproduces best the magnitudes m_s. The input data are m_s, Φ, and r, and the resulting parameters can be seen in Table 3, for all satellites, on-, and off-station samples. Figure 4 shows the observed stationary magnitudes and the ones computed by the model, as well as the residuals between them, for all satellites. The same parameters results (within 10⁻³) are obtained when the fitting is performed using scaled magnitudes H¹²⁰⁰. Parameters results for the V filter and R filter samples can be seen in Tables A.2 and A.3, respectively, in Appendix A.

6.2 Lambertian diffuse specular spherical model

Likewise, in the previous section, the DSS model in equations (4) and (6) describes the stationary magnitude m_s as a function of Φ and fixed parameters ζ and β. The DEMC algorithm finds the values of ζ and β for which the DSS model reproduces best the magnitudes m_s. The input data are the same than in the DS model, and the resulting parameters can be seen in Table 3, for all satellites, on-, and off-station samples. Figure 5 shows the observed stationary magnitudes, the model-computed magnitudes, and the residuals between them, for all satellites. Again, the same parameters results (within 10⁻³) are obtained when the fitting is performed using scaled magnitudes H¹²⁰⁰. Parameters results for the V filter and R filter samples can be seen in Tables A.2 and A.3, respectively, in Appendix A.

6.3 DS versus DSS results

For each stationary magnitude observed and plotted in the left panels of Figures 4 and 5, the results from the DS and DSS models, respectively, are very similar. Both plots look almost the same; however, point by point, they differ typically within ±0.051 magnitudes.

Fig. 4 Results of the DS model applied to all the ONEWEB LEOsats of this study. Left: observed ms (gray points) and DS model estimates ms,ds (purple points) vs. Φ. Right: histogram of residuals ms − ms,ds.

Fig. 5 Results of the DSS model applied to all the ONEWEB LEOsats of this study. Left: observed ms (gray points) and DSS model estimates ms,ds (purple points) vs. Φ. Right: histogram of residuals ms − ms,dss.

Fig. 6 Effective albedo for the diffuse and half-specular spheres. Left: ρeff,ds/ρeff,dss vs. Φ. Center and right: H1200 vs. Φ, color-coded by ρeff,ds and ρeff,dss, respectively. For the same effective albedo, H1200 follows the shape for each model (see Appendix C).

6.4 Effective albedo for the DS and DSS models

Some algebra proves that the ratio ρ_eff,ds/ρ_eff,dss follows a defined curve below 1 for Φ ∼ 84^◦ and greater than 1 for larger values of Φ. This curve is plotted in the left panel of Figure 6. In other words, each observation can be associated to either model’s effective albedo but a diffuse sphere must be less reflective for Φ ≲ 84^◦ and more reflective for Φ ≳ 84^◦ than the half-specular one. What seems at first counterintuitive has an explanation. An observer receives flux from the whole visible illuminated surface of a diffuse sphere because all parts of the illuminated surface reflect in all directions. For the half-specular version, the light it reflects - which for any differential surface area is more than its diffuse counterpart - will reach the observer only if it is properly aligned. The total flux received, of course, also depends on the orientation, as parameterized by Φ. When we observe a fully illuminated satellite face-on (i.e., at full phase or small Φ), the diffuse sphere does not need to be so reflective to explain the observed brightness. When we observe only a small illuminated portion of the satellite (i.e., close to new phase or large Φ), the diffuse sphere must be more reflective to fit the observed magnitude as compared to the half-specular version, for which a smaller illuminated area but precisely oriented reflects as much light.

From Equations (9) and (10), we computed for every single observation its corresponding effective albedo, one for each model. Figure 6 shows the distribution of H¹²⁰⁰ versus Φ, color-coded by the effective albedo ρ_eff,ds (middle panel) and ρ_eff,dss (right panel), respectively. It is easy to notice that the observations with the same effective albedo follow the expected shape for each model (Figures C.2 and C.3.). Also, as expected, at each Sun phase angle Φ, larger effective albedo yield brighter magnitudes. The overall observed trend for most of the range-scale magnitudes H¹²⁰⁰ in Figure 6 is as follows:

• Small Φ (large γ): the satellite is close to full phase; solar panels point to the Sun and also roughly toward the observer. Therefore the satellite is brighter. We observe a lot of reflection from the solar panels, and they have a larger area than the bus.

• Large Φ (small γ): the satellite is close to new phase; solar panels point to the Sun and away from the observer. If it is close enough to the Sun, the satellite is not actually visible. But, when the Sun is below the horizon, and the satellite is at an angular distance γ large enough from the Sun, its bus side toward Earth and its antennas (also pointing toward Earth) may act as a mirror, making the satellite visible and bright. This would explain the small increase in brightness observed for Φ ≳ 100^◦.

This is readily seen in the polar plots of the horizontal coordinates of the satellites and the Sun shown in Figure 7. The satellites’ symbols are colored according to the corresponding Sun phase angle Φ. The left panel shows how the satellites’ range-corrected magnitude is generally brighter when they are opposite the Sun (full phase) and fainter with higher dispersion for the rest of the positions. This happens for both off- and on-station satellites. The increased brightness at full phase is most probably due to specular reflection by the LEOsats solar panels oriented consistently toward the Sun. The center and right panels show the effective albedos, indicated by the symbol size, for the diffuse and the half-specular models, respectively. At full phase or small Φ, the diffuse sphere yields slightly lower effective albedos than the half-specular one. At new phase or large Φ, the opposite occurs, the diffuse sphere yields in average higher effective albedos than the half-specular one. Overall, the diffuse sphere effective albedo shows more visible mean changes in horizontal coordinates than the half-specular one. As proposed by Krantz et al. (2023), the diffuse sphere effective albedo is a measure of the reflection’s specularity, thereby possibly tracing the different attitudes of the ONEWEB satellites. On the other hand, the half-specular sphere effective albedo appears less structured in horizontal coordinates. It can either be closer to the real mean physical albedo of the ONEWEB satellites or the half-specular component absorbs or partially accounts for the specularity changes.

In Appendix C, Figure C.4, we show additional plots with respect to the azimuth coordinates, for comparison with work such as Mallama (2021a,b) and Mróz et al. (2022). The observed trends of H¹²⁰⁰ with ∆Az are explained by the strong correlations in specific intervals of ∆Az with Φ. In our data, we find that the range-corrected magnitude follows a tighter zig-zag pattern with ∆Az for the 2021 observations (October to November, the Sun at sunset to the south), while the 2022 observations (April to September, the Sun at sunset to the north) follow a smoother more disperse oscillation. These two groups correspond to the two clusters of gray symbols in Figure 7. No difference is found in the plot H¹²⁰⁰ versus Φ, although at the low end of Φ, the 2021 observations are, on average, slightly brighter, up to about half a magnitude for Φ < 25^◦ (see Figure C.4). More extended observations are needed to discard a bias effect caused by the Sun’s azimuth interval of the 2021 observations being half the size of that of the 2022 observations. Another possibility is a seasonal and/or geographical variation of the Earth’s albedo playing a role, as discussed in Section 7.

We do notice a few observations in Figure 6 that deviate from what we just described. The faintest observations for which the effective albedo is the lowest, and a nonnegligible number of the brightest observations and largest effective albedos (at Φ ∼ 80^◦ and Φ ∼ 100^◦). The latter appear in Figure 7 as the largest size symbols, more or less following a line from direction NNE to WSW. We checked these subsamples in detail but find nothing remarkable that could explain their extreme values.

Fig. 7 Horizontal coordinates of the ONEWEB LEOsats and the Sun projected on a polar plane. North and east are at 0◦ and 90◦, respectively. Both visible and non-visible sky are projected onto the same plane. The Sun is always below the horizon. Left panel: LEOsats range-corrected magnitude H1200 distribution. The larger the symbol, the brighter the satellite. Center panel: LEOsats diffuse sphere effective albedo distribution. The larger the symbol, the larger the effective albedo ρeff,ds. Right panel: same as center panel, but for the half-specular sphere effective albedo ρeff,dss.

6.5 Relative reflective model – RR

For the RR model, Equation (11) describes the range-corrected magnitude H¹²⁰⁰ as a function of (θ₀, θ₁) and the fixed parameters k, and $H_{(70,0)}^{1200}$ which measures the constant brightness of the satellites at 1 200 km range and fixed orientation (70^◦, 0^◦). The DEMC algorithm finds the values of k and $H_{(70,0)}^{1200}$ for which the RR model reproduces best the magnitudes H¹²⁰⁰. The input data are H¹²⁰⁰ and (θ₀, θ₁), and the resulting parameters can be seen in Table 3, for all satellites, on-, and off-station samples. Figure 8 shows the observed range-corrected magnitudes the model-computed magnitudes, the residuals between them, for all satellites. Parameters results for the V filter and R filter samples can be seen in Tables A.2 and A.3, respectively, in Appendix A.

7 Discussion and conclusion

In all models, we find that the residuals have a dispersion of about 0.6–0.7 magnitudes when having a large sample of data points (σ_demc in Table 2). That amount cannot be explained by the input data errors. The observed m_s are affected only by photometric errors, which are typically about 0.2 magnitudes. Errors in range r are negligible, since the TLE-SGP4-ranges have estimated errors below 1 km (see e.g., Figures 3 to 5 in Vallado & Cefola (2012) and Figure 4 in Geul et al. (2017)).

Therefore, although all the brightness models used can fit the mean value of the observed stationary magnitudes, they still lack other details that account for the rest of the observed residuals. More importantly, we noticed that the neither model reproduces the overall shape of m_s versus Φ (DS, DSS) or H¹²⁰⁰ versus (θ₀, θ₁) (RR) plots. The fact is that both Lambertian models yield fainter magnitudes as Φ increases, but we observe a slight increase in H¹²⁰⁰ for Φ ≳ 100^◦ (see Figure 6).

We therefore conclude that a spherical shape is a first approach to modeling the observed brightness of the ONEWEB LEOsats. However, there is a need for improvement to reproduce the brightness observations within the limits of their errors and the observed trend that seems to exist for larger Φ. In this work, we proved that spherical and RR models, whose magnitude modulation depends fundamentally on the relative position Sun-satellite-observer as parameterized by Φ (or (θ₀, θ₁)) but are independent of the orientation of the body itself (attitude) cannot explain the observations well enough. This means that more information is needed regarding the specific orientation (and size and albedo) of the various surfaces, that are actually reflecting the Sun’s light onto Earth, with respect to the Sun and the observer. Since we know the reflecting surfaces of the ONEWEB satellites are close to flat (including the antennas), a flat reflecting surface model is being considered for future work by the authors. McCue et al. (1971) cites flat surface brightness models from Liemohn (1968) (specular) and Giese (1963) (diffuse). In Maley (1987), the specular flat plate model is cited as the analytical justification for the observed flashes or outburst of brightness for several Cosmos satellites. Those presumed to be tumbling could generate random flashes. ONEWEB satellites examined by this work do not, however, show such behavior.

Another factor that has not been considered so far is Earth’s albedo. Our planet reflects sunlight to space and some of it reaches the satellites, LEOsats being closer are therefore more affected. In other words, the flux received by the satellite is higher than F_⊙. If we apply the formalism explained in Appendix B to Earth as the reflecting sphere, i.e., we model Earth as a Lambertian sphere and express the flux from Earth that reaches the satellite F_⊕ = κF_⊙ with κ > 0, then F = (1 + κ)F_⊙ and the new magnitude will be brighter by −2.5 log(1 + κ). More specifically $$\kappa =A_{B,\oplus }\pi F_{0,\oplus }(z_{\odot ,\text{sat}}){\left(\frac{R_{\odot }}{H_{\text{sat}}}\right)}^2,$$(13)

with A_B,⊕ being the Earth’s Bond albedo and F_0,⊕ the corresponding phase function in terms of the phase angle which in this case is z_⊙,sat. The last term in (13) between ∼25 and ∼120 for our sample of satellites. If Earth were a Lambertian diffuse sphere then its phase function F_0,⊕ would be the DS model in equation (5). For a fixed H_sat (e.g., the on-station sample), such F_0,⊕ decreases with z_⊙,sat, similar to the associated flux on the satellite. Earth flux is greater when the Sun is closer to the satellite’s zenith (smaller z_⊙,sat) and the Earth looks fuller from it.

In our sample, there is anticorrelation between z_⊙,sat and Φ, although with a high dispersion: 125^◦ ≳ z_⊙,sat ≳ 90^◦ when 15^◦ ≲ Φ ≲ 130^◦. The observed mean increase in H¹²⁰⁰ versus Φ for Φ ≳ 100^◦ (z_⊙,sat ≲ 110^◦) could be explained if for these data points, Earth’s albedo becomes important enough to increase the range-corrected magnitudes. In the absence of earthshine H¹²⁰⁰ should keep decreasing. Work studying earth-shine from observations of the Moon at various phases have found that the Earth’s phase function F_0,⊕ is observed to be roughly Lambertian around crescent Moon phases (or gibbous Earth, as seen from the Moon) (Qiu et al. 2003; Robinson et al. 2010; Goode et al. 2021). In such a case, the Earth’s apparent albedo is roughly comparable to the Earth’s visible light albedo (∼0.3) (Goode et al. 2021). It is important to note that Earth appears Lambertian at gibbous phases not because it truly is Lambertian, but because thick cloud cover and geometric averaging suppress anisotropic scattering, making the disk-integrated reflectance mimic that of a diffuse sphere. On the other hand, crescent and full Earth disk-integrated brightnesses deviate significantly from the Lambertian model. Forward scattering by clouds and aerosols, as well as ocean specular glint can dominate total flux for crescent phases, while coherent backscatter and shadow hiding (especially over clouds and ice) are important in the full phase (Robinson et al. 2010; Goode et al. 2021). Moreover, the long-term monitoring of earthshine (e.g., Goode et al. 2021) has found temporal variation of the Earth albedo within the Lambertian behavior tied to the seasons and climate cycles.

Another important factor is that this approximation holds for distant observers, the Moon is far enough from Earth but LEOsats are not. For a distant observer all the surface elements of the Earth’s sphere can be assumed to be at the same phase angle, while for a LEOsat, different surface elements on Earth have measurably different phase angles. In any case, a Lambertian diffuse and specular Earth may be considered as a first approach to model the Earth albedo effect on the observed brightness of LEOsats.

Fig. 8 Results of the RR model applied to all the ONEWEB LEOsats of this study. Left and center: observed H1200 (gray points) and DSS model estimates $H_{\text{rr}}^{1200}$ (purple points) vs. θ0 and θ1, respectively. Right: histogram of residuals $H^{1200}-H_{\text{rr}}^{1200}$.

Data availability

The full content of Table 1 can be found at https://zenodo.org/records/19204696.

Acknowledgements

Author M.R.C. acknowledges the Universidad de Atacama PhD Scholarship; and the Air Force Office of Scientific Research (AFORS) grant # FA9550-22-1-0292.

References

Appendix A Observations and results tables for V filter and R filter samples

Appendix B Magnitudes and phase function

From the definitions of magnitudes and fluxes, we want to compare the flux we receive from a sphere F_sph to the flux we receive from the Sun, then $$m_{\text{sph}}-m_{\odot }=-2.5\log \left(\frac{F_{\text{sph}}}{F_{\odot }}\right).$$

We start with the definitions of geometric albedo p and phase function F(Φ), as in Shepard (2017, Section 5.5): $$\begin{array}{cc}p& =\frac{I(0)}{F{R}^2},\\ F(\text{Φ})& =\frac{I(\text{Φ})}{I(0)},\end{array}$$

where F is the flux received by the sphere, R is its radius, I is the radiant intensity from the sphere and Φ is the phase angle of that flux. As noted by Shepard (2017, Section 5.5.2), the geometric albedo compares the reflectance of the sphere to that of a flat disk made of Lambertian material with the same radius R, both seen face-on. I(Φ) is the amount of energy we receive from (reflected by) the sphere per unit time and unit solid angle. Therefore: $$F_{\text{sph}}(\text{Φ})=\frac{I(\text{Φ})}{r^2}=\frac{I(0)F(\text{Φ})}{r^2}=\frac{pF{R}^{2}F(\text{Φ})}{r^2}.$$

From Shepard (2017, Section 5.5.4), we have the definitions of phase integral q and Bond (or spherical) albedo A_B: $$\begin{array}{cc}q& =\frac{{{\displaystyle \int }}_0^{4\pi }F(\text{Φ})d\text{Ω}}{\pi },\\ A_B& =p\cdot q,\end{array}$$

from which we obtain that $p=\frac{A_B\pi }{{{\displaystyle \int }}_0^{4\pi }F(\text{Φ})d\text{Ω}}.$. Substituting in the equation for the sphere flux above, we have $$\begin{array}{l}{F}_{\text{sph}}(\text{Φ})=\frac{A_B\pi }{{{\displaystyle \int }}_0^{4\pi }F(\text{Φ})d\text{Ω}}\frac{F{R}^{2}F(\text{Φ})}{r^2}\hfill \\ =\frac{A_{B}F\pi R^2}{r^2}\frac{F(\text{Φ})}{{{\displaystyle \int }}_0^{4\pi }F(\text{Φ})d\text{Ω}}\hfill \\ =\frac{A_{B}F\pi R^2}{r^2}{F}_0(\text{Φ}).\hfill \end{array}$$

F₀(Φ) is referred to by Williams & McCue (1966) as the phase function normalized by the integral over a solid angle of 4π stereoradians. Substituting F = F_⊙, i.e., the sphere is illuminated by the Sun at 1 AU, we obtain $$\begin{array}{lll}\frac{F_{\text{sph}}(\text{Φ})}{F_{\odot }}\hfill & =\hfill & \frac{A_B\pi R^2{F}_0(\text{Φ})}{r^2}\Rightarrow \hfill \\ m_{\text{sph}}-m_{\odot }\hfill & =\hfill & -2.5{\log }_{10}\left(\frac{F_{\text{sph}}(\text{Φ})}{F_{\odot }}\right)\Rightarrow \hfill \\ m_{\text{sph}}\hfill & =\hfill & m_{\odot }-2.5{\log }_{10}\left[A_B\pi R^2{F}_0(\text{Φ})\right]+5{\log }_{10}(r).\hfill \end{array}$$

In equation (4): m_s = m_sph, ρ = A_B, A = πR² is the cross section of the sphere, and ζ = ρ A. For a Lambertian sphere, the Bond Albedo is equal to the albedo (Shepard 2017, equation 5.71). Keep in mind that the units of A have to be compatible with the units of r. In other words, if the sphere cross section is given in square meters, the range has to be given in meters (Karttunen et al. 2003).

Appendix C Additional figures

Figures C.1, C.2 and C.3 contain complementary information.

Fig. C.1 ONEWEB LEOsats of this investigation, separated by their on-station (orange) and off-station (blue) orbital stages based on their altitude. The vertical axis shows their range, i.e., distance between the observer and the satellite. Off- and on-station satellites have different altitudes, but they may have similar ranges.

Fig. C.2 Plots of ms for the DS model, considering different values of p, ρ, and A, from left to right panels, respectively, as indicated by the legends. All plots consider r = 1 200 km. Left: A = 1 and ρ = 0.5. Middle: p = ρ = 0.5. Right: p = 0.5 and A = 1.

Fig. C.3 Plots of ms for the DSS model, considering different values of β, ρ, and A, from left to right panels, respectively, as indicated by the legends. All plots consider r = 1 200 km. Left: A = 1 and ρ = 0.5. Middle: β = ρ = 0.5. Right: β = 0.5 and A = 1.

Fig. C.4 Variation of H1200 with azimuth. In all figures, points are color-coded by the azimuth of the Sun: north is red, south is green. Upper left: H1200 vs. Φ, as in Figure 6, for reference. Lower left: ∆Az = Azsat − Az⊙ vs. Φ. Lower right: ∆Az vs. H1200. Upper right: Az⊙ vs. Azsat.

Appendix D Glossary of terms

Glossary

Nomenclature

m_s: stationary magnitude.

m₀: zero-point magnitude of the stars.

m_ins: instrumental aperture magnitude.

ϕ: angular velocity satellite in the sky.

L: length of the line left by the LEOsat.

t_exp: exposure time.

m_⊙: Sun’s apparent magnitude.

H¹²⁰⁰: scaled magnitude, range-corrected magnitude.

DS: Diffuse Spherical (model).

DSS: Diffuse Specular Spherical (model).

RR: Relative Reflectance (model).

ℛ: ratio of fluxes in RR model.

F₀: phase function.

ζ: effective area of the sphere reflecting light.

r: range of the satellite (satellite-observer distance).

ρ: albedo.

ρ_eff,ds: effective albedo of DS model.

ρ_eff,dss: effective albedo of DSS model.

A: πR² = sphere cross-area.

β: mixing coefficient of diffuse-specular reflection.

$F_{{\theta }_0,{\theta }_1}(r)$: flux coming from direction θ₀, θ₁ at a distance r.

k: edge darkening parameter.

Φ: Sun phase angle.

θ₀: Sun incidence angle or solar incidence angle.

θ₁: observer angle.

γ: elongation angle.

$H_{({\theta }_0,{\theta }_1)}^{1200}$: scaled magnitude in orientation (θ₀, θ₁).

$H_{(70,0)}^{1200}$: scaled magnitude in orientation (70^◦, 0^◦).

∆Geog: angular distance between observer’s longitude and latitude and the satellite’s longitude and latitude.

∆Az: Az_sat − Az_⊙ = azimuth of the satellite minus the azimuth of the Sun.

z_⊙,obs: Sun’s zenith distance for the observer.

z_⊙,sat: Sun’s zenith distance for the satellite.

O_⊕: Sun’s Position.

O_sat: satellite’s position.

O_obs: observer’s position.

R_⊕: Earth’s radius.

H_sat: perpendicular height above Earth of the satellite.

σ_tle: angular distance between TLE-predicted and observed sky position of the satellite.

All Tables

All Figures

	Fig. 1 Distribution of stationary magnitudes of all ONEWEB satellites studied. The sample is divided by filter and orbit type. There are 451 on-station satellites, of which 410 were observed with the V filter and 41 with the R filter. There are 175 off-station satellites observed 158 with the V filter and 17 with the R filter. The average µ and standard deviation σ of ms are shown in the legend.
In the text

	Fig. 2 Graphical representation of the position of the satellite and definitions of the Sun incidence angle θ0 (in orange), the observer angle θ1 (in dark blue), the Sun phase angle Φ (in green), and the solar elongation angle γ (in brown). Notice that θ0, θ1 and Φ are not contained in the same plane but form a triangular pyramid. The distance is not to scale.
In the text

	Fig. 3 Histogram of σtle for all the ONEWEB satellites studied, showing the difference between the TLE-predicted and the observed equatorial coordinated of satellites.
In the text

	Fig. 4 Results of the DS model applied to all the ONEWEB LEOsats of this study. Left: observed ms (gray points) and DS model estimates ms,ds (purple points) vs. Φ. Right: histogram of residuals ms − ms,ds.
In the text

	Fig. 5 Results of the DSS model applied to all the ONEWEB LEOsats of this study. Left: observed ms (gray points) and DSS model estimates ms,ds (purple points) vs. Φ. Right: histogram of residuals ms − ms,dss.
In the text

	Fig. 6 Effective albedo for the diffuse and half-specular spheres. Left: ρeff,ds/ρeff,dss vs. Φ. Center and right: H1200 vs. Φ, color-coded by ρeff,ds and ρeff,dss, respectively. For the same effective albedo, H1200 follows the shape for each model (see Appendix C).
In the text

Fig. 7 Horizontal coordinates of the ONEWEB LEOsats and the Sun projected on a polar plane. North and east are at 0◦ and 90◦, respectively. Both visible and non-visible sky are projected onto the same plane. The Sun is always below the horizon. Left panel: LEOsats range-corrected magnitude H1200 distribution. The larger the symbol, the brighter the satellite. Center panel: LEOsats diffuse sphere effective albedo distribution. The larger the symbol, the larger the effective albedo ρeff,ds. Right panel: same as center panel, but for the half-specular sphere effective albedo ρeff,dss.

In the text

	Fig. 8 Results of the RR model applied to all the ONEWEB LEOsats of this study. Left and center: observed H1200 (gray points) and DSS model estimates $H_{\text{rr}}^{1200}$ (purple points) vs. θ0 and θ1, respectively. Right: histogram of residuals $H^{1200}-H_{\text{rr}}^{1200}$.
In the text

	Fig. C.1 ONEWEB LEOsats of this investigation, separated by their on-station (orange) and off-station (blue) orbital stages based on their altitude. The vertical axis shows their range, i.e., distance between the observer and the satellite. Off- and on-station satellites have different altitudes, but they may have similar ranges.
In the text

	Fig. C.2 Plots of ms for the DS model, considering different values of p, ρ, and A, from left to right panels, respectively, as indicated by the legends. All plots consider r = 1 200 km. Left: A = 1 and ρ = 0.5. Middle: p = ρ = 0.5. Right: p = 0.5 and A = 1.
In the text

	Fig. C.3 Plots of ms for the DSS model, considering different values of β, ρ, and A, from left to right panels, respectively, as indicated by the legends. All plots consider r = 1 200 km. Left: A = 1 and ρ = 0.5. Middle: β = ρ = 0.5. Right: β = 0.5 and A = 1.
In the text

	Fig. C.4 Variation of H1200 with azimuth. In all figures, points are color-coded by the azimuth of the Sun: north is red, south is green. Upper left: H1200 vs. Φ, as in Figure 6, for reference. Lower left: ∆Az = Azsat − Az⊙ vs. Φ. Lower right: ∆Az vs. H1200. Upper right: Az⊙ vs. Azsat.
In the text