rssn 0.2.9

//! # Coordinate System Transformations
//!
//! This module provides functions for transforming points, expressions, and vector/tensor
//! fields between different coordinate systems (Cartesian, Cylindrical, Spherical).
//! It includes utilities for computing Jacobian matrices and metric tensors, which are
//! fundamental for transformations in curvilinear coordinate systems.

use std::sync::Arc;

use crate::symbolic::calculus::differentiate;
use crate::symbolic::calculus::substitute;
use crate::symbolic::core::Expr;
use crate::symbolic::matrix;
use crate::symbolic::matrix::inverse_matrix;
use crate::symbolic::matrix::mul_matrices;
use crate::symbolic::matrix::transpose_matrix;
use crate::symbolic::simplify_dag::simplify;

/// A tuple representing coordinate transformation rules: (`source_vars`, `target_vars`, formulas).
pub type TransformationRules = (Vec<String>, Vec<String>, Vec<Expr>);

use serde::Deserialize;
use serde::Serialize;

/// Supported coordinate systems for symbolic transformations.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]
#[repr(C)]
pub enum CoordinateSystem {
    /// Standard 3D Cartesian coordinates (x, y, z).
    Cartesian,
    /// Cylindrical coordinates (r, theta, z).
    Cylindrical,
    /// Spherical coordinates (rho, theta, phi).
    Spherical,
}

/// Transforms a point from one coordinate system to another.
///
/// This function acts as a dispatcher, converting the point first to Cartesian coordinates
/// (if not already in Cartesian) and then from Cartesian to the target system.
///
/// # Arguments
/// * `point` - A slice of `Expr` representing the coordinates of the point.
/// * `from` - The `CoordinateSystem` of the input point.
/// * `to` - The target `CoordinateSystem`.
///
/// # Returns
/// A `Result` containing a `Vec<Expr>` of the transformed coordinates.
///
/// # Errors
///
/// This function will return an error if the input `point` has an invalid number
/// of components for the given `from` coordinate system, or if an underlying
/// transformation (e.g., `to_cartesian`, `from_cartesian`) fails.
pub fn transform_point(
    point: &[Expr],
    from: CoordinateSystem,
    to: CoordinateSystem,
) -> Result<Vec<Expr>, String> {
    if from == to {
        return Ok(point.to_vec());
    }

    let cartesian_point = to_cartesian(point, from)?;

    from_cartesian(&cartesian_point, to)
}

/// Converts a point from any system to Cartesian coordinates.
pub(crate) fn to_cartesian(
    point: &[Expr],
    from: CoordinateSystem,
) -> Result<Vec<Expr>, String> {
    match from {
        | CoordinateSystem::Cartesian => Ok(point.to_vec()),
        | CoordinateSystem::Cylindrical => {
            if point.len() != 3 {
                return Err("Cylindrical point must have 3 components (r, theta, z)".to_string());
            }

            let r = &point[0];

            let theta = &point[1];

            let z = &point[2];

            let x = simplify(&Expr::new_mul(r.clone(), Expr::new_cos(theta.clone())));

            let y = simplify(&Expr::new_mul(r.clone(), Expr::new_sin(theta.clone())));

            Ok(vec![x, y, z.clone()])
        },
        | CoordinateSystem::Spherical => {
            if point.len() != 3 {
                return Err("Spherical point must have 3 components (rho, theta, phi)".to_string());
            }

            let rho = &point[0];

            let theta = &point[1];

            let phi = &point[2];

            let x = simplify(&Expr::new_mul(
                rho.clone(),
                Expr::new_mul(Expr::new_sin(phi.clone()), Expr::new_cos(theta.clone())),
            ));

            let y = simplify(&Expr::new_mul(
                rho.clone(),
                Expr::new_mul(Expr::new_sin(phi.clone()), Expr::new_sin(theta.clone())),
            ));

            let z = simplify(&Expr::new_mul(rho.clone(), Expr::new_cos(phi.clone())));

            Ok(vec![x, y, z])
        },
    }
}

/// Converts a point from Cartesian coordinates to any other system.
pub(crate) fn from_cartesian(
    point: &[Expr],
    to: CoordinateSystem,
) -> Result<Vec<Expr>, String> {
    match to {
        | CoordinateSystem::Cartesian => Ok(point.to_vec()),
        | CoordinateSystem::Cylindrical => {
            if point.len() < 2 {
                return Err("Cartesian point must have at least 2 components (x, y)".to_string());
            }

            let x = &point[0];

            let y = &point[1];

            let r = simplify(&Expr::new_sqrt(Expr::new_add(
                Expr::new_pow(x.clone(), Expr::Constant(2.0)),
                Expr::new_pow(y.clone(), Expr::Constant(2.0)),
            )));

            let theta = simplify(&Expr::new_atan2(y.clone(), x.clone()));

            let mut result = vec![r, theta];

            if point.len() > 2 {
                result.push(point[2].clone());
            }

            Ok(result)
        },
        | CoordinateSystem::Spherical => {
            if point.len() != 3 {
                return Err("Cartesian point must have 3 components (x, y, z)".to_string());
            }

            let x = &point[0];

            let y = &point[1];

            let z = &point[2];

            let rho = simplify(&Expr::new_sqrt(Expr::new_add(
                Expr::new_add(
                    Expr::new_pow(x.clone(), Expr::Constant(2.0)),
                    Expr::new_pow(y.clone(), Expr::Constant(2.0)),
                ),
                Expr::new_pow(z.clone(), Expr::Constant(2.0)),
            )));

            let theta = simplify(&Expr::new_atan2(y.clone(), x.clone()));

            let phi = simplify(&Expr::new_arccos(Expr::new_div(z.clone(), rho.clone())));

            Ok(vec![rho, theta, phi])
        },
    }
}

/// Transforms a symbolic expression from one coordinate system to another.
///
/// This function substitutes the variables of the `from` coordinate system with their
/// equivalent expressions in the `to` coordinate system. It uses Cartesian coordinates
/// as an intermediate pivot for transformations between non-Cartesian systems.
///
/// # Arguments
/// * `expr` - The symbolic expression to transform.
/// * `from` - The `CoordinateSystem` of the input expression.
/// * `to` - The target `CoordinateSystem`.
///
/// # Returns
/// A `Result` containing the transformed `Expr`.
///
/// # Errors
///
/// This function will return an error if the transformation rules for the specified
/// `from` and `to` coordinate systems cannot be found or are not supported.
pub fn transform_expression(
    expr: &Expr,
    from: CoordinateSystem,
    to: CoordinateSystem,
) -> Result<Expr, String> {
    if from == to {
        return Ok(expr.clone());
    }

    let (from_vars, _to_vars, rules) = get_transform_rules(from, to)?;

    let mut current_expr = expr.clone();

    for (from_var, rule) in from_vars.iter().zip(rules.iter()) {
        current_expr = substitute(&current_expr, from_var, rule);
    }

    Ok(simplify(&current_expr))
}

/// Helper function to get the variables and transformation rules between two coordinate systems.
///
/// This function provides the formulas to express the `from` system's coordinates
/// in terms of the `to` system's coordinates, typically by pivoting through Cartesian.
///
/// # Arguments
/// * `from` - The source `CoordinateSystem`.
/// * `to` - The target `CoordinateSystem`.
///
/// # Returns
/// A `Result` containing a tuple `(from_vars, to_vars, rules)`.
/// `from_vars` are the variable names of the source system.
/// `to_vars` are the variable names of the target system.
/// `rules` are the expressions for `from_vars` in terms of `to_vars`.
///
/// # Errors
///
/// This function will return an error if a direct transformation between two
/// non-Cartesian systems is requested, as this is not yet supported.
pub fn get_transform_rules(
    from: CoordinateSystem,
    to: CoordinateSystem,
) -> Result<TransformationRules, String> {
    let (_to_vars, cartesian_vars, _to_cart_rules) = get_to_cartesian_rules(to)?;

    let (from_vars, _, from_cart_rules) = get_to_cartesian_rules(from)?;

    if from == CoordinateSystem::Cartesian {
        let (res_from, res_to, res_rules) = get_from_cartesian_rules(to);

        Ok((res_from, res_to, res_rules))
    } else if to == CoordinateSystem::Cartesian {
        Ok((from_vars, cartesian_vars, from_cart_rules))
    } else {
        Err("Direct transformation \
             between two \
             non-Cartesian systems is \
             not yet supported."
            .to_string())
    }
}

/// Provides the transformation rules from a given coordinate system to Cartesian coordinates.
///
/// # Arguments
/// * `from` - The source `CoordinateSystem`.
///
/// # Returns
/// A `Result` containing a tuple `(from_vars, cartesian_vars, rules)`.
/// `from_vars` are the variable names of the source system.
/// `cartesian_vars` are the variable names of the Cartesian system.
/// `rules` are the expressions for Cartesian coordinates in terms of `from_vars`.
///
/// # Errors
///
/// This function will return an error if the transformation rules for the
/// requested coordinate system are not implemented.
pub fn get_to_cartesian_rules(from: CoordinateSystem) -> Result<TransformationRules, String> {
    let cartesian_vars = vec!["x".to_string(), "y".to_string(), "z".to_string()];

    match from {
        | CoordinateSystem::Cartesian => {
            Ok((
                cartesian_vars.clone(),
                cartesian_vars,
                vec![
                    Expr::Variable("x".to_string()),
                    Expr::Variable("y".to_string()),
                    Expr::Variable("z".to_string()),
                ],
            ))
        },
        | CoordinateSystem::Cylindrical => {
            let cyl_vars = vec!["r".to_string(), "theta".to_string(), "z_cyl".to_string()];

            let r = Expr::Variable("r".to_string());

            let theta = Expr::Variable("theta".to_string());

            let rules = vec![
                simplify(&Expr::Mul(
                    Arc::new(r.clone()),
                    Arc::new(Expr::Cos(Arc::new(theta.clone()))),
                )),
                simplify(&Expr::Mul(
                    Arc::new(r),
                    Arc::new(Expr::Sin(Arc::new(theta))),
                )),
                Expr::Variable("z_cyl".to_string()),
            ];

            Ok((cyl_vars, cartesian_vars, rules))
        },
        | CoordinateSystem::Spherical => {
            let sph_vars = vec![
                "rho".to_string(),
                "theta_sph".to_string(),
                "phi".to_string(),
            ];

            let rho = Expr::Variable("rho".to_string());

            let theta = Expr::Variable("theta_sph".to_string());

            let phi = Expr::Variable("phi".to_string());

            let rules = vec![
                simplify(&Expr::Mul(
                    Arc::new(rho.clone()),
                    Arc::new(Expr::Mul(
                        Arc::new(Expr::Sin(Arc::new(phi.clone()))),
                        Arc::new(Expr::Cos(Arc::new(theta.clone()))),
                    )),
                )),
                simplify(&Expr::Mul(
                    Arc::new(rho.clone()),
                    Arc::new(Expr::Mul(
                        Arc::new(Expr::Sin(Arc::new(phi.clone()))),
                        Arc::new(Expr::Sin(Arc::new(theta))),
                    )),
                )),
                simplify(&Expr::Mul(
                    Arc::new(rho),
                    Arc::new(Expr::Cos(Arc::new(phi))),
                )),
            ];

            Ok((sph_vars, cartesian_vars, rules))
        },
    }
}

/// Provides the transformation rules from Cartesian coordinates to a given coordinate system.
///
/// # Arguments
/// * `to` - The target `CoordinateSystem`.
///
/// # Returns
/// A `Result` containing a tuple `(cartesian_vars, to_vars, rules)`.
/// `cartesian_vars` are the variable names of the Cartesian system.
/// `to_vars` are the variable names of the target system.
/// `rules` are the expressions for `to_vars` in terms of Cartesian coordinates.
pub(crate) fn get_from_cartesian_rules(to: CoordinateSystem) -> TransformationRules {
    let cartesian_vars = vec!["x".to_string(), "y".to_string(), "z".to_string()];

    let x = Expr::Variable("x".to_string());

    let y = Expr::Variable("y".to_string());

    let z = Expr::Variable("z".to_string());

    match to {
        | CoordinateSystem::Cartesian => (cartesian_vars.clone(), cartesian_vars, vec![x, y, z]),
        | CoordinateSystem::Cylindrical => {
            let cyl_vars = vec!["r".to_string(), "theta".to_string(), "z_cyl".to_string()];

            let rules = vec![
                simplify(&Expr::Sqrt(Arc::new(Expr::Add(
                    Arc::new(Expr::Power(
                        Arc::new(x.clone()),
                        Arc::new(Expr::Constant(2.0)),
                    )),
                    Arc::new(Expr::Power(
                        Arc::new(y.clone()),
                        Arc::new(Expr::Constant(2.0)),
                    )),
                )))),
                simplify(&Expr::Atan2(Arc::new(y), Arc::new(x))),
                z,
            ];

            (cartesian_vars, cyl_vars, rules)
        },
        | CoordinateSystem::Spherical => {
            let sph_vars = vec![
                "rho".to_string(),
                "theta_sph".to_string(),
                "phi".to_string(),
            ];

            let rho_rule = simplify(&Expr::new_sqrt(Expr::new_add(
                Expr::new_add(
                    Expr::new_pow(x.clone(), Expr::Constant(2.0)),
                    Expr::new_pow(y.clone(), Expr::Constant(2.0)),
                ),
                Expr::new_pow(z.clone(), Expr::Constant(2.0)),
            )));

            let rules = vec![
                rho_rule.clone(),
                simplify(&Expr::Atan2(Arc::new(y), Arc::new(x))),
                simplify(&Expr::ArcCos(Arc::new(Expr::Div(
                    Arc::new(z),
                    Arc::new(rho_rule),
                )))),
            ];

            (cartesian_vars, sph_vars, rules)
        },
    }
}

/// Transforms a contravariant vector field (e.g., velocity) from one coordinate system to another.
///
/// Contravariant vectors transform with the Jacobian matrix of the coordinate transformation.
///
/// # Returns
/// A `Result` containing a `Vec<Expr>` of the transformed components.
///
/// # Errors
///
/// This function will return an error if:
/// - The `get_from_cartesian_rules` function fails to provide transformation rules.
/// - The symbolic matrix-vector multiplication fails due to incompatible dimensions.
pub fn transform_contravariant_vector(
    components: &[Expr],
    from: CoordinateSystem,
    to: CoordinateSystem,
) -> Result<Vec<Expr>, String> {
    if from == to {
        return Ok(components.to_vec());
    }

    let (vars_from, _, rules_to) = get_from_cartesian_rules(from);

    let jacobian = compute_jacobian(&rules_to, &vars_from);

    let new_comps_old_coords = symbolic_mat_vec_mul(&jacobian, components)?;

    let (_, _, rules_from) = get_to_cartesian_rules(from)?;

    let mut final_comps = Vec::new();

    for comp in new_comps_old_coords {
        let mut final_comp = comp;

        for (i, var) in vars_from.iter().enumerate() {
            final_comp = substitute(&final_comp, var, &rules_from[i]);
        }

        final_comps.push(simplify(&final_comp));
    }

    Ok(final_comps)
}

/// Transforms a covariant vector field (e.g., gradient) from one coordinate system to another.
///
/// Covariant vectors transform with the inverse transpose of the Jacobian matrix of the
/// coordinate transformation. `V'_i = (∂x_j / ∂x'_i) * V_j`.
///
/// # Arguments
/// * `components` - A slice of `Expr` representing the components of the vector field.
/// * `from` - The source `CoordinateSystem`.
/// * `to` - The target `CoordinateSystem`.
///
/// # Returns
/// A `Result` containing a `Vec<Expr>` of the transformed components.
///
/// # Errors
///
/// This function will return an error if:
/// - The `get_transform_rules` function fails to provide transformation rules.
/// - The Jacobian matrix is singular and cannot be inverted.
/// - The symbolic matrix multiplication results in a non-vector expression,
///   indicating an issue with the transformation or simplification.
pub fn transform_covariant_vector(
    components: &[Expr],
    from: CoordinateSystem,
    to: CoordinateSystem,
) -> Result<Vec<Expr>, String> {
    if from == to {
        return Ok(components.to_vec());
    }

    let (vars_from, _, rules_to) = get_transform_rules(from, to)?;

    let jacobian_vec = compute_jacobian(&rules_to, &vars_from);

    let jacobian = Expr::Matrix(jacobian_vec);

    let jacobian_inv = inverse_matrix(&jacobian);

    let jacobian_inv_t = transpose_matrix(&jacobian_inv);

    let old_vec = Expr::Matrix(components.iter().map(|c| vec![c.clone()]).collect());

    let new_vec_expr = mul_matrices(&jacobian_inv_t, &old_vec);

    let (from_vars, _, rules) = get_to_cartesian_rules(from)?;

    let mut final_comps_expr = new_vec_expr;

    for (i, var) in from_vars.iter().enumerate() {
        final_comps_expr = substitute(&final_comps_expr, var, &rules[i]);
    }

    match simplify(&final_comps_expr) {
        | Expr::Matrix(rows) => Ok(rows.into_iter().map(|row| row[0].clone()).collect()),
        | _ => {
            Err("Transformation \
                 resulted in a \
                 non-vector expression"
                .to_string())
        },
    }
}

/// Computes the Jacobian matrix for a set of transformation rules.
pub(crate) fn compute_jacobian(
    rules: &[Expr],
    vars: &[String],
) -> Vec<Vec<Expr>> {
    let mut jacobian = Vec::new();

    for rule in rules {
        let mut row = Vec::new();

        for var in vars {
            row.push(differentiate(rule, var));
        }

        jacobian.push(row);
    }

    jacobian
}

/// Performs symbolic matrix-vector multiplication.
///
/// # Errors
///
/// This function will return an error if the dimensions of the `matrix` and `vector`
/// are incompatible for multiplication.
pub(crate) fn symbolic_mat_vec_mul(
    matrix: &[Vec<Expr>],
    vector: &[Expr],
) -> Result<Vec<Expr>, String> {
    if matrix.is_empty() || (!matrix.is_empty() && matrix[0].len() != vector.len()) {
        return Err("Matrix and \
                    vector dimensions \
                    are incompatible.\
                    "
        .to_string());
    }

    let mut result = Vec::new();

    for row in matrix {
        let mut sum = Expr::Constant(0.0);

        for (i, val) in row.iter().enumerate() {
            sum = simplify(&Expr::new_add(
                sum,
                Expr::new_mul(val.clone(), vector[i].clone()),
            ));
        }

        result.push(sum);
    }

    Ok(result)
}

/// Classification of tensor transformation behavior.
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub enum TensorType {
    /// Transforms with the Jacobian matrix.
    Contravariant,
    /// Transforms with the inverse transpose of the Jacobian.
    Covariant,
    /// Transforms with both Jacobian and its inverse.
    Mixed,
}

/// Transforms a rank-2 tensor field from one coordinate system to another.
///
/// The transformation rule depends on the tensor's type (contravariant, covariant, or mixed).
/// - **Contravariant**: `T'^ij = (∂x'^i/∂x^a)(∂x'^j/∂x^b) T^ab`
/// - **Covariant**: `T'_ij = (∂x^a/∂x'^i)(∂x^b/∂x'^j) T_ab`
/// - **Mixed**: `T'^i_j = (∂x'^i/∂x^a)(∂x^b/∂x'^j) T^a_b`
///
/// # Arguments
/// * `tensor` - The rank-2 tensor as an `Expr::Matrix`.
/// * `from` - The source `CoordinateSystem`.
/// * `to` - The target `CoordinateSystem`.
/// * `tensor_type` - The `TensorType` (Contravariant, Covariant, Mixed).
///
/// # Returns
/// A `Result` containing the transformed `Expr::Matrix`.
///
/// # Errors
///
/// This function will return an error if:
/// - The `get_transform_rules` function fails to provide transformation rules.
/// - The Jacobian matrix is singular and cannot be inverted.
/// - An underlying symbolic matrix multiplication operation fails.
pub fn transform_tensor2(
    tensor: &Expr,
    from: CoordinateSystem,
    to: CoordinateSystem,
    tensor_type: TensorType,
) -> Result<Expr, String> {
    if from == to {
        return Ok(tensor.clone());
    }

    let (from_vars, _, rules) = get_transform_rules(from, to)?;

    let jacobian_vec = compute_jacobian(&rules, &from_vars);

    let jacobian = Expr::Matrix(jacobian_vec);

    let jacobian_inv = inverse_matrix(&jacobian);

    let transformed_tensor = match tensor_type {
        | TensorType::Contravariant => {
            mul_matrices(
                &jacobian,
                &mul_matrices(tensor, &transpose_matrix(&jacobian)),
            )
        },
        | TensorType::Covariant => {
            mul_matrices(
                &transpose_matrix(&jacobian_inv),
                &mul_matrices(tensor, &jacobian_inv),
            )
        },
        | TensorType::Mixed => mul_matrices(&jacobian, &mul_matrices(tensor, &jacobian_inv)),
    };

    Ok(transformed_tensor)
}

/// Performs symbolic matrix-matrix multiplication.
///
/// # Errors
///
/// This function will return an error if the dimensions of `m1` and `m2` are
/// incompatible for multiplication (i.e., number of columns in `m1` does not
/// equal number of rows in `m2`).
pub fn symbolic_mat_mat_mul(
    m1: &[Vec<Expr>],
    m2: &[Vec<Expr>],
) -> Result<Vec<Vec<Expr>>, String> {
    let m1_rows = m1.len();

    if m1_rows == 0 {
        return Ok(vec![]);
    }

    let m1_cols = m1[0].len();

    let m2_rows = m2.len();

    if m2_rows == 0 {
        return Ok(vec![]);
    }

    let m2_cols = m2[0].len();

    if m1_cols != m2_rows {
        return Err("Matrix dimensions are \
             incompatible for \
             multiplication."
            .to_string());
    }

    let mut result = vec![vec![Expr::Constant(0.0); m2_cols]; m1_rows];

    for (i, row) in m1.iter().enumerate() {
        for j in 0..m2_cols {
            let mut sum = Expr::Constant(0.0);

            for (k, val) in row.iter().enumerate() {
                sum = simplify(&Expr::new_add(
                    sum,
                    Expr::new_mul(val.clone(), m2[k][j].clone()),
                ));
            }

            result[i][j] = sum;
        }
    }

    Ok(result)
}

/// Computes and returns the metric tensor for a given orthogonal coordinate system.
///
/// The metric tensor `g_ij` defines the inner product in a coordinate system and is
/// crucial for calculating distances, angles, and volumes. For orthogonal coordinates,
/// it is a diagonal matrix where `g_ii = (∂x/∂u_i)^2 + (∂y/∂u_i)^2 + (∂z/∂u_i)^2`.
///
/// # Arguments
/// * `system` - The `CoordinateSystem` for which to compute the metric tensor.
///
/// # Returns
/// A `Result` containing an `Expr::Matrix` representing the metric tensor.
///
/// # Errors
///
/// This function will return an error if `get_to_cartesian_rules` fails to provide
/// transformation rules for the given coordinate system.
pub fn get_metric_tensor(system: CoordinateSystem) -> Result<Expr, String> {
    let rules = match system {
        | CoordinateSystem::Cartesian => return Ok(matrix::identity_matrix(3)),
        | _ => get_to_cartesian_rules(system)?.2,
    };

    let vars = match system {
        | CoordinateSystem::Cylindrical => {
            vec!["r".to_string(), "theta".to_string(), "z_cyl".to_string()]
        },
        | CoordinateSystem::Spherical => {
            vec![
                "rho".to_string(),
                "theta_sph".to_string(),
                "phi".to_string(),
            ]
        },
        | CoordinateSystem::Cartesian => unreachable!(),
    };

    let jacobian_vec = compute_jacobian(&rules, &vars);

    let jacobian = Expr::Matrix(jacobian_vec);

    Ok(mul_matrices(&transpose_matrix(&jacobian), &jacobian))
}

/// Computes the divergence of a contravariant vector field in any orthogonal coordinate system.
///
/// The divergence of a vector field `V` in orthogonal curvilinear coordinates `u_1, u_2, u_3` is given by:
/// `div(V) = (1/sqrt(g)) * [∂/∂u_1(V^1 * sqrt(g)/h_1) + ∂/∂u_2(V^2 * sqrt(g)/h_2) + ∂/∂u_3(V^3 * sqrt(g)/h_3)]`
/// where `g` is the determinant of the metric tensor and `h_i` are scale factors.
/// This simplified implementation assumes `h_i` are implicitly handled by `sqrt(g)`.
///
/// # Arguments
/// * `vector_comps` - A slice of `Expr` representing the contravariant components of the vector field.
/// * `from` - The `CoordinateSystem` of the vector field.
///
/// # Returns
/// A `Result` containing an `Expr` representing the divergence.
///
/// # Errors
///
/// This function will return an error if `get_metric_tensor` fails to retrieve
/// the metric tensor for the given coordinate system.
pub fn transform_divergence(
    vector_comps: &[Expr],
    from: CoordinateSystem,
) -> Result<Expr, String> {
    let g_matrix = get_metric_tensor(from)?;

    let g = matrix::determinant(&g_matrix);

    let sqrt_g = simplify(&Expr::new_sqrt(g));

    let (vars, _, _) = get_to_cartesian_rules(from)?;

    let mut total_divergence = Expr::Constant(0.0);

    for i in 0..vector_comps.len() {
        let term_to_diff = simplify(&Expr::new_mul(sqrt_g.clone(), vector_comps[i].clone()));

        let partial_deriv = differentiate(&term_to_diff, &vars[i]);

        total_divergence = simplify(&Expr::new_add(total_divergence, partial_deriv));
    }

    Ok(simplify(&Expr::new_div(total_divergence, sqrt_g)))
}

/// Computes the curl of a covariant vector field in any orthogonal coordinate system.
///
/// The curl of a vector field `V` in orthogonal curvilinear coordinates `u_1, u_2, u_3` is given by:
/// `curl(V) = (1/(h_1*h_2*h_3)) * det(h_1*e_1, h_2*e_2, h_3*e_3; ∂/∂u_1, ∂/∂u_2, ∂/∂u_3; h_1*V_1, h_2*V_2, h_3*V_3)`
/// where `h_i` are the scale factors.
///
/// # Arguments
/// * `vector_comps` - A slice of `Expr` representing the covariant components of the vector field.
/// * `from` - The `CoordinateSystem` of the vector field.
///
/// # Returns
/// A `Result` containing a `Vec<Expr>` of the transformed components.
///
/// # Errors
///
/// This function will return an error if:
/// - The `vector_comps` length is not 3 (curl is only defined for 3D vectors).
/// - `get_metric_tensor` fails to retrieve the metric tensor.
/// - The metric tensor is not an `Expr::Matrix`.
pub fn transform_curl(
    vector_comps: &[Expr],
    from: CoordinateSystem,
) -> Result<Vec<Expr>, String> {
    if vector_comps.len() != 3 {
        return Err("Curl is only \
                    defined for 3D \
                    vectors."
            .to_string());
    }

    let g_matrix = get_metric_tensor(from)?;

    let g_rows = if let Expr::Matrix(rows) = g_matrix {
        rows
    } else {
        return Err("Metric tensor is not \
                 a matrix"
            .to_string());
    };

    let h1 = simplify(&Expr::new_sqrt(g_rows[0][0].clone()));

    let h2 = simplify(&Expr::new_sqrt(g_rows[1][1].clone()));

    let h3 = simplify(&Expr::new_sqrt(g_rows[2][2].clone()));

    let (vars, _, _) = get_to_cartesian_rules(from)?;

    let u1 = &vars[0];

    let u2 = &vars[1];

    let u3 = &vars[2];

    let v1 = &vector_comps[0];

    let v2 = &vector_comps[1];

    let v3 = &vector_comps[2];

    let curl_1 = simplify(&Expr::new_div(
        Expr::new_sub(
            differentiate(&simplify(&Expr::new_mul(h3.clone(), v3.clone())), u2),
            differentiate(&simplify(&Expr::new_mul(h2.clone(), v2.clone())), u3),
        ),
        simplify(&Expr::new_mul(h2.clone(), h3.clone())),
    ));

    let curl_2 = simplify(&Expr::new_div(
        Expr::new_sub(
            differentiate(&simplify(&Expr::new_mul(h1.clone(), v1.clone())), u3),
            differentiate(&simplify(&Expr::new_mul(h3.clone(), v3.clone())), u1),
        ),
        simplify(&Expr::new_mul(h3, h1.clone())),
    ));

    let curl_3 = simplify(&Expr::new_div(
        Expr::new_sub(
            differentiate(&simplify(&Expr::new_mul(h2.clone(), v2.clone())), u1),
            differentiate(&simplify(&Expr::new_mul(h1.clone(), v1.clone())), u2),
        ),
        simplify(&Expr::new_mul(h1, h2)),
    ));

    Ok(vec![curl_1, curl_2, curl_3])
}

/// Transforms the gradient of a scalar field from one coordinate system to another.
///
/// The gradient of a scalar field `f` is a covariant vector field. This function
/// computes the gradient in the source system, then transforms it to Cartesian coordinates,
/// and finally to the target coordinate system.
///
/// # Arguments
/// * `scalar_field` - The scalar field as an `Expr`.
/// * `from_vars` - A slice of strings representing the variables of the source system.
/// * `from` - The source `CoordinateSystem`.
/// * `to` - The target `CoordinateSystem`.
///
/// # Returns
/// A `Result` containing a `Vec<Expr>` of the transformed gradient components.
///
/// # Errors
///
/// This function will return an error if:
/// - The transformation rules for the `from` system cannot be retrieved.
/// - The underlying covariant vector transformation fails (e.g., if the Jacobian
///   matrix is singular or dimensions are incompatible).
pub fn transform_gradient(
    scalar_field: &Expr,
    from_vars: &[String],
    from: CoordinateSystem,
    to: CoordinateSystem,
) -> Result<Vec<Expr>, String> {
    if from == to {
        let mut grad_comps = Vec::new();

        for var in from_vars {
            grad_comps.push(differentiate(scalar_field, var));
        }

        return Ok(grad_comps);
    }

    let (_, _, rules) = get_to_cartesian_rules(from)?;

    let mut field_cart = scalar_field.clone();

    for (i, var) in from_vars.iter().enumerate() {
        field_cart = substitute(&field_cart, var, &rules[i]);
    }

    field_cart = simplify(&field_cart);

    let cartesian_vars = vec!["x".to_string(), "y".to_string(), "z".to_string()];

    let mut grad_cart_comps = Vec::new();

    for var in &cartesian_vars {
        grad_cart_comps.push(differentiate(&field_cart, var));
    }

    transform_covariant_vector(&grad_cart_comps, CoordinateSystem::Cartesian, to)
}