#![allow(clippy::excessive_precision)]
#![cfg(feature = "dft")]
use approx::assert_relative_eq;
use feos::gc_pcsaft::{GcPcSaft, GcPcSaftFunctional, GcPcSaftParameters};
use feos_core::parameter::{ChemicalRecord, Identifier, IdentifierOption, SegmentRecord};
use feos_core::{PhaseEquilibrium, State, StateBuilder, Verbosity};
use feos_dft::adsorption::{ExternalPotential, Pore1D, PoreSpecification};
use feos_dft::interface::PlanarInterface;
use feos_dft::{DFTSolver, Geometry};
use nalgebra::dvector;
use quantity::*;
use std::error::Error;
use typenum::P3;

#[test]
#[allow(non_snake_case)]
fn test_bulk_implementation() -> Result<(), Box<dyn Error>> {
    // correct for different k_B in old code
    let KB_old = 1.38064852e-23 * JOULE / KELVIN;
    let NAV_old = 6.022140857e23 / MOL;

    let parameters = GcPcSaftParameters::from_json_segments_hetero(
        &["propane"],
        "../../parameters/pcsaft/gc_substances.json",
        "../../parameters/pcsaft/sauer2014_hetero.json",
        None,
        IdentifierOption::Name,
    )
    .unwrap();

    let parameters_func = GcPcSaftParameters::from_json_segments_hetero(
        &["propane"],
        "../../parameters/pcsaft/gc_substances.json",
        "../../parameters/pcsaft/sauer2014_hetero.json",
        None,
        IdentifierOption::Name,
    )
    .unwrap();

    let eos = GcPcSaft::new(parameters);
    let func = GcPcSaftFunctional::new(parameters_func);
    let t = 200.0 * KELVIN;
    let v = 0.002 * METER.powi::<P3>() * NAV / NAV_old;
    let n = dvector![1.5] * MOL;
    let state_eos = State::new_nvt(&&eos, t, v, &n)?;
    let state_func = State::new_nvt(&&func, t, v, &n)?;
    let p_eos = state_eos.pressure_contributions();
    let p_func = state_func.pressure_contributions();

    println!("{:29}: {}", p_eos[0].0, p_eos[0].1);
    println!("{:29}: {}", p_func[0].0, p_func[0].1);
    println!();
    println!("{:29}: {}", p_eos[1].0, p_eos[1].1);
    println!("{:29}: {}", p_func[1].0, p_func[1].1);
    println!();
    println!("{:29}: {}", p_eos[2].0, p_eos[2].1);
    println!("{:29}: {}", p_func[2].0, p_func[2].1);
    println!();
    println!("{:29}: {}", p_eos[3].0, p_eos[3].1);
    println!("{:29}: {}", p_func[3].0, p_func[3].1);

    assert_relative_eq!(
        p_eos[0].1,
        1.2471689792172869 * MEGA * PASCAL * KB / KB_old,
        max_relative = 1e-14,
    );
    assert_relative_eq!(
        p_eos[1].1,
        280.0635060891395938 * KILO * PASCAL * KB / KB_old,
        max_relative = 1e-14,
    );
    assert_relative_eq!(
        p_eos[2].1,
        -141.9023918353318550 * KILO * PASCAL * KB / KB_old,
        max_relative = 1e-14,
    );
    assert_relative_eq!(
        p_eos[3].1,
        -763.2289230004602132 * KILO * PASCAL * KB / KB_old,
        max_relative = 1e-14,
    );

    assert_relative_eq!(
        p_func[0].1,
        1.2471689792172869 * MEGA * PASCAL * KB / KB_old,
        max_relative = 1e-14,
    );
    assert_relative_eq!(
        p_func[1].1,
        280.0635060891395938 * KILO * PASCAL * KB / KB_old,
        max_relative = 1e-14,
    );
    assert_relative_eq!(
        p_func[2].1,
        -141.9023918353318550 * KILO * PASCAL * KB / KB_old,
        max_relative = 1e-14,
    );
    assert_relative_eq!(
        p_func[3].1,
        -763.2289230004602132 * KILO * PASCAL * KB / KB_old,
        max_relative = 1e-14,
    );
    Ok(())
}

#[test]
fn test_bulk_association() -> Result<(), Box<dyn Error>> {
    let segment_records =
        SegmentRecord::from_json("../../parameters/pcsaft/sauer2014_hetero.json")?;
    let ethylene_glycol = ChemicalRecord::new(
        Identifier::default(),
        vec!["OH".into(), "CH2".into(), "CH2".into(), "OH".into()],
        None,
    );
    let eos_parameters = GcPcSaftParameters::from_segments_hetero(
        vec![ethylene_glycol.clone()],
        &segment_records,
        None,
    )?;
    let eos = GcPcSaft::new(eos_parameters);
    let func_parameters =
        GcPcSaftParameters::from_segments_hetero(vec![ethylene_glycol], &segment_records, None)?;
    let func = GcPcSaftFunctional::new(func_parameters);

    let t = 200.0 * KELVIN;
    let v = 0.002 * METER.powi::<P3>();
    let n = dvector![1.5] * MOL;
    let state_eos = State::new_nvt(&&eos, t, v, &n)?;
    let state_func = State::new_nvt(&&func, t, v, &n)?;
    let p_eos = state_eos.pressure_contributions();
    let p_func = state_func.pressure_contributions();
    for (s, x) in &p_eos {
        println!("{s:18}: {x:21.16}");
    }
    for (s, x) in &p_func {
        println!("{s:26}: {x:21.16}");
    }
    assert_relative_eq!(p_eos[4].1, p_func[4].1, max_relative = 1e-14);
    Ok(())
}

#[test]
#[allow(non_snake_case)]
fn test_dft() -> Result<(), Box<dyn Error>> {
    // correct for different k_B in old code
    let KB_old = 1.38064852e-23 * JOULE / KELVIN;
    let NAV_old = 6.022140857e23 / MOL;

    let parameters = GcPcSaftParameters::from_json_segments_hetero(
        &["propane"],
        "../../parameters/pcsaft/gc_substances.json",
        "../../parameters/pcsaft/sauer2014_hetero.json",
        None,
        IdentifierOption::Name,
    )
    .unwrap();

    let func = GcPcSaftFunctional::new(parameters);
    let t = 200.0 * KELVIN;
    let w = 150.0 * ANGSTROM;
    let points = 2048;
    let tc = State::critical_point(&&func, None, None, None, Default::default())?.temperature;
    let vle = PhaseEquilibrium::pure(&&func, t, None, Default::default())?;
    let profile = PlanarInterface::from_tanh(&vle, points, w, tc, false).solve(None)?;
    println!(
        "hetero {} {} {}",
        profile.surface_tension.unwrap(),
        vle.vapor().density,
        vle.liquid().density,
    );

    assert_relative_eq!(
        vle.vapor().density,
        12.8820179191167643 * MOL / METER.powi::<P3>() * NAV_old / NAV,
        max_relative = 1e-13,
    );

    assert_relative_eq!(
        vle.liquid().density,
        13.2705903446123212 * KILO * MOL / METER.powi::<P3>() * NAV_old / NAV,
        max_relative = 1e-13,
    );

    assert_relative_eq!(
        profile.surface_tension.unwrap(),
        21.1478901897016272 * MILLI * NEWTON / METER * KB / KB_old,
        max_relative = 6e-5,
    );
    Ok(())
}

#[test]
#[allow(non_snake_case)]
fn test_dft_assoc() -> Result<(), Box<dyn Error>> {
    let parameters = GcPcSaftParameters::from_json_segments_hetero(
        &["1-pentanol"],
        "../../parameters/pcsaft/gc_substances.json",
        "../../parameters/pcsaft/sauer2014_hetero.json",
        None,
        IdentifierOption::Name,
    )
    .unwrap();

    let func = &GcPcSaftFunctional::new(parameters);
    let t = 300.0 * KELVIN;
    let w = 100.0 * ANGSTROM;
    let points = 4096;
    let vle = PhaseEquilibrium::pure(&func, t, None, Default::default())?;
    let profile = PlanarInterface::from_tanh(&vle, points, w, 600.0 * KELVIN, false).solve(None)?;
    println!(
        "hetero {} {} {}",
        profile.surface_tension.unwrap(),
        vle.vapor().density,
        vle.liquid().density,
    );

    let solver = DFTSolver::new(Some(Verbosity::Iter))
        .picard_iteration(None, None, Some(1e-5), Some(0.05))
        .anderson_mixing(None, None, None, None, None);
    let bulk = StateBuilder::new(&func)
        .temperature(t)
        .pressure(5.0 * BAR)
        .build()?;
    Pore1D::new(
        Geometry::Cartesian,
        20.0 * ANGSTROM,
        ExternalPotential::LJ93 {
            epsilon_k_ss: 10.0,
            sigma_ss: 3.0,
            rho_s: 0.08,
        },
        None,
        None,
    )
    .initialize(&bulk, None, None)
    .unwrap()
    .solve(Some(&solver))?;
    Ok(())
}

#[test]
#[allow(non_snake_case)]
fn test_dft_newton() -> Result<(), Box<dyn Error>> {
    let params = GcPcSaftParameters::from_json_segments_hetero(
        &["propane"],
        "../../parameters/pcsaft/gc_substances.json",
        "../../parameters/pcsaft/sauer2014_hetero.json",
        None,
        IdentifierOption::Name,
    )
    .unwrap();
    let func = GcPcSaftFunctional::new(params);
    let t = 200.0 * KELVIN;
    let w = 150.0 * ANGSTROM;
    let points = 512;
    let tc = State::critical_point(&&func, None, None, None, Default::default())?.temperature;
    let vle = PhaseEquilibrium::pure(&&func, t, None, Default::default())?;
    let solver = DFTSolver::new(Some(Verbosity::Iter))
        .picard_iteration(None, Some(10), None, None)
        .newton(None, None, None, None);
    PlanarInterface::from_tanh(&vle, points, w, tc, false).solve(Some(&solver))?;
    Ok(())
}
