Current densities

def EIS_current(t, parameters):
    """
    Represents a current density used for creating an EIS curve and Bode diagrams.
    The current density is first equilibrated at i_EIS A.m-2 from 0 to t0_EIS seconds using a step increase.
    Then, a sinusoidal perturbation is added to the current density. This perturbation has an amplitude of
    (ratio_EIS * i_EIS) A.m-2 and a frequency of f[n_inf] Hz.

    Parameters:
    ----------
    t : float
        Time in seconds.
    parameters : dict
        A dictionary containing the parameters for the current density function.

    Returns:
    -------
    i_fc : float
        The polarization current density at time t.
    """

    # Initialisation
    i_EIS, ratio_EIS = parameters['i_EIS'], parameters['ratio_EIS']  # (A/m², ). i_EIS is the current for which a
    #                                                                  ratio_EIS perturbation is added.
    t0_EIS, t_new_start_EIS, tf_EIS, delta_t_break_EIS, delta_t_measurement_EIS = parameters['t_EIS']  # It is the initial
    #         EIS time after stack equilibrium, a list of time parameters which gives the beginning of each frequency
    #         change, the final time, a list of time parameters which gives the estimated time for reaching equilibrium
    #         at each frequency, and a list of time parameters which gives the estimated time for measuring the voltage
    #         response at each frequency.
    f_power_min_EIS, f_power_max_EIS, nb_f_EIS, nb_points_EIS = parameters['f_EIS']  # It is the power of the initial
    #         frequency: f_min_EIS = 10**f_power_min_EIS, the power of the final frequency, the number of frequencies
    #         tested and the number of points calculated per specific period.
    f = np.logspace(f_power_min_EIS, f_power_max_EIS, num=nb_f_EIS)  # It is a list of all the frequency tested,
    #                                                      ranged logarithmically.

    # Current density for the EIS curve
    if t < t0_EIS:
        delta_t_ini = t0_EIS / 4  # s. It is the required time for elevating i_fc from 0 to i_EIS without starving the
        #                              cell.
        i_fc = i_EIS * (1.0 + np.tanh(4 * (t - 2 * (delta_t_ini / 2)) / delta_t_ini)) / 2
    else:
        n_inf = np.where(t_new_start_EIS <= t)[0][-1]  # It is the number of frequency changes which has been made so far.
        i_disruption = (ratio_EIS * i_EIS) * np.cos(2 * np.pi * f[n_inf] * t)
        i_fc = i_EIS + i_disruption

    return i_fc

def polarization_current(t, parameters):
    """Represents a current density used for creating a polarization curve.
    Starting from 0, the current density increases by the value of delta_i_pola every delta_t, following C∞ step current
    increments, until it reaches i_max_pola. Each increment lasts for delta_t_load_pola seconds. After each increment,
    there is a pause of delta_t_break_pola seconds to allow the stack to reach equilibrium.

    Parameters:
    ----------
    t : float
        Time in seconds.
    parameters : dict
        A dictionary containing the parameters for the current density function.

    Returns:
    -------
    i_fc : float
        The polarization current density at time t.
    """

    # Initialisation
    delta_t_load_pola, delta_t_break_pola, delta_i_pola, delta_t_ini_pola = parameters['delta_pola']
    #                                 (s, s, A.m-2, s). It is the loading time, the breaking time, the current density
    #                                 step, and the initial breaking time.
    i_max_pola = parameters['i_max_pola']  # A.m-2. It is the maximum current density for the polarization curve.
    delta_t = delta_t_load_pola + delta_t_break_pola  # s. It is the time of one load.
    tf = delta_t_ini_pola + int(i_max_pola / delta_i_pola + 1) * delta_t  # s. It is the duration of this polarization current.
    n = int(tf / delta_t)  # . It is the number of loads made for this polarization current.

    # Current density for the polarization curve
    i_fc = 0  # A.m-2. Initialisation of the current density.
    if t < delta_t_ini_pola:  # It is the initial break for having homogeneity inside the cell
        i_fc = 0  # before starting the measurements.
    else:
        for i in range(n):
            t_switch = delta_t * i  # The current density value changes around this time.
            i_fc += delta_i_pola * (1.0 + np.tanh(4 * (t - delta_t_ini_pola - delta_t - t_switch - (delta_t_load_pola / 2)) /
                                             delta_t_load_pola)) / 2

    return i_fc

def step_current(t, parameters):
    """Represents a step change in current density.
    The current starts at 0 and smoothly stabilizes at i_ini_step A.m-2 in delta_t_load seconds.
    Around t_switch seconds, the current increases smoothly and stabilizes at i_final_step A.m-2 in delta_t_load seconds.
    This is a C∞ function, which is advantageous for enhancing the overall stability of the results.

    Parameters:
    ----------
    t : float
        Time in seconds.
    parameters : dict
        A dictionary containing the parameters for the current density function.

    Returns:
    -------
    i_fc : float
        The step current density at time t.
    """

    # Initialisation
    t0_step, tf_step, delta_t_load_step, delta_t_dyn_step = parameters['t_step']  # (s, s, s, s). It is the initial, final,
    #                                                                     loading and dynamic time for display.
    i_ini_step, i_final_step = parameters['i_step']  # (A.m-2, A.m-2). It is the initial and final current density values.
    t_switch = tf_step // 2  # The current density value changes around this time.

    # Step current density
    i_fc = i_ini_step * (1.0 + np.tanh(4 * (t - 2 * (delta_t_load_step / 2)) / delta_t_load_step)) / 2 + \
           + (i_final_step - i_ini_step) * (1.0 + np.tanh(4 * (t - t_switch - (delta_t_load_step / 2)) / delta_t_load_step)) / 2

    return i_fc