te1d/teg/

mod.rs

//! Simulates steady-state thermoelectric power generation.

use std::fmt;
use std::rc::Rc;

use ndarray::{Array, Array2, array, s};

use crate::calc::{ColVec, integrate_simpson, is_close};
use crate::func::{RealValuedFunction};
use crate::ode::{explicit_rk4, explicit_rk4_fvp};
use crate::opt::{OptimizeResult, MinimizerKind, MinimizerParam, minimize, minimize_brent};
use crate::tep::{Tep};
use crate::bc::{BoundaryCondition, OptVarConverter, FullOptVar};
use crate::simul::{SimulationErrorKind, convert_opt_error_to_simul_error};
use crate::interp::{PiecewiseCubicPolynomialConstExt};

pub mod prelude;


/// Parameters for energy generation module simulation
/// 
/// # Tips
/// * Setting `dx: length/22.0` and `grad_tol: 1e-07` is a good choice for high accuracy.
/// * Setting `dx: length/10.0` and `grad_tol: 1e-06` is a good choice for fast computation.
pub struct SimulationParam {
    /// A suggested size of `x`-meshes to discretize the space in legs. Real size of `x`-meshes can be different.
    pub dx: f64,
    /// A suggested size of `T`-meshes to discretize the temperature intervals. Real size of `T`-meshes can be different.
    pub dtemp: f64,
    /// Minimum number of `x`-meshes. If `dx` is too large, `min_num_x_mesh` is used to prevent a coarse mesh.
    pub min_num_x_mesh: usize,
    /// Minimum number of `T`-meshes. If `dtemp` is too large, `min_num_temp_mesh` is used to prevent a coarse mesh.
    pub min_num_temp_mesh: usize,
    /// Maximum number of `x`-meshes. If `dx` is too small, `max_num_x_mesh` is used to prevent a excessively fine mesh.
    pub max_num_x_mesh: usize,
    /// Maximum number of `T`-meshes. If `dtemp` is too small, `max_num_temp_mesh` is used to prevent a excessively fine mesh.
    pub max_num_temp_mesh: usize,
    /// The kind of minimizer that is used for matching the boundary conditions
    pub minimizer_kind: MinimizerKind,
    /// The parameters of the minimizer
    pub minimizer_param: MinimizerParam,
    /// Relative tolerance in temperature (K) matching
    pub rel_tol_temp: f64,
    /// Absolute tolerance in temperature (K) matching
    pub abs_tol_temp: f64,
    /// Relative tolerance in heat rate (W) matching
    pub rel_tol_heat_rate: f64,
    /// Absolute tolerance in heat rate (W) matching
    pub abs_tol_heat_rate: f64,
    /// Maximum number of iterations to perform Brent-Dekker minimization
    pub max_iter_brent: usize,
    /// Absolute tolerance in electrical current (A), used in Brent-Dekker minimization algorithm
    pub abs_tol_elec_cur_brent: f64,
}

impl SimulationParam {
    /// Returns the default parameters.
    /// * **Modify** `dx` and `dtemp` according to the length of the leg and the size of temperature interval.
    pub fn default() -> Self {
        SimulationParam {
            dx: 1e-4,      // (m) (=0.1 mm)
            dtemp: 10.0,   // (K)
            min_num_x_mesh: 11,
            min_num_temp_mesh: 11,
            max_num_x_mesh: 101,
            max_num_temp_mesh: 101,
            minimizer_kind: MinimizerKind::Bfgs,
            minimizer_param: MinimizerParam::default(MinimizerKind::Bfgs),
            rel_tol_temp: 1e-4,
            abs_tol_temp: 1e-2,
            rel_tol_heat_rate: 1e-4,
            abs_tol_heat_rate: 1e-4,
            max_iter_brent: 40,
            abs_tol_elec_cur_brent: 1e-6,
        }
    }
}


/// Results of leg simulation for energy generation
pub struct LegSimulationResult {
    /// The temperature distribution inside a leg
    pub temp: Rc<dyn RealValuedFunction>,
    /// The heat flux distribution inside a leg
    pub heat_flux: Rc<dyn RealValuedFunction>,
    /// Electrical Current (A)
    pub elec_cur: f64,
    /// Voltage (V) induced from the Seebeck effect
    pub voltage: f64,
    /// Resistance (Ohm) inside a leg
    pub internal_resistance: f64,
    /// Temperature (K) at the left end
    pub templ: f64,
    /// Temperature (K) at the right end
    pub tempr: f64,
    /// Heat rate (W) at the left end
    pub heat_ratel: f64,
    /// Heat rate (W) at the right end
    pub heat_rater: f64,
    /// Thermoelectric power (W)
    pub power: f64,
    /// Energy conversion efficiency (1)
    pub efficiency: f64,
    /// A boolean flag indicating if the simulation exited successfully
    pub success: bool,
    /// Error from the simulation, if nay
    pub error: Option<SimulationErrorKind>,
    /// Number of iterations in optimization process
    pub num_iter: usize,
}

impl fmt::Display for LegSimulationResult {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "LegSimulationResult {{ success: {}, elec_cur: {}, voltage: {}, internal_resistance: {}, templ: {}, tempr: {}, heat_ratel: {}, heat_rater: {}, power: {}, efficiency: {}, .. }}",
            self.success, self.elec_cur, self.voltage, self.internal_resistance, self.templ, self.tempr, self.heat_ratel, self.heat_rater, self.power, self.efficiency)
    }
}

impl LegSimulationResult {
    /// Returns the sum of internal and external load resistances (Ohm)
    #[inline]
    pub fn total_resistance(&self) -> f64 {
        self.voltage / self.elec_cur
    }

    /// Returns the external load resistance (Ohm); if the `elec_cur==0`, then `NAN` is returned.
    #[inline]
    pub fn load_resistance(&self) -> f64 {
        self.voltage / self.elec_cur - self.internal_resistance
    }

    // Returns the external-to-internal reistance ratio (`=load_resistance/internal_reistance`)
    #[inline]
    pub fn load_ratio(&self) -> f64 {
        self.voltage / (self.elec_cur * self.internal_resistance) - 1.0
    }
}


/// Simulator of single-leg thermoelectric generators
/// 
/// # Examples
/// The following is a simple example simulating a constant property model under a fixed electrical current.
/// ```
/// use std::rc::Rc;
/// use ndarray::prelude::*;
/// use te1d::teg::prelude::*;
/// use te1d::func::ConstantFunction;
///
/// // set tep and leg parameters
/// let alpha0: f64 = 1.2;
/// let rho0: f64 = 3.0;
/// let kappa0: f64 = 1.5;
/// let tep = Tep::new(Rc::new(ConstantFunction::new(alpha0)),
///                    Rc::new(ConstantFunction::new(rho0)),
///                    Rc::new(ConstantFunction::new(kappa0)));
/// let length = 0.8;
/// let area = 0.3;
/// 
/// // set simulation parameters
/// let minimizer_kind = MinimizerKind::Bfgs;
/// let minimizer_param = MinimizerParam {
///     grad_tol: 1e-05,
///     a_max: 2.0,
///     max_iter: 200,
///     ..MinimizerParam::default(minimizer_kind)
/// };
/// let simul_param = SimulationParam {
///     dx: length/10.0,
///     minimizer_kind,
///     minimizer_param,
///     ..SimulationParam::default()
/// };
/// 
/// // create a leg
/// let leg = SingleLeg::new(&tep, length, area, simul_param);
/// let elec_cur: f64 = 1.3;
/// let j: f64 = elec_cur / area;
/// let hot_temp: f64 = 600.0;
/// let cold_temp: f64 = 300.0;
/// let left_bc = FixedTempBc::new(hot_temp);
/// let right_bc = FixedTempBc::new(cold_temp);
/// 
/// // simulate
/// let simul_result = leg.simulate_fixed_elec_cur(elec_cur, &left_bc, &right_bc, hot_temp, cold_temp);
/// assert!(simul_result.success);
/// 
/// // check the solution
/// let x_mesh: ColVec = Array::linspace(0.0, length, 31);
/// let temp_numerical: ColVec = simul_result.temp.ats_incr(&x_mesh);
/// let heat_flux_numerical: ColVec = simul_result.heat_flux.ats_incr(&x_mesh);
/// let temp_exact: ColVec = hot_temp + (cold_temp-hot_temp)*&x_mesh/length + j.powi(2)*0.5*rho0/kappa0*&x_mesh*(length-&x_mesh);
/// let dtempdx_exact: ColVec = (cold_temp-hot_temp)/length + j.powi(2)*0.5*rho0/kappa0*(length-2.0*&x_mesh);
/// let heat_flux_exact: ColVec = alpha0*&temp_exact*j - kappa0*&dtempdx_exact;
/// 
/// assert!(all_close(&temp_numerical, &temp_exact, 1e-05, 0.0));
/// assert!(all_close(&heat_flux_numerical, &heat_flux_exact, 1e-05, 0.0));
/// ```
pub struct SingleLeg {
    /// Thermoelectric material properties
    pub tep: Tep,
    /// Length (m) of a leg
    pub length: f64,
    /// Cross-sectional area (m^2) of a leg
    pub area: f64,
    /// Simulation parameters
    pub simulation_param: SimulationParam,
    /// The mesh points discretizing the sptial coordinate inside the leg. Must be equidistant.
    pub xs: ColVec,
    /// The left-half mesh points.
    xs_left: ColVec,
    /// The right-half mesh points.
    xs_right: ColVec,
}

impl SingleLeg {
    /// Create a new simulator
    /// 
    /// # Panics
    /// * when `length` is not positive.
    /// * when `area` is not positive.
    pub fn new(tep: &Tep, length: f64, area: f64, simulation_param: SimulationParam) -> Self {
        if length <= 0.0 { panic!("`length` must be a positive number!"); }
        if area <= 0.0 { panic!("`area` must be a positive number!"); }

        // generate mesh points
        let num_x_mesh: usize = get_num_mesh(length, simulation_param.dx, simulation_param.min_num_x_mesh, simulation_param.max_num_x_mesh);
        let xs: ColVec = Array::linspace(0.0, length, num_x_mesh);
        let half_num_x_mesh: usize = (num_x_mesh+1)/2;  // `num_x_mesh` must be odd
        let xs_left: ColVec = xs.slice(s![0..half_num_x_mesh]).to_owned();
        let xs_right: ColVec = xs.slice(s![half_num_x_mesh-1..num_x_mesh]).to_owned();

        // generate a thermoelectric leg
        SingleLeg { tep: tep.clone(), length, area, simulation_param, xs, xs_left, xs_right }
    }

    /// Find the steady state under the given boundary condition and electrical current (A).
    /// * You need to provide hints for electrical current (A), left-end temperature (K) and right-end temperature (K).
    /// * If this does not converge, use [`SingleLeg::simulate_fixed_elec_cur_with_manual_init()`] to provide initial values manually.
    pub fn simulate_fixed_elec_cur(&self, fixed_elec_cur: f64, left_bc: &(dyn BoundaryCondition), right_bc: &(dyn BoundaryCondition), templ_hint: f64, tempr_hint: f64) -> LegSimulationResult {
        let fixed_elec_cur_density: f64 = fixed_elec_cur / self.area;
        let heat_rate_hint: f64 = self.estimate_heat_flux(templ_hint, tempr_hint, fixed_elec_cur_density) * self.area;
        let init_var: FullOptVar = FullOptVar { templ: templ_hint, heat_ratel: heat_rate_hint, tempr: tempr_hint, heat_rater: heat_rate_hint, elec_cur: fixed_elec_cur };

        self.simulate_fixed_elec_cur_with_manual_init(fixed_elec_cur, left_bc, right_bc, &init_var)
    }

    /// Find the steady state under the given boundary condition and electrical current (A).
    /// * You need to provide initial values for optimization by completing [`FullOptVar`](crate::bc::FullOptVar)
    pub fn simulate_fixed_elec_cur_with_manual_init(&self, fixed_elec_cur: f64, left_bc: &(dyn BoundaryCondition), right_bc: &(dyn BoundaryCondition), init_var: &FullOptVar ) -> LegSimulationResult {
        let converter: OptVarConverter = OptVarConverter::new(left_bc, right_bc, Some(fixed_elec_cur), self.area);
        let n_half_cols: usize = self.xs_left.len();

        let cost_f = |var: &ColVec| -> f64 {
            let full_var: FullOptVar = converter.to_full_var(left_bc, right_bc, var);
            let elec_cur_density: f64 = full_var.elec_cur / self.area;
            let ivp_ys = self.solve_ivp(elec_cur_density, full_var.templ, full_var.heat_ratel/self.area, &self.xs_left);
            let fvp_ys = self.solve_fvp(elec_cur_density, full_var.tempr, full_var.heat_rater/self.area, &self.xs_right);

            let ivp_temp_mid = ivp_ys[[0, n_half_cols-1]];
            let ivp_heat_rate_mid = ivp_ys[[1, n_half_cols-1]] * self.area;
            let fvp_temp_mid = fvp_ys[[0, 0]];
            let fvp_heat_rate_mid = fvp_ys[[1, 0]] * self.area;

            let mut cost = (ivp_temp_mid - fvp_temp_mid).powi(2) + (ivp_heat_rate_mid - fvp_heat_rate_mid).powi(2);
            cost *= 0.5;

            cost
        };

        self.minimize_cost_by_shooting(&cost_f, left_bc, right_bc, &converter, init_var)
    }

    /// Find the steady state under the given boundary condition and external resistance (Ohm).
    /// * You need to provide hints for electrical current (A), left-end temperature (K) and right-end temperature (K).
    /// * If this does not converge, use [`SingleLeg::simulate_fixed_load_resistance_with_manual_init()`] to provide initial values manually.
    pub fn simulate_fixed_load_resistance(&self, fixed_load_resistance: f64, left_bc: &(dyn BoundaryCondition), right_bc: &(dyn BoundaryCondition), elec_cur_hint: f64, templ_hint: f64, tempr_hint: f64) -> LegSimulationResult {
        let elec_cur_density_hint: f64 = elec_cur_hint / self.area;
        let heat_rate_hint: f64 = self.estimate_heat_flux(templ_hint, tempr_hint, elec_cur_density_hint) * self.area;
        let init_var: FullOptVar = FullOptVar { templ: templ_hint, heat_ratel: heat_rate_hint, tempr: tempr_hint, heat_rater: heat_rate_hint, elec_cur: elec_cur_hint };

        self.simulate_fixed_load_resistance_with_manual_init(fixed_load_resistance, left_bc, right_bc, &init_var)
    }

    /// Find the steady state under the given boundary condition and external resistance (Ohm).
    /// * You need to provide initial values for optimization by completing [`FullOptVar`](crate::bc::FullOptVar)
    pub fn simulate_fixed_load_resistance_with_manual_init(&self, fixed_load_resistance: f64, left_bc: &(dyn BoundaryCondition), right_bc: &(dyn BoundaryCondition), init_var: &FullOptVar) -> LegSimulationResult {
        let converter: OptVarConverter = OptVarConverter::new(left_bc, right_bc, None, self.area);
        let n_half_cols: usize = self.xs_left.len();
        let n_cols: usize = self.xs.len();

        let cost_f = |var: &ColVec| -> f64 {
            let full_var = converter.to_full_var(left_bc, right_bc, var);    
            let elec_cur_density: f64 = full_var.elec_cur / self.area;
            let ivp_ys = self.solve_ivp(elec_cur_density, full_var.templ, full_var.heat_ratel / self.area, &self.xs);
            let fvp_ys = self.solve_fvp(elec_cur_density, full_var.tempr, full_var.heat_rater / self.area, &self.xs);
            let ivp_temp_mid = ivp_ys[[0, n_half_cols-1]];
            let ivp_heat_rate_mid = ivp_ys[[1, n_half_cols-1]] * self.area;
            let fvp_temp_mid = fvp_ys[[0, n_half_cols-1]];
            let fvp_heat_rate_mid = fvp_ys[[1, n_half_cols-1]] * self.area;
            let ys = (1.0-&self.xs/self.length) * &ivp_ys + &self.xs/self.length * &fvp_ys;  // this enforces the correct BCs

            let voltage = self.compute_voltage(ys[[0, 0]], ys[[0, n_cols-1]]);
            let internal_resistance = self.compute_internal_resistance(&ys.row(0).to_owned());
            let elec_cur = full_var.elec_cur;

            let mut cost = (ivp_temp_mid - fvp_temp_mid).powi(2) + (ivp_heat_rate_mid - fvp_heat_rate_mid).powi(2);
            cost *= 0.5;
            cost += (voltage - elec_cur*(internal_resistance + fixed_load_resistance)).powi(2) * 1e3;

            cost
        };

        self.minimize_cost_by_shooting(&cost_f, left_bc, right_bc, &converter, init_var)
    }

    /// Find the steady state under the given boundary condition and external-to-internal load ratio (1).
    /// * You need to provide hints for electrical current (A), left-end temperature (K) and right-end temperature (K).
    /// * If this does not converge, use [`SingleLeg::simulate_fixed_load_ratio_with_manual_init()`] to provide initial values manually.
    pub fn simulate_fixed_load_ratio(&self, fixed_load_ratio: f64, left_bc: &(dyn BoundaryCondition), right_bc: &(dyn BoundaryCondition), elec_cur_hint: f64, templ_hint: f64, tempr_hint: f64) -> LegSimulationResult {
        let elec_cur_density_hint: f64 = elec_cur_hint / self.area;
        let heat_rate_hint: f64 = self.estimate_heat_flux(templ_hint, tempr_hint, elec_cur_density_hint) * self.area;
        let init_var: FullOptVar = FullOptVar { templ: templ_hint, heat_ratel: heat_rate_hint, tempr: tempr_hint, heat_rater: heat_rate_hint, elec_cur: elec_cur_hint };

        self.simulate_fixed_load_ratio_with_manual_init(fixed_load_ratio, left_bc, right_bc, &init_var)
    }

    /// Find the steady state under the given boundary condition and external-to-internal load ratio (1).
    /// * You need to provide initial values for optimization by completing [`FullOptVar`](crate::bc::FullOptVar)
    /// 
    /// # Panics
    /// * when `fixed_load_ratio` is infinite
    /// * when `fixed_load_ratio` is negative
    pub fn simulate_fixed_load_ratio_with_manual_init(&self, fixed_load_ratio: f64, left_bc: &(dyn BoundaryCondition), right_bc: &(dyn BoundaryCondition), init_var: &FullOptVar) -> LegSimulationResult {
        if fixed_load_ratio.is_infinite() {
            panic!("`fixed_load_ratio` must be a finite value!");
        }
        if fixed_load_ratio < 0.0 {
            panic!("`fixed_load_ratio` must be a nonngeative value!");
        }

        let converter: OptVarConverter = OptVarConverter::new(left_bc, right_bc, None, self.area);
        let n_half_cols: usize = self.xs_left.len();
        let n_cols: usize = self.xs.len();

        let cost_f = |var: &ColVec| -> f64 {
            let full_var = converter.to_full_var(left_bc, right_bc, var);    
            let elec_cur_density: f64 = full_var.elec_cur / self.area;
            let ivp_ys = self.solve_ivp(elec_cur_density, full_var.templ, full_var.heat_ratel / self.area, &self.xs);
            let fvp_ys = self.solve_fvp(elec_cur_density, full_var.tempr, full_var.heat_rater / self.area, &self.xs);
            let ivp_temp_mid = ivp_ys[[0, n_half_cols-1]];
            let ivp_heat_rate_mid = ivp_ys[[1, n_half_cols-1]] * self.area;
            let fvp_temp_mid = fvp_ys[[0, n_half_cols-1]];
            let fvp_heat_rate_mid = fvp_ys[[1, n_half_cols-1]] * self.area;
            let ys = (1.0-&self.xs/self.length) * &ivp_ys + &self.xs/self.length * &fvp_ys;  // this enforces the correct BCs

            let voltage = self.compute_voltage(ys[[0, 0]], ys[[0, n_cols-1]]);
            let internal_resistance = self.compute_internal_resistance(&ys.row(0).to_owned());
            let elec_cur = full_var.elec_cur;

            let mut cost = (ivp_temp_mid - fvp_temp_mid).powi(2) + (ivp_heat_rate_mid - fvp_heat_rate_mid).powi(2);
            cost *= 0.5;
            cost += (voltage - elec_cur*(1.0 + fixed_load_ratio)*internal_resistance).powi(2) * 1e3;

            //println!("var={}, cost={}", var, cost);

            cost
        };

        self.minimize_cost_by_shooting(&cost_f, left_bc, right_bc, &converter, init_var)
    }

    /// Find the steady state minimizing the cost function. A shooting on electrical current, boundary temperature and boundary heat rate is performed.
    /// * You need to provide initial values for optimization by completing [`FullOptVar`](crate::bc::FullOptVar)
    fn minimize_cost_by_shooting(&self, cost_f: &(dyn Fn(&ColVec) -> f64), left_bc: &(dyn BoundaryCondition), right_bc: &(dyn BoundaryCondition), converter: &OptVarConverter, init_var: &FullOptVar) -> LegSimulationResult {
        let var0: ColVec = converter.to_opt_var(init_var);
        let n_half_cols: usize = self.xs_left.len();
        let optimize_result: OptimizeResult = minimize(&cost_f, &var0, self.simulation_param.minimizer_kind, &self.simulation_param.minimizer_param);
        let full_var: FullOptVar = converter.to_full_var(left_bc, right_bc, &optimize_result.x);

        // Report the result
        let mut success: bool = true;
        let mut error: Option<SimulationErrorKind> = None;
        // generate the solution by taking the average of two-side shooting solutions.
        let elec_cur_density: f64 = full_var.elec_cur / self.area;
        let ivp_ys = self.solve_ivp(elec_cur_density, full_var.templ, full_var.heat_ratel / self.area, &self.xs);
        let fvp_ys = self.solve_fvp(elec_cur_density, full_var.tempr, full_var.heat_rater / self.area, &self.xs);
        let ys = (1.0-&self.xs/self.length) * &ivp_ys + &self.xs/self.length * &fvp_ys;  // this enforces the correct BCs
        // if the mid-point values are different, it is a failure
        if !(is_close(ivp_ys[[0, n_half_cols]], fvp_ys[[0, n_half_cols]], self.simulation_param.rel_tol_temp, self.simulation_param.abs_tol_temp))
           || !(is_close(ivp_ys[[1, n_half_cols]]/self.area, fvp_ys[[1, n_half_cols]]/self.area, self.simulation_param.rel_tol_heat_rate, self.simulation_param.abs_tol_heat_rate)) {
            success = false;
            error = Some(SimulationErrorKind::BcNotMatchError);
        }
        // if the optimization is incomplete, it is a failure
        if !(optimize_result.success) {
            success = false;
            error = Some(convert_opt_error_to_simul_error(optimize_result.error.unwrap()));
        }

        let temps: ColVec = ys.row(0).to_owned();
        let heat_fluxs: ColVec = ys.row(1).to_owned();
        let heat_fluxl: f64 = full_var.heat_ratel/self.area;
        let heat_fluxr: f64 = full_var.heat_rater/self.area;

        LegSimulationResult {
            temp: Rc::new(PiecewiseCubicPolynomialConstExt::new(&self.xs, &temps, full_var.templ, full_var.tempr)),
            heat_flux: Rc::new(PiecewiseCubicPolynomialConstExt::new(&self.xs, &heat_fluxs, heat_fluxl, heat_fluxr)),
            elec_cur: full_var.elec_cur,
            voltage: self.compute_voltage(full_var.templ, full_var.tempr),
            internal_resistance: self.compute_internal_resistance(&temps),
            templ: full_var.templ,
            tempr: full_var.tempr,
            heat_ratel: full_var.heat_ratel,
            heat_rater: full_var.heat_rater,
            power: self.get_power(heat_fluxl, heat_fluxr),
            efficiency: get_efficiency(heat_fluxl, heat_fluxr),
            success,
            error,
            num_iter: optimize_result.num_iter,
        }
    }

    /// Find the steady state under the given boundary condition and external-to-internal load ratio (1).
    /// * You need to provide hints search range for electrical current (A), left-end temperature (K) and right-end temperature (K).
    /// * If this does not converge, use [`SingleLeg::simulate_fixed_load_ratio_in_range_with_manual_init()`] to provide initial values manually.
    pub fn simulate_fixed_load_ratio_in_range(&self, fixed_load_ratio: f64, left_bc: &(dyn BoundaryCondition), right_bc: &(dyn BoundaryCondition), min_elec_cur: f64, max_elec_cur: f64, templ_hint: f64, tempr_hint: f64) -> LegSimulationResult {
        let elec_cur_hint: f64 = (min_elec_cur + max_elec_cur) * 0.5;
        let elec_cur_density_hint: f64 = elec_cur_hint / self.area;
        let heat_rate_hint: f64 = self.estimate_heat_flux(templ_hint, tempr_hint, elec_cur_density_hint) * self.area;
        let init_var: FullOptVar = FullOptVar { templ: templ_hint, heat_ratel: heat_rate_hint, tempr: tempr_hint, heat_rater: heat_rate_hint, elec_cur: elec_cur_hint };

        self.simulate_fixed_load_ratio_in_range_with_manual_init(fixed_load_ratio, left_bc, right_bc, min_elec_cur, max_elec_cur, &init_var)
    }

    /// Find the steady state under the given boundary condition and external-to-internal load ratio (1).
    /// * You need to provide the search range for electrical current (A).
    /// * You need to provide initial values for optimization by completing [`FullOptVar`](crate::bc::FullOptVar).
    pub fn simulate_fixed_load_ratio_in_range_with_manual_init(&self, fixed_load_ratio: f64, left_bc: &(dyn BoundaryCondition), right_bc: &(dyn BoundaryCondition), min_elec_cur: f64, max_elec_cur: f64, init_var: &FullOptVar) -> LegSimulationResult {
        let cost_f = |elec_cur: f64| -> f64 {
            let result_bc = self.simulate_fixed_elec_cur_with_manual_init(elec_cur, left_bc, right_bc, init_var);

            let cost: f64 = (result_bc.load_ratio() - fixed_load_ratio).powi(2);

            cost
        };
        let optimize_result: OptimizeResult = minimize_brent(&cost_f, min_elec_cur, max_elec_cur, self.simulation_param.max_iter_brent, self.simulation_param.abs_tol_elec_cur_brent);
        let elec_cur = optimize_result.x[0];
        let result_bc = self.simulate_fixed_elec_cur_with_manual_init(elec_cur, left_bc, right_bc, init_var);

        result_bc
    }

    /// Find the steady state where the power (W) is maximized.
    /// * You need to provide hints for left-end temperature (K) and right-end temperature (K).
    /// * If this does not converge, use [`SingleLeg::maximize_power_in_range()`] to provide a current range manually.
    pub fn maximize_power_auto_range(&self, left_bc: &(dyn BoundaryCondition), right_bc: &(dyn BoundaryCondition), templ_hint: f64, tempr_hint: f64) -> LegSimulationResult {
        let elec_cur_hint: f64 = self.estimate_elec_cur_for_max_power(templ_hint, tempr_hint);
        let min_elec_cur: f64;
        let max_elec_cur: f64;
        if elec_cur_hint > 0.0 {
            min_elec_cur = elec_cur_hint * 0.5;
            max_elec_cur = elec_cur_hint * 1.5;
        } else {
            min_elec_cur = elec_cur_hint * 1.5;
            max_elec_cur = elec_cur_hint * 0.5;
        }

        self.maximize_power_in_range(left_bc, right_bc, min_elec_cur, max_elec_cur, templ_hint, tempr_hint)
    }

    /// Find the steady state where the power (W) is maximized.
    /// * You need to provide hints for search range for electrical current (A) and hints for left-end temperature (K) and right-end temperature (K).
    /// * If this does not converge, use [`SingleLeg::maximize_power_in_range_with_manual_init()`] to provide initial values manually.
    pub fn maximize_power_in_range(&self, left_bc: &(dyn BoundaryCondition), right_bc: &(dyn BoundaryCondition), min_elec_cur: f64, max_elec_cur: f64, templ_hint: f64, tempr_hint: f64) -> LegSimulationResult {
        let elec_cur_hint: f64 = (min_elec_cur + max_elec_cur) * 0.5;
        let elec_cur_density_hint: f64 = elec_cur_hint / self.area;
        let heat_rate_hint: f64 = self.estimate_heat_flux(templ_hint, tempr_hint, elec_cur_density_hint) * self.area;
        let init_var: FullOptVar = FullOptVar { templ: templ_hint, heat_ratel: heat_rate_hint, tempr: tempr_hint, heat_rater: heat_rate_hint, elec_cur: elec_cur_hint };

        self.maximize_power_in_range_with_manual_init(left_bc, right_bc, min_elec_cur, max_elec_cur, &init_var)
    }

    /// Find the steady state where the power (W) is maximized.
    /// * You need to provide the search range for electrical current (A).
    /// * You need to provide initial values for optimization by completing [`FullOptVar`](crate::bc::FullOptVar).
    pub fn maximize_power_in_range_with_manual_init(&self, left_bc: &(dyn BoundaryCondition), right_bc: &(dyn BoundaryCondition), min_elec_cur: f64, max_elec_cur: f64, init_var: &FullOptVar) -> LegSimulationResult {
        let cost_f = |elec_cur: f64| -> f64 {
            let result_bc = self.simulate_fixed_elec_cur_with_manual_init(elec_cur, left_bc, right_bc, init_var);

            let cost: f64 = -result_bc.power;

            cost
        };
        let optimize_result: OptimizeResult = minimize_brent(&cost_f, min_elec_cur, max_elec_cur, self.simulation_param.max_iter_brent, self.simulation_param.abs_tol_elec_cur_brent);
        let elec_cur = optimize_result.x[0];
        let result_bc = self.simulate_fixed_elec_cur_with_manual_init(elec_cur, left_bc, right_bc, init_var);

        result_bc
    }

    /// Find the steady state where the efficiency (1) is maximized.
    /// * You need to provide hints for left-end temperature (K) and right-end temperature (K).
    /// * If this does not converge, use [`SingleLeg::maximize_efficiency_in_range()`] to provide a current range manually.
    pub fn maximize_efficiency_auto_range(&self, left_bc: &(dyn BoundaryCondition), right_bc: &(dyn BoundaryCondition), templ_hint: f64, tempr_hint: f64) -> LegSimulationResult {
        let elec_cur_hint: f64 = self.estimate_elec_cur_for_max_efficiency(templ_hint, tempr_hint);
        let min_elec_cur: f64;
        let max_elec_cur: f64;
        if elec_cur_hint > 0.0 {
            min_elec_cur = elec_cur_hint * 0.5;
            max_elec_cur = elec_cur_hint * 1.5;
        } else {
            min_elec_cur = elec_cur_hint * 1.5;
            max_elec_cur = elec_cur_hint * 0.5;
        }

        self.maximize_efficiency_in_range(left_bc, right_bc, min_elec_cur, max_elec_cur, templ_hint, tempr_hint)
    }

    /// Find the steady state where the efficiency (1) is maximized.
    /// * You need to provide hints for search range for electrical current (A) and hints for left-end temperature (K) and right-end temperature (K).
    /// * If this does not converge, use [`SingleLeg::maximize_efficiency_in_range_with_manual_init()`] to provide initial values manually.
    pub fn maximize_efficiency_in_range(&self, left_bc: &(dyn BoundaryCondition), right_bc: &(dyn BoundaryCondition), min_elec_cur: f64, max_elec_cur: f64, templ_hint: f64, tempr_hint: f64) -> LegSimulationResult {
        let elec_cur_hint: f64 = (min_elec_cur + max_elec_cur) * 0.5;
        let elec_cur_density_hint: f64 = elec_cur_hint / self.area;
        let heat_rate_hint: f64 = self.estimate_heat_flux(templ_hint, tempr_hint, elec_cur_density_hint) * self.area;
        let init_var: FullOptVar = FullOptVar { templ: templ_hint, heat_ratel: heat_rate_hint, tempr: tempr_hint, heat_rater: heat_rate_hint, elec_cur: elec_cur_hint };

        self.maximize_efficiency_in_range_with_manual_init(left_bc, right_bc, min_elec_cur, max_elec_cur, &init_var)
    }

    /// Find the steady state where the efficiency (1) is maximized.
    /// * You need to provide the search range for electrical current (A).
    /// * You need to provide initial values for optimization by completing [`FullOptVar`](crate::bc::FullOptVar).
    pub fn maximize_efficiency_in_range_with_manual_init(&self, left_bc: &(dyn BoundaryCondition), right_bc: &(dyn BoundaryCondition), min_elec_cur: f64, max_elec_cur: f64, init_var: &FullOptVar) -> LegSimulationResult {
        let cost_f = |elec_cur: f64| -> f64 {
            let result_bc = self.simulate_fixed_elec_cur_with_manual_init(elec_cur, left_bc, right_bc, init_var);

            let cost: f64 = -result_bc.efficiency;

            cost
        };
        let optimize_result: OptimizeResult = minimize_brent(&cost_f, min_elec_cur, max_elec_cur, self.simulation_param.max_iter_brent, self.simulation_param.abs_tol_elec_cur_brent);
        let elec_cur = optimize_result.x[0];
        let result_bc = self.simulate_fixed_elec_cur_with_manual_init(elec_cur, left_bc, right_bc, init_var);

        result_bc
    }

    /// Solve the initial-value problem and return the solution matrix (`T`, `q`).
    /// * The result is a 2x`N` matrix, where `N` is the length of `xs`.
    fn solve_ivp(&self, elec_cur_density: f64, templ: f64, heat_fluxl: f64, xs: &ColVec) -> Array2<f64> {
        let f = |_x: f64, y: &ColVec| -> ColVec {
            teg_rhs(&self.tep, y[0], y[1], elec_cur_density)
        };

        explicit_rk4(&f, xs, &array![templ, heat_fluxl])
    }

    /// Solve the final-value problem and return the solution matrix (`T`, `q`).
    /// * The result is a 2x`N` matrix, where `N` is the length of `xs`.
    fn solve_fvp(&self, elec_cur_density: f64, tempr: f64, heat_fluxr: f64, xs: &ColVec) -> Array2<f64> {
        let f = |_x: f64, y: &ColVec| -> ColVec {
            teg_rhs(&self.tep, y[0], y[1], elec_cur_density)
        };

        explicit_rk4_fvp(&f, xs, &array![tempr, heat_fluxr])
    }

    /// Returns the voltage (V) induced from the Seebeck effect
    fn compute_voltage(&self, templ: f64, tempr: f64) -> f64 {
        let temp_diff = templ - tempr;
        if temp_diff == 0.0 {
            return 0.0;
        }

        let num_temp_mesh: usize = get_num_mesh(temp_diff.abs(), self.simulation_param.dtemp, self.simulation_param.min_num_temp_mesh, self.simulation_param.max_num_temp_mesh);
        let temp_mesh: ColVec;
        let h: f64;

        if tempr < templ {
            temp_mesh = Array::linspace(tempr, templ, num_temp_mesh);
            h = temp_mesh[1] - temp_mesh[0];
        } else {
            temp_mesh = Array::linspace(templ, tempr, num_temp_mesh);
            h = temp_mesh[0] - temp_mesh[1];  // Here the mesh size `h` is negative.
        }
        
        integrate_simpson(&self.tep.seebeck_ats_incr(&temp_mesh), h)
    }

    /// Returns the resistance (Ohm) inside the leg.
    /// * The argument `temps` must be `T(xs)` where `xs` is a member variable of `SingleLeg` that is a equidistant spatial mesh grid.
    fn compute_internal_resistance(&self, temps: &ColVec) -> f64 {
        let rho_of_x: ColVec = self.tep.elec_resi_ats(temps);
        let h = self.xs[1] - self.xs[0];

        integrate_simpson(&rho_of_x, h) / self.area
    }

    /// Returns the power (W).
    #[inline]
    fn get_power(&self, heat_fluxl: f64, heat_fluxr: f64) -> f64 {
        (heat_fluxl - heat_fluxr) * self.area
    }

    /// Returns an estimation of the heat flux
    pub fn estimate_heat_flux(&self, templ: f64, tempr: f64, elec_cur_density: f64) -> f64 {
        let alpha = self.tep.seebeck_at(templ);
        let kappa = self.tep.thrm_cond_at(templ);

        alpha*templ*elec_cur_density - kappa*(tempr-templ)/self.length
    }

    /// Returns an approximate value of the maximum efficiency.
    /// * The approximation is performed using the Equation 2--4 in [B. Ryu, J. Chung and S. Park (2021)](https://doi.org/10.1016/j.isci.2021.102934).
    /// * The approximation is based on the theory of three thermoelectric degrees of freedom suggested by the authors; the theory uses three figures of merit (`Z_gen`, `tau` and `beta`).
    pub fn estimate_max_efficiency(&self, templ: f64, tempr: f64) -> f64 {
        let hot_temp: f64;
        let cold_temp: f64;

        if templ < tempr {
            cold_temp = templ;
            hot_temp = tempr;
        } else {
            cold_temp = tempr;
            hot_temp = templ;
        }

        let delta_temp: f64 = hot_temp - cold_temp;
        let num_temp_mesh: usize = get_num_mesh(delta_temp, self.simulation_param.dtemp, self.simulation_param.min_num_temp_mesh, self.simulation_param.max_num_temp_mesh);
        let temp_mesh: ColVec = Array::linspace(cold_temp, hot_temp, num_temp_mesh);
        let h: f64 = temp_mesh[1] - temp_mesh[0];

        let alpha_bar: f64 = integrate_simpson(&self.tep.seebeck_ats_incr(&temp_mesh), h) / delta_temp;
        let rho_kappa_bar: f64 = integrate_simpson(&(self.tep.elec_resi_ats_incr(&temp_mesh)*self.tep.thrm_cond_ats_incr(&temp_mesh)), h) / delta_temp;

        // approximate the three degrees of freedom
        let z_gen: f64 = alpha_bar.powi(2) / rho_kappa_bar;
        let tau: f64 = (-1.0/6.0) * (self.tep.seebeck_at(hot_temp) - self.tep.seebeck_at(cold_temp)) / alpha_bar;
        let beta: f64 = (1.0/6.0) * (self.tep.elec_resi_at(hot_temp)*self.tep.thrm_cond_at(hot_temp) - self.tep.elec_resi_at(cold_temp)*self.tep.thrm_cond_at(cold_temp)) / rho_kappa_bar;

        // perturb the temperature range
        let temp_hp: f64 = hot_temp - tau * delta_temp;
        let temp_cp: f64 = cold_temp - (tau + beta)*delta_temp;
        let temp_midp: f64 = (temp_hp+temp_cp)*0.5;

        // approximate the max. efficiency
        let approx_max_eff = (delta_temp/temp_hp) * ((1.0+z_gen*temp_midp).sqrt()-1.0) / ((1.0+z_gen*temp_midp).sqrt()+(temp_cp/temp_hp));

        approx_max_eff
    }

    /// Returns an approximate value of the maximum efficiency (1).
    /// * The approximation is performed using the only `Z_gen` in the Equation 2--4 in [B. Ryu, J. Chung and S. Park (2021)](https://doi.org/10.1016/j.isci.2021.102934); that is, we assume `tau=0` and `beta=0`.
    /// * This approximation is exact if the Seebeck coefficient is temperature-independent; for its proof, refer to [J. Chung, B. Ryu and H. Seo](https://arxiv.org/abs/2106.10957).
    pub fn estimate_max_efficiency_by_z(&self, templ: f64, tempr: f64) -> f64 {
        let hot_temp: f64;
        let cold_temp: f64;

        if templ < tempr {
            cold_temp = templ;
            hot_temp = tempr;
        } else {
            cold_temp = tempr;
            hot_temp = templ;
        }

        let delta_temp: f64 = hot_temp - cold_temp;
        let num_temp_mesh: usize = get_num_mesh(delta_temp, self.simulation_param.dtemp, self.simulation_param.min_num_temp_mesh, self.simulation_param.max_num_temp_mesh);
        let temp_mesh: ColVec = Array::linspace(cold_temp, hot_temp, num_temp_mesh);
        let h: f64 = temp_mesh[1] - temp_mesh[0];

        let alpha_bar: f64 = integrate_simpson(&self.tep.seebeck_ats_incr(&temp_mesh), h) / delta_temp;
        let rho_kappa_bar: f64 = integrate_simpson(&(self.tep.elec_resi_ats_incr(&temp_mesh)*self.tep.thrm_cond_ats_incr(&temp_mesh)), h) / delta_temp;

        // approximate the `Z_gen`
        let z_gen: f64 = alpha_bar.powi(2) / rho_kappa_bar;

        // perturb the temperature range
        let mid_temp: f64 = (hot_temp + cold_temp) * 0.5;

        // approximate the max. efficiency (`tau=0` and `beta=0` case)
        let approx_max_eff = (delta_temp/hot_temp) * ((1.0+z_gen*mid_temp).sqrt()-1.0) / ((1.0+z_gen*mid_temp).sqrt()+(cold_temp/hot_temp));

        approx_max_eff
    }

    /// Returns an approximate value of the maximum power (W).
    /// * The approximation is performed using the formula `P_max = V^2/(4*R)` from the constant property model, where averaged material properties are used in computation of voltage `V` and internal resistance `R`.
    /// * Seebeck coefficient is averaged over temperature range.
    /// * Electrical resistivity is computed by the average of `rho*kappa` divided by the average of `kappa` (`rho` is electrical resistivity and `kappa` is thermal conductivity).
    pub fn estimate_max_power(&self, templ: f64, tempr: f64) -> f64 {
        let hot_temp: f64;
        let cold_temp: f64;

        if templ < tempr {
            cold_temp = templ;
            hot_temp = tempr;
        } else {
            cold_temp = tempr;
            hot_temp = templ;
        }

        let delta_temp: f64 = hot_temp - cold_temp;
        let num_temp_mesh: usize = get_num_mesh(delta_temp, self.simulation_param.dtemp, self.simulation_param.min_num_temp_mesh, self.simulation_param.max_num_temp_mesh);
        let temp_mesh: ColVec = Array::linspace(cold_temp, hot_temp, num_temp_mesh);
        let h: f64 = temp_mesh[1] - temp_mesh[0];

        let alpha_bar: f64 = integrate_simpson(&self.tep.seebeck_ats_incr(&temp_mesh), h) / delta_temp;
        let rho_kappa_bar: f64 = integrate_simpson(&(self.tep.elec_resi_ats_incr(&temp_mesh)*self.tep.thrm_cond_ats_incr(&temp_mesh)), h) / delta_temp;
        let kappa_bar: f64 = integrate_simpson(&self.tep.thrm_cond_ats_incr(&temp_mesh), h) / delta_temp;

        // averaged material properties
        let alpha0: f64 = alpha_bar;
        let rho0: f64 = rho_kappa_bar / kappa_bar;

        // device parameters
        let voltage: f64 = alpha0 * delta_temp;
        let internal_resistance: f64 = rho0 * self.length / self.area;

        // approximate the max. power
        let approx_max_power = voltage.powi(2) / (4.0 * internal_resistance);

        approx_max_power
    }

    /// Returns an approximate value of the electrical current (A) attaing the maximum power. If the material is n-type, the return value is negative, representing the reversed direction of the electrical current.
    /// * The approximation is performed using the formula `I_(P_max) = V/(2*R)`, where averaged material properties are used in computation of voltage `V` and internal resistance `R`.
    /// * Seebeck coefficient is averaged over temperature range.
    /// * Electrical resistivity is computed by the average of `rho*kappa` divided by the average of `kappa` (`rho` is electrical resistivity and `kappa` is thermal conductivity).
    pub fn estimate_elec_cur_for_max_power(&self, templ: f64, tempr: f64) -> f64 {
        let hot_temp: f64;
        let cold_temp: f64;

        if templ < tempr {
            cold_temp = templ;
            hot_temp = tempr;
        } else {
            cold_temp = tempr;
            hot_temp = templ;
        }

        let delta_temp: f64 = hot_temp - cold_temp;
        let num_temp_mesh: usize = get_num_mesh(delta_temp, self.simulation_param.dtemp, self.simulation_param.min_num_temp_mesh, self.simulation_param.max_num_temp_mesh);
        let temp_mesh: ColVec = Array::linspace(cold_temp, hot_temp, num_temp_mesh);
        let h: f64 = temp_mesh[1] - temp_mesh[0];

        let alpha_bar: f64 = integrate_simpson(&self.tep.seebeck_ats_incr(&temp_mesh), h) / delta_temp;
        let rho_kappa_bar: f64 = integrate_simpson(&(self.tep.elec_resi_ats_incr(&temp_mesh)*self.tep.thrm_cond_ats_incr(&temp_mesh)), h) / delta_temp;
        let kappa_bar: f64 = integrate_simpson(&self.tep.thrm_cond_ats_incr(&temp_mesh), h) / delta_temp;

        // averaged material properties
        let alpha0: f64 = alpha_bar;
        let rho0: f64 = rho_kappa_bar / kappa_bar;

        // device parameters
        let voltage: f64 = alpha0 * delta_temp;
        let internal_resistance: f64 = rho0 * self.length / self.area;

        // approximate the electrical current
        let load_ratio = 1.0;
        let approx_elec_cur_for_max_power = voltage / ((1.0 + load_ratio) * internal_resistance);

        approx_elec_cur_for_max_power
    }

    /// Returns an approximate value of the electrical current (A) attaing the maximum efficiency. If the material is n-type, the return value is negative, representing the reversed direction of the electrical current.
    /// * The approximation is performed using the formula `I_(eff_max) = V/((1+sqrt(1+z_gen*Tm))*R)`, where averaged material properties are used in computation of voltage `V` and internal resistance `R`.
    /// * Seebeck coefficient is averaged over temperature range.
    /// * Electrical resistivity is computed by the average of `rho*kappa` divided by the average of `kappa` (`rho` is electrical resistivity and `kappa` is thermal conductivity).
    pub fn estimate_elec_cur_for_max_efficiency(&self, templ: f64, tempr: f64) -> f64 {
        let hot_temp: f64;
        let cold_temp: f64;

        if templ < tempr {
            cold_temp = templ;
            hot_temp = tempr;
        } else {
            cold_temp = tempr;
            hot_temp = templ;
        }

        let delta_temp: f64 = hot_temp - cold_temp;
        let mid_temp: f64 = (hot_temp + cold_temp) * 0.5;
        let num_temp_mesh: usize = get_num_mesh(delta_temp, self.simulation_param.dtemp, self.simulation_param.min_num_temp_mesh, self.simulation_param.max_num_temp_mesh);
        let temp_mesh: ColVec = Array::linspace(cold_temp, hot_temp, num_temp_mesh);
        let h: f64 = temp_mesh[1] - temp_mesh[0];

        let alpha_bar: f64 = integrate_simpson(&self.tep.seebeck_ats_incr(&temp_mesh), h) / delta_temp;
        let rho_kappa_bar: f64 = integrate_simpson(&(self.tep.elec_resi_ats_incr(&temp_mesh)*self.tep.thrm_cond_ats_incr(&temp_mesh)), h) / delta_temp;
        let kappa_bar: f64 = integrate_simpson(&self.tep.thrm_cond_ats_incr(&temp_mesh), h) / delta_temp;

        // averaged material properties
        let alpha0: f64 = alpha_bar;
        let rho0: f64 = rho_kappa_bar / kappa_bar;

        // device parameters
        let voltage: f64 = alpha0 * delta_temp;
        let internal_resistance: f64 = rho0 * self.length / self.area;

        // approximate the `Z_gen`
        let z_gen: f64 = alpha_bar.powi(2) / rho_kappa_bar;

        // approximate the max. electrical current
        let load_ratio = (1.0 + z_gen*mid_temp).sqrt();
        let approx_elec_cur_for_max_efficiency = voltage / ((1.0+load_ratio)*internal_resistance);

        approx_elec_cur_for_max_efficiency
    }
}


/// Returns the right-hand side of `[T',q']=f(T,q,J)` where `T` is temperature, `q` is heat flux and `J` is electrical current density.
/// * The derivative is a spatial derivative (position inside a leg).
/// * For details, refer to p.43 and p.45 in [C. Goupil (ed.)](https://onlinelibrary.wiley.com/doi/book/10.1002/9783527338405).
/// 
/// # Arguments
/// `tep`: Thermoelectric material properties
/// `temp`: temperature
/// `heat_flux`: heat flux
/// `elec_cur_density`: electrical current density flowing through a leg
fn teg_rhs(tep: &Tep, temp: f64, heat_flux: f64, elec_cur_density: f64) -> ColVec {
    let alpha = tep.seebeck_at(temp);
    let rho = tep.elec_resi_at(temp);
    let kappa = tep.thrm_cond_at(temp);
    let q = heat_flux;
    let j = elec_cur_density;

    let dtempdx = (alpha*temp*j - q)/kappa;
    let dqdx = (alpha*dtempdx + rho*j) * j;

    array![dtempdx, dqdx]
}


/// Returns the number of mesh points satisfying the following conditions:
/// * If possible, return a close number of meshes having the length of `suggested_mesh_size`.
/// * The number must be larger than or equal to `min_num_mesh`.
/// * The number must be larger than or equal to four. This is becuase we use [`PiecewiseCubicPolynomial`](crate::interp::PiecewiseCubicPolynomial).
/// * The number must be odd. This is because we use [`integrate_simpson`](create::calc::integrate_simpson) for integration.
/// 
/// # Panics
/// * when `total_size` is not positive.
/// * when `suggested_mesh_size` is not positive.
fn get_num_mesh(total_size: f64, suggested_mesh_size: f64, min_num_mesh: usize, max_num_mesh: usize) -> usize {
    if total_size <= 0.0 {
        panic!("`total_size` must be a positive number!");
    }
    if suggested_mesh_size <= 0.0 {
        panic!("`suggested_mesh_size` must be a positive number!");
    }
    let mut num_mesh: usize = (total_size / suggested_mesh_size).ceil() as usize + 1;
    if num_mesh > max_num_mesh {
        num_mesh = max_num_mesh - 1;
    }
    if num_mesh < min_num_mesh {
        num_mesh = min_num_mesh;
    }
    if num_mesh < 4 {
        num_mesh = 4;
    }
    if num_mesh > max_num_mesh {
        num_mesh = max_num_mesh - 1;
    }
    // if the number is even, make it odd
    if num_mesh & 1 == 0 {
        num_mesh += 1;
    }

    num_mesh
}

/// Returns the efficiency
#[inline]
fn get_efficiency(heat_fluxl: f64, heat_fluxr: f64) -> f64 {
    1.0 - heat_fluxr/heat_fluxl
}


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

    use crate::calc::{all_close, is_close, max};
    use crate::func::{ConstantFunction, ExplicitFunction};
    use crate::interp::PiecewiseLinearConstExt;
    use crate::bc::{FixedTempBc, FixedHeatFluxBc, ConvectionBc};

    /// Test using linearly increasing thermal conductivity.
    #[test]
    #[ignore]
    fn test_estimate_max_power_for_linear_thrm_cond() {
        // environment condition
        let hot_temp: f64 = 500.0;
        let cold_temp: f64 = 300.0;

        // set tep and leg parameters
        let alpha0: f64 = 150e-6; // [V/K]
        let rho0: f64 = 1e-5; // [Ohm m]
        let kappa0: f64 = 2.0; // [W/m/K]
        let temp0: f64 = 300.0; // [K]
        let seebeck = ConstantFunction::new(alpha0);
        let elec_resi = ConstantFunction::new(rho0);
        let thrm_cond = ExplicitFunction::new(Rc::new(move |temp: f64| { max(kappa0*temp/temp0, kappa0*cold_temp/temp0) }));   // `thrm_cond` must be positive.
        let tep = Tep::new(Rc::new(seebeck), Rc::new(elec_resi), Rc::new(thrm_cond));
        let length = 1e-3; // [m]
        let area = 1e-6;  // m^2]

        // set simulation parameters
        let minimizer_kind = MinimizerKind::Bfgs;
        let minimizer_param = MinimizerParam {
            grad_tol: 1e-5,
            max_iter: 100,
            ..MinimizerParam::default(minimizer_kind)
        };
        let simul_param = SimulationParam {
            dx: length/10.0,
            abs_tol_elec_cur_brent: 1e-10,
            minimizer_kind,
            minimizer_param,
            ..SimulationParam::default()
        };

        // create a leg
        let leg = SingleLeg::new(&tep, length, area, simul_param);

        // check the temperature
        let now = Instant::now();
        let simul_result = leg.maximize_power_auto_range(&FixedTempBc::new(500.0), &FixedTempBc::new(300.0), hot_temp, cold_temp);
        println!("test_estimate_max_power_for_linear_thrm_cond: Elapsed {:.2?}", now.elapsed());
        if !simul_result.success {
            println!("error: {:?}", simul_result.error);
        }
        assert!(simul_result.success);
        let approx_power: f64 = leg.estimate_max_power(hot_temp, cold_temp);
        let approx_elec_cur: f64 = leg.estimate_elec_cur_for_max_power(hot_temp, cold_temp);
        println!("estimated power={}, numerical power={}", approx_power, simul_result.power);
        assert!(is_close(simul_result.power, approx_power, 5e-06, 0.0));  // 3.0e-6 rel. error, 6.7e-8 abs. error
        println!("estimated elec_cur={}, numerical elec_cur={}", approx_elec_cur, simul_result.elec_cur);
        assert!(is_close(approx_elec_cur, simul_result.elec_cur, 5e-06, 0.0));  // 3.6e-6 rel. error, 5.4e-6 abs. error
    }

    /// Test using linearly increasing thermal conductivity.
    #[test]
    #[ignore]
    fn test_estimate_max_efficiency_for_linear_thrm_cond() {
        // environment condition
        let hot_temp: f64 = 500.0;
        let cold_temp: f64 = 300.0;

        // set tep and leg parameters
        let alpha0: f64 = 150e-6; // [V/K]
        let rho0: f64 = 1e-5; // [Ohm m]
        let kappa0: f64 = 2.0; // [W/m/K]
        let temp0: f64 = 300.0; // [K]
        let seebeck = ConstantFunction::new(alpha0);
        let elec_resi = ConstantFunction::new(rho0);
        let thrm_cond = ExplicitFunction::new(Rc::new(move |temp: f64| { max(kappa0*temp/temp0, kappa0*cold_temp/temp0) }));   // `thrm_cond` must be positive.
        let tep = Tep::new(Rc::new(seebeck), Rc::new(elec_resi), Rc::new(thrm_cond));
        let length = 1e-3; // [m]
        let area = 1e-6;  // m^2]

        // set simulation parameters
        let minimizer_kind = MinimizerKind::Bfgs;
        let minimizer_param = MinimizerParam {
            grad_tol: 1e-5,
            max_iter: 100,
            ..MinimizerParam::default(minimizer_kind)
        };
        let simul_param = SimulationParam {
            dx: length/10.0,
            abs_tol_elec_cur_brent: 1e-6,
            minimizer_kind,
            minimizer_param,
            ..SimulationParam::default()
        };

        // create a leg
        let leg = SingleLeg::new(&tep, length, area, simul_param);

        // check the temperature
        let now = Instant::now();
        let simul_result = leg.maximize_efficiency_auto_range(&FixedTempBc::new(500.0), &FixedTempBc::new(300.0), hot_temp, cold_temp);
        println!("test_estimate_max_efficiency_for_linear_thrm_cond: Elapsed {:.2?}", now.elapsed());
        if !simul_result.success {
            println!("error: {:?}", simul_result.error);
        }
        assert!(simul_result.success);
        let approx_eff: f64 = leg.estimate_max_efficiency(hot_temp, cold_temp);
        let approx_eff_by_z: f64 = leg.estimate_max_efficiency_by_z(hot_temp, cold_temp);
        let approx_elec_cur: f64 = leg.estimate_elec_cur_for_max_efficiency(hot_temp, cold_temp);
        println!("estimated eff by Ztb={}, estimated eff by CSM={}, numerical eff={}", approx_eff, approx_eff_by_z, simul_result.efficiency);
        assert!(is_close(simul_result.efficiency, approx_eff, 2e-03, 0.0));  // 1.3e-3 rel. error, 4.6e-5 abs. error
        assert!(is_close(simul_result.efficiency, approx_eff_by_z, 5e-06, 0.0));  // 3.6e-6 rel. error, 1.3e-7 abs. error
        println!("estimated elec_cur={}, numerical elec_cur={}", approx_elec_cur, simul_result.elec_cur);
        assert!(is_close(approx_elec_cur, simul_result.elec_cur, 5e-05, 0.0));  // 3.0e-5 rel. error, 4.2e-5 abs. error
    }

    /// Test using constant property model.
    /// The example is from p.89 in [C. Goupil (ed.)](https://onlinelibrary.wiley.com/doi/book/10.1002/9783527338405).
    #[test]
    #[ignore]
    fn test_maximize_power_for_cpm() {
        // environment condition
        let hot_temp: f64 = 500.0;
        let cold_temp: f64 = 300.0;
        let delta_temp = hot_temp - cold_temp;
        
        // set tep and leg parameters
        let alpha0: f64 = 150e-6;
        let rho0: f64 = 1e-5;
        let kappa0: f64 = 2.0;
        let seebeck = ConstantFunction::new(alpha0);
        let elec_resi = ConstantFunction::new(rho0);
        let thrm_cond = ConstantFunction::new(kappa0);
        let tep = Tep::new(Rc::new(seebeck), Rc::new(elec_resi), Rc::new(thrm_cond));
        let length = 1e-3;
        let area = 1e-6;

        // set simulation parameters
        let minimizer_kind = MinimizerKind::Bfgs;
        let minimizer_param = MinimizerParam {
            grad_tol: 1e-5,
            a_max: 2.0,
            max_iter: 200,
            ..MinimizerParam::default(minimizer_kind)
        };
        let simul_param = SimulationParam {
            dx: length/10.0,
            dtemp: (hot_temp-cold_temp)/10.0,
            abs_tol_elec_cur_brent: 1e-6,
            minimizer_kind,
            minimizer_param,
            ..SimulationParam::default()
        };

        // create a leg and simulate
        let leg = SingleLeg::new(&tep, length, area, simul_param);
        let left_bc = FixedTempBc::new(hot_temp);
        let right_bc = FixedTempBc::new(cold_temp);

        // exact solution
        let internal_resistance = rho0*length/area;
        let voltage = alpha0*delta_temp;
        let max_power = voltage.powi(2)/(4.0 * internal_resistance);
        let estimated_elec_cur = voltage/(2.0 * internal_resistance);
        println!("estimated_elec_cur={}", estimated_elec_cur);

        let now = Instant::now();
        // let result_fixed_temp = leg.maximize_power_in_range(&left_bc, &right_bc, estimated_elec_cur*0.8, estimated_elec_cur*1.2, hot_temp, cold_temp);  // 370.3089ms
        let result_fixed_temp = leg.maximize_power_auto_range(&left_bc, &right_bc, hot_temp, cold_temp);  // 370.3089ms
        println!("test_maximize_power_for_cpm(): elapsed: {:?}", now.elapsed());
        assert!(result_fixed_temp.success);

        println!("{}", result_fixed_temp);
        println!("max_power={}, internal_resistance={}", max_power, internal_resistance);

        assert!(is_close(result_fixed_temp.power, max_power, 1e-04, 0.0));  // 0.01% rel. error
        assert!(is_close(result_fixed_temp.internal_resistance, internal_resistance, 1e-04, 0.0));
    }

    /// Test using linearly increasing electrical resistivity.
    #[test]
    #[ignore]
    fn test_maximize_efficiency_for_linear_elec_resi() {
        // environment condition
        let hot_temp: f64 = 500.0;
        let cold_temp: f64 = 300.0;
        let delta_temp = hot_temp - cold_temp;
        let mid_temp: f64 = (hot_temp+cold_temp) * 0.5;

        // set tep and leg parameters
        let alpha0: f64 = 150e-6;
        let rho0: f64 = 1e-5;
        let kappa0: f64 = 2.0;
        let temp0: f64 = 400.0;
        let seebeck = ConstantFunction::new(alpha0);
        let elec_resi = ExplicitFunction::new(Rc::new(move |temp: f64| { max(rho0*temp/temp0, 1e-08) })); // `elec_resi` must be positive.
        let thrm_cond = ConstantFunction::new(kappa0);
        let tep = Tep::new(Rc::new(seebeck), Rc::new(elec_resi), Rc::new(thrm_cond));
        let length = 1e-3;
        let area = 1e-6;

        // set simulation parameters
        let minimizer_kind = MinimizerKind::Bfgs;
        let minimizer_param = MinimizerParam {
            grad_tol: 1e-7,
            a_max: 2.0,
            max_iter: 200,
            ..MinimizerParam::default(minimizer_kind)
        };
        let simul_param = SimulationParam {
            dx: length/10.0,
            dtemp: (hot_temp-cold_temp)/10.0,
            abs_tol_elec_cur_brent: 1e-6,
            minimizer_kind,
            minimizer_param,
            ..SimulationParam::default()
        };

        // create a leg
        let leg = SingleLeg::new(&tep, length, area, simul_param);

        // simulate the leg
        let left_bc = FixedTempBc::new(hot_temp);
        let right_bc = FixedTempBc::new(cold_temp);
        // let result_fixed_temp = leg.maximize_efficiency_in_range(&left_bc, &right_bc, elec_cur*0.5, elec_cur*1.5, hot_temp, cold_temp);
        let result_fixed_temp = leg.maximize_efficiency_auto_range(&left_bc, &right_bc, hot_temp, cold_temp);
        assert!(result_fixed_temp.success);
        println!("result = {}", result_fixed_temp);
        // test helper functions
        let load_resistance = result_fixed_temp.load_resistance();
        assert!(is_close(result_fixed_temp.total_resistance(), result_fixed_temp.internal_resistance + load_resistance, 1e-05, 0.0));
        assert!(is_close(result_fixed_temp.load_ratio(), load_resistance / result_fixed_temp.internal_resistance, 1e-05, 0.0));
        // compute exact max. efficiency
        let z = alpha0.powi(2)/(rho0*kappa0);
        let max_efficiency = (delta_temp/hot_temp) * ((1.0+z*mid_temp).sqrt()-1.0) / ((1.0+z*mid_temp).sqrt()+(cold_temp/hot_temp));
        println!("computed eff={}, exact max eff={}", result_fixed_temp.efficiency, max_efficiency);
        assert!(is_close(result_fixed_temp.efficiency, max_efficiency, 1e-04, 0.0));  // 0.01% rel. error
    }

    /// Test using linearly increasing thermal conductivity.
    #[test]
    #[ignore]
    fn test_maximize_efficiency_for_cpm() {
        // environment condition
        let hot_temp: f64 = 500.0;
        let cold_temp: f64 = 300.0;
        let delta_temp = hot_temp-cold_temp;
        let mid_temp: f64 = (hot_temp+cold_temp) * 0.5;

        // set tep and leg parameters
        let alpha0: f64 = 150e-6; // [V/K]
        let rho0: f64 = 1e-5; // [Ohm m]
        let kappa0: f64 = 2.0; // [W/m/K]
        let seebeck = ConstantFunction::new(alpha0);
        let elec_resi = ConstantFunction::new(rho0);
        let thrm_cond = ConstantFunction::new(kappa0);
        let tep = Tep::new(Rc::new(seebeck), Rc::new(elec_resi), Rc::new(thrm_cond));
        let length = 1e-3; // [m]
        let area = 1e-6;  // m^2]

        // compute exact max. efficiency
        let z = alpha0.powi(2)/(rho0*kappa0);
        let max_efficiency = (delta_temp/hot_temp) * ((1.0+z*mid_temp).sqrt()-1.0) / ((1.0+z*mid_temp).sqrt()+(cold_temp/hot_temp));

        // set simulation parameters
        let minimizer_kind = MinimizerKind::Bfgs;
        let minimizer_param = MinimizerParam {
            grad_tol: 1e-5,
            a_max: 2.0,
            max_iter: 200,
            ..MinimizerParam::default(minimizer_kind)
        };
        let simul_param = SimulationParam {
            dx: length/10.0,
            dtemp: (hot_temp-cold_temp)/10.0,
            abs_tol_elec_cur_brent: 1e-6,
            minimizer_kind,
            minimizer_param,
            ..SimulationParam::default()
        };

        // create a leg
        let leg = SingleLeg::new(&tep, length, area, simul_param);

        // simulate the leg
        let left_bc = FixedTempBc::new(hot_temp);
        let right_bc = FixedTempBc::new(cold_temp);
        // let result_fixed_temp = leg.maximize_efficiency_in_range(&left_bc, &right_bc, elec_cur*0.8, elec_cur*1.2, hot_temp, cold_temp);
        let result_fixed_temp = leg.maximize_efficiency_auto_range(&left_bc, &right_bc, hot_temp, cold_temp);
        println!("result={}", result_fixed_temp);
        println!("numerical eff={}, exact max eff={}", result_fixed_temp.efficiency, max_efficiency);
        assert!(result_fixed_temp.success);
        assert!(is_close(result_fixed_temp.efficiency, max_efficiency, 1e-04, 0.0));  // 0.01% rel. error
    }

    /// Test using constant property model.
    #[test]
    #[ignore]
    fn test_single_leg_simulate_fixed_load_ratio_for_cpm() {
        // environment condition
        let hot_temp: f64 = 500.0;
        let cold_temp: f64 = 300.0;
        let delta_temp: f64 = hot_temp - cold_temp;

        // set tep and leg parameters
        let alpha0: f64 = 150e-6;
        let rho0: f64 = 1e-5;
        let kappa0: f64 = 2.0;
        let seebeck = ConstantFunction::new(alpha0);
        let elec_resi = ConstantFunction::new(rho0);
        let thrm_cond = ConstantFunction::new(kappa0);
        let tep = Tep::new(Rc::new(seebeck), Rc::new(elec_resi), Rc::new(thrm_cond));
        let length = 1e-3;
        let area = 1e-6;
        let elec_cur: f64 = 1.0;

        // set simulation parameters
        let minimizer_kind = MinimizerKind::DampedBfgs;
        let minimizer_param = MinimizerParam {
            grad_tol: 1e-7,
            a_max: 2.0,
            max_iter: 200,
            ..MinimizerParam::default(minimizer_kind)
        };
        let simul_param = SimulationParam {
            dx: length/10.0,
            dtemp: delta_temp/10.0,
            minimizer_kind,
            minimizer_param,
            ..SimulationParam::default()
        };

        // create a leg
        let leg = SingleLeg::new(&tep, length, area, simul_param);

        // case 1: fixed temperatures
        let left_bc = FixedTempBc::new(hot_temp);
        let right_bc = FixedTempBc::new(cold_temp);
        let simul_result_fixed_elec_cur = leg.simulate_fixed_elec_cur(elec_cur, &left_bc, &right_bc, hot_temp, cold_temp);
        assert!(simul_result_fixed_elec_cur.success);
        let fixed_load_ratio = simul_result_fixed_elec_cur.load_ratio();

        // test a time-consuming method
        let now = Instant::now();
        // let simul_result_fixed_load_ratio = leg.simulate_fixed_load_ratio(simul_result_fixed_elec_cur.load_ratio(), &left_bc, &right_bc, elec_cur, hot_temp, cold_temp);
        let simul_result_fixed_load_ratio = leg.simulate_fixed_load_ratio(fixed_load_ratio, &left_bc, &right_bc, elec_cur/2.0, hot_temp, cold_temp);
        println!("For fixed BC and fixed load ratio: Elapsed {:.2?}", now.elapsed());
        assert!(simul_result_fixed_load_ratio.success);
        assert!(is_close(simul_result_fixed_load_ratio.load_ratio(), fixed_load_ratio, 2e-04, 0.0));
        assert!(is_close(simul_result_fixed_load_ratio.elec_cur, elec_cur, 2e-04, 0.0));
        println!("result={}", simul_result_fixed_load_ratio);
        println!("computed load_ratio={}", simul_result_fixed_load_ratio.load_ratio());

        // test rtg
        let left_bc = FixedHeatFluxBc::new(100.0);
        let right_bc = ConvectionBc::new(30.0, 8.0);
        // let fixed_load_ratio = 0.2;
        let now = Instant::now();
        let simul_result_fixed_load_ratio = leg.simulate_fixed_load_ratio(fixed_load_ratio, &left_bc, &right_bc, elec_cur, hot_temp, cold_temp);
        println!("For RTG BC and fixed load ratio: Elapsed {:.2?}", now.elapsed());
        println!("result={}", simul_result_fixed_load_ratio);
        println!("computed load_ratio={}", simul_result_fixed_load_ratio.load_ratio());
        println!("templ={}, tempr={}", simul_result_fixed_load_ratio.temp.at(0.0), simul_result_fixed_load_ratio.temp.at(leg.length));
        if !simul_result_fixed_load_ratio.success {
            println!("rtg error={:?}", simul_result_fixed_load_ratio.error);
        }
        //assert!(simul_result_fixed_load_ratio.success);
        assert!(is_close(simul_result_fixed_load_ratio.load_ratio(), fixed_load_ratio, 2e-04, 0.0));
    }

    /// Test using temperature-dependent electrical resistivity.
    /// The example is from p.123 in [C. Goupil (ed.)](https://onlinelibrary.wiley.com/doi/book/10.1002/9783527338405).
    #[test]
    #[ignore]
    fn test_single_leg_simulate_fixed_load_resistance_for_linear_elec_resi() {
        // environment condition
        let hot_temp: f64 = 500.0;
        let cold_temp: f64 = 300.0;

        // set tep and leg parameters
        let alpha0: f64 = 150e-6;
        let rho0: f64 = 1e-5;
        let temp0: f64 = 400.0;
        let kappa0: f64 = 2.0;
        let seebeck = PiecewiseLinearConstExt::new(&array![0.0, cold_temp, hot_temp], &array![0.0, alpha0, alpha0], 0.0, alpha0);  // cut off below `cold_temp`
        let elec_resi = ExplicitFunction::new(Rc::new(move |temp: f64| {
            if temp > hot_temp {
                return rho0*hot_temp/temp0;  // avoid infinitely growing value
            }
            max(rho0*temp/temp0, rho0*cold_temp/temp0)   // `elec_resi` must be positive.
        }));
        let thrm_cond = ConstantFunction::new(kappa0);
        let tep = Tep::new(Rc::new(seebeck), Rc::new(elec_resi), Rc::new(thrm_cond));
        let length = 1e-3;
        let area = 1e-6;
        let elec_cur: f64 = 1.0;
        let j: f64 = elec_cur / area;  // electrical current density (W/m^2)

        // set simulation parameters
        let minimizer_kind = MinimizerKind::Bfgs;
        let minimizer_param = MinimizerParam {
            grad_tol: 1e-5,
            a_max: 2.0,
            max_iter: 200,
            ..MinimizerParam::default(minimizer_kind)
        };
        let simul_param = SimulationParam {
            dx: length/10.0,
            dtemp: (hot_temp-cold_temp)/10.0,
            minimizer_kind,
            minimizer_param,
            ..SimulationParam::default()
        };

        // fix an exact solution
        let num_x_mesh: usize = 31;
        let x_mesh: ColVec = Array::linspace(0.0, length, num_x_mesh);
        let omega: f64 = (rho0/(kappa0*temp0)).sqrt() * j;
        let cos_omega_x: ColVec = x_mesh.mapv(|x| (omega*x).cos());
        let sin_omega_x: ColVec = x_mesh.mapv(|x| (omega*x).sin());
        let temp_exact: ColVec = hot_temp * &cos_omega_x + (cold_temp - hot_temp*(omega*length).cos())/(omega*length).sin() * &sin_omega_x;
        let dtempdx_exact: ColVec = -hot_temp * omega * &sin_omega_x + (cold_temp - hot_temp*(omega*length).cos())/(omega*length).sin() * omega * &cos_omega_x;
        let heat_flux_exact: ColVec = alpha0*&temp_exact*j - kappa0*&dtempdx_exact;

        // create a leg
        let leg = SingleLeg::new(&tep, length, area, simul_param);

        // case 1: fixed temperatures
        let left_bc = FixedTempBc::new(hot_temp);
        let right_bc = FixedTempBc::new(cold_temp);
        let simul_result_fixed_elec_cur = leg.simulate_fixed_elec_cur(elec_cur, &left_bc, &right_bc, hot_temp, cold_temp);
        assert!(simul_result_fixed_elec_cur.success);
        // test helper functions
        let load_resistance = simul_result_fixed_elec_cur.load_resistance();
        assert!(is_close(simul_result_fixed_elec_cur.total_resistance(), simul_result_fixed_elec_cur.internal_resistance + load_resistance, 1e-04, 0.0));
        assert!(is_close(simul_result_fixed_elec_cur.load_ratio(), load_resistance / simul_result_fixed_elec_cur.internal_resistance, 1e-04, 1.0));
        // check if `simulate_fixed_load_resistance` works
        // let simul_result_fixed_load_resistance = leg.simulate_fixed_load_resistance_with_manual_init(load_resistance, &left_bc, &right_bc, 0.0, 2.0, hot_temp, 1000.0);
        let simul_result_fixed_load_resistance = leg.simulate_fixed_load_resistance(load_resistance, &left_bc, &right_bc, 0.5, hot_temp, cold_temp);
        println!("dirichlet; numerical elec_cur={}, exact elec_cur={}", simul_result_fixed_load_resistance.elec_cur, elec_cur);
        println!("dirichlet; numerical load resistance={}, exact load resistance={}", simul_result_fixed_load_resistance.load_resistance(), load_resistance);
        println!("dirichlet; result={}", simul_result_fixed_load_resistance);
        assert!(simul_result_fixed_load_resistance.success);
        assert!(is_close(simul_result_fixed_load_resistance.elec_cur, elec_cur, 2e-03, 0.0));  // 0.2% rel. error
        assert!(is_close(simul_result_fixed_load_resistance.load_resistance(), load_resistance, 2e-03, 0.0));

        // case 2: RTG BC
        let heat_fluxl: f64 = heat_flux_exact[0];
        let fluid_temp: f64 = cold_temp - 100.0;
        let h: f64 = heat_flux_exact[num_x_mesh-1] / (temp_exact[num_x_mesh-1] - fluid_temp);
        let left_bc = FixedHeatFluxBc::new(heat_fluxl);
        let right_bc = ConvectionBc::new(h, fluid_temp);
        let elec_cur: f64 = 1.2;
        let simul_result_fixed_elec_cur = leg.simulate_fixed_elec_cur(elec_cur, &left_bc, &right_bc, hot_temp, cold_temp);
        assert!(simul_result_fixed_elec_cur.success);
        // println!("rtg templ={}, tempr={}, heat_ratel={}, heat_rater={}", simul_result_fixed_elec_cur.templ, simul_result_fixed_elec_cur.tempr, simul_result_fixed_elec_cur.heat_fluxl*leg.area, simul_result_fixed_elec_cur.heat_fluxr*leg.area);
        println!("rtg result={}", simul_result_fixed_elec_cur);
        let load_resistance = simul_result_fixed_elec_cur.load_resistance();
        // let simul_result_fixed_load_resistance = leg.simulate_fixed_load_resistance_with_manual_init(load_resistance, &left_bc, &right_bc, 0.0, 2.0, 700.0, 3600.0);
        // let simul_result_fixed_load_resistance = leg.simulate_fixed_load_resistance(load_resistance, &left_bc, &right_bc, 1.0, hot_temp, cold_temp);
        let simul_result_fixed_load_resistance = leg.simulate_fixed_load_resistance(load_resistance, &left_bc, &right_bc, elec_cur*0.9, hot_temp, cold_temp);
        // let simul_result_fixed_load_resistance = leg.obsolete_simulate_fixed_load_resistance(load_resistance, &left_bc, &right_bc, elec_cur*2.0, hot_temp, cold_temp);
        println!("rtg; numerical elec_cur={}, exact elec_cur={}", simul_result_fixed_load_resistance.elec_cur, elec_cur);
        println!("rtg; numerical load resistance={}, exact load resistance={}", simul_result_fixed_load_resistance.load_resistance(), load_resistance);
        println!("rtg; result={}", simul_result_fixed_load_resistance);
        assert!(simul_result_fixed_load_resistance.success);
        assert!(is_close(simul_result_fixed_load_resistance.elec_cur, elec_cur, 2e-03, 0.0));  // 0.2% rel. error
        assert!(is_close(simul_result_fixed_load_resistance.load_resistance(), load_resistance, 2e-03, 0.0));
    }

    /// Test using logarithmic Seebeck coefficient.
    /// The example is from p.112 in [C. Goupil (ed.)](https://onlinelibrary.wiley.com/doi/book/10.1002/9783527338405).
    #[test]
    #[ignore]
    fn test_single_leg_simulate_fixed_elec_cur_for_log_seebeck() {
        // environment condition
        let hot_temp: f64 = 500.0;
        let cold_temp: f64 = 300.0;

        // set tep and leg parameters
        let tau0: f64 = 150e-6;
        let rho0: f64 = 1e-5;
        let kappa0: f64 = 2.0;
        let seebeck = ExplicitFunction::new(Rc::new(move |temp: f64| { tau0*(max(temp/cold_temp, 1.0)).ln() }));   // constant Thomson
        let elec_resi = ConstantFunction::new(rho0);
        let thrm_cond = ConstantFunction::new(kappa0);
        let tep = Tep::new(Rc::new(seebeck), Rc::new(elec_resi), Rc::new(thrm_cond));
        let length = 1e-3;
        let area = 1e-6;
        let elec_cur: f64 = 1.0;
        let j: f64 = elec_cur / area;  // electrical current density (W/m^2)

        // set simulation parameters
        let minimizer_kind = MinimizerKind::Bfgs;
        let minimizer_param = MinimizerParam {
            grad_tol: 1e-5,
            max_iter: 100,
            ..MinimizerParam::default(minimizer_kind)
        };
        let simul_param = SimulationParam {
            dx: length/10.0,
            minimizer_kind,
            minimizer_param,
            ..SimulationParam::default()
        };

        // fix an exact solution
        let num_x_mesh: usize = 31;
        let x_mesh: ColVec = Array::linspace(0.0, length, num_x_mesh);
        let temp_exact: ColVec = hot_temp + ((cold_temp-hot_temp) - (rho0*j/tau0)*length)/(1.0-(tau0*j*length/kappa0).exp()) * (1.0 - x_mesh.mapv(|x| { (tau0*j*x/kappa0).exp() })) + rho0*j/tau0*&x_mesh;
        assert!(all_close(&array![temp_exact[0], temp_exact[num_x_mesh-1]], &array![hot_temp, cold_temp], 1e-04, 0.0));

        // create a leg
        let leg = SingleLeg::new(&tep, length, area, simul_param);

        // check the temperature
        let simul_result = leg.simulate_fixed_elec_cur(elec_cur, &FixedTempBc::new(hot_temp), &FixedTempBc::new(cold_temp), hot_temp, cold_temp);
        assert!(simul_result.success);
        assert!(all_close(&simul_result.temp.ats_incr(&x_mesh), &temp_exact, 1e-04, 0.0));
    }

    /// Test using linearly increasing electrical resistivity.
    /// The example is from p.123 in [C. Goupil (ed.)](https://onlinelibrary.wiley.com/doi/book/10.1002/9783527338405).
    #[test]
    #[ignore]
    fn test_single_leg_simulate_fixed_elec_cur_for_linear_elec_resi() {
        // environment condition
        let hot_temp: f64 = 500.0;
        let cold_temp: f64 = 300.0;

        // set tep and leg parameters
        let alpha0: f64 = 150e-6;
        let rho0: f64 = 1e-5;
        let kappa0: f64 = 2.0;
        let temp0: f64 = 400.0;
        let seebeck = PiecewiseLinearConstExt::new(&array![0.0, cold_temp, hot_temp], &array![0.0, alpha0, alpha0], 0.0, alpha0); // cut off below `cold_temp`
        let elec_resi = ExplicitFunction::new(Rc::new(move |temp: f64| { max(rho0*temp/temp0, rho0*cold_temp/temp0) }));
        let thrm_cond = ConstantFunction::new(kappa0);
        let tep = Tep::new(Rc::new(seebeck), Rc::new(elec_resi), Rc::new(thrm_cond));
        let length = 1e-3;
        let area = 1e-6;
        let elec_cur: f64 = 1.0;
        let j: f64 = elec_cur / area;  // electrical current density (W/m^2)

        // set simulation parameters
        // let minimizer_kind = MinimizerKind::DampedBfgs;
        let minimizer_kind = MinimizerKind::Bfgs;
        let minimizer_param = MinimizerParam {
            grad_tol: 1e-5,
            max_iter: 100,
            ..MinimizerParam::default(minimizer_kind)
        };
        let simul_param = SimulationParam {
            dx: length/10.0,
            minimizer_kind,
            minimizer_param,
            ..SimulationParam::default()
        };

        // fix an exact solution
        let num_x_mesh: usize = 31;
        let x_mesh: ColVec = Array::linspace(0.0, length, num_x_mesh);
        let omega: f64 = (rho0/(kappa0*temp0)).sqrt() * j;
        let cos_omega_x: ColVec = x_mesh.mapv(|x| (omega*x).cos());
        let sin_omega_x: ColVec = x_mesh.mapv(|x| (omega*x).sin());
        let temp_exact: ColVec = hot_temp * &cos_omega_x + (cold_temp - hot_temp*(omega*length).cos())/(omega*length).sin() * &sin_omega_x;
        let dtempdx_exact: ColVec = -hot_temp * omega * &sin_omega_x + (cold_temp - hot_temp*(omega*length).cos())/(omega*length).sin() * omega * &cos_omega_x;
        let heat_flux_exact: ColVec = alpha0*&temp_exact*j - kappa0*&dtempdx_exact;

        // create a leg
        let leg = SingleLeg::new(&tep, length, area, simul_param);

        // case 1: fixed temperatures
        let now = Instant::now();
        let dirichlet_simul_result = leg.simulate_fixed_elec_cur(elec_cur, &FixedTempBc::new(hot_temp), &FixedTempBc::new(cold_temp), hot_temp, cold_temp);
        println!("test_single_leg_simulate_fixed_elec_cur_for_linear_elec_resi--case 1: Elapsed {:.2?}", now.elapsed());
        if !dirichlet_simul_result.success {
            println!("dirichlet failed: error={:?}", dirichlet_simul_result.error);
        }
        assert!(dirichlet_simul_result.success);
        assert!(all_close(&dirichlet_simul_result.temp.ats_incr(&x_mesh), &temp_exact, 1e-03, 0.0));
        assert!(all_close(&dirichlet_simul_result.heat_flux.ats_incr(&x_mesh), &heat_flux_exact, 1e-03, 0.0));

        // case 2: RTG BC
        let heat_fluxl: f64 = heat_flux_exact[0];
        let fluid_temp: f64 = cold_temp - 100.0;
        let h: f64 = heat_flux_exact[num_x_mesh-1] / (temp_exact[num_x_mesh-1] - fluid_temp);
        let now = Instant::now();
        let rtg_simul_result = leg.simulate_fixed_elec_cur(elec_cur, &FixedHeatFluxBc::new(heat_fluxl), &ConvectionBc::new(h, fluid_temp), hot_temp, cold_temp);
        println!("test_single_leg_simulate_fixed_elec_cur_for_linear_elec_resi--case 2: Elapsed {:.2?}", now.elapsed());
        if !rtg_simul_result.success {
            println!("rtg failed: error={:?}", rtg_simul_result.error);
        }
        assert!(rtg_simul_result.success);
        println!("rtg numerical temp = {}", rtg_simul_result.temp.ats_incr(&x_mesh));
        println!("rtg exact temp = {}", temp_exact);
        println!("rtg numerical heat_flux = {}", rtg_simul_result.heat_flux.ats_incr(&x_mesh));
        println!("rtg exact heat_flux = {}", heat_flux_exact);
        assert!(all_close(&rtg_simul_result.temp.ats_incr(&x_mesh), &temp_exact, 1e-03, 0.0));
        assert!(all_close(&rtg_simul_result.heat_flux.ats_incr(&x_mesh), &heat_flux_exact, 1e-03, 0.0));
    }

    /// Test using linearly increasing thermal conductivity.
    /// The example is from p.125 in [C. Goupil (ed.)](https://onlinelibrary.wiley.com/doi/book/10.1002/9783527338405).
    #[test]
    #[ignore]
    fn test_single_leg_simulate_fixed_elec_cur_for_linear_thrm_cond() {
        // environment condition
        let hot_temp: f64 = 500.0;
        let cold_temp: f64 = 300.0;

        // set tep and leg parameters
        let alpha0: f64 = 150e-6; // [V/K]
        let rho0: f64 = 1e-5; // [Ohm m]
        let kappa0: f64 = 2.0; // [W/m/K]
        let temp0: f64 = 300.0; // [K]
        let seebeck = ConstantFunction::new(alpha0);
        let elec_resi = ConstantFunction::new(rho0);
        let thrm_cond = ExplicitFunction::new(Rc::new(move |temp: f64| { max(kappa0*temp/temp0, kappa0*cold_temp/temp0) }));   // `thrm_cond` must be positive.
        let tep = Tep::new(Rc::new(seebeck), Rc::new(elec_resi), Rc::new(thrm_cond));
        let length = 1e-3; // [m]
        let area = 1e-6;  // m^2]
        let elec_cur: f64 = 1.0;
        let j: f64 = elec_cur / area;  // electrical current density (W/m^2)

        // set simulation parameters
        let minimizer_kind = MinimizerKind::Bfgs;
        let minimizer_param = MinimizerParam {
            grad_tol: 1e-5,
            max_iter: 100,
            ..MinimizerParam::default(minimizer_kind)
        };
        let simul_param = SimulationParam {
            dx: length/10.0,
            minimizer_kind,
            minimizer_param,
            ..SimulationParam::default()
        };

        // fix an exact solution
        let num_x_mesh: usize = 31;
        let x_mesh: ColVec = Array::linspace(0.0, length, num_x_mesh);
        let temp_exact: ColVec = ( hot_temp.powi(2) + (cold_temp.powi(2)-hot_temp.powi(2))*&x_mesh/length + j.powi(2)*temp0*rho0/kappa0*&x_mesh*(length-&x_mesh) ).mapv(|e| {e.sqrt()});
        assert!(all_close(&array![temp_exact[0], temp_exact[num_x_mesh-1]], &array![hot_temp, cold_temp], 0.0, 0.0)); // the BC is exact

        // create a leg
        let leg = SingleLeg::new(&tep, length, area, simul_param);

        // check the temperature
        let now = Instant::now();
        let simul_result = leg.simulate_fixed_elec_cur(elec_cur, &FixedTempBc::new(hot_temp), &FixedTempBc::new(cold_temp), hot_temp, cold_temp);
        println!("test_single_leg_simulate_fixed_elec_cur_for_linear_thrm_cond: Elapsed {:.2?}", now.elapsed());
        if !simul_result.success {
            println!("error: {:?}", simul_result.error);
        }
        assert!(simul_result.success);
        println!("result = {}", simul_result);
        println!("num_iter = {}", simul_result.num_iter);
        println!("fixed cur, linear_thrm_cond={}", simul_result);
        println!("temps={}", simul_result.temp.ats_incr(&x_mesh));
        assert!(all_close(&simul_result.temp.ats_incr(&x_mesh), &temp_exact, 2e-04, 0.0));  // 16 iterations
    }

    #[test]
    fn test_get_num_mesh() {
        assert_eq!(get_num_mesh(10.0, 1.0, 5, 11), 11);
        assert_eq!(get_num_mesh(10.0, 10.0, 3, 11), 5);
    }
}