function all_stable = verifyStabilization(sol_matrix, t_array, time_fraction)
%VERIFYSTABILIZATION - Checks if the solution reached a stabilized status comparing the final time point with a manually specified time point
% Beware that this can have false positives when the analyzed solution had
% a tmax so small that the ionic movement did not even start.
% A false negative can be obtained when the relative tolerance of the pdepe
% solver is larger than the default 1e-3.
%
% Syntax:  all_stable = verifyStabilization(sol_matrix, t_array, time_fraction)
%
% Inputs:
%   SOL_MATRIX - a solution matrix, not the full struct, as contained in
%     structs created by DF
%   T_ARRAY - the time mesh array, as generated by meshgen_t and included
%     in structs created by DF
%   TIME_FRACTION - numeric between 0 and 1. The smaller the stricter the
%     check. When 1: disable the stability check and simply return true.
%     When zero or strictly included between 0 and 1: check if the
%     solution reached a stabilization. Suggested values: 1e-3 for
%     normal cases, 0 for cases in which the stabilization is very
%     important, 1 for cases in which the first solution is needed.
%
% Outputs:
%   ALL_STABLE - a logic asserting if the solution in input reached its
%     stability
%
% Example:
%   verifyStabilization(soleq.ion.u, soleq.ion.t, 1e-3)
%     check if the solution did not change much from a time which is one
%     thousandth of the final time point and the final time point
%
% Other m-files required: none
% Subfunctions: none
% MAT-files required: none
%
% See also DF.
%
%
%% LICENSE
% Copyright (C) 2020  Philip Calado, Ilario Gelmetti, and Piers R. F. Barnes
% Imperial College London
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU Affero General Public License as published
% by the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Author: Ilario Gelmetti, Ph.D. student, perovskite photovoltaics
% Institute of Chemical Research of Catalonia (ICIQ)
% Research Group Prof. Emilio Palomares
% email address: iochesonome@gmail.com
% October 2017; Last revision: May 2020
%
%------------- BEGIN CODE --------------

all_stable = true;

if time_fraction >= 1 || time_fraction < 0 || ~isnumeric(time_fraction)
    warning('Driftfusion:verifyStabilizationSkip', [mfilename...
        ' - time_fraction is out of [0,1) range: skipping stabilization check']);
    return
end

% name of the variables
names = ["potential", "electrons", "holes", "cations", "anions"];
% which values have to be considered in a linear or in a log10 scale
compare_log = [false, true, true, false, false];

% no need to calculate end_time for each of the 4 solutions: if they
% break they break at the same time
% using t(end) could get a time out of the solution when the computation broke before reaching the final time
end_time = t_array(length(sol_matrix(:, 1, 1)));
[~, time_index] = min(abs(t_array - time_fraction * end_time)); % get time mesh index at a specified percentage of maximum time reached

% for each component of the solution verify that didn't change significantly over the requested time
for i = 1:length(sol_matrix(1, 1, :))
    profile_at_time = sol_matrix(time_index, :, i); % take profile of values at a certain time of evolution
    profile_end = sol_matrix(end, :, i); % take profile of values at the end of time
    if compare_log(i) % for variables ranging on huge scales comparing the log values makes more sense
        profile_at_time = log10(profile_at_time);
        profile_end = log10(profile_end);
    end

    difference = sum(abs(profile_end - profile_at_time)); % sum up all the differences between the profiles
    % if the values change more than one tenthousandth, the solution is
    % considered not stable
    threshold = 1e-5 * sum(abs(profile_end - mean(profile_end))); % sum up absolute values, ignore constant bias
    stable = difference <= threshold;

    if stable
        %display(['variable ', num2str(i), ' stabilisation verified']);        
    else
        warning('Driftfusion:verifyStabilization',...
            'Comparing final solutions at %s s and %s s showed that the %s distribution did not reach stability. Consider trying with a greater tmax.',...
            num2str(t_array(time_index)), num2str(t_array(end)), names(i));
    end
    % true just if all the variables are stabilized
    all_stable = all_stable && stable;
end

if all_stable
    disp("Stabilisation verified");
end

%------------- END OF CODE --------------