#[cfg(feature = "rayon")]
use crate::adsorption::fea_potential::calculate_fea_potential;
use crate::functional::HelmholtzEnergyFunctional;
#[cfg(feature = "rayon")]
use crate::geometry::Geometry;
use libm::tgamma;
use nalgebra::DVector;
use ndarray::{Array1, Array2, Axis as Axis_nd};
#[cfg(feature = "rayon")]
use quantity::Length;
use std::f64::consts::PI;
use std::ops::Deref;

const DELTA_STEELE: f64 = 3.35;

/// A collection of external potentials.
#[derive(Clone)]
pub enum ExternalPotential {
    /// Hard wall potential: $V_i^\mathrm{ext}(z)=\begin{cases}\infty&z\leq\sigma_{si}\\\\0&z>\sigma_{si}\end{cases},~~~~\sigma_{si}=\frac{1}{2}\left(\sigma_{ss}+\sigma_{ii}\right)$
    HardWall { sigma_ss: f64 },
    /// 9-3 Lennard-Jones potential: $V_i^\mathrm{ext}(z)=\frac{2\pi}{45} m_i\varepsilon_{si}\sigma_{si}^3\rho_s\left(2\left(\frac{\sigma_{si}}{z}\right)^9-15\left(\frac{\sigma_{si}}{z}\right)^3\right),~~~~\varepsilon_{si}=\sqrt{\varepsilon_{ss}\varepsilon_{ii}},~~~~\sigma_{si}=\frac{1}{2}\left(\sigma_{ss}+\sigma_{ii}\right)$
    LJ93 {
        sigma_ss: f64,
        epsilon_k_ss: f64,
        rho_s: f64,
    },
    /// Simple 9-3 Lennard-Jones potential: $V_i^\mathrm{ext}(z)=\varepsilon_{si}\left(\left(\frac{\sigma_{si}}{z}\right)^9-\left(\frac{\sigma_{si}}{z}\right)^3\right),~~~~\varepsilon_{si}=\sqrt{\varepsilon_{ss}\varepsilon_{ii}},~~~~\sigma_{si}=\frac{1}{2}\left(\sigma_{ss}+\sigma_{ii}\right)$
    SimpleLJ93 { sigma_ss: f64, epsilon_k_ss: f64 },
    /// Custom 9-3 Lennard-Jones potential: $V_i^\mathrm{ext}(z)=\varepsilon_{si}\left(\left(\frac{\sigma_{si}}{z}\right)^9-\left(\frac{\sigma_{si}}{z}\right)^3\right)$
    CustomLJ93 {
        sigma_sf: Array1<f64>,
        epsilon_k_sf: Array1<f64>,
    },
    /// Steele potential: $V_i^\mathrm{ext}(z)=2\pi m_i\xi\varepsilon_{si}\sigma_{si}^2\Delta\rho_s\left(0.4\left(\frac{\sigma_{si}}{z}\right)^{10}-\left(\frac{\sigma_{si}}{z}\right)^4-\frac{\sigma_{si}^4}{3\Delta\left(z+0.61\Delta\right)^3}\right),~~~~\varepsilon_{si}=\sqrt{\varepsilon_{ss}\varepsilon_{ii}},~~~~\sigma_{si}=\frac{1}{2}\left(\sigma_{ss}+\sigma_{ii}\right),~~~~\Delta=3.35$
    Steele {
        sigma_ss: f64,
        epsilon_k_ss: f64,
        rho_s: f64,
        xi: Option<f64>,
    },
    /// Steele potential with custom combining rules: $V_i^\mathrm{ext}(z)=2\pi m_i\xi\varepsilon_{si}\sigma_{si}^2\Delta\rho_s\left(0.4\left(\frac{\sigma_{si}}{z}\right)^{10}-\left(\frac{\sigma_{si}}{z}\right)^4-\frac{\sigma_{si}^4}{3\Delta\left(z+0.61\Delta\right)^3}\right),~~~~\Delta=3.35$
    CustomSteele {
        sigma_sf: Array1<f64>,
        epsilon_k_sf: Array1<f64>,
        rho_s: f64,
        xi: Option<f64>,
    },
    /// Double well potential: $V_i^\mathrm{ext}(z)=\mathrm{min}\left(\frac{2\pi}{45} m_i\varepsilon_{2si}\sigma_{si}^3\rho_s\left(2\left(\frac{2\sigma_{si}}{z}\right)^9-15\left(\frac{2\sigma_{si}}{z}\right)^3\right),0\right)+\frac{2\pi}{45} m_i\varepsilon_{1si}\sigma_{si}^3\rho_s\left(2\left(\frac{\sigma_{si}}{z}\right)^9-15\left(\frac{\sigma_{si}}{z}\right)^3\right),~~~~\varepsilon_{1si}=\sqrt{\varepsilon_{1ss}\varepsilon_{ii}},~~~~\varepsilon_{2si}=\sqrt{\varepsilon_{2ss}\varepsilon_{ii}},~~~~\sigma_{si}=\frac{1}{2}\left(\sigma_{ss}+\sigma_{ii}\right)$
    DoubleWell {
        sigma_ss: f64,
        epsilon1_k_ss: f64,
        epsilon2_k_ss: f64,
        rho_s: f64,
    },
    /// Free-energy averaged potential:
    #[cfg(feature = "rayon")]
    FreeEnergyAveraged {
        coordinates: Length<Array2<f64>>,
        sigma_ss: Array1<f64>,
        epsilon_k_ss: Array1<f64>,
        pore_center: [f64; 3],
        system_size: [Length; 3],
        n_grid: [usize; 2],
        cutoff_radius: Option<f64>,
    },

    /// Custom potential
    Custom(Array2<f64>),
}

/// Parameters of the fluid required to evaluate the external potential.
pub trait FluidParameters {
    fn epsilon_k_ff(&self) -> DVector<f64>;
    fn sigma_ff(&self) -> DVector<f64>;
}

impl<C: Deref<Target = T>, T: FluidParameters> FluidParameters for C {
    fn epsilon_k_ff(&self) -> DVector<f64> {
        T::epsilon_k_ff(self)
    }
    fn sigma_ff(&self) -> DVector<f64> {
        T::sigma_ff(self)
    }
}

impl ExternalPotential {
    // Evaluate the external potential in cartesian coordinates for a given grid and fluid parameters.
    pub fn calculate_cartesian_potential<P: HelmholtzEnergyFunctional + FluidParameters>(
        &self,
        z_grid: &Array1<f64>,
        fluid_parameters: &P,
        #[cfg_attr(not(feature = "rayon"), expect(unused_variables))] temperature: f64,
    ) -> Array2<f64> {
        if let ExternalPotential::Custom(potential) = self {
            return potential.clone();
        }

        // Allocate external potential
        let m = fluid_parameters.m();
        let mut ext_pot = Array2::zeros((m.len(), z_grid.len()));

        for (i, &mi) in m.iter().enumerate() {
            ext_pot.index_axis_mut(Axis_nd(0), i).assign(&match self {
                Self::HardWall { sigma_ss } => {
                    let sigma_sf = (fluid_parameters.sigma_ff()[i] + *sigma_ss) * 0.5;
                    z_grid.mapv(|z| if z < sigma_sf { f64::INFINITY } else { 0.0 })
                }
                Self::LJ93 {
                    sigma_ss,
                    epsilon_k_ss,
                    rho_s,
                } => {
                    // combining rules
                    let epsilon_k_sf =
                        (fluid_parameters.epsilon_k_ff() * *epsilon_k_ss).map(|e| e.sqrt());
                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;

                    2.0 * PI * mi * epsilon_k_sf[i] * sigma_sf[i].powi(3) * rho_s / 45.0
                        * (2.0 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(9))
                            - 15.0 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(3)))
                }
                Self::SimpleLJ93 {
                    sigma_ss,
                    epsilon_k_ss,
                } => {
                    // combining rules
                    let epsilon_k_sf =
                        (fluid_parameters.epsilon_k_ff() * *epsilon_k_ss).map(|e| e.sqrt());
                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;

                    epsilon_k_sf[i]
                        * ((sigma_sf[i] / z_grid).mapv(|x| x.powi(9))
                            - (sigma_sf[i] / z_grid).mapv(|x| x.powi(3)))
                }
                Self::CustomLJ93 {
                    sigma_sf,
                    epsilon_k_sf,
                } => {
                    epsilon_k_sf[i]
                        * ((sigma_sf[i] / z_grid).mapv(|x| x.powi(9))
                            - (sigma_sf[i] / z_grid).mapv(|x| x.powi(3)))
                }
                Self::Steele {
                    sigma_ss,
                    epsilon_k_ss,
                    rho_s,
                    xi,
                } => {
                    // combining rules
                    let epsilon_k_sf =
                        (fluid_parameters.epsilon_k_ff() * *epsilon_k_ss).map(|e| e.sqrt());
                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;

                    (2.0 * PI * mi * xi.unwrap_or(1.0) * epsilon_k_sf[i])
                        * (sigma_sf[i].powi(2) * DELTA_STEELE * rho_s)
                        * (0.4 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(10))
                            - (sigma_sf[i] / z_grid).mapv(|x| x.powi(4))
                            - sigma_sf[i].powi(4)
                                / ((3.0 * DELTA_STEELE)
                                    * (z_grid + 0.61 * DELTA_STEELE).mapv(|x| x.powi(3))))
                }
                Self::CustomSteele {
                    sigma_sf,
                    epsilon_k_sf,
                    rho_s,
                    xi,
                } => {
                    (2.0 * PI * mi * xi.unwrap_or(1.0) * epsilon_k_sf[i])
                        * (sigma_sf[i].powi(2) * DELTA_STEELE * rho_s)
                        * (0.4 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(10))
                            - (sigma_sf[i] / z_grid).mapv(|x| x.powi(4))
                            - sigma_sf[i].powi(4)
                                / ((3.0 * DELTA_STEELE)
                                    * (z_grid + 0.61 * DELTA_STEELE).mapv(|x| x.powi(3))))
                }
                Self::DoubleWell {
                    sigma_ss,
                    epsilon1_k_ss,
                    epsilon2_k_ss,
                    rho_s,
                } => {
                    // combining rules
                    let epsilon1_k_sf =
                        (fluid_parameters.epsilon_k_ff() * *epsilon1_k_ss).map(|e| e.sqrt());
                    let epsilon2_k_sf =
                        (fluid_parameters.epsilon_k_ff() * *epsilon2_k_ss).map(|e| e.sqrt());
                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;

                    let bh_tail = (2.0 * PI * mi * epsilon2_k_sf[i] * sigma_sf[i].powi(3) * rho_s
                        / 45.0
                        * (2.0 * (2.0 * sigma_sf[i] / z_grid).mapv(|x| x.powi(9))
                            - 15.0 * (2.0 * sigma_sf[i] / z_grid).mapv(|x| x.powi(3))))
                    .mapv(|x| x.min(0.0));

                    bh_tail
                        + &(2.0 * PI * mi * epsilon1_k_sf[i] * sigma_sf[i].powi(3) * rho_s / 45.0
                            * (2.0 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(9))
                                - 15.0 * (sigma_sf[i] / z_grid).mapv(|x| x.powi(3))))
                }
                #[cfg(feature = "rayon")]
                Self::FreeEnergyAveraged {
                    coordinates,
                    sigma_ss,
                    epsilon_k_ss,
                    pore_center,
                    system_size,
                    n_grid,
                    cutoff_radius,
                } => {
                    // combining rules
                    let epsilon_k_sf =
                        (fluid_parameters.epsilon_k_ff()[i] * epsilon_k_ss).map(|e| e.sqrt());
                    let sigma_sf = (fluid_parameters.sigma_ff()[i] + sigma_ss) * 0.5;

                    calculate_fea_potential(
                        z_grid,
                        mi,
                        coordinates,
                        sigma_sf,
                        epsilon_k_sf,
                        pore_center,
                        system_size,
                        n_grid,
                        temperature,
                        Geometry::Cartesian,
                        *cutoff_radius,
                    )
                }
                _ => unreachable!(),
            });
        }
        ext_pot
    }

    // Evaluate the external potential in cylindrical coordinates for a given grid and fluid parameters.
    pub fn calculate_cylindrical_potential<P: HelmholtzEnergyFunctional + FluidParameters>(
        &self,
        r_grid: &Array1<f64>,
        pore_size: f64,
        fluid_parameters: &P,
        #[cfg_attr(not(feature = "rayon"), expect(unused_variables))] temperature: f64,
    ) -> Array2<f64> {
        if let ExternalPotential::Custom(potential) = self {
            return potential.clone();
        }

        // Allocate external potential
        let m = fluid_parameters.m();
        let mut ext_pot = Array2::zeros((m.len(), r_grid.len()));

        for (i, &mi) in m.iter().enumerate() {
            ext_pot.index_axis_mut(Axis_nd(0), i).assign(&match self {
                Self::HardWall { sigma_ss } => {
                    let sigma_sf = (fluid_parameters.sigma_ff()[i] + *sigma_ss) * 0.5;
                    r_grid.mapv(|r| {
                        if r > pore_size - sigma_sf {
                            f64::INFINITY
                        } else {
                            0.0
                        }
                    })
                }
                Self::LJ93 {
                    sigma_ss,
                    epsilon_k_ss,
                    rho_s,
                } => {
                    // combining rules
                    let epsilon_k_sf =
                        (fluid_parameters.epsilon_k_ff() * *epsilon_k_ss).map(|e| e.sqrt());
                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;

                    (phi(6, &(r_grid / pore_size), sigma_sf[i] / pore_size)
                        - phi(3, &(r_grid / pore_size), sigma_sf[i] / pore_size))
                        * 2.0
                        * PI
                        * mi
                        * epsilon_k_sf[i]
                        * sigma_sf[i].powi(3)
                        * *rho_s
                }
                Self::SimpleLJ93 {
                    sigma_ss: _,
                    epsilon_k_ss: _,
                } => {
                    unimplemented!()
                }
                Self::CustomLJ93 {
                    sigma_sf: _,
                    epsilon_k_sf: _,
                } => {
                    unimplemented!()
                }
                Self::Steele {
                    sigma_ss,
                    epsilon_k_ss,
                    rho_s,
                    xi,
                } => {
                    // combining rules
                    let epsilon_k_sf =
                        (fluid_parameters.epsilon_k_ff() * *epsilon_k_ss).map(|e| e.sqrt());
                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;

                    (2.0 * PI * mi * xi.unwrap_or(1.0) * epsilon_k_sf[i])
                        * (sigma_sf[i].powi(2) * DELTA_STEELE * rho_s)
                        * (psi(6, &(r_grid / pore_size), sigma_sf[i] / pore_size)
                            - psi(3, &(r_grid / pore_size), sigma_sf[i] / pore_size)
                            - sigma_sf[i] / DELTA_STEELE
                                * phi(
                                    3,
                                    &(r_grid / (pore_size + DELTA_STEELE * 0.61)),
                                    sigma_sf[i] / (pore_size + DELTA_STEELE * 0.61),
                                ))
                }
                Self::CustomSteele {
                    sigma_sf,
                    epsilon_k_sf,
                    rho_s,
                    xi,
                } => {
                    (2.0 * PI * mi * xi.unwrap_or(1.0) * epsilon_k_sf[i])
                        * (sigma_sf[i].powi(2) * DELTA_STEELE * rho_s)
                        * (psi(6, &(r_grid / pore_size), sigma_sf[i] / pore_size)
                            - psi(3, &(r_grid / pore_size), sigma_sf[i] / pore_size)
                            - sigma_sf[i] / DELTA_STEELE
                                * phi(
                                    3,
                                    &(r_grid / (pore_size + DELTA_STEELE * 0.61)),
                                    sigma_sf[i] / (pore_size + DELTA_STEELE * 0.61),
                                ))
                }
                Self::DoubleWell {
                    sigma_ss,
                    epsilon1_k_ss,
                    epsilon2_k_ss,
                    rho_s,
                } => {
                    // combining rules
                    let epsilon1_k_sf =
                        (fluid_parameters.epsilon_k_ff() * *epsilon1_k_ss).map(|e| e.sqrt());
                    let epsilon2_k_sf =
                        (fluid_parameters.epsilon_k_ff() * *epsilon2_k_ss).map(|e| e.sqrt());
                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;

                    let bh_tail = ((phi(6, &(r_grid / pore_size), 2.0 * sigma_sf[i] / pore_size)
                        - phi(3, &(r_grid / pore_size), 2.0 * sigma_sf[i] / pore_size))
                        * 2.0
                        * PI
                        * mi
                        * epsilon2_k_sf[i]
                        * sigma_sf[i].powi(3)
                        * *rho_s)
                        .mapv(|x| x.min(0.0));

                    bh_tail
                        + &((phi(6, &(r_grid / pore_size), sigma_sf[i] / pore_size)
                            - phi(3, &(r_grid / pore_size), sigma_sf[i] / pore_size))
                            * 2.0
                            * PI
                            * mi
                            * epsilon1_k_sf[i]
                            * sigma_sf[i].powi(3)
                            * *rho_s)
                }
                #[cfg(feature = "rayon")]
                Self::FreeEnergyAveraged {
                    coordinates,
                    sigma_ss,
                    epsilon_k_ss,
                    pore_center,
                    system_size,
                    n_grid,
                    cutoff_radius,
                } => {
                    // combining rules
                    let epsilon_k_sf =
                        (fluid_parameters.epsilon_k_ff()[i] * epsilon_k_ss).map(|e| e.sqrt());
                    let sigma_sf = (fluid_parameters.sigma_ff()[i] + sigma_ss) * 0.5;

                    calculate_fea_potential(
                        r_grid,
                        mi,
                        coordinates,
                        sigma_sf,
                        epsilon_k_sf,
                        pore_center,
                        system_size,
                        n_grid,
                        temperature,
                        Geometry::Cylindrical,
                        *cutoff_radius,
                    )
                }
                _ => unreachable!(),
            });
        }
        ext_pot
    }

    // Evaluate the external potential in spherical coordinates for a given grid and fluid parameters.
    pub fn calculate_spherical_potential<P: HelmholtzEnergyFunctional + FluidParameters>(
        &self,
        r_grid: &Array1<f64>,
        pore_size: f64,
        fluid_parameters: &P,
        #[cfg_attr(not(feature = "rayon"), expect(unused_variables))] temperature: f64,
    ) -> Array2<f64> {
        if let ExternalPotential::Custom(potential) = self {
            return potential.clone();
        }

        // Allocate external potential
        let m = fluid_parameters.m();
        let mut ext_pot = Array2::zeros((m.len(), r_grid.len()));

        for (i, &mi) in m.iter().enumerate() {
            ext_pot.index_axis_mut(Axis_nd(0), i).assign(&match self {
                Self::HardWall { sigma_ss } => {
                    let sigma_sf = (fluid_parameters.sigma_ff()[i] + *sigma_ss) * 0.5;
                    r_grid.mapv(|r| {
                        if r > pore_size - sigma_sf {
                            f64::INFINITY
                        } else {
                            0.0
                        }
                    })
                }
                Self::LJ93 {
                    sigma_ss,
                    epsilon_k_ss,
                    rho_s,
                } => {
                    // combining rules
                    let epsilon_k_sf =
                        (fluid_parameters.epsilon_k_ff() * *epsilon_k_ss).map(|e| e.sqrt());
                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;

                    PI * mi
                        * epsilon_k_sf[i]
                        * rho_s
                        * (sigma_sf[i].powi(12) / 90.
                            * ((r_grid - 9.0 * pore_size)
                                / (r_grid - pore_size).mapv(|x| x.powi(9))
                                - (r_grid + 9.0 * pore_size)
                                    / (r_grid + pore_size).mapv(|x| x.powi(9)))
                            - sigma_sf[i].powi(6) / 3.
                                * ((r_grid - 3.0 * pore_size)
                                    / (r_grid - pore_size).mapv(|x| x.powi(3))
                                    - (r_grid + 3.0 * pore_size)
                                        / (r_grid + pore_size).mapv(|x| x.powi(3))))
                        / r_grid
                }
                Self::SimpleLJ93 {
                    sigma_ss: _,
                    epsilon_k_ss: _,
                } => {
                    unimplemented!()
                }
                Self::CustomLJ93 {
                    sigma_sf: _,
                    epsilon_k_sf: _,
                } => {
                    unimplemented!()
                }
                Self::Steele {
                    sigma_ss,
                    epsilon_k_ss,
                    rho_s,
                    xi,
                } => {
                    // combining rules
                    let epsilon_k_sf =
                        (fluid_parameters.epsilon_k_ff() * *epsilon_k_ss).map(|e| e.sqrt());
                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;

                    (2.0 * PI * mi * xi.unwrap_or(1.0) * epsilon_k_sf[i])
                        * (sigma_sf[i].powi(2) * DELTA_STEELE * rho_s)
                        * (2.0 / 5.0 * sum_n(10, r_grid, sigma_sf[i], pore_size)
                            - sum_n(4, r_grid, sigma_sf[i], pore_size)
                            - sigma_sf[i] / (3.0 * DELTA_STEELE)
                                * (sigma_sf[i].powi(3)
                                    / r_grid
                                        .mapv(|r| (pore_size + 0.61 * DELTA_STEELE - r).powi(3))
                                    + sigma_sf[i].powi(3)
                                        / r_grid.mapv(|r| {
                                            (pore_size + 0.61 * DELTA_STEELE + r).powi(3)
                                        })
                                    + 1.5
                                        * sum_n(
                                            3,
                                            r_grid,
                                            sigma_sf[i],
                                            pore_size + 0.61 * DELTA_STEELE,
                                        )))
                }
                Self::CustomSteele {
                    sigma_sf,
                    epsilon_k_sf,
                    rho_s,
                    xi,
                } => {
                    (2.0 * PI * mi * xi.unwrap_or(1.0) * epsilon_k_sf[i])
                        * (sigma_sf[i].powi(2) * DELTA_STEELE * rho_s)
                        * (2.0 / 5.0 * sum_n(10, r_grid, sigma_sf[i], pore_size)
                            - sum_n(4, r_grid, sigma_sf[i], pore_size)
                            - sigma_sf[i] / (3.0 * DELTA_STEELE)
                                * (sigma_sf[i].powi(3)
                                    / r_grid
                                        .mapv(|r| (pore_size + 0.61 * DELTA_STEELE - r).powi(3))
                                    + sigma_sf[i].powi(3)
                                        / r_grid.mapv(|r| {
                                            (pore_size + 0.61 * DELTA_STEELE + r).powi(3)
                                        })
                                    + 1.5
                                        * sum_n(
                                            3,
                                            r_grid,
                                            sigma_sf[i],
                                            pore_size + 0.61 * DELTA_STEELE,
                                        )))
                }
                Self::DoubleWell {
                    sigma_ss,
                    epsilon1_k_ss,
                    epsilon2_k_ss,
                    rho_s,
                } => {
                    // combining rules
                    let epsilon1_k_sf =
                        (fluid_parameters.epsilon_k_ff() * *epsilon1_k_ss).map(|e| e.sqrt());
                    let epsilon2_k_sf =
                        (fluid_parameters.epsilon_k_ff() * *epsilon2_k_ss).map(|e| e.sqrt());
                    let sigma_sf = fluid_parameters.sigma_ff().add_scalar(*sigma_ss) * 0.5;

                    let bh_tail = (2.0
                        * PI
                        * mi
                        * epsilon2_k_sf[i]
                        * sigma_sf[i].powi(2)
                        * rho_s
                        * (2.0 / 5.0 * sum_n(10, r_grid, 2.0 * sigma_sf[i], pore_size)
                            - sum_n(4, r_grid, 2.0 * sigma_sf[i], pore_size)))
                    .mapv(|x| x.min(0.0));

                    bh_tail
                        + &(2.0
                            * PI
                            * mi
                            * epsilon1_k_sf[i]
                            * sigma_sf[i].powi(2)
                            * rho_s
                            * (2.0 / 5.0 * sum_n(10, r_grid, sigma_sf[i], pore_size)
                                - sum_n(4, r_grid, sigma_sf[i], pore_size)))
                }
                #[cfg(feature = "rayon")]
                Self::FreeEnergyAveraged {
                    coordinates,
                    sigma_ss,
                    epsilon_k_ss,
                    pore_center,
                    system_size,
                    n_grid,
                    cutoff_radius,
                } => {
                    // combining rules
                    let epsilon_k_sf =
                        (fluid_parameters.epsilon_k_ff()[i] * epsilon_k_ss).map(|e| e.sqrt());
                    let sigma_sf = (fluid_parameters.sigma_ff()[i] + sigma_ss) * 0.5;

                    calculate_fea_potential(
                        r_grid,
                        mi,
                        coordinates,
                        sigma_sf,
                        epsilon_k_sf,
                        pore_center,
                        system_size,
                        n_grid,
                        temperature,
                        Geometry::Spherical,
                        *cutoff_radius,
                    )
                }
                _ => unreachable!(),
            });
        }
        ext_pot
    }
}

fn phi(n: i32, r_r: &Array1<f64>, sigma_r: f64) -> Array1<f64> {
    let m3n2 = 3.0 - 2.0 * n as f64;
    let n2m3 = 2.0 * n as f64 - 3.0;

    (1.0 - &r_r.mapv(|r| r.powi(2))).mapv(|r| r.powf(m3n2)) * 4.0 * PI.sqrt() / n2m3
        * sigma_r.powf(n2m3)
        * tgamma(n as f64 - 0.5)
        / tgamma(n as f64)
        * taylor_2f1_phi(r_r, n)
}

fn psi(n: i32, r_r: &Array1<f64>, sigma_r: f64) -> Array1<f64> {
    (1.0 - &r_r.mapv(|r| r.powi(2))).mapv(|r| r.powf(2.0 - 2.0 * n as f64))
        * 4.0
        * PI.sqrt()
        * tgamma(n as f64 - 0.5)
        / tgamma(n as f64)
        * sigma_r.powf(2.0 * n as f64 - 2.0)
        * taylor_2f1_psi(r_r, n)
}

fn sum_n(n: i32, r_grid: &Array1<f64>, sigma: f64, pore_size: f64) -> Array1<f64> {
    let mut result = Array1::zeros(r_grid.len());
    for i in 0..n {
        result = result
            + sigma.powi(n) / (pore_size.powi(i) * r_grid.mapv(|r| (pore_size - r).powi(n - i)))
            + sigma.powi(n) / (pore_size.powi(i) * r_grid.mapv(|r| (pore_size + r).powi(n - i)));
    }
    result
}

fn taylor_2f1_phi(x: &Array1<f64>, n: i32) -> Array1<f64> {
    match n {
        3 => {
            1.0 + 3.0 / 4.0 * x.mapv(|x| x.powi(2)) + 3.0 / 64.0 * x.mapv(|x| x.powi(4))
                - 1.0 / 256.0 * x.mapv(|x| x.powi(6))
                - 15.0 / 16384.0 * x.mapv(|x| x.powi(8))
        }
        6 => {
            1.0 + 63.0 / 4.0 * x.mapv(|x| x.powi(2))
                + 2205.0 / 64.0 * x.mapv(|x| x.powi(4))
                + 3675.0 / 256.0 * x.mapv(|x| x.powi(6))
                + 11025.0 / 16384.0 * x.mapv(|x| x.powi(8))
        }
        _ => unreachable!(),
    }
}

fn taylor_2f1_psi(x: &Array1<f64>, n: i32) -> Array1<f64> {
    match n {
        3 => {
            1.0 + 9.0 / 4.0 * x.mapv(|x| x.powi(2))
                + 9.0 / 64.0 * x.mapv(|x| x.powi(4))
                + 1.0 / 256.0 * x.mapv(|x| x.powi(6))
                + 9.0 / 16384.0 * x.mapv(|x| x.powi(8))
        }
        6 => {
            1.0 + 81.0 / 4.0 * x.mapv(|x| x.powi(2))
                + 3969.0 / 64.0 * x.mapv(|x| x.powi(4))
                + 11025.0 / 256.0 * x.mapv(|x| x.powi(6))
                + 99125.0 / 16384.0 * x.mapv(|x| x.powi(8))
        }
        _ => unreachable!(),
    }
}
