use crate::core::constants::{
    ASTRONAUMICAL_UNIT, MINIMUM_PRECISION, NUMBER_THERMAL_SKIN_DEPTH, NUMERICAL_STABILITY,
    SOLAR_CONSTANT, TAU,
};
use crate::core::object_3d::Object3D;
use crate::toolbox::matrix::{clip, dot_vectors, linspace};
use crate::toolbox::physics::{
    compute_conductivity, compute_diffusivity, compute_ground_step_from_time_and_stability,
};
use crate::{M3, V3, V3X, VX};
use itertools::zip;
use na::{DMatrix, Dynamic, MatrixSliceMutXx1, MatrixSliceXx1, Rotation3, Unit, U1, U3};
use std::path::Path;

/// Message for expected error when calling a field of Properties that has not been computed.
const EXPECT_TIME_STEP_NOT_SET: &str = "The time step of the simulation has not been shared yet with the properties of the celestial body. Once the World is initialized, the time step is shared and you can call this variable.";

/// A celestial body representation.
///
/// # Definition
///
/// [`Body`] is defined by its name, its position with respect to the Sun, the shape model of its object file, the physical properties of its ground layers and its surface temperatures
///
/// # Example
///
/// ```no_run
/// use kalast::{Properties, Body, ASTRONAUMICAL_UNIT, HOUR, V3};
///
/// let mut body = Body::new(
///     "Dimorphos", // name
///     V3::new(0., 1., 0.) * ASTRONAUMICAL_UNIT, // distance to Sun
///     "path/to/dimorphos.obj", // 3D object path
///     Properties::new(
///         11.92 * HOUR, // rotation period
///         11.92 * HOUR, // revolution period
///         0.,           // obliquity
///         500.,         // thermal inertia
///         2146.,        // density
///         600.,         // heat capacity
///         0.07,         // albedo
///         0.9,          // emissivity
///     ),
/// );
/// ```
pub struct Body {
    /// Name of the celestial body.
    pub name: String,
    /// Position of the celestial body with respect to the Sun (in meters).
    pub position: V3<f64>,
    /// Path to the 3D object file for the shape model of the celestial body.
    pub object: Object3D,
    /// Ground properties of the celestial body (completed with time).
    pub properties: Properties,
    /// Ground temperatures of the body. Is initialized when time step is given from [`World`].
    /// Rows are faces, columns are ground layers.
    _ground_temperatures: Option<DMatrix<f64>>,
    /// Heat flux at the surface layer.
    _surface_heat_flux: VX<f64>,
    /// Heat flux at the deepest ground layer.
    _ground_limit_heat_flux: VX<f64>,
    /// Whether the body is fixed or rotating.
    _is_fixed: bool,
    /// Daily surface temperatures of a single facet during the last revolution, it is called the
    /// hour angle.
    _daily_surface_temperatures: Option<VX<f64>>,
    /// Index of the filling of the daily surface temperature vector.
    _index_current_daily: usize,
    /// Index of the temperature to pick for the hour angle.
    _index_daily_surface_temperature: usize,
    /// Whether the temperatures at the surface should be fixed (e.g. Dirichlet problem).
    _fixed_surface_temperatures: bool,
    /// Whether the temperatures at the ground limit should be fixed (e.g. Dirichlet problem).
    _fixed_ground_limit_temperatures: bool,
    /// Self view factor.
    _self_view_factor: Option<DMatrix<f64>>,
}

impl Body {
    /// Celestial body constructor.
    pub fn new<P: AsRef<Path>>(
        name: &str,
        position: V3<f64>,
        path: P,
        properties: Properties,
    ) -> Self {
        let object = Object3D::new(path);
        let number_faces_raw = object.number_faces_raw();
        let mut body = Self {
            name: String::from(name),
            position,
            object,
            properties,
            _ground_temperatures: None,
            _surface_heat_flux: VX::zeros(number_faces_raw),
            _ground_limit_heat_flux: VX::zeros(number_faces_raw),
            _is_fixed: false,
            _daily_surface_temperatures: None,
            _index_current_daily: 0,
            _index_daily_surface_temperature: 0,
            _fixed_surface_temperatures: false,
            _fixed_ground_limit_temperatures: false,
            _self_view_factor: None,
        };
        body.initiate_obliquity();
        body.recompute();
        body
    }

    pub fn empty(name: &str, position: V3<f64>, properties: Properties) -> Self {
        let mut body = Self {
            name: String::from(name),
            position,
            object: Object3D::empty(),
            properties,
            _ground_temperatures: None,
            _surface_heat_flux: VX::from_element(0, f64::NAN),
            _ground_limit_heat_flux: VX::from_element(0, f64::NAN),
            _is_fixed: false,
            _daily_surface_temperatures: None,
            _index_current_daily: 0,
            _index_daily_surface_temperature: 0,
            _fixed_surface_temperatures: false,
            _fixed_ground_limit_temperatures: false,
            _self_view_factor: None,
        };
        body.initiate_obliquity();
        body.recompute();
        body
    }

    /// Share time with Body to initialize ground vector, rotation angle, and more.
    pub fn share_time_step(&mut self, time_step: f64) {
        // Share time step with properties.
        let modified = self.properties.set_time_step(time_step);
        if modified {
            // Force recomputation of ground vector.
            self.properties.recompute_ground_vector();
            // Compute ground temperatures if not computed already.
            self.recompute_ground_temperatures(None);
        }
    }

    /// Create a new matrix for ground temperatures initialized with given value.
    pub fn new_ground_temperatures(&self, value: f64) -> DMatrix<f64> {
        DMatrix::<f64>::from_element(
            self.object.number_faces(),
            self.properties.ground_vector().len(),
            value,
        )
    }

    /// Recompute ground temperatures to match 3D object faces' mask.
    /// If the ground temperatures is already set, make it match faces' mask. But if the argument
    /// `temperature` is some, force reset of ground temperatures matrix to this value. If ground
    /// temperatures unset, initialize it with `temperatures` value or 0.
    pub fn recompute_ground_temperatures(&mut self, temperature: Option<f64>) {
        let (ground_temperatures_is_some, ground_size_match) =
            match self._ground_temperatures.as_mut() {
                Some(ground_temperatures) => (
                    true,
                    ground_temperatures.ncols() == self.properties.ground_vector().len(),
                ),
                None => (false, false),
            };
        if temperature.is_none() && ground_temperatures_is_some && ground_size_match {
            // Copy previous ground temperatures matrix to be able to copy values at faces' mask.
            let old_ground_temperatures = self._ground_temperatures.as_mut().unwrap().clone();
            // Create new ground temperatures matrix matching faces' mask and ground vector.
            let mut ground_temperatures = self.new_ground_temperatures(0.0);
            let mask = self.object.faces_mask();
            // Match 3D object faces' mask.
            for (mut ground_temperature, (_, old_ground_temperature)) in zip(
                ground_temperatures.row_iter_mut(),
                old_ground_temperatures
                    .row_iter()
                    .enumerate()
                    .filter(|&(i, _)| mask[i]),
            ) {
                ground_temperature.copy_from(&old_ground_temperature);
            }
            // Apply ground temperatures matching 3D object faces' mask.
            self._ground_temperatures = Some(ground_temperatures);
        } else {
            let initial_value = match temperature {
                Some(value) => value,
                None => 0.0,
            };
            self._ground_temperatures = Some(self.new_ground_temperatures(initial_value));
        }
    }

    /// Ground temperatures getter.
    pub fn ground_temperatures(&self) -> &DMatrix<f64> {
        self._ground_temperatures.as_ref().unwrap()
    }

    /// Ground temperatures getter as mutable.
    pub fn ground_temperatures_mut(&mut self) -> &mut DMatrix<f64> {
        self._ground_temperatures.as_mut().unwrap()
    }

    /// Ground temperatures setter.
    pub fn set_ground_temperatures(&mut self, ground_temperatures: DMatrix<f64>) {
        self._ground_temperatures = Some(ground_temperatures);
    }

    /// Surface temperatures getter.
    pub fn surface_temperatures(&self) -> MatrixSliceXx1<f64, U1, Dynamic> {
        self._ground_temperatures.as_ref().unwrap().column(0)
    }

    /// Surface temperatures getter as mutable.
    pub fn surface_temperatures_mut(&mut self) -> MatrixSliceMutXx1<f64, U1, Dynamic> {
        self._ground_temperatures.as_mut().unwrap().column_mut(0)
    }

    /// Surface temperatures setter.
    pub fn set_surface_temperatures(&mut self, surface_temperatures: VX<f64>) {
        self._ground_temperatures
            .as_mut()
            .unwrap()
            .column_mut(0)
            .copy_from(&surface_temperatures);
    }

    /// Surface temperatures getter as vector (cloned version).
    pub fn surface_temperatures_as_vec(&self) -> VX<f64> {
        VX::from_iterator(
            self.object.number_faces(),
            self._ground_temperatures
                .as_ref()
                .unwrap()
                .column(0)
                .iter()
                .cloned(),
        )
    }

    /// Surface temperatures getter.
    pub fn deepest_temperatures(&self) -> MatrixSliceXx1<f64, U1, Dynamic> {
        let ground_size = self.properties.ground_vector().len();
        self._ground_temperatures
            .as_ref()
            .unwrap()
            .column(ground_size - 1)
    }

    /// Surface heat flux getter.
    pub fn surface_heat_flux(&self) -> &VX<f64> {
        &self._surface_heat_flux
    }

    /// Set the surface heat flux.
    pub fn set_surface_heat_flux(&mut self, heat_flux: VX<f64>) {
        self._surface_heat_flux = heat_flux;
    }

    /// Set the index of the facet that record the daily surface temperatures.
    pub fn set_index_daily_surface_temperature(&mut self, index: usize) {
        self._index_daily_surface_temperature = index;
    }

    /// Check if the record of the last revolution has started, i.e. the index of the filled got
    /// updated to more than 0.
    pub fn has_last_revolution_started(&self) -> bool {
        self._index_current_daily > 0
    }

    /// Daily surface temperatures getter.
    pub fn daily_surface_temperatures(&self) -> &VX<f64> {
        self._daily_surface_temperatures.as_ref().unwrap()
    }

    /// Daily surface temperatures mutable getter.
    pub fn daily_surface_temperatures_mut(&mut self) -> &mut VX<f64> {
        self._daily_surface_temperatures.as_mut().unwrap()
    }

    /// Initiate the record of the daily surface temperatures.
    pub fn initiate_daily_record(&mut self, size: usize) {
        self._daily_surface_temperatures = Some(VX::zeros(size));
        self.set_daily_surface_temperatures();
    }

    /// Save daily surface temperature.
    pub fn set_daily_surface_temperatures(&mut self) {
        let surface_temperatures = self.surface_temperatures_as_vec();
        if let Some(daily_surface_temperatures) = &mut self._daily_surface_temperatures {
            daily_surface_temperatures[self._index_current_daily] =
                surface_temperatures[self._index_daily_surface_temperature];
            self._index_current_daily += 1;
        }
    }

    /// Iterate in the revolution of the celestial body around its spin axis over a time step.
    pub fn revolution_iteration(&mut self) {
        if self._is_fixed {
            return;
        }
        let new_normals = self.properties.rotation_matrix() * self.object.normals();
        self.object.normals_mut().copy_from(&new_normals);
    }

    /// Compute the direct solar flux at the surface of the celestial body.
    ///
    /// ## Expression
    ///
    /// U=\frac{F_s\left(1-A\right)\cos\varsigma\left(t\right)}{r_H^2\left(t\right)}
    pub fn direct_solar_flux(&self, position_sun: &V3<f64>) -> VX<f64> {
        clip(
            &(SOLAR_CONSTANT
                * (1. - self.properties.albedo())
                * self.illumination_cosine(position_sun)
                / (self.distance_to(position_sun) / ASTRONAUMICAL_UNIT).powi(2)),
            Some(0.),
            None,
        )
    }

    /// Compute the illumation angles for each face at the surface of the celestial body (angle
    /// between the normal of each face and the direction of the Sun).
    pub fn illumination_cosine(&self, position_sun: &V3<f64>) -> VX<f64> {
        let facets_directions_to_sun = self.facets_directions_to_sun(position_sun);
        dot_vectors(
            &-&facets_directions_to_sun,
            &self.object.normals().into_owned(),
        )
    }

    /// Compute the direction of the Sun for each face at the surface of the celestial body.
    pub fn facets_directions_to_sun(&self, position_sun: &V3<f64>) -> V3X<f64> {
        let mut directions = self.object.centers().into_owned();
        for mut face in directions.column_iter_mut() {
            face.copy_from(&(self.position + &face - position_sun).normalize());
        }
        directions
    }

    /// Compute the distance of the celestial body to the given position.
    pub fn distance_to(&self, other_position: &V3<f64>) -> f64 {
        (self.position - other_position).norm()
    }

    /// Update the view with the indices of the faces that are located at the equator.
    pub fn set_faces_mask_equator(&mut self) {
        self.object.set_equator_mask(None);
        self.recompute_ground_temperatures(None);
    }

    /// Initiate the obliquity by tilting the revolution axis
    fn initiate_obliquity(&mut self) {
        let old_vertices = self.object.vertices().clone();
        self.object
            .set_vertices(self.properties.obliquity_initiation_rotation_matrix() * old_vertices);
    }

    /// Add a face to the 3D object.
    pub fn new_face(&mut self, vertices: M3<f64>) {
        self.object.new_face(vertices);
        self.recompute();
    }

    /// Recompute data on the body after some changes (like number of facets or sublayers).
    pub fn recompute(&mut self) {
        self.object.recompute_faces();
        self.recompute_ground_temperatures(None);
    }

    /// Change whether the body is fixed or rotating.
    pub fn fixed(&mut self, is_fixed: bool) {
        self._is_fixed = is_fixed;
    }

    /// Get whether the temperatures at the surface are fixed (e.g. Dirichlet problem).
    pub fn fixed_surface_temperatures(&self) -> bool {
        self._fixed_surface_temperatures
    }

    /// Change whether the temperatures at the surface are fixed (e.g. Dirichlet problem).
    pub fn set_fixed_surface_temperatures(&mut self, is_fixed: bool) {
        self._fixed_surface_temperatures = is_fixed;
    }

    /// Get whether the temperatures at the ground limit is be fixed (e.g. Dirichlet problem).
    pub fn fixed_ground_limit_temperatures(&self) -> bool {
        self._fixed_ground_limit_temperatures
    }

    /// Change whether the temperatures at the ground limit is be fixed (e.g. Dirichlet problem).
    pub fn set_fixed_ground_limit_temperatures(&mut self, is_fixed: bool) {
        self._fixed_ground_limit_temperatures = is_fixed;
    }

    /// Compute the view factor of the target body with itself.
    pub fn compute_self_view_factor(&mut self) {
        let mut view_factor =
            DMatrix::<f64>::zeros(self.object.number_faces(), self.object.number_faces());
        // Cartesian product of secondary's and primary's faces.
        for ((face_target, face_other), view) in izip!(
            iproduct!(self.object.face_iter_all(), self.object.face_iter_all()),
            view_factor.iter_mut()
        ) {
            let normal_target = face_target.fixed_rows::<U3>(3);
            let normal_other = face_other.fixed_rows::<U3>(3);
            let center_target = face_target.fixed_rows::<U3>(0);
            let center_other = face_other.fixed_rows::<U3>(0);
            if center_target == center_other {
                continue;
            }
            let area_other = face_other[9];
            let relative_position_faces = center_other - center_target;
            let direction_faces = relative_position_faces.normalize();
            let distance_faces = relative_position_faces.norm();
            if distance_faces < area_other.sqrt() {
                continue;
            }
            let angle_target = normal_target.angle(&direction_faces);
            let angle_other = normal_other.angle(&-direction_faces);
            if angle_target > TAU / 4.0 || angle_other > TAU / 4.0 {
                continue;
            }
            *view = 2. * area_other * angle_target.cos() * angle_other.cos()
                / (TAU * distance_faces.powi(2));
        }
        self._self_view_factor = Some(view_factor);
    }

    /// Get the self view factor.
    pub fn self_view_factor(&mut self) -> &DMatrix<f64> {
        self._self_view_factor.as_ref().unwrap()
    }
}

/// Ground properties of celestial bodies
///
/// # Definition
///
/// [`Properties`] are initialized
///
/// # Example
///
/// ```
/// use kalast::{Properties, HOUR};
///
/// let properties = Properties::new(
///     11.92 * HOUR, // rotation period
///     11.92 * HOUR, // revolution period
///     0.,           // obliquity
///     500.,         // thermal inertia
///     2146.,        // density
///     600.,         // heat capacity
///     0.07,         // albedo
///     0.9,          // emissivity
/// );
/// ```
#[derive(Clone)]
pub struct Properties {
    /// Time to complete a rotation of the celestial body around its primary body (corresponds to a body year), given in seconds.
    _rotation_period: f64,
    /// Time to complete a rotation of the celestial around its spin axis (corresponds to a day), given in seconds.
    _revolution_period: f64,
    /// The obliquity defines the inclination of the spin axis, given in radians.
    _obliquity: f64,
    /// The thermal inertia characterize the sensitivity to temperature changes (dimensions: M.T^{-5/2}.Θ^{-1}).
    _thermal_inertia: f64,
    /// Material density of the surface (dimensions: M.L^{-3}).
    _density: f64,
    /// Heat capacity of the surface (dimensions: L^{2}.T^{-2}.Θ^{-1}).
    _heat_capacity: f64,
    /// The surface albedo defines the capacity to reflect the light.
    _albedo: f64,
    /// Surface emissivity.
    _emissivity: f64,
    /// The number of thermal skin depth to consider before adiabatic condition in the ground.
    _number_thermal_skin_depth: usize,
    /// Time step of the simulation.
    _time_step: Option<f64>,
    /// The surface conductivity (dimensions: M.L.T^{-3}.Θ^{-1}).
    _conductivity: Option<f64>,
    /// The surface diffusivity, given in (dimensions: L^{2}.T^{-1}).
    _diffusivity: Option<f64>,
    /// The thermal skin depth (dimensions: L).
    _thermal_skin_depth: Option<f64>,
    /// The ground depth is the number of thermal skin depth time the thermal skin depth
    /// (dimensions: L).
    _ground_depth: Option<f64>,
    /// The ground step (dimensions: L).
    _ground_step: Option<f64>,
    /// The vector of 1D grounds layers defined by the ground steps.
    _ground_vector: Option<VX<f64>>,
    /// The surface diffusivity times the space and time steps (without dimension).
    _diffusivity_field: Option<f64>,
    /// The axis of the revolution rotation.
    _rotation_axis: Option<Unit<V3<f64>>>,
    /// The angle for a revolution of one time step, given in radians.
    _rotation_angle: Option<f64>,
    /// The rotation matrix for one time step built from the axis and angle of rotation.
    _rotation_matrix: Option<Rotation3<f64>>,
    /// The axis for initial obliquity incination.
    _obliquity_initiation_axis: Option<Unit<V3<f64>>>,
    /// The rotation matrix for initial obliquity incination.
    _obliquity_initiation_rotation_matrix: Option<Rotation3<f64>>,
}

impl Properties {
    /// Basic properties constructor.
    pub fn new(
        rotation_period: f64,
        revolution_period: f64,
        obliquity: f64,
        thermal_inertia: f64,
        density: f64,
        heat_capacity: f64,
        albedo: f64,
        emissivity: f64,
    ) -> Self {
        let mut prop = Self {
            _rotation_period: rotation_period,
            _revolution_period: revolution_period,
            _obliquity: obliquity,
            _thermal_inertia: thermal_inertia,
            _density: density,
            _heat_capacity: heat_capacity,
            _albedo: albedo,
            _emissivity: emissivity,
            _number_thermal_skin_depth: NUMBER_THERMAL_SKIN_DEPTH,
            _time_step: None,
            _conductivity: None,
            _diffusivity: None,
            _thermal_skin_depth: None,
            _ground_depth: None,
            _ground_step: None,
            _ground_vector: None,
            _diffusivity_field: None,
            _rotation_axis: None,
            _rotation_angle: None,
            _rotation_matrix: None,
            _obliquity_initiation_axis: None,
            _obliquity_initiation_rotation_matrix: None,
        };
        // Initialize the obliquity tilt axis.
        prop.compute_obliquity_initiation_axis();
        prop.compute_obliquity_initiation_rotation_matrix();
        // Initialize the ground vector to be the surface.
        prop.compute_ground_vector();
        prop
    }

    /// Rotation period getter.
    pub fn rotation_period(&self) -> f64 {
        self._rotation_period
    }

    /// Revolution period getter.
    pub fn revolution_period(&self) -> f64 {
        self._revolution_period
    }

    /// Obliquity getter.
    pub fn obliquity(&self) -> f64 {
        self._obliquity
    }

    /// Thermal inertia getter.
    pub fn thermal_inertia(&self) -> f64 {
        self._thermal_inertia
    }

    /// Density getter.
    pub fn density(&self) -> f64 {
        self._density
    }

    /// Heat capacity getter.
    pub fn heat_capacity(&self) -> f64 {
        self._heat_capacity
    }

    /// Albedo getter.
    pub fn albedo(&self) -> f64 {
        self._albedo
    }

    /// Emissivity getter.
    pub fn emissivity(&self) -> f64 {
        self._emissivity
    }

    /// Number of thermal skin depth getter.
    pub fn number_thermal_skin_depth(&self) -> usize {
        self._number_thermal_skin_depth
    }

    /// Number of thermal skin depth setter.
    pub fn set_number_thermal_skin_depth(&mut self, number_thermal_skin_depth: usize) {
        self._number_thermal_skin_depth = number_thermal_skin_depth;
        self._ground_depth = None;
        self._ground_vector = None;
        self.compute_properties();
    }

    /// Time step getter. Fails if the time step of the simulation has not been shared with this
    /// struct.
    pub fn time_step(&self) -> f64 {
        self._time_step.expect(EXPECT_TIME_STEP_NOT_SET)
    }

    /// Time step setter. Return `true` if the time step has been modified, else `false`.
    pub fn set_time_step(&mut self, time_step: f64) -> bool {
        if self._time_step.is_some() && (self.time_step() - time_step).abs() < MINIMUM_PRECISION {
            return false;
        }
        // crate::trace!("old time step: {:?} => new: {}", self._time_step, time_step);
        self._time_step = Some(time_step);
        self._diffusivity_field = None;
        self._rotation_angle = None;
        self.compute_properties();
        true
    }

    /// Compute properties not initialized.
    pub fn compute_properties(&mut self) {
        self.compute_conductivity();
        self.compute_diffusivity();
        self.compute_thermal_skin_depth();
        self.compute_ground_depth();
        self.compute_ground_step();
        self.compute_ground_vector();
        self.compute_diffusivity_field();
        self.compute_rotation_axis();
        self.compute_rotation_angle();
        self.compute_rotation_matrix();
    }

    /// Conductivity getter.
    pub fn conductivity(&self) -> f64 {
        self._conductivity.unwrap()
    }

    /// Compute conductivity.
    pub fn compute_conductivity(&mut self) {
        if self._conductivity.is_none() {
            self._conductivity = Some(compute_conductivity(
                self.thermal_inertia(),
                self.density(),
                self.heat_capacity(),
            ));
        }
    }

    /// Diffusivity getter.
    pub fn diffusivity(&self) -> f64 {
        self._diffusivity.expect(EXPECT_TIME_STEP_NOT_SET)
    }

    /// Compute diffusivity.
    pub fn compute_diffusivity(&mut self) {
        if self._diffusivity.is_none() {
            self._diffusivity = Some(compute_diffusivity(
                self.conductivity(),
                self.density(),
                self.heat_capacity(),
            ));
        }
    }

    /// Thermal skin depth getter.
    pub fn thermal_skin_depth(&self) -> f64 {
        self._thermal_skin_depth.expect(EXPECT_TIME_STEP_NOT_SET)
    }

    /// Compute thermal skin depth.
    pub fn compute_thermal_skin_depth(&mut self) {
        if self._thermal_skin_depth.is_none() {
            self._thermal_skin_depth =
                Some((self.diffusivity() * TAU / 2. * self.revolution_period()).sqrt());
        }
    }

    /// Ground depth getter.
    pub fn ground_depth(&self) -> f64 {
        self._ground_depth.expect(EXPECT_TIME_STEP_NOT_SET)
    }

    /// Ground depth setter.
    pub fn set_ground_depth(&mut self, ground_depth: f64) {
        self._ground_depth = Some(ground_depth);
        self._ground_vector = None;
        self.compute_ground_vector();
    }

    /// Compute depth memoized.
    pub fn compute_ground_depth(&mut self) {
        if self._ground_depth.is_none() {
            self._ground_depth =
                Some(self.thermal_skin_depth() * self.number_thermal_skin_depth() as f64);
        }
    }

    /// Ground step getter.
    pub fn ground_step(&self) -> f64 {
        self._ground_step.expect(EXPECT_TIME_STEP_NOT_SET)
    }

    /// Compute ground step.
    pub fn compute_ground_step(&mut self) {
        if self._ground_step.is_none() {
            self._ground_step = Some(compute_ground_step_from_time_and_stability(
                self.time_step(),
                self.diffusivity(),
                NUMERICAL_STABILITY,
            ));
        }
    }

    /// Ground step setter.
    pub fn set_ground_step(&mut self, ground_step: f64) {
        self._ground_step = Some(ground_step);
    }

    /// Ground vector getter as reference.
    pub fn ground_vector(&self) -> &VX<f64> {
        self._ground_vector
            .as_ref()
            .expect(EXPECT_TIME_STEP_NOT_SET)
    }

    /// Force ground vector recomputation.
    pub fn recompute_ground_vector(&mut self) {
        self._ground_vector = None;
        self.compute_ground_vector();
    }

    /// Compute ground vector. Initialize to be at least the surface if no ground step or depth are
    /// set already.
    pub fn compute_ground_vector(&mut self) {
        if self._ground_vector.is_none() {
            if self._ground_depth.is_some() && self._ground_step.is_some() {
                self._ground_vector = Some(linspace(0., self.ground_depth(), self.ground_step()));
            }
            // If ground step and depth are not set yet, the ground vector is juste the surface.
            else {
                self._ground_vector = Some(VX::zeros(1));
            }
        }
    }

    /// Diffusivity field getter.
    pub fn diffusivity_field(&self) -> &f64 {
        self._diffusivity_field
            .as_ref()
            .expect(EXPECT_TIME_STEP_NOT_SET)
    }

    /// Compute diffusivity field.
    pub fn compute_diffusivity_field(&mut self) {
        if self._diffusivity_field.is_none() {
            self._diffusivity_field =
                Some(self.diffusivity() * self.time_step() / self.ground_step().powi(2));
        }
    }

    /// Rotation axis getter as reference.
    pub fn rotation_axis(&self) -> &Unit<V3<f64>> {
        self._rotation_axis
            .as_ref()
            .expect(EXPECT_TIME_STEP_NOT_SET)
    }

    /// Compute rotation axis.
    pub fn compute_rotation_axis(&mut self) {
        if self._rotation_axis.is_none() {
            self._rotation_axis = Some(Unit::new_normalize(V3::new(
                0.,
                -self.obliquity().sin(),
                self.obliquity().cos(),
            )));
            self._rotation_matrix = None;
        }
    }

    /// Rotation angle getter.
    pub fn rotation_angle(&self) -> f64 {
        self._rotation_angle.expect(EXPECT_TIME_STEP_NOT_SET)
    }

    /// Compute rotation angle.
    pub fn compute_rotation_angle(&mut self) {
        if self._rotation_angle.is_none() {
            if self.revolution_period() > 0. {
                self._rotation_angle = Some(TAU * self.time_step() / self.revolution_period());
            } else {
                self._rotation_angle = Some(0.);
            }
            self._rotation_matrix = None;
        }
    }

    /// Rotation matrix getter as reference.
    pub fn rotation_matrix(&self) -> &Rotation3<f64> {
        self._rotation_matrix
            .as_ref()
            .expect(EXPECT_TIME_STEP_NOT_SET)
    }

    /// Compute rotation matrix.
    pub fn compute_rotation_matrix(&mut self) {
        if self._rotation_matrix.is_none() {
            self._rotation_matrix = Some(Rotation3::from_axis_angle(
                self.rotation_axis(),
                self.rotation_angle(),
            ));
        }
    }

    /// Obliquity initiation axis getter as reference.
    pub fn obliquity_initiation_axis(&self) -> &Unit<V3<f64>> {
        self._obliquity_initiation_axis
            .as_ref()
            .expect("Obliquity initiation axis not set.")
    }

    /// Compute obliquity initiation axis.
    pub fn compute_obliquity_initiation_axis(&mut self) {
        if self._obliquity_initiation_axis.is_none() {
            self._obliquity_initiation_axis = Some(Unit::new_normalize(V3::new(1., 0., 0.)));
        }
    }

    /// Obliquity initiation rotation matrix getter as reference.
    pub fn obliquity_initiation_rotation_matrix(&self) -> &Rotation3<f64> {
        self._obliquity_initiation_rotation_matrix
            .as_ref()
            .expect("Obliquity initiation rotation matrix no set.")
    }

    /// Compute obliquity initiation rotation matrix.
    pub fn compute_obliquity_initiation_rotation_matrix(&mut self) {
        if self._obliquity_initiation_rotation_matrix.is_none() {
            self._obliquity_initiation_rotation_matrix = Some(Rotation3::from_axis_angle(
                self.obliquity_initiation_axis(),
                self.obliquity(),
            ));
        }
    }
}