//! Forward modeling for MEG/EEG source localization.
//!
//! This crate provides forward models that compute the magnetic field (MEG)
//! or electric potential (EEG) at sensor locations given source dipoles
//! in the brain.
//!
//! ## Spherical Head Models
//!
//! - **MEG**: Sarvas formula for magnetic field in a spherically symmetric conductor
//! - **EEG**: Berg parameters for a three-shell spherical model
//!
//! ## Usage
//!
//! ```rust,ignore
//! use rtx_neuro_forward::{SphericalMeg, SphericalEeg, SourceSpace, SensorArray};
//!
//! // Create MEG forward model
//! let meg_model = SphericalMeg::new([0.0, 0.0, 0.04], 0.08);
//!
//! // Compute gain matrix
//! let gain = meg_model.compute_gain(&sources, &sensors)?;
//! ```

#![warn(missing_docs)]

pub mod bem;
pub mod gain;
pub mod overlapping_spheres;
pub mod sensors;
pub mod source_space;
pub mod sphere_eeg;
pub mod sphere_meg;

pub use bem::{BemConfig, BemModel, BemSurface};
pub use gain::GainMatrix;
pub use overlapping_spheres::{
    LocalSphere, OverlappingSpheres, OverlappingSpheresConfig, SurfaceMesh,
};
pub use sensors::{CoilType, MegCoil, Sensor, SensorArray, SensorType};
pub use source_space::{SourceOrientation, SourcePoint, SourceSpace};
pub use sphere_eeg::SphericalEeg;
pub use sphere_meg::SphericalMeg;

use nalgebra::Vector3;

/// Errors in forward modeling
#[derive(Debug, thiserror::Error)]
pub enum ForwardError {
    /// Invalid geometry parameters
    #[error("Invalid geometry: {0}")]
    InvalidGeometry(String),

    /// Source outside valid region
    #[error("Source outside head model: {0}")]
    SourceOutsideHead(String),

    /// Sensor configuration error
    #[error("Sensor error: {0}")]
    SensorError(String),

    /// Computation error
    #[error("Computation error: {0}")]
    ComputationError(String),

    /// Dimension mismatch
    #[error("Dimension mismatch: {0}")]
    DimensionMismatch(String),
}

/// Result type for forward modeling
pub type ForwardResult<T> = Result<T, ForwardError>;

/// 3D position in meters
pub type Position = Vector3<f64>;

/// 3D orientation (unit vector)
pub type Orientation = Vector3<f64>;

/// Cross product helper
#[inline]
fn cross(a: &Vector3<f64>, b: &Vector3<f64>) -> Vector3<f64> {
    Vector3::new(
        a.y * b.z - a.z * b.y,
        a.z * b.x - a.x * b.z,
        a.x * b.y - a.y * b.x,
    )
}

/// Dot product helper
#[inline]
fn dot(a: &Vector3<f64>, b: &Vector3<f64>) -> f64 {
    a.x * b.x + a.y * b.y + a.z * b.z
}

/// Compute the norm of a vector
#[inline]
fn norm(v: &Vector3<f64>) -> f64 {
    (v.x * v.x + v.y * v.y + v.z * v.z).sqrt()
}

/// Normalize a vector
#[inline]
fn normalize(v: &Vector3<f64>) -> Vector3<f64> {
    let n = norm(v);
    if n > 1e-15 { v / n } else { Vector3::zeros() }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_cross_product() {
        let a = Vector3::new(1.0, 0.0, 0.0);
        let b = Vector3::new(0.0, 1.0, 0.0);
        let c = cross(&a, &b);
        assert!((c.x - 0.0).abs() < 1e-10);
        assert!((c.y - 0.0).abs() < 1e-10);
        assert!((c.z - 1.0).abs() < 1e-10);
    }

    #[test]
    fn test_normalize() {
        let v = Vector3::new(3.0, 4.0, 0.0);
        let n = normalize(&v);
        assert!((norm(&n) - 1.0).abs() < 1e-10);
    }
}