rssn 0.2.9

//! # Numerical Coordinate Transformations
//!
//! This module provides numerical tools for coordinate transformations.
//! It supports converting points between various coordinate systems (Cartesian, Cylindrical, Spherical)
//! and computing numerical Jacobians of these transformations.

use std::collections::HashMap;

use crate::numerical::calculus::gradient;
use crate::numerical::matrix::Matrix;
use crate::symbolic::coordinates::CoordinateSystem;
use crate::symbolic::coordinates::{
    self,
};
use crate::symbolic::core::Expr;

/// Transforms a numerical point from one coordinate system to another.
///
/// This function leverages symbolic transformation rules and then numerically evaluates
/// the resulting expressions. For performance-critical applications, consider using
/// `transform_point_pure`.
///
/// # Arguments
/// * `point` - The point to transform as a slice of `f64` coordinates.
/// * `from` - The `CoordinateSystem` of the input point.
/// * `to` - The target `CoordinateSystem`.
///
/// # Returns
/// A `Result` containing a `Vec<f64>` of the transformed coordinates.
///
/// # Errors
///
/// Returns an error if the transformation rules cannot be found or if expression evaluation fails.
pub fn transform_point(
    point: &[f64],
    from: CoordinateSystem,
    to: CoordinateSystem,
) -> Result<Vec<f64>, String> {
    let point_expr: Vec<Expr> = point.iter().map(|&v| Expr::Constant(v)).collect();

    let transformed_expr = coordinates::transform_point(&point_expr, from, to)?;

    let mut result = Vec::new();

    for expr in transformed_expr {
        result.push(crate::numerical::elementary::eval_expr(
            &expr,
            &HashMap::new(),
        )?);
    }

    Ok(result)
}

/// Computes the numerical Jacobian matrix of a coordinate transformation at a specific point.
///
/// The Jacobian matrix `J` contains the partial derivatives of the new coordinates with
/// respect to the old coordinates: `J_ij = ∂(new_coord_i) / ∂(old_coord_j)`.
/// This is computed using symbolic differentiation followed by numerical evaluation.
///
/// # Arguments
/// * `from` - The source `CoordinateSystem`.
/// * `to` - The target `CoordinateSystem`.
/// * `at_point` - The point at which to evaluate the Jacobian.
///
/// # Returns
/// A `Result` containing a `Matrix<f64>` representing the Jacobian matrix.
///
/// # Errors
///
/// Returns an error if the transformation rules cannot be found or if numerical gradient computation fails.
pub fn numerical_jacobian(
    from: CoordinateSystem,
    to: CoordinateSystem,
    at_point: &[f64],
) -> Result<Matrix<f64>, String> {
    let (from_vars, _, rules) = coordinates::get_transform_rules(from, to)?;

    let mut jacobian_rows = Vec::new();

    for rule in &rules {
        let grad = gradient(
            rule,
            &from_vars
                .iter()
                .map(std::string::String::as_str)
                .collect::<Vec<_>>(),
            at_point,
        )?;

        jacobian_rows.push(grad);
    }

    let rows = jacobian_rows.len();

    let cols = if rows > 0 {
        jacobian_rows[0].len()
    } else {
        0
    };

    Ok(Matrix::new(rows, cols, jacobian_rows.concat()))
}

/// Transforms a numerical point using direct `f64` calculations for high performance.
///
/// This function provides a more optimized approach for coordinate transformations
/// by directly applying numerical formulas without intermediate symbolic manipulation.
///
/// # Arguments
/// * `point` - The point to transform as a slice of `f64` coordinates.
/// * `from` - The `CoordinateSystem` of the input point.
/// * `to` - The target `CoordinateSystem`.
///
/// # Returns
/// A `Result` containing a `Vec<f64>` of the transformed coordinates.
///
/// # Errors
///
/// Returns an error if the input point has an incorrect number of components for the specified coordinate systems.
pub fn transform_point_pure(
    point: &[f64],
    from: CoordinateSystem,
    to: CoordinateSystem,
) -> Result<Vec<f64>, String> {
    if from == to {
        return Ok(point.to_vec());
    }

    let cartesian_point = to_cartesian_pure(point, from)?;

    from_cartesian_pure(&cartesian_point, to)
}

pub(crate) fn to_cartesian_pure(
    point: &[f64],
    from: CoordinateSystem,
) -> Result<Vec<f64>, String> {
    /// Converts a numerical point from a given coordinate system to Cartesian coordinates.
    ///
    /// This is a helper function for `transform_point_pure`.
    ///
    /// # Arguments
    /// * `point` - The point to convert.
    /// * `from` - The source `CoordinateSystem`.
    ///
    /// # Returns
    /// A `Result` containing a `Vec<f64>` of the Cartesian coordinates, or an error 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 = r * theta.cos();

            let y = r * theta.sin();

            Ok(vec![x, y, z])
        },
        | 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 = rho * phi.sin() * theta.cos();

            let y = rho * phi.sin() * theta.sin();

            let z = rho * phi.cos();

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

pub(crate) fn from_cartesian_pure(
    point: &[f64],
    to: CoordinateSystem,
) -> Result<Vec<f64>, String> {
    /// Converts a numerical point from Cartesian coordinates to a given target coordinate system.
    ///
    /// This is a helper function for `transform_point_pure`.
    ///
    /// # Arguments
    /// * `point` - The Cartesian point to convert.
    /// * `to` - The target `CoordinateSystem`.
    ///
    /// # Returns
    /// A `Result` containing a `Vec<f64>` of the transformed coordinates, or an error 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 = x.hypot(y);

            let theta = y.atan2(x);

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

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

            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 = z.mul_add(z, y.mul_add(y, x.powi(2))).sqrt();

            let theta = y.atan2(x);

            let phi = (z / rho).acos();

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