//! # 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 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;
use std::sync::Arc;
pub type TransformationRules = (Vec<String>, Vec<String>, Vec<Expr>);
use serde::{Deserialize, Serialize};

#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Serialize, Deserialize)]
#[repr(C)]
pub enum CoordinateSystem {
    Cartesian,
    Cylindrical,
    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, or an error string
/// if the input is invalid or transformation is not supported.
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`, or an error string if the transformation
/// is not supported or rules cannot be found.
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`.
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`.
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.
/// `V'_i = (∂x'_i / ∂x_j) * 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, or an error string
/// if the transformation is not supported or computation fails.
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, or an error string
/// if the transformation is not supported or computation fails.
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]);
    }
    if let Expr::Matrix(rows) = simplify(&final_comps_expr) {
        Ok(rows.into_iter().map(|row| row[0].clone()).collect())
    } else {
        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.
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)
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, Hash)]
pub enum TensorType {
    Contravariant,
    Covariant,
    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`, or an error string
/// if the transformation is not supported or computation 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.
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,
/// or an error string if the system is not supported or computation fails.
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, or an error string
/// if the system is not supported or computation fails.
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, or an error string
/// if the system is not supported or computation fails.
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,
/// or an error string if the transformation is not supported or computation fails.
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)
}
