© The Authors 2025

1 Introduction

Inferring the shape of a Solar System body from its rotational light curve is a well-established technique (e.g., Ostro & Connelly 1984; Ostro et al. 1988; Kaasalainen & Torppa 2001; Kaasalainen et al. 2001; Durech et al. 2010 - DAMIT¹), typically relying on observations obtained over a range of viewing geometries and with sufficient signal-to-noise ratio. However, most inversion tools developed for this purpose are based on the assumption that the target body is convex, meaning that its surface lacks significant concavities or selfshadowing features. This assumption is clearly violated in the case of irregular or bilobate bodies such as 67P/Churyumov-Gerasimenko, (216) Kleopatra, or Arrokoth (Preusker et al. 2017; Shepard et al. 2018; Keane et al. 2022). Moreover, accurate modeling often depends on the availability of high-quality data collected from multiple observing geometries to constrain the 3D shape effectively. However, these conditions are rarely met for distant bodies such as Centaurs or trans-Neptunian objects (TNOs). Because of their great distance, the aspect angle remains nearly constant over time. In addition, their faintness often leads to low signal-to-noise ratios in the light curves. As a result, traditional shape-inversion techniques become inapplicable or unreliable in this context.

In this work, we adopt a different approach. Instead of inferring the shape from the light curve, we start from a predefined 3D shape model and compute its projection onto the sky plane and the rotational light curve, as would be observed from Earth, using its known orbital parameters and spin axis orientation (when available).

This tool was primarily designed to support the analysis of data obtained from stellar occultations, where the projected silhouette of a Solar System body on the sky plane can be directly constrained. To this end, the application incorporates a dedicated module that computes both the rotational phase and the projected silhouette at the epoch of the occultation, enabling direct comparison with observed chords. This functionality proves especially useful for visualizing and analyzing the set of rotational solutions that are consistent with the occultation profile. A recent application of this approach can be found in Rizos et al. (2024). Nevertheless, the tool can also be employed in other planetary science contexts, such as studying the rotational properties and photometric behavior of minor bodies through light-curve observations.

While similar efforts have been presented in previous studies - see DAMIT (Durech et al. 2010) or Showalter et al. (2021) - our implementation introduces several novel features. First, it supports arbitrary non-convex 3D shapes with a high degree of topographic complexity, including user-supplied shape models. Second, our tool allows for the use of multiple photometric scattering models to compute the reflected light. Users can choose between empirical approaches, such as the Lommel-Seeliger or Minnaert laws, and more physically grounded models, such as the Hapke formulation. Additionally, our application includes advanced options, such as albedo variegation, axial precession, and self-shadowing, which further increase the physical realism of the simulations. The user can also upload their own observational data for direct comparison with simulated outputs. Another significant advantage of our tool is its accessibility: it is entirely web-based and requires no local installation. The user only needs to enter the desired parameters, with the option to upload their own 3D shape model and/or a rotational light curve for comparison. The server then performs the calculations, which are displayed in the interface and can be downloaded. Hosted on the servers of the Instituto de Astrofísica de Andalucía (IAA-CSIC), the application can be run remotely through an intuitive browser interface². To guide new users, a brief online manual is available under the title How does this work? in the web interface, which provides practical instructions for basic usage without requiring a full reading of this article.

The structure of this paper is as follows. In Section 2, we describe the overall system architecture, including the main components of the web application. Section 3 outlines the functionalities offered by the tool and explains how users can interact with the different simulation options. In Section 4, we present validation cases, comparing simulations with real observational data. Finally, Section 5 summarizes the conclusions.

2 System architecture and implementation

The tool is implemented as a web application built using the Django framework³, a high-level Python-based web platform designed for rapid development and clean design. The architecture follows the Model-View-Template (MVT) pattern, which allows for a clear separation between the application logic, data structure, and presentation layer.

The interface is based on a series of HTML templates, which are dynamically rendered through Django views. These views serve as controllers that handle user input, validate parameters, and initiate the corresponding simulation processes on the server side. Once a simulation is executed, the results are returned to the front-end and displayed to the user.

The backend computation is entirely written in Python. The user provides input parameters such as the shape model, spin axis orientation, or scattering law through the web interface, and all operations are handled server-side, ensuring consistency and reproducibility.

The application is deployed on a dedicated Linux server hosted at the Instituto de Astrofísica de Andalucía (IAA-CSIC). The system is optimized for consistent performance and direct access to hardware resources, without the use of virtualization or containerization layers.

The interaction between the user and the web application follows a well-defined workflow, as illustrated in Figure 1. The process begins when the user accesses the web interface and enters the required parameters. Upon submission, the data are passed to a Django view, which first performs a validation step (e.g., verifying that the object identifier, latitude, and longitude are within acceptable ranges). If any inconsistency is detected, the view returns an error message to the user through the same web interface, highlighting the problematic field(s).

If the input passes validation, the view performs automated queries to the JPL Horizons system⁴ through its public batch interface. The tool retrieves the geometric parameters for the specified date without accounting for the duration of the rotation period. For most small bodies, this simplification is valid, but for very slow rotators (e.g., Apophis or Toutatis), users should simulate multiple dates to capture geometry changes. Once these data are retrieved, the server-side Python routines perform the relevant calculations: projection of the 3D model onto the sky plane, application of the selected photometric scattering law, and generation of the expected light curve.

The results - such as plots of the simulated light curve, skyplane projection snapshots, and downloadable data files - are then passed back to the HTML template. They are subsequently rendered in the user’s browser.

Fig. 1 Workflow of the web application from user input to result rendering.

3 Tool features and user configuration

The web application provides two operational modes specifically designed for photometric studies of rotational light curves (the Known and Generic body modes). In addition, it includes an Occultation mode, which projects the silhouette of a selected body at a given epoch, enabling the analysis and planning of stellar occultations.

3.1 Rotational light curves

• Known body mode: This mode is designed for the simulation of Solar System objects with a known rotation pole⁵ and orbital parameters. The user enters the Minor Planet Center (MPC) designation and the observation date. The tool generates both the rotational light curve and the projected view sky-plane as seen from Earth using the ephemeris data from the JPL Horizons system.

• Generic body mode: This mode enables simulations under user-defined geometric conditions. The observing and illumination geometries are specified directly by the user, including the phase angle, and the spin axis is defined within a custom reference frame. This is particularly useful for testing hypotheses or exploring different subsolar and subobserver configurations. It offers full flexibility for cases where the viewing geometry changes significantly over the course of a single rotation.

In both modes, users can provide a custom shape model in .obj format (with triangulated facets) or define a triaxial ellipsoid by specifying axes a, b, and c. Shape models should be oriented such that the forward axis is Y and the up axis is Z. This ensures that the semiaxes a, b, and c align with the X, Y, and Z axes, respectively. Users are encouraged to limit models (when possible) to fewer than 3000 facets for optimal performance. The internal ellipsoid is generated with 642 facets by default.

Several additional simulation features are available to users:

• Albedo spot: enables the simulation of localized surface heterogeneities by allowing the user to define a specific surface region—based on the body’s geocentric latitude and longitude—whose radiance factor (RADF) is modified by a user-defined percentage.

• Precession: simulates the motion of the rotation pole as it precesses over time, instead of remaining fixed. The pole (λ and β_c represent longitude and co-latitude, respectively) follows a conical trajectory, varying sinusoidally as a function of the rotational phase. This behavior is controlled by user-defined parameters as follows: $$\lambda (t)={\lambda }_0+A_{\lambda }\cos (2\pi {\omega }_{p}t),$$(1a) $${\beta }_c(t)={\beta }_{c0}+A_{{\beta }_c}\sin (2\pi {\omega }_{p}t),$$(1b)

  where λ₀ and β_c0 are the initial orientation of the rotation pole, A_λ and Aβ_c are the amplitudes of the precessional motion, ω_p is the precession frequency (in cycles per rotation), and t is the rotational phase (ranging from 0 to 1). This option enables the simulation of dynamic spin states and allows for the exploration of more realistic or perturbed rotational scenarios.

• Self-shadowing: enables the simulation of shadows cast by the object’s own surface features, enhancing the physical realism of the reflected light distribution. For each facet, the algorithm checks whether the direction of incoming sunlight intersects any other facet of the shape model, thus determining whether the facet is illuminated or in shadow. This is implemented using the Möller-Trumbore ray-triangle intersection algorithm Möller & Trumbore (1997), which efficiently tests whether a ray defined by the facet centroid and the Sun’s direction vector intersects any of the mesh triangles. While this process increases computational time, it is essential for accurately modeling objects with significant concavities or complex topography.

Photometric modeling

In this tool, the RADF - representing the ratio of the observed bidirectional reflectance to that of a perfectly diffuse surface illuminated at i = 0° (Hapke 2012) - can be computed in several ways. We provide both empirical formulations and physically based models, allowing users to choose the level of complexity and realism appropriate for their application. In all cases, RADF depends on three photometric angles: the phase angle α (Sun-object-observer), the incidence angle i (Sun-object surface normal), and the emission angle e (observer-object surface normal).

Notice that this tool does not account for the object’s distance to the Sun; instead, the solar constant is normalized to S = 1. Therefore, if absolute RADF values are needed, users must rescale the output accordingly. The value output by the server corresponds to the integrated RADF as seen from Earth, computed by summing the contribution from each individual facet, as detailed below.

Empirical models. For empirical models, the RADF is approximated as the product of a phase function f (α) and a disk function d(α, e, i) as follows: $$RADF=f(\alpha )d(\alpha ,e,i).$$(2)

Phase functions

• Constant phase function: We provide this simplest case for users, who do not need to account for any dependence on the phase angle, f (α) = 1.

• Minnaert phase function: This is a generalization of Lambert’s law, suggested by Minnaert (1941) with the modification suggested by Takir et al. (2015): $$f(\alpha )=A_{\text{min}}\pi {10}^{-0.4(\beta \alpha +\gamma {\alpha }^2+\delta {\alpha }^3)},$$(3)

  where A_min, β, γ, and δ are user-defined parameters.

Disk functions

• Lambert: This model describes a perfectly diffuse (Lambertian) surface, which reflects incident light equally in all directions (Lambert 1760). The reflected intensity is proportional to the cosine of the incidence angle, representing the projection of the incident flux onto the surface. The division by π ensures energy conservation when integrating the radiance over the entire hemisphere. Although idealized, the Lambertian model is commonly used due to its simplicity and serves as a useful reference in photometric analyses: $$d(i,e)=\frac{\cos i}{\pi }.$$(4)

• Lommel-Seeliger: the Lommel-Seeliger law is based on a simple physical model of diffuse reflection under the assumption of single, isotropic scattering. It remains a valuable first-order approximation in planetary photometry and is still widely used in applications such as light curve inversion and shape modeling of minor bodies (Fairbairn 2005). Unlike the Lambertian model, this formulation introduces an angular dependence on both incidence (i) and emission (e) angles, and it is particularly useful for modeling dark planetary surfaces or regolith-covered bodies such as asteroids or moons: $$d(i,e)=\frac{\cos i}{\cos i+\cos e}.$$(5)

• Minnaert: this disk function generalizes Lambert’s law by introducing the empirical exponent k, which accounts for nonuniform scattering characteristics and modulates the dependence on incidence (i) and emission (e) angles. The formulation adopted here is consistent with that introduced by Minnaert (1941): $$d(i,e)={\cos }^{k}i{\cos }^{k-1}e.$$(6)

Hapke model. In addition to empirical formulations, we implement the Hapke model (Hapke 2012), which accounts for multiple scattering, shadow hiding, coherent backscatter, and macroscopic roughness. This model depends on nine physical parameters:

• w: single scattering albedo.

• b_HG, c_HG: Henyey-Greenstein parameters for the particle phase function.

• B_C0, h_c: coherent backscatter amplitude and width.

• B_s0, h_s: shadow hiding amplitude and width.

• θ_p : macroscopic roughness slope angle.

• φ: porosity factor.

The RADF is computed internally using the full formulation. For a complete detailed description of all the equations see Appendix A.

3.2 Occultation module

This mode builds directly on the functionalities described above but is specifically designed for a direct application in the context of stellar occultations. Given the MPC designation of a target body, its sidereal rotation period (in hours), and the orientation of its spin axis, the tool computes the rotational phase for a user-specified date at which a stellar occultation is predicted, using the reference epoch t₀ (in Julian Date) corresponding to a photometric rotational light curve maximum. As in the other modes, users may also upload a custom shape model in .obj format or define the principal axes of an ellipsoid with 642 facets by default.

A major advantage of this tool is that it enables the exploration of viewing-geometry variations across different dates. Since the on-screen output is rendered to scale, users can empirically determine the expected dimensions along any given direction and directly compare them with the observed occultation chords.

3.3 Output

The web interface provides a complete summary of all the parameters used in each simulation. If the Generate images checkbox is selected for Known and Generic modes, an animated GIF is displayed showing the projected shape of the object on the sky plane. For the Occultation mode, the interface also displays the rotational phase of the object at the inserted occultation date and the corresponding projected figure. The colors of the surface facets reflect the solar incidence angle, with red indicating lower values (more irradiated areas) and blue indicating higher ones. The white arrow corresponds to the orientation of the rotation pole at the time of the observation. A scale bar is shown in the bottom-right corner. The default units of the tool are kilometers.

In the Known and Generic modes, we can see the synthetic light curve in relative magnitude. This is computed as $$\text{Relative magnitude}=-2.5{\log }_{10}\left(\text{RADF}/\overline{\text{RADF}}\right),$$(7)

where RADF is the sum of theradiance factor from all facets at a given rotational phase and RADF is the average RADF.

By clicking the Download files button, a ZIP archive is generated and downloaded to the user’s local machine. This archive contains:

• All generated images, each named with the corresponding rotational phase.

• The synthetic light curve as synthetic_LC.png.

• The shape model in .obj format.

• The animated projection GIF (when Generate images checkbox was selected).

• A text file named output.dat with the computed numerical results.

In the output.dat file, the first column lists the rotational phase, the second the RADF, and the third the relative magnitude as defined in Eq. (7).

Similarly, in the case of the Occultation mode, the user can download a ZIP file containing the projected image for the occultation date and the shape model in .obj format. The package also includes the output.dat file with all the data corresponding to this case.

4 Validation

4.1 Projected shape

First, we validated that the projected shape generated by the tool reproduces the actual silhouette of an object at a specific date using the Occultation mode. For this test, we selected the trans-Neptunian dwarf planet (136108) Haumea, since projected occultation chords on the sky plane are available for January 21, 2017, when the body was observed near the absolute brightness minimum of its rotational phase (Ortiz et al. 2017). Haumea’s triaxial shape is well constrained (a = 1161 ± 30 km, b = 852 ± 4 km, c = 513 ± 16 km), as well as its spin-axis orientation (α,δ) = (285.1° ± 0.5°, −10.6° ± 1.2°, J2000 equatorial), and rotation period (3.915341 ± 0.000005 h) (Ortiz et al. 2017). From our light-curve database, we find a brightness maximum (t₀) at 2016-12-19 20:54:04.554 UTC, corresponding to JD 2457742.37088604.

The spin-axis orientation was then converted into ecliptic coordinates, as required by the tool (lon = 285.17°, lat = 12.05°), and we generated the projection for the 2017 occultation date. It yields Fig. 2, with a rotational phase of 0.75 from the light curve maximum. We then extracted the stellar occultation timings from Ortiz et al. (2017), projected the corresponding chords onto the sky plane, and overplotted them on the figure (pink and yellow lines; see also Fig. 2 in Ortiz et al. 2017). This comparison shows the good agreement between the projected shape model and the observational data, supporting the consistency of our results.

4.2 Rotational light curve

To ensure that this validation is meaningful, it is necessary to have a highly accurate 3D shape model of the object, a well-characterized photometric model, and observed light curves from Earth for comparison. For these reasons, we first selected asteroid (101955) Bennu as our validation target. This asteroid was orbited by the OSIRIS-REx spacecraft (Lauretta 2015), which provided a high-resolution shape model. For this analysis, we adopt Bennu’s shape model version 20, which features a spatial resolution of 75 cm and comprises approximately 200 000 triangular facets. The model is distributed in .obj format⁶.

The mission also enabled a well-constrained characterization of Bennu’s photometric properties for the four color filters of the OCAMS MapCam instrument: b′, v, w, and x, centered at 473, 550, 698, and 847 nm, respectively (Rizk et al. 2018). These filters were designed to match the passbands of the EightColor Asteroid Survey (ECAS) system (Tedesco et al. 1982). Conveniently, Bennu was also observed with the ECAS w filter using the Kuiper 1.54 m reflector telescope in Arizona, USA (Hergenrother et al. 2013). A rotational light curve was constructed by (Hergenrother et al. 2013) from observations conducted between 14 and 17 September 2005. Accordingly, we simulated a synthetic rotational light curve for the same filter to ensure consistency.

The rotation pole of Bennu is located at α = 85.65°, δ = −60.17° (Barnouin et al. 2019), which corresponds to ecliptic coordinates of lon = 71.11°, lat = −83.32°. Since Bennu is a near-Earth asteroid (NEA), its viewing geometry changes significantly over the course of the four-day observational campaign. The solar phase angle for each observed day was 58.8°, 61.9°, 65.4°, and 69.4°, respectively. Therefore, it is necessary to take into account not only rotational variations but also phase angle changes. For this reason, we adopted the Minnaert photometric model, which incorporates both a phase function and a disk function, and whose parameters have been well constrained by Golish et al. (2021), as summarized in Table 1.

To reduce computation time, we downsampled the number of facets in the shape model by 10%. Figure 3 shows the projected shape of Bennu on 15 September 2025.

In Figure 4, we present the synthetic light curves generated for each day, with a sampling resolution of one data point per degree of rotational phase. To facilitate comparison, we extracted the observational data from the light curve published by Hergenrother et al. (2013) using the Automeris⁷ tool and plotted them together. An offset was applied to the synthetic curves for clarity.

The first result that highlights the usefulness of this tool is the identification of scatter in the observational light curve caused by the combination of multiple light curves affected by changes in viewing geometry. This dispersion could otherwise be mistakenly attributed to photometric noise if not properly accounted for.

It should be noted that, due to the complexity of this case, matching the observed photometric range required scaling the synthetic light curve amplitude by a factor of 1.5 - a correction that is justified for several reasons:

First, the photometric models provide a simplified description of surface scattering. It is well-known that photometric models tend to perform poorly for large values of phase, emission, or incidence angles (Rizos et al. 2021). In particular, regions near the limb, combined with limitations in the shape model, are expected to introduce biases. Furthermore, Bennu’s surface is covered with numerous boulders, some of which are of exogenous origin (Tatsumi et al. 2021), exhibiting albedo variations of up to 15% (Dellagiustina et al. 2019). Each of these surface elements may contribute to distinct photometric behaviors, which are not captured by the simulation. Finally, the shape model is inherently limited and does not fully capture the complexity of Bennu’s surface (see Fig. B.1). To assess the sensitivity of the results to the shape model resolution, we repeated the same experiment by further downgrading the shape model by another 10% (using instead only 2000 facets). We observed a peak-to-peak amplitude variation of approximately 20%, underscoring the limitations imposed by model resolution.

We selected a second, less demanding example: asteroid (433) Eros. This object was observed by Aznar Macias et al. (2018) on July 17, 18, 19, and 24, 2016, from the Aznar Observatory, located at the Centro Astronómico del Alto Turia in Aras de los Olmos, Spain. During this period, the solar phase angle varied from 24.0° to 20.7°. Observations were conducted using a 0.36 m Schmidt-Cassegrain telescope equipped with a CCD camera and broad-band Astrodon Sloan r and J-V filters. The light curves were obtained over multiple nights and reduced to a common reference phase angle of 24.0°. For comparison, we simulated a rotational light curve of Eros on July 17, 2016, when the observed phase angle matched 24.0° (the spin axis of Eros in ecliptic coordinates is located at lon = 17.23°, lat = 11.34° (Yeomans et al. 2000). Since no specific photometric model has been derived for Eros using these particular filters (it requires multi-epoch observations covering a wide range of phase and illumination angles with a consistent filter set, which is not the case here), we applied only disk function corrections using the Minnaert model, adopting a coefficient of k = 0.7, a value lying within the typical range reported for Solar System bodies (e.g., Li 2005; Jehl et al. 2008; Li et al. 2013; Masoumzadeh et al. 2015). The shape model of Eros⁸ (Gaskell 2021) was downsampled to approximately 20000 facets to match the resolution used in the Bennu test case and reduce computation time. Fig 5 shows the orientation of Eros as seen from Earth during the simulation, and Fig. 6 displays the synthetic light curve together with the observational data extracted using the Automeris tool. In this simpler case, the comparison did not require any scaling of the light curve amplitude.

Fig. 2 Haumea on January 21, 2017, as seen from Earth. On this date, Ortiz et al. (2017) reported observations from multiple Earth-based observatories of Haumea passing in front of a distant star (a multi-chord stellar occultation). These observations allowed them to constrain the shape of Haumea’s projected ellipse on the sky. The pink lines represent the projected occultation chords derived from the reported timings, while the yellow endpoints indicate the associated uncertainties in the ingress and egress times. The facet colors on the projected body indicate the relative level of irradiance received by each facet, with red corresponding to lower incidence angles and therefore higher irradiation. The overall agreement between the projected shape and the occultation chords supports the validity of our simulation.

Fig. 3 Projected view of Bennu, generated with an approximately 20000-facet downsampled shape model, as seen from Earth on 15 September 2005 (phase angle of 61.9°).

Fig. 4 Observational data from the light curve published by Hergenrother et al. (2013) shown with black dots. Colored lines show the synthetic rotational light curves of Bennu as seen from Earth on 14-17 September 2005, sampled at one point per degree of rotation. An offset is applied to the synthetic light curves for clarity.

Fig. 5 Projected view of Eros, generated with an approximately 20 000-facet downsampled shape model, as seen from Earth on July 17, 2016 (phase angle of 24.07°).

Fig. 6 Observational data from the light curve published by Aznar Macias et al. (2018) represented by black dots. The green line shows the synthetic rotational light curve of Eros as seen from Earth on July 17, 2016, sampled at one point per degree of rotation.

5 Conclusions

We presented a web-based simulation tool for generating projected shapes and synthetic rotational light curves of small Solar System bodies. By leveraging high-resolution shape models, validated photometric laws, and realistic viewing geometries, the tool enables users to explore the physical properties and photometric behavior of asteroids, moons, and comets in a fully customizable framework.

As a validation case, we reproduced the silhouette of (136108) Haumea, whose shape is well constrained by multichord occultations, obtaining very good agreement. Then, we modeled the rotational light curves of (101955) Bennu and (433) Eros. Bennu is an almost spherical body (more precisely, diamond-shaped), with a light curve amplitude Δm below 0.16 magnitudes. Therefore, this object constitutes a particularly stringent benchmark, pushing the tool to its sensitivity limits. The results demonstrate that the tool can reproduce subtle photometric variations under varying geometries, while also highlighting the limitations imposed by surface complexity and photometric modeling assumptions.

In summary, this tool provides a flexible and extensible framework for photometric modeling, with a focus on accessibility, scientific rigor, and future adaptability. It is our hope that it will contribute to the analysis and interpretation of small body observations across a wide range of planetary science applications.

Acknowledgements

We sincerely thank the anonymous referee for their constructive suggestions, which helped improve the clarity and quality of this manuscript. J.L. Rizos, P.J. Gutierrez, I.M. Navajas and L.M. Lara acknowledge financial support from project PID2021-126365NB-C21 (MCI/AEI/FEDER, UE) and from the Severo Ochoa grant CEX2021-001131-S funded by MCIN/AEI/10.13039/501100011033. J.L. Rizos acknowledges support from the Ministry of Science and Innovation under the funding of the European Union NextGeneration EU/PRTR. J.L. Ortiz acknowledges the Spanish projects PID2020-112789GB-I00 (AEI) and Proyecto de Excelencia de la Junta de Andalucía PY20-01309.

References

Appendix A Hapke photometric model equations

The Hapke model (Hapke 2012) is based on radiative transfer theory and provides a more physically grounded framework for describing surface scattering. It accounts for single and multiple scattering, surface roughness, and opposition effects. The RADF in our implementation is computed as $$RADF=LS(i,e)K\frac{w}{4}[p(\alpha )(1+B_{S0}{B}_S(\alpha ))+M(i,e)][1+B_{C0}{B}_C(\alpha )]S(i,e,\psi )$$(A.1)

Here:

• LS (i, e) is the Lommel-Seeliger function, $\frac{\cos i}{\cos i+\cos e}$

• K is the porosity factor, $\frac{-\ln (1-1.209{\varphi }^{2/3})}{1.209{\varphi }^{2/3}}$, with φ the filling factor (1 - porosity)

• p(α) is the phase function, for which we use the Henyey-Greenstein form: $$\begin{array}{rlr}p(\alpha )& =\frac{(1+c_{HG})}{2}\frac{(1-b_{HG}^2)}{(1-2b_{HG}\cos \alpha +b_{HG}^2{)}^{3/2}}& +\frac{(1-c_{HG})}{2}\frac{(1-b_{HG}^2)}{(1+2b_{HG}\cos \alpha +b_{HG}^2{)}^{3/2}}\end{array}$$(A.2)

• B_S (α) is the Shadow Hiding Opposition Effect (SHOE), $\frac{1}{1+\tan (\alpha /2)/h_S}$

• M(i, e) is the Isotropic Multiple-Scattering Approximation (IMSA): $$M(i,e)=H(\frac{\cos (i)}{K},w)H(\frac{\cos (e)}{K},w)-1$$(A.3)

  where H(x, w) is the Ambartsumian-Chandrasekhar function: $$H(x,w)\approx {[1-wx(r_0+\frac{1-2r_{0}x}{2}\ln \left(\frac{1+x}{x}\right))]}^{-1}$$(A.4)

  and r₀ is the diffusive reflectance: $$r_0=\frac{1-\sqrt{1-w}}{1+\sqrt{1-w}}$$(A.5)

• B_C(α) is the Coherent Backscatter Opposition Effect (CBOE): $$B_C(\alpha )=\frac{1+\frac{1-\exp (-\tan (\alpha /2)/h_C)}{\tan (\alpha /2)/h_C}}{2(1+\tan (\alpha /2)/h_C{)}^2}$$(A.6)

• The shadowing function S (i, e, ψ) accounts for the effect of macroscopic surface roughness. Its formulation varies depending on the relative values of the incidence angle i and the emission angle e:

  Case 1: i ≤ e $$S(i,e,\psi )=\frac{\mu }{\eta (e)}\frac{{\mu }_0}{\eta (i)}\frac{\chi ({\theta }_p)}{1-f(\psi )+f(\psi )\chi ({\theta }_p)\left(\frac{{\mu }_0}{\eta (i)}\right)}$$(A.7)

  with $${\mu }_0=\chi ({\theta }_p)[\cos (i)+\sin (i)\tan ({\theta }_p)\frac{\cos (\psi )E_2(e)+{\sin }^2(\psi /2)E_2(i)}{2-E_1(e)-(\psi /\pi )E_1(i)}]$$(A.8) $$\mu =\chi ({\theta }_p)[\cos (e)+\sin (e)\tan ({\theta }_p)\frac{E_2(e)-{\sin }^2(\psi /2)E_2(i)}{2-E_1(i)-(\psi /\pi )E_1(i)}]$$(A.9)

  Case 2: e < i $$S(i,e,\psi )=\frac{\mu }{\eta (e)}\frac{{\mu }_0}{\eta (i)}\frac{\chi ({\theta }_p)}{1-f(\psi )+f(\psi )\chi ({\theta }_p)\left(\frac{\mu }{\eta (e)}\right)}$$(A.10)

  with $${\mu }_0=\chi ({\theta }_p)[\cos (i)+\sin (i)\tan ({\theta }_p)\frac{E_2(i)-{\sin }^2(\psi /2)E_2(e)}{2-E_1(i)-(\psi /\pi )E_1(e)}]$$(A.11) $$\mu =\chi ({\theta }_p)[\cos (e)+\sin (e)\tan ({\theta }_p)\frac{\cos (\psi )E_2(i)+{\sin }^2(\psi /2)E_2(e)}{2-E_1(i)-(\psi /\pi )E_1(e)}]$$(A.12)

  The terms involved in the expressions above are defined as follows: $$\chi ({\theta }_p)=\frac{1}{\sqrt{1+\pi {\tan }^2({\theta }_p)}}$$(A.13) $$\eta (y)=\chi ({\theta }_p)[\cos (y)+\sin (y)\tan ({\theta }_p)\frac{E_2(y)}{2-E_1(y)}]$$(A.14) $$E_1(y)=\exp [-\frac{2}{\pi }\cot ({\theta }_p)\cot (y)]$$(A.15) $$E_2(y)=\exp [-\frac{1}{\pi }{\cot }^2({\theta }_p){\cot }^2(y)]$$(A.16) $$f(\psi )=\exp [-2\tan \left(\frac{\psi }{2}\right)]$$(A.17)

Appendix B OSIRIS-REx Bennu’s image

To provide visual context, we include in Figure B.1 a global view of Bennu’s surface acquired by the OSIRIS-REx spacecraft using the PolyCam imager during its approach. This mosaic offers a comprehensive look at Bennu’s topography and albedo variations across its entire surface. The image was rotated to approximately match the viewing geometry from Earth on 14-17 September 2005, the period for which we generated synthetic light curves for validation against ground-based data. This orientation helps illustrate the correspondence between surface features and their potential photometric signature in rotational brightness variations.

Although the orientation does not exactly reproduce the sub-Earth latitude or illumination conditions of the 2005 observations, it serves as a valuable qualitative reference for assessing discrepancies between the synthetic simulations and the actual observational data.

Fig. B.1 Bennu imaged by OSIRIS-REx (PolyCam) during the approach phase. The image is rotated so that the spin axis orientation approximates the aspect as seen from Earth in September 2005. Credit: NASA/Goddard/University of Arizona.

All Tables

All Figures

	Fig. 1 Workflow of the web application from user input to result rendering.
In the text

Fig. 2 Haumea on January 21, 2017, as seen from Earth. On this date, Ortiz et al. (2017) reported observations from multiple Earth-based observatories of Haumea passing in front of a distant star (a multi-chord stellar occultation). These observations allowed them to constrain the shape of Haumea’s projected ellipse on the sky. The pink lines represent the projected occultation chords derived from the reported timings, while the yellow endpoints indicate the associated uncertainties in the ingress and egress times. The facet colors on the projected body indicate the relative level of irradiance received by each facet, with red corresponding to lower incidence angles and therefore higher irradiation. The overall agreement between the projected shape and the occultation chords supports the validity of our simulation.

In the text

	Fig. 3 Projected view of Bennu, generated with an approximately 20000-facet downsampled shape model, as seen from Earth on 15 September 2005 (phase angle of 61.9°).
In the text

	Fig. 4 Observational data from the light curve published by Hergenrother et al. (2013) shown with black dots. Colored lines show the synthetic rotational light curves of Bennu as seen from Earth on 14-17 September 2005, sampled at one point per degree of rotation. An offset is applied to the synthetic light curves for clarity.
In the text

	Fig. 5 Projected view of Eros, generated with an approximately 20 000-facet downsampled shape model, as seen from Earth on July 17, 2016 (phase angle of 24.07°).
In the text

	Fig. 6 Observational data from the light curve published by Aznar Macias et al. (2018) represented by black dots. The green line shows the synthetic rotational light curve of Eros as seen from Earth on July 17, 2016, sampled at one point per degree of rotation.
In the text

	Fig. B.1 Bennu imaged by OSIRIS-REx (PolyCam) during the approach phase. The image is rotated so that the spin axis orientation approximates the aspect as seen from Earth in September 2005. Credit: NASA/Goddard/University of Arizona.
In the text