Light Curve Tool for Small Bodies

Introduction

This application is specifically designed to generate light curves for small, atmosphere-less bodies in the Solar System, with the primary objective of supporting research in this domain. While alternative tools exist, this application offers a distinctive combination of features: simulations are executed directly on the server, eliminating the need for local installation, and it generates projected images of objects in the sky plane. These functionalities are particularly valuable for validating research in fields such as stellar occultations. Developed in Python with Django as its framework, the tool currently provides three operational modes:

• Rotational light curves
• Known Body Mode: Allows users to generate the light curve and projected shape of a known object by entering its Minor Planet designation and observation date.
• Generic Body Mode: Enables users to perform generic simulations of an object by specifying its aspect angle.
• Phase Curve: It is used to represent the phase curve of the object. In this case, the object will rotate as its phase angle moves between the values entered by the user.
• Phase curves
• Known Body Mode: This tool allows generating light curves of rotating objects at different phase values defined by the user.

In all modes, users have the option to upload a custom shape model in .obj format or define the triaxial axes (a, b, and c) of an ellipsoid. The following sections provide a detailed explanation of each tool's features. For a demonstration of the rotational light curve functionality, watch this video. To explore the phase curve feature, refer to this video.

Known Body Mode

Object: In this mode, the user is required to provide the Minor Planet number assigned to the object. This information can be retrieved from the JPL Small-Body Database.

Date: The user must also specify the observation date from which the light curve should be simulated. The required format is yyyy-mm-dd hh:mm:ss.ssss (e.g., 2001-01-01 12:30:12.4320). If the date provided is in the yyyy-mm-dd format, the tool will automatically assume the time to be 00:00:00.0000.

Generic Body Mode

Phase Angle: The phase angle is the angle formed between the vector pointing toward the Sun and the vector pointing toward the observer. In this mode, Earth is fixed along the X-axis of the ecliptic coordinate system, the object is located at the origin of the coordinate system, and the Sun moves according to the angle specified by the user within the XY plane (from the X-axis toward the Y-axis, see Fig. 1). In case the check box Generate Images is selected, the positive Z-axis points upward, while the positive Y-axis is oriented to the right in the image. The user can specify the object’s rotation pole in ecliptic coordinates by providing values for longitude and latitude as explained below.

Common Input Parameters for the Rotational Light Curves

Rotation Pole: The rotation pole must be specified in ecliptic coordinates, using decimal degrees. The longitude value is defined such that it starts at x = 0 and increase in the direction of the Y-axis. For instance, a pole with a longitude of 90 degrees will be aligned with the Y-axis, pointing in its positive direction. If the rotation pole is known in another reference frame, it can be converted using this external tool.

LC number of points: This parameter defines the temporal resolution of the rotational light curve. An integer value is required to determine the number of points sampled from the rotational phase.

Rotational phase offset: A rotational phase offset, ranging between 0 and 1, can be applied to adjust the starting point of the light curve. This adjustment ensures that the curve begins at a specific phase, such as a maximum, minimum, or any other desired point.

Input Parameters for the Phase Curves

This mode allows generating phase curves of rotating objects at different phase values. To achieve this, the body is positioned at the center of the ecliptic reference system, while the observer is located along the X-axis. The Sun will rotate along the object's equatorial plane, sweeping through the entire phase range specified by the user. The object has a rotation pole at latitude +90° (Pole-off), and it completes a full rotation for each phase value. As a result, a graph is generated with brightness normalized to unity at the maximum phase curve value, where error bars represent the standard deviation of the reflectance values obtained for each rotation.

Phase value range: Specifies the minimum and maximum phase values (in degrees) at which the object's computation is performed.

LC number of points: Determines the resolution of the phase curve by setting the number of data points used in the calculation.

Shape Model

Custom Shape Model: A custom shape model can be uploaded in .obj format, with facets required to be triangular for compatibility. It is recommended to use Blender, a free and versatile software, to create and export models in the appropriate format for this tool. Models with a higher number of facets will result in longer computation times, and it is recommended not to exceed 3000 facets. In this video, the user can find a basic example of how to generate a shape model. It is important to specify that the forward axis is set to Y and the up axis is set to Z before exporting. Otherwise, an undesired orientation may occur. An example of a shape model in the required format can be downloaded at this link.

Triaxial Ellipsoid: A triaxial ellipsoid composed of 642 triangular facets can be generated by specifying the values of the axes a, b, and c. The shape model generated through this process will be included in the download folder.

Photometric Model

Radiance represents the intrinsic brightness of an object, independent of the observation distance or measurement instrument, but dependent on the photometric angles of phase (α), emission (e), and incidence (i). It is quantified as energy per unit area, solid angle, and time. However, for small bodies reflecting sunlight, it is often more practical to consider the ratio between incident and reflected energy. A commonly used term for this is I/F, which represents the ratio of reflected light to that of a perfectly Lambertian surface under the same illumination conditions. In this tool, several photometric models are implemented. On the one hand, we use empirical models, assuming that radiance can be expressed as the product of two components: a phase function (f(α)) and a disk function (d(α,e,i)). This approach allows the simulation to separately account for the effects of the Sun-object-observer geometry (phase angle) and the surface scattering properties of the object.

$$\frac{I}{F} = f(\alpha) \cdot d(\alpha, e, i)$$

On the other hand, we implement the Hapke model (Hapke 2012). This model is significantly more complex than empirical models but has the advantage that its parameters have a physical meaning (or at least that's what it's supposed to be), making them more intuitive and applicable in scientific studies.

For simplicity, the solar irradiance (S) received by the object is normalized to 1 W/m². This assumption facilitates rescaling of the results if needed.

The currently implemented and functional models are as follows:

Empirical models:

Phase functions

• No phase function: Assumes no dependence on the phase angle:
  
  $$f(\alpha) = 1$$
  
  
• Minnaert: Assumes a dependence as follows:
  
  $$ f(\alpha) = A_{\text{min}} \cdot \pi \cdot 10^{-0.4 (\beta \alpha + \gamma \alpha^2 + \delta \alpha^3)} $$
  
  where Amin, β, γ, and δ are fitting parameters.

Disk functions

• Lambertian: Assumes isotropic scattering given by:
  
  $$d(i) = \frac{\cos i}{\pi}$$
  
  
• Lommel-Seeliger: Models as:
  
  $$ d(e, i) = \frac{\cos i}{\cos i + \cos e} $$
  
  
• Minnaert: Dependence on incidence and emission, parameterized through the parameter k:
  
  $$ d(e, i) = \cos^k (i) \cdot \cos^{k-1} (e) $$
  
  

Hapke model:

The Hapke model is based on radiative transfer theory and provides a more intuitive and interpretable approach to surface scattering. The model accounts for multiple scattering, macroscopic roughness, and opposition effects, which improve the accuracy of brightness predictions for irregular and granular surfaces. The reflectance in the Hapke model used in this tool, which depends on 9 parameters, is given by:

$$ \frac{I}{F} = LS(i, e) \cdot K \cdot \frac{w}{4} \left[ p(\alpha)(1 + B_{S0} B_{S}(\alpha)) + M(i, e) \right] \left[ 1 + B_{C0} B_{C}(\alpha) \right] S(i, e, \psi) $$

where:

• LS(e, i) is the Lommel-Seeliger function.
• K is the porosity factor, accounting for voids between particles in the surface.
• w is the single-scattering albedo, representing the fraction of incident light scattered by a single particle.
• p(α) is the phase function.
• BS0 is the amplitude of the Shadow Hiding Opposition Effect (SHOE).
• BS(g) is the SHOE function.
• M(i,e) is the Isotropic Multiple-Scattering Approximation (IMSA).
• BC0 is the amplitude of the Coherent Backscatter Opposition Effect (CBOE).
• BC(α) is the CBOE function.
• S(e, i, ψ) represents the macroscopic roughness correction, accounting for large-scale variations in the surface slope, where ψ is the azimuthal angle between the planes of incidence and emission:

$$ \psi = \arccos \left( \frac{\cos(α) - \cos(e) \cos(i)}{\sin(e) \sin(i)} \right) $$

Each term of the above equation is calculated as:

• Lommel-Seeliger function (LS(e, i)):

  $$LS(e, i) = \frac{\cos i}{\cos i + \cos e}$$

• Porosity factor (K):

  $$ K = \frac{-\ln (1 - 1.209\phi^{2/3})}{1.209\phi^{2/3}} $$

  where Φ is the filling factor, equivalent to 1 - porosity.


• Phase function (p(α)), we employ the Henyey-Greenstein double-lobed single particle phase function:

  $$ p(α) = \frac{(1 + c_{HG})}{2} \cdot \frac{(1 - b_{HG}^2)}{(1 - 2b_{HG} \cos α + b_{HG}^2)^{3/2}} + \frac{(1 - c_{HG})}{2} \cdot \frac{(1 - b_{HG}^2)}{(1 + 2b_{HG} \cos α + b_{HG}^2)^{3/2}} $$

  
  where bHG (0 ≤ bHG ≤ 1) is the shape-controlling parameter, and cHG (−1 ≤ cHG ≤ 1)is the relative strength of backward and forward lobes.


• SHOE function, or Shadow Hiding Opposition Effect:

  $$ B_S(α) = \frac{1}{1 + \tan(α/2)/h_S} $$

  where hS is the angular width parameter of SHOE.


• Isotropic Multiple-Scattering Approximation (IMSA):

  $$ M(i, e) = H\left(\frac{\cos(i)}{K}, w\right) H\left(\frac{\cos(e)}{K}, w\right) - 1 $$

  where H(x,w) is the Ambartsumian-Chandrasekhar H function, which can be approximated as:

  $$ H(x, w) \approx \left\{ 1 - w x \left[ r_0 + \frac{1 - 2 r_0 x}{2} \ln \left(\frac{1 + x}{x}\right) \right] \right\}^{-1} $$

  where r0 is the diffusive reflectance, given by

  $$ r_0 = \frac{1 - \sqrt{1 - w}}{1 + \sqrt{1 - w}} $$

• Coherent Backscatter Opposition Effect function:

  $$ B_C(\alpha) = \frac{1 + \frac{1 - \exp(-\tan(\alpha/2) / h_C)}{\tan(\alpha/2) / h_C}}{2 \left[ 1 + \tan(\alpha/2) / h_C \right]^2} $$

  where hC is the angular width of CBOE.


• Shadowing function \( S(i, e, \psi) \), which is function of θp, the effective value of the photometric roughness, and given by:
  if 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) [\mu_0 / \eta(i)]} $$

  with

  $$ \mu_{0} = \chi({\theta}_p) \left[ \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)} \right] $$

  $$ \mu = \chi({\theta}_p) \left[ \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)} \right] $$

  if 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) [\mu / \eta(e)]} $$

  with

  $$ \mu_{0} = \chi({\theta}_p) \left[ \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)} \right] $$

  $$ \mu = \chi({\theta}_p) \left[ \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)} \right] $$

  With the following expressions, which are independent of e and i values, as follows:

  $$ \chi({\theta}_p) = \frac{1}{\sqrt{1 + \pi \tan^2({\theta}_p)}} $$

  $$ \eta(y) = \chi({\theta}_p) \left[ \cos(y) + \sin(y) \tan({\theta}_p) \frac{E_2(y)}{2 - E_1(y)} \right] $$

  $$ E_1(y) = \exp \left[ -\frac{2}{\pi} \cot({\theta}_p) \cot(y) \right] $$

  $$ E_2(y) = \exp \left[ -\frac{1}{\pi} \cot^2({\theta}_p) \cot^2(y) \right] $$

  $$ f(\psi) = \exp \left[ -2 \tan \left( \frac{\psi}{2} \right) \right] $$

To use the Hapke model in this tool, the following 9 parameters are required:

Albedo Spot

The Albedo Spot checkbox allows users to simulate surface features with different reflectance properties (albedo) on the surface. This feature is particularly useful for studying the effects of localized surface variations on the rotational light curve.

When this checkbox is enabled, the user can define the following parameters:

• Albedo variation (%): Specifies the percentage reduction in reflectance for the selected spot. For example, a value of 20% reduces the reflectance by 20% compared to the surrounding surface.
• Latitude range: Defines the minimum and maximum latitudes (in decimal degrees) where the albedo spot is located.
• Longitude range: Defines the minimum and maximum longitudes (in decimal degrees) where the albedo spot is located.

By adjusting these parameters, users can model and analyze how surface heterogeneities influence the rotational light curve of the object.

Precession

The Precession checkbox allows users to simulate the motion of the rotation pole as it precesses over time. Instead of remaining fixed, the rotation pole follows a conical trajectory, modifying its longitude and latitude dynamically during the simulation.

When this checkbox is enabled, the user can define the following parameters:

• Amplitude Lat (°): The maximum variation of the rotation pole's latitude due to precession. A higher value results in a more pronounced inclination change.
• Amplitude Lon (°): The maximum variation of the rotation pole's longitude due to precession. This determines how much the pole drifts in the longitude direction.
• Precession freq. (cycles per body rotation): Defines how many full precession cycles occur per complete rotation of the body. A frequency of 1 means the pole completes one full precession cycle for every full rotation.

The precession is modeled using sinusoidal variations for both the longitude and co-latitude of the pole:

β = β0 + Aβ * cos(2π * ωp * t)

λ = λ0 + Aλ * sin(2π * ωp * t)

Where:

• β0 and λ0 are the initial latitude and longitude of the pole.
• Aβ and Aλ define the amplitude of the precession in latitude and longitude .
• ωp is the precession frequency in cycles per body rotation.
• t represents the rotational phase, varying from 0 to 1 over a full rotation.

By adjusting these parameters, users can simulate different precession behaviors, from small oscillations to large variations where the rotation pole traces a wide circular path over time.

Check Boxes

The tool includes the following checkboxes to customize the output:

• Generate Images: Enables the generation of projected images for the object in addition to the light curve. These images can be used for visualization or further validation purposes.
• Self-shadowing: Includes the effect of self-shadowing in the simulation, which considers how different parts of the object's surface may block sunlight. Be aware that enabling this option can significantly increase computation time.

Output

After pressing the Submit button, the server will perform the necessary calculations and display a summary of the input parameters on the screen. The rotational light curve, representing relative magnitude as a function of rotational phase, will be displayed below.

If the Generate Images checkbox is selected, an animation of the object's projections in the skyplane will also be displayed alongside the light curve.

If the result meets expectations, the generated data can be downloaded by clicking the Download Data button. This action will prompt the user to save a compressed .zip file containing:

• Individual projection images and rotational light curves. Each image is named according to the corresponding rotational phase value (e.g., 0.250_synthetic_LC_g indicates a rotational phase of 0.250).
• Animations (GIFs) of the body projections and the rotational light curve.
• A plain text file summarizing the input parameters and including three columns of data: rotational phase, total reflectance, and relative magnitude. The total reflectance is the sum of all I/F values calculated for each facet and is converted to relative magnitude using the following equation:
  
  Relative magnitude = -2.5 · np.log10(total_reflectance / np.average(total_reflectance))
  
  
• If the user provides axes a, b, and c, the generated shape model in .obj format will also be included.