# Getting Started

## Installation
You can install this library directly using pip:
```bash
pip install phenigraph
```
Alternatively, you can download it from GitLab, after which you will need to compile it.
In the main folder, run the following command:
```bash
python setup.py bdist_wheel
```
The installation wheel file will be created in the newly generated `dist` folder.
To install the library in your local environment, use the following command:
```bash
pip install /path/to/wheelfile.whl
```

### Directory and Data Required Before Starting
The library assumes that the working directory has a specific structure. You can change the reference directory and the location of necessary data files by creating a `paths.txt` file in the working directory.

The required data are:
- **The water cross-section**, which can be found on [this page](https://home.strw.leidenuniv.nl/~ewine/photo/display_h2o_b89e91ccf14bfd7f485dd7be7d789b0a.html).
  The file should be saved as `H2O.txt` in a directory called `Leiden_database`.
  The name and location of this directory can be changed in the `paths.txt` file.
- **The interstellar spectrum**, which can be found on [this page](https://home.strw.leidenuniv.nl/~ewine/photo/radiation_fields.html).
  The file should be saved as `ISRF.dat` in the `Leiden_database` directory.
  Several radiation files are available on this page; the standard ISRF works for most cases.
- **The luminosity data**. These data are more difficult to obtain, as they must be produced by stellar simulation codes such as [cesam2k20](https://www.ias.u-psud.fr/cesam2k20/).
  The structure of this file will be described in a dedicated section. Without such files, the luminosity will be assumed to be a constant `Ls`, which can be set in a configuration dictionary (see below), and the water photodissociation due to the central star will not be modelled.

All theses datas including luminosity datas (generated by Louis Manchon using [cesam2k20](https://www.ias.u-psud.fr/cesam2k20/) ) can be download [here](https://github.com/HuetPaul/Diffenix/blob/c2659196cdd5b8e544ca5dc6f411e9babdc53dfa/datas.zip)

#### Directory Structure
The main directory, referred to as `racine` in the paths file, is assumed to be the current working directory.
It is assumed that this directory contains a subdirectory named `Leiden_database`, which holds the `H2O.txt` and `ISRF.dat` files.
If the user has data modelling the central star, it is assumed to be in a `data_star` directory inside the main `racine` directory, in a `.npz` file called `data_fus.npz`.
Finally, a `Results` directory is required to save the `Sol` objects.

#### The `paths.txt` File
If you do not want to work inside the `racine` directory or wish to organise the data differently, you can create a `paths.txt` file.
This file is structured as follows:
- One line represents one directory.
- Each line starts with the name of the directory to be relocated; this can be `racine`, `Results`, `data_leiden`, or `data_stars`.
- A comma separates this name from the path to the new directory.

Example: Reproducing the default structure.
```
racine,.
results,./Results
data_leiden,./Leiden_database
data_stars,./data_stars/data_fus.npz
```

---

## A First Example
In the following sections, we will assume that these libraries are imported with the following names:
```python
import numpy as np
from scipy.integrate import cumulative_trapezoid
import phenigraph as g  # To generate plots
from diffenix import time_integration as ti  # Only the time_integration submodule is required
```
All information characterising the system to be modelled, along with details of the models used and numerical settings, should be included in a single dictionary. Examples of such dictionaries are provided with the library and can be adapted as needed.

### Setting Up a New Simulation
Let’s start with a disc that initially contains one Earth mass of hydrogen, with no asteroid belts or Kuiper belt, in a solar system-like configuration.
This system contains one planet with the same mass as Earth, located one astronomical unit from its central star.
To achieve this, let’s load the `sys_sol_ppd` configuration dictionary:
```python
conf_dic: dict = ti.sys_sol_ppd.copy()
```
`sys_sol_ppd` is defined as follows in the library:
```python
sys_sol_ppd: dict = {
    "racine": racine, "results": results,
    "C_L": C_L, "C_M": C_M, "C_t": C_t,
    "Ms": Msun, "Ts": np.double(5500.), "Rs": Rsun, "Ls": Lsun,
    "a_in": 0.05 * au, "a_out": 100 * au, "a0": 2.65 * au,
    "T0": T0, "pls_temp": np.double(-1 / 2), "mH_init": Mearth,
    "alpha": np.double(1e-3), "gamma": np.double(1.),
    "AB_mdot_st_init": np.double(1e-10) * Mearth / Myr,
    "AB_Mbelt": np.double(0.12) * Mearth,
    "AB_Abelt": np.double(0.06),
    "AB_a0": np.double(2.65) * au, "AB_delta_r": np.double(2.65 / 2.) * au,
    "AB_K": np.double(1e-5 * C_L * C_L / C_t), "AB_Phi": np.double(0.6),
    "AB_rp": np.double(1e-6 * C_L),
    "AB_f_ice": np.double(0.2), "AB_rho_refr": np.double(2.5e3 * C_M / (C_L * C_L * C_L)),
    "AB_rho_ice": np.double(1.e3 * C_M / (C_L * C_L * C_L)), "AB_sublimation_model": "none",
    "AB_const_mdot": np.double(1e-1) * Mearth / Myr,
    "AB_SD_amax": np.double(1.e6) * C_L, "AB_SD_abig": np.double(120e3) * C_L,
    "AB_SD_amed": np.double(20.e3) * C_L, "AB_SD_amin": C_L,
    "AB_SD_qh": np.double(-4.5), "AB_SD_qm": np.double(-1.2), "AB_SD_ql": np.double(-3.6),
    "AB_sb": np.double(316) * C_L, "AB_e": np.double(0.075), "AB_i": np.double(0.075 / 2.),
    "AB_As": np.double(5.) * C_E / C_M, "AB_bs": np.double(-0.1), "AB_bg": np.double(0.5),
    "AB_t0_diss": 150 * Myr, "AB_t1_diss": Myr,
    "KB_mdot_st_init": zero,
    "KB_Mbelt": np.double(0.01 * Mearth), "KB_Abelt": np.double(0.06),
    "KB_a0": np.double(44) * au, "KB_delta_r": np.double(8) * au,
    "KB_K": np.double(1e-10 * C_L * C_L / C_t), "KB_Phi": np.double(0.6),
    "KB_rp": np.double(1e-6 * C_L),
    "KB_f_ice": np.double(0.2), "KB_rho_refr": np.double(2.5e3 * C_M / (C_L * C_L * C_L)),
    "KB_rho_ice": np.double(1.e3 * C_M / (C_L * C_L * C_L)), "KB_sublimation_model": "none",
    "KB_const_mdot": zero,
    "KB_SD_amax": np.double(1.e6) * C_L, "KB_SD_abig": np.double(120e3) * C_L,
    "KB_SD_amed": np.double(20.e3) * C_L, "KB_SD_amin": C_L,
    "KB_SD_qh": np.double(-4.5), "KB_SD_qm": np.double(-1.2), "KB_SD_ql": np.double(-3.6),
    "KB_sb": np.double(316) * C_L, "KB_e": np.double(0.075), "KB_i": np.double(0.075 / 2.),
    "KB_As": np.double(5.) * C_E / C_M, "KB_bs": np.double(-0.1), "KB_bg": np.double(0.5),
    "KB_t0_diss": np.double(10.) * Gyr, "KB_t1_diss": np.double(10.) * Gyr,
    "a_planets": np.array([au]), "m_planets": np.array([Mearth]), "f_accr": np.double(0.5),
    "t0": np.double(5. * Myr),
    "IT_s": 3000, "dip": 10, "rtol": np.double(1e-6), "atol": np.double(1e-250), "method_root": 'RK45',
}
```

### Time Integration
Let’s start our first simulation:
```python
sol: ti.Sol = ti.time_integration(tf=ti.Myr, save_time_step=0.1*ti.Myr, consts_params=configs)
```
The evolution of the disc is modelled from an initial time of zero to a final time of one million years (`tf`).
The initial time is defined as the point at which the protoplanetary disc dissipates.
This can be set relative to the star’s lifetime in the configuration dictionary via the `t0` parameter.
This is only useful if luminosity data are provided for the central star.
The result is saved in a dedicated `Sol` object, `sol`.
If the integration stops before the end of the integration time `tf`, the data will not be lost.
This object is automatically saved to the `Results` directory at each 0.1 Myr time step.
As there is no `Simulation` entry in the configuration dictionary, the simulation’s name is generated automatically and displayed when the integration begins.
If the code is run for the first time, it should be `Simulation_1`. The number will then increase for each new simulation.
A new directory with the simulation’s name should have been created inside the `Results` directory.
The `.npz` file containing the `Sol` object is located inside this directory and is named `Simulation_diffusion.npz`.
To reload this simulation, follow these steps:
```python
sol = ti.Sol(filename="Simulation_diffusion", directory=ti.results+"/Simulation_1")
```

### Diagnostics

#### Surface Density
All diagnostics can be run for this simulation.
Let’s begin with the surface density. The `Sol` object has a specific method that generates a [Graphique](https://huetpaul.github.io/PheniGraph/) representing the surface density at different times as a function of the distance to the central star:
```python
gr = sol.graph_surface_density(save=True)
```
![result](Simulation_ppd_surface_density.png)
The argument `save=True` means that the graph and the PNG image will be saved in the simulation directory.
If you want to save the [Graphique](https://huetpaul.github.io/PheniGraph/) under a different name or in a different location, you can do so by specifying the desired name and location:
```python
gr = sol.graph_surface_density()
gr.filename = "your_new_name"
gr.directory = "your_directory"  # To save the plot in the current working directory, you can comment out this line
gr.save_figure(dpi=100)  # A 100 dpi PNG image will be saved
gr.save()  # The Graphique object will be saved for later reuse (and modification if needed)
```

#### Total Mass
It is also possible to access the total disc mass as a function of time:
```python
gr = g.loglog(sol.tps / ti.Myr, sol.total_mass() / ti.Mearth, show=False)  # Create the Graphique
gr.config_ax(xlabel="Time [Myr]", ylabel="Total mass [M$_\oplus$]")  # Set labels for the x and y axes
gr.show()  # Show the plot
gr.filename = "Total_mass"
gr.directory = sol.directory  # To save the plot in the current working directory, you can comment out this line
gr.save_figure(dpi=100)  # A 100 dpi PNG image will be saved
gr.save()  # The Graphique object will be saved for later reuse (and modification if needed)
```
![result](Simulation_ppd_total_mass.png)

#### Mass and Momentum Conservation
To save memory, the saved time steps can be larger than the real time steps.
However, the values used for diagnostics are saved at each real time step to monitor mass and momentum conservation:
```python
gr = g.loglog(sol.diag["time"] / ti.Myr, [sol.mass_conservation(), sol.momentum_conservation()],
              label=["Mass conservation", "Momentum conservation"], show=False)  # Create the Graphique
gr.config_ax(xlabel="Time [Myr]", ylabel="Relative variations")  # Set labels for the x and y axes
gr.show()  # Show the plot
gr.filename = "Mass_conservation"
gr.directory = sol.directory  # To save the plot in the current working directory, you can comment out this line
gr.save_figure(dpi=100)  # A 100 dpi PNG image will be saved
gr.save()  # The Graphique object will be saved for later reuse (and modification if needed)
```
![result](Simulation_ppd_mass_momentum_conservation.png)

#### Mass Flux
It is possible to obtain the mass flux for an index `i` of `sol.tps`. Let’s plot the last one at one million years:
```python
gr = g.loglog(sol.r, [sol.mass_flux(-1), -sol.mass_flux(-1)], label=[">0", "<0"], show=False)  # Create the Graphique
gr.config_ax(xlabel="Distance to the central star [au]", ylabel="Mass flux [M$_\oplus$ / Myr]")
gr.show()  # Show the plot
gr.filename = "mass_flux"
gr.directory = sol.directory  # To save the plot in the current working directory, you can comment out this line
gr.save_figure(dpi=100)  # A 100 dpi PNG image will be saved
gr.save()  # The Graphique object will be saved for later reuse (and modification if needed)
```
![result](Simulation_ppd_mass_flux.png)
As can be seen, there is no outer flux. The boundary condition implies that the outer flux is zero.
There is a sharp decrease in the mass flux in the inner part of the system due to the planet’s accretion of gas.
This can be confirmed by plotting a grey-shaded area on the planet’s sink cells:
```python
sol.add_planets_to_graph(gr, show=True)
gr.save_figure(dpi=100)  # A 100 dpi PNG image will be saved
gr.save()  # The Graphique object will be saved for later reuse (and modification if needed)
```
![result](Simulation_ppd_mass_flux2.png)
Similar plots can be created for the radial velocity.

#### The Planet
In addition to the diagnostic values saved at each time step, the mass flux at the outer and inner boundaries of the planet’s sink cells is saved. The accretion rate and the final accreted mass can be deduced from these values.
```python
gr = g.loglog(sol.diag["time"], [-np.array(sol.flux_planets)[:,1,0], -np.array(sol.flux_planets)[:, 0,0]],
             label=["Inner flux", "Outer flux"], show=False)
gr.config_ax(xlabel="Time [Myr]", ylabel="Mass flux [M$_\oplus$ / Myr]")
gr.filename = "planet_mass_flux"
gr.save_figure(dpi=100)  # A 100 dpi PNG image will be saved
gr.save()  # The Graphique object will be saved for later reuse (and modification if needed)
```
![result](Simulation_ppd_planet_mass_flux.png)
The shape of `sol.flux_planets` is `(len(sol.tps), 2, number_of_planets)`.
We can then estimate the accreted mass onto the planet as a function of time:
```python
gr = g.loglog(sol.diag["time"],
              ti.cumulative_trapezoid(np.array(sol.flux_planets)[:,1,0] - np.array(sol.flux_planets)[:, 0,0],
                                      sol.diag["time"], initial=0.) / ti.Mearth, show=False)
gr.config_ax(xlabel="Time [Myr]", ylabel="Accreted mass [M$_\oplus$]")
gr.filename = "Simulation_ppd_planet_mass"
gr.save()
gr.save_figure(dpi=400)
```
![result](Simulation_ppd_planet_mass.png)

---

## An Asteroid Belt with a Constant Gas Mass Generation Rate
The previous simulation was initialised with a given mass of hydrogen, and no further mass was produced. There was also only one molecular species.
However, this code can model the evolution of multiple chemical species, from their production at the asteroid belt to their photodissociation.
Let’s start with a simple case involving a constant gas production rate at the asteroid belt level, using the same system as before.
After the evolution of the total surface density, the evolution of the mass fraction of each chemical species is integrated and saved at different times.
To avoid divergence in the numerical integration of this second set of equations, it is recommended that the simulation is initialised with a small amount of water.
This is controlled by the `AB_mdot_st_init` parameter in the configuration dictionary.
The system is initialised with a steady-state water gas disc and a gas generation rate of `AB_mdot_st_init`.
In this simulation, there is no primordial hydrogen remaining.
```python
configs = ti.sys_sol_ast.copy()
configs["AB_sublimation_model"] = "constant_rate"
configs["AB_mdot_st_init"] = np.double(1e-10) * ti.Mearth / ti.Myr  # Set the initial state
sol = ti.time_integration(tf=ti.Myr, save_time_step=0.1*ti.Myr, consts_params=configs)
```

### Surface Density
Since there are multiple species, the `graph_surface_density` is a bit more difficult to read:
```python
gr = sol.graph_surface_density(save=True)
```
![result](Simulation_ab1_surface_density.png)
But it is possible to plot only the water and neutral oxygen (the behaviour of neutral hydrogen is very similar to that of neutral oxygen):
```python
gr = sol.graph_surface_density(species=["H2O","O"], show=False)
gr.config_ax(ylim=[1e-15, 0.1])
gr.filename = "Surface_density"
gr.directory = sol.directory
gr.save_figure(dpi=100)
gr.save()
gr.show()
```
![result](Simulation_ab1_surface_density2.png)
The species names are all in capital letters and can be found in `sol.nom_mu`.
To gain a better understanding of how the surface density has changed over time, you can also use a second type of plot:
```python
gr_tot = sol.image_surface_density(save=True, show=True)  # For the total surface density
gr_tot = sol.image_surface_density("H2O", save=True, show=True, cmap="default", color_min=g.C1,
                                   color_max=g.C2)  # For the water surface density
gr_tot = sol.image_surface_density("O", save=True, show=True, plot_sigma_crit=False, cmap="default",
                                   color_min=g.C17, color_max=g.C13)  # For the oxygen surface density
```
![result](Simulation_ab1_Image_total_surface_density.png)
![result](Simulation_ab1_Image_H2O_surface_density.png)
![result](Simulation_ab1_Image_O_surface_density.png)
These plots can be used to evaluate changes in surface density over time, such as those in the neutral oxygen levels.
For example, it is produced at the beginning of the simulation until the critical surface density is reached. Then, the surface density of neutral oxygen decreases.
You can choose a Matplotlib colour bar; the default is `inferno`. You can also build a custom colour bar from the colour minimum to the colour maximum (the `cmap` parameter should be set to `default`).
See the [PheniGraph](https://huetpaul.github.io/PheniGraph/) documentation for more information.

---

## Thermal Sublimation Models
There are different models for the gas generation rate. Two are based on gas sublimation by thermal heating.
The first and simplest is referred to in the configuration dictionary as `"thermal_diff"`.
This model assumes that the sublimation timescale and gas diffusion into the asteroid or Kuiper Belt Object (KBO) are negligible.
It only evaluates the mass generation rate by the thermal diffusion timescale, set by the thermal diffusion coefficient `K` in the configuration dictionary (`AB_K` and `KB_K` for the asteroid belt and the Kuiper belt, respectively).
Let’s model a young Kuiper belt similar to that of [Huet et al. 2025](https://ui.adsabs.harvard.edu/abs/2025arXiv250603412H/abstract):
In this simulation, only CO is produced at the Kuiper belt level.
```python
configs = ti.sys_sol_kuip.copy()
configs["KB_sublimation_model"] = "thermal_diff"
sol = ti.time_integration(tf=50 * ti.Myr, save_time_step=ti.Myr, consts_params=configs)
```

---

## The Configuration Dictionary
Here, we introduce each value needed to build the configuration dictionary. Best practice is to adapt one of the default configuration dictionaries, such as `sys_sol_ast`, `sys_sol_ast_kuip`, or `sys_sol_appd`, rather than building a new one from scratch.

- `racine` → Path to the root directory.
- `results` → Path to the results subdirectory.
- `C_L`: `np.float64` → Length conversion unit (cannot differ from the global value; see the [Default Unit System](#default-unit-system) subsection to modify the global unit system).
- `C_M`: `np.float64` → Mass conversion unit.
- `C_t`: `np.float64` → Time conversion unit.
- `Ms`: `np.float64` → Central star’s mass.
- `Ts`: `np.float64` → Central star’s surface temperature (ignored if luminosity data are provided).
- `Rs`: `np.float64` → Central star’s radius (ignored if luminosity data are provided).
- `Ls`: `np.float64` → Central star’s luminosity (ignored if luminosity data are provided).
- `a_in`: `np.float64` → Inner boundary of the integration domain.
- `a_out`: `np.float64` → Outer boundary of the integration domain.
- `a0`: `np.float64` → Radius used for the primordial hydrogen distribution (assumed to be the steady state centred on `a0`).
- `T0`: `np.float64` → Temperature at `a0` (ignored if luminosity data are provided).
- `pls_temp`: `np.float64` → Radial power slope of the temperature.
- `mH_init`: `np.float64` → Initial mass of primordial hydrogen.
- `alpha`: `np.float64` → Viscous parameter.
- `gamma`: `np.float64` → Adiabatic constant.
- `AB_mdot_st_init`: `np.float64` → Initial steady-state mass generation rate for the asteroid belt.
- `AB_Mbelt`: `np.float64` → Total mass of the asteroid belt.
- `AB_Abelt`: `np.float64` → Asteroids’ albedo.
- `AB_a0`: `np.float64` → Mean semi-major axis of the asteroid belt.
- `AB_delta_r`: `np.float64` → Half radial width of the asteroid belt (the belt extends from `AB_a0 - AB_delta_r` to `AB_a0 + AB_delta_r`; ensure that the inner and outer boundaries of the integration domain are smaller/larger than the inner/outer limits of the belt, i.e., `a_in < AB_a0 - AB_delta_r` and `AB_a0 + AB_delta_r < a_out`).
- `AB_K`: `np.float64` → Thermal diffusivity in asteroids (dimension: $[L]^2 . [T]^{-1}$).
- `AB_Phi`: `np.float64` → Asteroids’ porosity (dimensionless).
- `AB_rp`: `np.float64` → Pore radius for asteroids.
- `AB_f_ice`: `np.float64` → Ice-to-refractory mass ratio for asteroids (dimensionless).
- `AB_rho_refr`: `np.float64` → Density of refractory materials for asteroids (dimension: $[M] . [L]{-3}$).
- `AB_rho_ice`: `np.float64` → Ice density for asteroids (dimension: $[M] . [L]{-3}$).
- `AB_sublimation_model`: `str` → Sublimation model; can be `none`, `constant_rate`, `thermal_diffusion`, or `thermal_full`.
- `AB_const_mdot`: `np.float64` → Mass generation rate for the `constant_rate` sublimation model.
- `AB_SD_amax`: `np.float64` → Maximum size used for the size distribution of the asteroid belt.
- `AB_SD_abig`: `np.float64` → Size used for the size distribution of the asteroid belt.
- `AB_SD_amed`: `np.float64` → Medium size used for the size distribution of the asteroid belt.
- `AB_SD_amin`: `np.float64` → Minimum size used for the size distribution of the asteroid belt.
- `AB_SD_qh`: `np.float64` → Power slope of the size distribution between `AB_SD_amax` and `AB_SD_abig`.
- `AB_SD_qm`: `np.float64` → Power slope of the size distribution between `AB_SD_abig` and `AB_SD_amed`.
- `AB_SD_ql`: `np.float64` → Power slope of the size distribution between `AB_SD_amed` and `AB_SD_amin`.
- `AB_e`: `np.float64` → Asteroids’ mean eccentricity (used by the `collisions` sublimation model).
- `AB_i`: `np.float64` → Asteroids’ mean inclination (used by the `collisions` sublimation model).
- `AB_As`: `np.float64` → Asteroids’ minimum mass disruption energy (used by the `collisions` sublimation model; dimension $[L]^2 . [T]^{-2}$); see [Löhne et al. 2008](https://ui.adsabs.harvard.edu/abs/2008ApJ...673.1123L/abstract).
- `AB_bs`: `np.float64` → Slope of the asteroids’ mass disruption energy (used by the `collisions` sublimation model; dimensionless); see [Löhne et al. 2008](https://ui.adsabs.harvard.edu/abs/2008ApJ...673.1123L/abstract).
- `AB_bg`: `np.float64` → Slope of the asteroids’ mass disruption energy (used by the `collisions` sublimation model; dimensionless); see [Löhne et al. 2008](https://ui.adsabs.harvard.edu/abs/2008ApJ...673.1123L/abstract).
- `AB_t0_diss`: `np.float64` → Asteroid belt lifetime before disruption begins.
- `AB_t1_diss`: `np.float64` → Asteroid belt disruption timescale.
- `KB_mdot_st_init`: `np.float64` → Initial steady-state mass generation rate for the Kuiper belt.
- `KB_Mbelt`: `np.float64` → Total mass of the Kuiper belt.
- `KB_Abelt`: `np.float64` → Kuiper Belt Objects’ albedo (dimensionless).
- `KB_a0`: `np.float64` → Mean semi-major axis of the Kuiper belt.
- `KB_delta_r`: `np.float64` → Half radial width of the Kuiper belt (the belt extends from `KB_a0 - KB_delta_r` to `KB_a0 + KB_delta_r`; ensure that the inner and outer boundaries of the integration domain are smaller/larger than the inner/outer limits of the belt, i.e., `a_in < KB_a0 - KB_delta_r` and `KB_a0 + KB_delta_r < a_out`).
- `KB_K`: `np.float64` → Thermal diffusivity in Kuiper Belt Objects (dimension: $[L]^2 . [T]^{-1}$).
- `KB_Phi`: `np.float64` → Kuiper Belt Objects’ porosity (dimensionless).
- `KB_rp`: `np.float64` → Pore radius for Kuiper Belt Objects.
- `KB_f_ice`: `np.float64` → Ice-to-refractory mass ratio for Kuiper Belt Objects (dimensionless).
- `KB_rho_refr`: `np.float64` → Density of refractory materials for Kuiper Belt Objects (dimension: $[M] . [L]{-3}$).
- `KB_rho_ice`: `np.float64` → Density of CO ice for Kuiper Belt Objects (dimension: $[M] . [L]{-3}$).
- `KB_sublimation_model`: `str` → Sublimation model; can be `none`, `constant_rate`, `thermal_diffusion`, or `thermal_full`.
- `KB_const_mdot`: `np.float64` → Mass generation rate for the `constant_rate` sublimation model.
- `KB_SD_amax`: `np.float64` → Maximum size used for the size distribution of the Kuiper belt.
- `KB_SD_abig`: `np.float64` → Size used for the size distribution of the Kuiper belt.
- `KB_SD_amed`: `np.float64` → Medium size used for the size distribution of the Kuiper belt.
- `KB_SD_amin`: `np.float64` → Minimum size used for the size distribution of the Kuiper belt.
- `KB_SD_qh`: `np.float64` → Power slope of the size distribution between `KB_SD_amax` and `KB_SD_abig`.
- `KB_SD_qm`: `np.float64` → Power slope of the size distribution between `KB_SD_abig` and `KB_SD_amed`.
- `KB_SD_ql`: `np.float64` → Power slope of the size distribution between `KB_SD_amed` and `KB_SD_amin`.
- `KB_e`: `np.float64` → Kuiper Belt Objects’ mean eccentricity (used by the `collisions` sublimation model).
- `KB_i`: `np.float64` → Kuiper Belt Objects’ mean inclination (used by the `collisions` sublimation model).
- `KB_As`: `np.float64` → Kuiper Belt Objects’ minimum mass disruption energy (used by the `collisions` sublimation model; dimension $[L]^2 . [T]^{-2}$); see [Löhne et al. 2008](https://ui.adsabs.harvard.edu/abs/2008ApJ...673.1123L/abstract).
- `KB_bs`: `np.float64` → Slope of the Kuiper Belt Objects’ mass disruption energy (used by the `collisions` sublimation model; dimensionless); see [Löhne et al. 2008](https://ui.adsabs.harvard.edu/abs/2008ApJ...673.1123L/abstract).
- `KB_bg`: `np.float64` → Slope of the Kuiper Belt Objects’ mass disruption energy (used by the `collisions` sublimation model; dimensionless); see [Löhne et al. 2008](https://ui.adsabs.harvard.edu/abs/2008ApJ...673.1123L/abstract).
- `KB_t0_diss`: `np.float64` → Kuiper belt lifetime before disruption begins.
- `KB_t1_diss`: `np.float64` → Kuiper belt disruption timescale.
- `a_planets`: `np.array[np.float64]` → Array of planets’ semi-major axes (one-dimensional array).
- `m_planets`: `np.array[np.float64]` → Array of planets’ masses (one-dimensional array; ensure it has the same size as `a_planets`).
- `f_accr`: `np.float64` → Planets’ accretion efficiency (dimensionless).
- `t0`: `np.float64` → Protoplanetary disc lifetime.
- `IT_s`: `int` → Number of spatial steps.
- `dip`: `int` → Half the number of sink cells around a planet.
- `rtol`: `np.float64` → Relative tolerance for time integration; see the [scipy.integrate.solve_ivp](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) documentation.
- `atol`: `np.float64` → Absolute tolerance for time integration; see the [scipy.integrate.solve_ivp](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) documentation.
- `method_root`: `str` → Time integration method: `RK45`, `RK23`, `DOP853`, `Radau`, `BDF`, or `LSODA`; see the [scipy.integrate.solve_ivp](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html) documentation.

---

## Default Unit System
The code uses its own unit system, fully described by the conversion factors from the universal MKSA system.
Only the temperature unit is set to Kelvin and cannot be modified.
There are three main conversion factors: `C_L`, `C_M`, and `C_t`, used to convert any length, mass, or time from metres, kilograms, and seconds to the code’s length, mass, or time units. Multiply the value by the correct conversion factor to perform the conversion.
For convenience, there are derived conversion factors, such as `C_F`, `C_rho`, `C_E`, and `C_Lum`, to convert force, density, energy, or luminosity from the MKSA system to the code’s unit system.
The default values are:
- For time: One code unit equals one million years (hence, $C_t = 1 / (3600 * 24 * 365.25 * 1e6)$).
- For mass: One code unit equals one Earth mass (hence, $C_M = 1 / M_\oplus$).
- For length: One code unit equals one astronomical unit.

To convert code values to MKSA units, divide them by the relevant conversion factor.
Example:
```python
mass_mksa = mass / C_M
length_mksa = length / C_L
time_mksa = time / C_t
surface_density_mksa = surface_density * C_L * C_L * C_L / C_M
```
Typical constants can also be used for unit conversion. For example, the constants `yr`, `kyr`, `Myr`, or `Gyr` can be used to convert times into years, kiloyears, megayears, or gigayears.
It is also possible to convert mass to solar mass by dividing the mass by one solar mass.
The default values were chosen to fit the typical values for our models. However, it is also possible to modify them if needed.
Since `Sol` objects retain their own conversion factors, they can be loaded with the unit system used to integrate them. All conversions will be performed automatically.
To do this, a new dictionary of constants must be created and saved in the working directory or the `racine` directory with the name `constantes`, using the [np.savez_compressed](https://numpy.org/doc/stable/reference/generated/numpy.savez_compressed.html) method:
```python
new_consts: dict = {
    "racine": new_racine, "results": new_results,
    "C_L": new_C_L, "C_M": new_C_M, "C_t": new_C_t,
    "Ms": new_Ms, "Ts": new_Ts, "Rs": new_Rs, "Ls": new_Ls,
    "a_in": new_a_in, "a_out": new_a_out, "a0": new_a0,
    "T0": new_T0, "pls_temp": new_pls_temp, "m_dot": new_m_dot,
    "alpha": new_alpha, "gamma": new_gamma,
    "AB_Mbelt": new_AB_Mbelt, "AB_Abelt": new_AB_Abelt,
    "AB_a0": new_AB_a0, "AB_delta_r": new_AB_delta_r,
    "AB_K": new_AB_K, "AB_Phi": new_AB_Phi, "AB_rp": new_AB_rp,
    "AB_f_ice": new_AB_f_ice, "AB_rho_refr": new_AB_rho_refr,
    "AB_rho_ice": new_AB_rho_ice, "AB_sublimation_model": new_AB_sublimation_model,
    "AB_SD_amax": new_AB_SD_amax, "AB_SD_abig": new_AB_SD_abig,
    "AB_SD_amed": new_AB_SD_amed, "AB_SD_amin": new_AB_SD_amin,
    "AB_SD_qh": new_AB_SD_qh, "AB_SD_qm": new_AB_SD_qm, "AB_SD_ql": new_AB_SD_ql,
    "AB_sb": new_AB_sb, "AB_e": new_AB_e, "AB_i": new_AB_i,
    "AB_As": new_AB_As, "AB_bs": new_AB_bs, "AB_bg": new_AB_bg, "AB_t0_diss": new_AB_t0_diss, "AB_t1_diss": new_AB_t1_diss,
    "KB_Mbelt": new_KB_Mbelt, "KB_Abelt": new_KB_Abelt,
    "KB_a0": new_KB_a0, "KB_delta_r": new_KB_delta_r,
    "KB_K": new_KB_K, "KB_Phi": new_KB_Phi, "KB_rp": new_KB_rp,
    "KB_f_ice": new_KB_f_ice, "KB_rho_refr": new_KB_rho_refr,
    "KB_rho_ice": new_KB_rho_ice, "KB_sublimation_model": new_KB_sublimation_model,
    "KB_SD_amax": new_KB_SD_amax, "KB_SD_abig": new_KB_SD_abig,
    "KB_SD_amed": new_KB_SD_amed, "KB_SD_amin": new_KB_SD_amin,
    "KB_SD_qh": new_KB_SD_qh, "KB_SD_qm": new_KB_SD_qm, "KB_SD_ql": new_KB_SD_ql,
    "KB_sb": new_KB_sb, "KB_e": new_KB_e, "KB_i": new_KB_i,
    "KB_As": new_KB_As, "KB_bs": new_KB_bs, "KB_bg": new_KB_bg, "KB_t0_diss": new_KB_t0_diss, "KB_t1_diss": new_KB_t1_diss,
    "a_planets": new_a_planets, "m_planets": new_m_planets, "f_accr": new_f_accr,
    "t0": new_t0,
    "IT_s": new_IT_s, "dip": new_dip, "rtol": new_rtol, "atol": new_atol, "method_root": new_method,
}
np.savez_compressed(racine + "/constantes.npz", **new_consts)
```
The current constants dictionary can be loaded using the `constantes()` method, which can then be modified.
When converting units, ensure that all resulting values are also converted, as every value in the constants dictionary is in code units.

---

## The Luminosity File
The central star’s luminosity is a key parameter for estimating the mass generation rate at the belt level in the `thermal_full` model.
It is also important for estimating the photodissociation rate.
However, the central star’s luminosity can vary significantly in the early stages of a system’s life, when secondary gas production is assumed to occur (see [Kral et al. 2024](https://ui.adsabs.harvard.edu/abs/2024A%26A...692A..70K/abstract) for the solar system and Figure 1).

![luminosity](stars_luminosity.png)
*Figure 1: Evolution of stellar luminosity as a function of time for different stellar masses.*

To provide the required data, you must first run stellar simulations (for example, using [cesam2k20](https://www.ias.u-psud.fr/cesam2k20/)).
These simulations must be performed for several masses (as many as possible to obtain reliable interpolation).
All data should then be merged into a single `.npz` file created using [np.savez_compressed](https://numpy.org/doc/stable/reference/generated/numpy.savez_compressed.html).
This file must contain the following keys associated with arrays of the same size:
- `Age`: Array containing the star’s age in Myr at each point.
- `M`: Star’s mass in solar masses (must be constant for each star) at each point.
- `Lbol`: Bolometric luminosity at each point in solar luminosity.
- `R`: Star’s radius in solar radii at each point.
- `Reuv_3692`: Luminosity fraction in the UV.
- `Reuv_lya`: Luminosity fraction in the UV.
```