function plot_irfs(M_, options_, oo_irfs, shock_list, var_list)
%PLOT_IRFS Plot IRF results using standard Dynare utilities
%
% SYNTAX:
%   heterogeneity.plot_irfs(M_, options_, oo_irfs, shock_list, var_list)
%
% INPUTS:
%   M_         [struct]  Dynare model structure
%   options_   [struct]  Dynare options structure
%   oo_irfs    [struct]  IRF structure with fields named 'var_shock' (from oo_.irfs)
%   shock_list [cell]    List of shock names to plot
%   var_list   [cell]    List of variable names to plot
%
% OUTPUTS:
%   none (creates and saves figures, optionally generates TeX file)
%
% DESCRIPTION:
%   Plots impulse response functions using standard Dynare utilities:
%   - pltorg() for automatic subplot layout
%   - dyn_figure() for figure creation respecting nodisplay option
%   - dyn_saveas() for saving in multiple formats
%   - Respects options_.impulse_responses.plot_threshold
%   - Supports TeX labels via options_.TeX
%   - Generates LaTeX figure loader when options_.TeX and EPS format requested

% Copyright © 2025 Dynare Team
%
% This file is part of Dynare.
%
% Dynare is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% Dynare is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with Dynare.  If not, see <https://www.gnu.org/licenses/>.
%
% Original author: Normann Rion <normann@dynare.org>

%% Input validation
if isempty(fieldnames(oo_irfs))
    return;
end

TeX = options_.TeX;

%% Create graphs directory if needed (for graphs or TeX output)
if ~options_.nograph || (TeX && any(strcmp('eps', cellstr(options_.graph_format))))
    if ~isfolder([M_.dname '/graphs'])
        mkdir(M_.dname, 'graphs');
    end
end

%% Open TeX file if needed
if TeX && any(strcmp('eps', cellstr(options_.graph_format)))
    fidTeX = fopen([M_.dname, '/graphs/' M_.fname '_IRF.tex'], 'w');
    fprintf(fidTeX, '%% TeX eps-loader file generated by heterogeneity.plot_irfs (Dynare).\n');
    fprintf(fidTeX, ['%% ' datestr(now, 0) '\n']);
    fprintf(fidTeX, ' \n');
else
    fidTeX = -1;
end

%% Get shock titles for captions
tit = M_.exo_names;
if TeX
    titTeX = M_.exo_names_tex;
end

%% Process each shock
for i_shock = 1:length(shock_list)
    shock_name = shock_list{i_shock};
    shock_idx = find(strcmp(shock_name, M_.exo_names));

    %% Build IRF matrix and variable list for this shock
    irfs = [];
    mylist = [];
    if TeX
        mylistTeX = [];
    end
    var_indices = [];

    for j = 1:length(var_list)
        var_name = var_list{j};
        field_name = [var_name '_' shock_name];

        if isfield(oo_irfs, field_name)
            irf_data = oo_irfs.(field_name);

            % Check threshold
            if max(abs(irf_data)) >= options_.impulse_responses.plot_threshold
                irfs = cat(1, irfs, irf_data(:)');
                if isempty(mylist)
                    mylist = var_name;
                else
                    mylist = char(mylist, var_name);
                end

                % Get variable index for TeX
                var_idx = find(strcmp(var_name, M_.endo_names), 1);
                var_indices = [var_indices; var_idx];

                if TeX
                    if ~isempty(var_idx)
                        if isempty(mylistTeX)
                            mylistTeX = M_.endo_names_tex{var_idx};
                        else
                            mylistTeX = char(mylistTeX, M_.endo_names_tex{var_idx});
                        end
                    else
                        if isempty(mylistTeX)
                            mylistTeX = var_name;
                        else
                            mylistTeX = char(mylistTeX, var_name);
                        end
                    end
                end
            elseif options_.debug
                fprintf('heterogeneity.plot_irfs: The IRF of %s to %s is smaller than the irf_plot_threshold of %4.3f and will not be displayed.\n', ...
                    var_name, shock_name, options_.impulse_responses.plot_threshold);
            end
        end
    end

    %% Skip if no IRFs to plot for this shock
    if isempty(irfs)
        continue;
    end

    %% Plot if requested
    if ~options_.nograph
        number_of_plots_to_draw = size(irfs, 1);
        [nbplt, nr, nc, lr, lc, nstar] = pltorg(number_of_plots_to_draw);

        if nbplt == 0
            continue;
        end

        T = size(irfs, 2);

        if nbplt == 1
            % Single figure
            if options_.relative_irf
                hh_fig = dyn_figure(options_.nodisplay, 'Name', ['Relative response to orthogonalized shock to ' tit{shock_idx}]);
            else
                hh_fig = dyn_figure(options_.nodisplay, 'Name', ['Orthogonalized shock to ' tit{shock_idx}]);
            end

            for j = 1:number_of_plots_to_draw
                subplot(nr, nc, j);
                plot(1:T, transpose(irfs(j,:)), '-k', 'linewidth', 1);
                hold on
                plot([1 T], [0 0], '-r', 'linewidth', 0.5);
                hold off
                xlim([1 T]);
                remove_fractional_xticks;
                if TeX
                    title(['$' deblank(mylistTeX(j,:)) '$'], 'Interpreter', 'latex');
                else
                    title(deblank(mylist(j,:)), 'Interpreter', 'none');
                end
            end

            dyn_saveas(hh_fig, [M_.dname, '/graphs/' M_.fname '_IRF_' tit{shock_idx}], options_.nodisplay, options_.graph_format);

            if fidTeX > 0
                fprintf(fidTeX, '\\begin{figure}[H]\n');
                fprintf(fidTeX, '\\centering \n');
                fprintf(fidTeX, '\\includegraphics[width=%2.2f\\textwidth]{%s_IRF_%s}\n', options_.figures.textwidth*min(j/nc, 1), [M_.dname, '/graphs/' M_.fname], tit{shock_idx});
                if options_.relative_irf
                    fprintf(fidTeX, '\\caption{Relative impulse response functions (orthogonalized shock to $%s$).}\n', titTeX{shock_idx});
                else
                    fprintf(fidTeX, '\\caption{Impulse response functions (orthogonalized shock to $%s$).}\n', titTeX{shock_idx});
                end
                fprintf(fidTeX, '\\label{Fig:IRF:%s}\n', tit{shock_idx});
                fprintf(fidTeX, '\\end{figure}\n');
                fprintf(fidTeX, ' \n');
            end
        else
            % Multiple figures
            for fig = 1:nbplt-1
                if options_.relative_irf
                    hh_fig = dyn_figure(options_.nodisplay, 'Name', ['Relative response to orthogonalized shock to ' tit{shock_idx} ' figure ' int2str(fig)]);
                else
                    hh_fig = dyn_figure(options_.nodisplay, 'Name', ['Orthogonalized shock to ' tit{shock_idx} ' figure ' int2str(fig)]);
                end

                for plt = 1:nstar
                    subplot(nr, nc, plt);
                    plot(1:T, transpose(irfs((fig-1)*nstar+plt,:)), '-k', 'linewidth', 1);
                    hold on
                    plot([1 T], [0 0], '-r', 'linewidth', 0.5);
                    hold off
                    xlim([1 T]);
                    remove_fractional_xticks;
                    if TeX
                        title(['$' deblank(mylistTeX((fig-1)*nstar+plt,:)) '$'], 'Interpreter', 'latex');
                    else
                        title(deblank(mylist((fig-1)*nstar+plt,:)), 'Interpreter', 'none');
                    end
                end

                dyn_saveas(hh_fig, [M_.dname, '/graphs/' M_.fname '_IRF_' tit{shock_idx} int2str(fig)], options_.nodisplay, options_.graph_format);

                if fidTeX > 0
                    fprintf(fidTeX, '\\begin{figure}[H]\n');
                    fprintf(fidTeX, '\\centering \n');
                    fprintf(fidTeX, '\\includegraphics[width=%2.2f\\textwidth]{%s_IRF_%s%s}\n', options_.figures.textwidth*min(plt/nc, 1), [M_.dname, '/graphs/' M_.fname], tit{shock_idx}, int2str(fig));
                    if options_.relative_irf
                        fprintf(fidTeX, '\\caption{Relative impulse response functions (orthogonalized shock to $%s$).}', titTeX{shock_idx});
                    else
                        fprintf(fidTeX, '\\caption{Impulse response functions (orthogonalized shock to $%s$).}', titTeX{shock_idx});
                    end
                    fprintf(fidTeX, '\\label{Fig:IRF:%s:%s}\n', tit{shock_idx}, int2str(fig));
                    fprintf(fidTeX, '\\end{figure}\n');
                    fprintf(fidTeX, ' \n');
                end
            end

            % Last figure (may have fewer plots)
            hh_fig = dyn_figure(options_.nodisplay, 'Name', ['Orthogonalized shock to ' tit{shock_idx} ' figure ' int2str(nbplt) '.']);
            m = 0;
            for plt = 1:number_of_plots_to_draw-(nbplt-1)*nstar
                m = m + 1;
                subplot(lr, lc, m);
                plot(1:T, transpose(irfs((nbplt-1)*nstar+plt,:)), '-k', 'linewidth', 1);
                hold on
                plot([1 T], [0 0], '-r', 'linewidth', 0.5);
                hold off
                xlim([1 T]);
                remove_fractional_xticks;
                if TeX
                    title(['$' deblank(mylistTeX((nbplt-1)*nstar+plt,:)) '$'], 'Interpreter', 'latex');
                else
                    title(deblank(mylist((nbplt-1)*nstar+plt,:)), 'Interpreter', 'none');
                end
            end

            dyn_saveas(hh_fig, [M_.dname, '/graphs/' M_.fname '_IRF_' tit{shock_idx} int2str(nbplt)], options_.nodisplay, options_.graph_format);

            if fidTeX > 0
                fprintf(fidTeX, '\\begin{figure}[H]\n');
                fprintf(fidTeX, '\\centering \n');
                fprintf(fidTeX, '\\includegraphics[width=%2.2f\\textwidth]{%s_IRF_%s%s}\n', options_.figures.textwidth*min(m/lc, 1), [M_.dname, '/graphs/' M_.fname], tit{shock_idx}, int2str(nbplt));
                if options_.relative_irf
                    fprintf(fidTeX, '\\caption{Relative impulse response functions (orthogonalized shock to $%s$).}', titTeX{shock_idx});
                else
                    fprintf(fidTeX, '\\caption{Impulse response functions (orthogonalized shock to $%s$).}', titTeX{shock_idx});
                end
                fprintf(fidTeX, '\\label{Fig:IRF:%s:%s}\n', tit{shock_idx}, int2str(nbplt));
                fprintf(fidTeX, '\\end{figure}\n');
                fprintf(fidTeX, ' \n');
            end
        end
    end
end

%% Close TeX file
if fidTeX > 0
    fprintf(fidTeX, ' \n');
    fprintf(fidTeX, '%% End Of TeX file. \n');
    fclose(fidTeX);
end

end
